MD

模拟过程中步骤和命令

蛋白和配体top文件的获取

蛋白质、DNA等大分子:

1
$ gmx pdb2gmx -f protein.pdb -o protein.gro -ignh

配体等小分子:(itp)

(1) GROMOS力场:ATB
{2} OPLS力场:TPPmktop
(3) 其他力场

组合gro和top文件

最后得到complex.gro和top文件

定义盒子并添加溶剂

定义方形盒子

1
$ gmx editconf -f complex.gro -o newbox.gro -bt cubic -d 1.0

加水模型

1
$ gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro

添加抗衡离子(成中性体系)

1
2
$ gmx grompp -f em.mdp -c solv.gro -p topol.top -o ions.tpr
$ gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral

能量最小化

1
2
$ gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr
$ gmx mdrun -v -deffnm em -nt 8

NVT和NPT平衡

平衡前需要给配体加位置限制

1
$ gmx genrestr -f lig.gro -o posre_lig.itp -fc 1000 1000 1000

并在top中加入限字段

1
2
3
4
;Ligand position restraints
#ifdef POSRES
#include "posre_lig.itp"
#endif

NVT:平衡温度

可能需要定义自己感兴趣的index文件

1
2
$ gmx grompp -f nvt.mdp -c em.gro -p topol.top -n index.ndx -o nvt.tpr
$ gmx mdrun -deffnm nvt -v -nt 8

NPT:压力平衡

1
2
$ gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -n index.ndx -o npt.tpr
$ gmx mdrun -deffnm npt -v -nt 8

成品MD

1
2
$ gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md.tpr
$ gmx mdrun -deffnm md -v -nt 8

分析

创建mdp文件脚本

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
#!/bin/bash
# a bash for generating mdp file for Protein systems
# custom main change parameters
Etmol=1000 # maximum force for EM
Rcoulomb=1.0 # cut-off of coulomb, generally different force field have different value.
Rvdw=1.0 # cut-off of vdW, generally different force field have different value.
Temp=298.0 # Temperature
NVT_time=200 # unit: ps
NPT_time=200 # unit: ps
MD_time=50000 # unit: ps
STEP_size=0.002 # unit: fs
NVT_nsteps=`awk -v NVT_time=$NVT_time -v STEP_size=$STEP_size 'BEGIN {print NVT_time/STEP_size}'`
NPT_nsteps=`awk -v NPT_time=$NPT_time -v STEP_size=$STEP_size 'BEGIN {print NPT_time/STEP_size}'`
MD_nsteps=`awk -v MD_time=$MD_time -v STEP_size=$STEP_size 'BEGIN {print MD_time/STEP_size}'`


# creat em.mdp
cat << EOF > em.mdp
; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = $Etmol ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Minimization step size
nsteps = 50000 ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; Method to determine neighbor list (simple, grid)
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = $Rcoulomb ; Short-range electrostatic cut-off
rvdw = $Rvdw ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions
EOF

# creat nvt.mdp
cat << EOF > nvt.mdp
title = NVT equilibration
define = -DPOSRES ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = $NVT_nsteps ;
dt = $STEP_size ; fs
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rcoulomb = $Rcoulomb ; short-range electrostatic cutoff (in nm)
rvdw = $Rvdw ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = $Temp $Temp ; reference temperature, one for each group, in K
; Pressure coupling is off
pcoupl = no ; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = $Temp ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed
EOF

# creat npt.mdp
cat << EOF > npt.mdp
title = NPT equilibration
define = -DPOSRES ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = $NPT_nsteps ;
dt = $STEP_size ; fs
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
; Bond parameters
continuation = yes ; Restarting after NVT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = $Rcoulomb ; short-range electrostatic cutoff (in nm)
rvdw = $Rvdw ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = $Temp $Temp ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = no ; Velocity generation is off
EOF

# creat md.mdp
cat << EOF > md.mdp
title = Product MD
; Run parameters
integrator = md ; leap-frog integrator
nsteps = $MD_nsteps ;
dt = $STEP_size ; fs
; Output control
nstxout = 500 ; suppress bulky .trr file by specifying
nstvout = 500 ; 0 for output frequency of nstxout,
nstfout = 500 ; nstvout, and nstfout
nstenergy = 500 ; save energies every 10.0 ps
nstlog = 500 ; update log file every 10.0 ps
; nstxout-compressed = 500 ; save compressed coordinates every 10.0 ps
; compressed-x-grps = System ; save the whole system
; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = $Rcoulomb ; short-range electrostatic cutoff (in nm)
rvdw = $Rvdw ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = $Temp $Temp ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
EOF