Back

gromacs下蛋白-配体复合物模拟[GAFF]

gromacs小分子的配体的电荷问题一直时gromacs模拟蛋白-配体复合物中的难点与重点~网上有许多相关方面的教程,但是都没有成为系统,应李继存老师之邀,把最近摸索的方法过程总结出来,分享给大家~感谢网络上的朋友对我的帮助~

我们以我现在在做的蛋白配体复合物进行模拟,蛋白质为AhR的PAS-B结构域,其是AhR转录因子的配体结合区~TCDD复合物是AhR的经典激动剂。

由于AhR晶体结构并未解析,建模的方法我在这里不累述

你可以在这里下载初始文件:tcddPhe289Ala.tar

mdp文件建议从压缩包中的进行修改,以免直接复制导致换行问题。

1.高斯输入文件建立

1
antechamber -i tcdd-h.mol2 -fi mol2 -o tcdd.gjf -fo gcrt -pf y -gm "%mem=1024MB" -gn "%nproc=2" -nc 0 -gk "#HF/6-31G*SCF=tight Test Pop=MK iop(6/33=2) iop(6/42=6) iop(6/50=1) opt nosymm" -ge tcdd_resp.gesp -gv 1

其中nosymm是是否保留坐标,可以不加。高斯windows版本内存最多使用1G

2.运行高斯

1
g09 tcdd.gjf

3.拟合成resp文件

1
antechamber -i tcdd_resp.gesp -fi gesp -o tcdd.mol2 -fo mol2 -pf y -c resp -rn TCD

4.生成缺失参数文件

1
parmchk2 -i tcdd.mol2 -f mol2 -o tcdd.frcmod

5.生成AMBER参数文件及坐标文件

1
2
3
4
5
6
7
source leaprc.ff99SB
 source leaprc.gaff
 loadamberparams tcdd.frcmod
 lig=loadmol2 tcdd.mol2
 check lig
 saveamberparm lig tcdd.prmtop tcdd.inpcrd
 quit

命名为leap.in

1
tleap -f leap.in

6.将AMBER文件转换为GROMACS文件

这里有两种办法amb2gmx.pl和acpype.py,一般使用后者

命令如下:

1
python2.7 acpype -p tcdd.prmtop -x tcdd.inpcrd -d

将会得到MOL_GMX.gro和MOL_GMX.top两个文件为下面做准备

7.准备蛋白质的拓扑

1
gmx pdb2gmx -f mAhR.pdb -o mAhR_processed.gro -water tip3p

注:群里说tip3p准确度比spec高

力场选择amber99SB 5

8.准备小分子拓扑结构

将得到的TCD_GMX.top中顶部的[defaults],[ atomtypes ],底部的[ system ],[ molecules ]删除,改名为TCD.itp。记得[ moleculetype ]需要被保留!!

复制mAhR_processed.gro,重命名为complex.gro 将TCD_GMX.gro中内容:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
1  TCD  CL1    1   1.163  -4.344   0.912
1  TCD  CL2    2   1.223  -4.088   0.726
1  TCD  CL3    3   2.103  -4.309   1.263
1  TCD  CL4    4   2.164  -4.054   1.077
1  TCD   O1    5   1.637  -4.310   1.076
1  TCD   O2    6   1.690  -4.088   0.913
1  TCD   C1    7   1.542  -4.258   0.994
1  TCD   C2    8   1.568  -4.148   0.914
1  TCD   C3    9   1.759  -4.250   1.075
1  TCD   C4   10   1.785  -4.140   0.995
1  TCD   C5   11   1.417  -4.317   0.992
1  TCD   C6   12   1.469  -4.098   0.832
1  TCD   C7   13   1.857  -4.300   1.157
1  TCD   C8   14   1.909  -4.081   0.997
1  TCD   C9   15   1.318  -4.266   0.911
1  TCD  C10   16   1.344  -4.156   0.830
1  TCD  C11   17   1.983  -4.242   1.159
1  TCD  C12   18   2.009  -4.132   1.079
1  TCD   H1   19   1.399  -4.402   1.055
1  TCD   H2   20   1.490  -4.013   0.771
1  TCD   H3   21   1.836  -4.385   1.218
1  TCD   H4   22   1.928  -3.996   0.934

添加到complex.gro底部数字行上面将顶部1741加22个原子改为1763,保存。

在topol文件中添加小分子的位置限制:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Include ligand topology 
#include "TCD.itp" 

; Include water topology
#include "amber99sb.ff/tip3p.itp"

[ molecules ]
; Compound        #mols
Protein             1
TCD                 1

[注:增加了itp文件,以及molecules增加了小分子,即:]

1
2
; Include ligand topology 
#include "TCD.itp"
1
TCD                 1

9.定义单元盒子

1
gmx editconf -f complex.gro -o newbox.gro -bt cubic -d 1.2
1
gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro

用vmd观察蛋白体系是否在盒子内部

vmdscene

个人感觉盒子可以调小一些~

10.添加离子

教程使用了能量最小化的文件,我这里稍微改了一下

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
; LINES STARTING WITH ';' ARE COMMENTS
title       = Minimization  ; 标题
define          = -DFLEXIBLE   ;使用柔性水模型而非刚性模型, 这样最陡下降法可进一步最小化能量

; Parameters describing what to do, when to stop and what to save
; 参数描述如何做,当模拟结束后如何保存
integrator  = steep     ; Algorithm (steep = steepest descent minimization) 最陡梯度法
emtol       = 1000.0    ; Stop minimization when the maximum force < 10.0 kJ/mol    当小于多少能量以后停止最小化
emstep      = 0.01      ; Energy step size  能量步大小
nsteps      = 50000     ; Maximum number of (minimization) steps to perform  最大最小化步骤
energygrps  = system    ; Which energy group(s) to write to disk   哪个组进行能量最小化

; 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 更新近邻列表的频率. 1表示每步都更新
cutoff-scheme   = Verlet           ; With cutoff-scheme=Verlet and verlet-buffer-tolerance set, nstlist is actually a minimum value and mdrun might increase it, unless it is set to 1 这样设置更加精确(看英文貌似是的)
ns_type         = grid      ; Method to determine neighbor list (simple, grid)
rlist           = 1.0       ; Cut-off for making neighbor list (short range forces) 短程力近邻列表的截断值
coulombtype     = PME       ; Treatment of long range electrostatic interactions 计算长程静电的方法. PME为粒子网格Ewald方法, 还可以使用cut-off
rcoulomb        = 1.0       ; long range electrostatic cut-off   长程库仑力截断值
rvdw            = 1.0       ; long range Van der Waals cut-off  范德华距离截断值
constraints         = none             ;无位置限制
pbc             = xyz       ; Periodic Boundary Conditions  周期性

我们这里还需要查看topol文件,看一下我们的体系电荷多少,需要补足电荷,可以看出我们需要补4个CL离子

topol

1
gmx grompp -f em.mdp -c solv.gro -p topol.top -o ions.tpr

做到这里会报错,我找了很久原因,后来查看amber力场文件,发现是大小写不同导致的 我们来看一下我们的TCD.itp文件

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
 1   Cl     1   TCD   CL1    1    -0.072622     35.45000 ; qtot -0.073
 2   Cl     1   TCD   CL2    2    -0.072622     35.45000 ; qtot -0.145
 3   Cl     1   TCD   CL3    3    -0.072622     35.45000 ; qtot -0.218
 4   Cl     1   TCD   CL4    4    -0.072622     35.45000 ; qtot -0.290
 5   OS     1   TCD    O1    5    -0.320484     16.00000 ; qtot -0.611
 6   OS     1   TCD    O2    6    -0.320484     16.00000 ; qtot -0.931
 7   CA     1   TCD    C1    7     0.247670     12.01000 ; qtot -0.684
 8   CA     1   TCD    C2    8     0.247670     12.01000 ; qtot -0.436
 9   CA     1   TCD    C3    9     0.247670     12.01000 ; qtot -0.188
10   CA     1   TCD    C4   10     0.247670     12.01000 ; qtot 0.059
11   CA     1   TCD    C5   11    -0.253268     12.01000 ; qtot -0.194
12   CA     1   TCD    C6   12    -0.253268     12.01000 ; qtot -0.447
13   CA     1   TCD    C7   13    -0.253268     12.01000 ; qtot -0.701
14   CA     1   TCD    C8   14    -0.253268     12.01000 ; qtot -0.954
15   CA     1   TCD    C9   15     0.041544     12.01000 ; qtot -0.912
16   CA     1   TCD   C10   16     0.041544     12.01000 ; qtot -0.871
17   CA     1   TCD   C11   17     0.041544     12.01000 ; qtot -0.829
18   CA     1   TCD   C12   18     0.041544     12.01000 ; qtot -0.788
19   HA     1   TCD    H1   19     0.196918      1.00800 ; qtot -0.591
20   HA     1   TCD    H2   20     0.196918      1.00800 ; qtot -0.394
21   HA     1   TCD    H3   21     0.196918      1.00800 ; qtot -0.197
22   HA     1   TCD    H4   22     0.196918      1.00800 ; qtot -0.000

然而amber99SB力场下如是写的:

1
Cl          17      35.45    0.0000  A   4.40104e-01  4.18400e-01

以下有两种方法解决: 方法一:

我们需要将top文件的[ atomtypes ]可以写入相应力场的ffnonbonded.itp中,这个方法更加的实用。此方法不建议再使用

个人认为比较稳妥且不会污染力场的方法为将力场文件复制出来,然后在当前目录中进行修改,如果采用apt-get安装,可以在usr/share/gromacs/top下查找对应力场文件,例如拷贝amber99sb-ildn-ligand.ff力场文件,重命名为amber99sb-ildn-ligand.ff,将; Include forcefield parameters下面一行修改为:

1
#include "amber99sb-ildn-ligand.ff/forcefield.itp"

topol文件其它引用力场的位置可修改可不修改,个人觉得最好修改

方法二:(建议使用该方法)

将[ atomtypes ]写入拓扑中,位于

1
2
; Include forcefield parameters
#include "amber99sb-ildn.ff/forcefield.itp"

1
2
3
[ moleculetype ]
; Name            nrexcl
Protein             3

之间,例如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
; Include forcefield parameters
#include "amber99sb-ildn.ff/forcefield.itp"

[ atomtypes ]
;name   bond_type     mass     charge   ptype   sigma         epsilon       Amb
 c3       c3          0.00000  0.00000   A     3.39967e-01   4.57730e-01 ; 1.91  0.1094
 nh       nh          0.00000  0.00000   A     3.25000e-01   7.11280e-01 ; 1.82  0.1700
 hn       hn          0.00000  0.00000   A     1.06908e-01   6.56888e-02 ; 0.60  0.0157
 cc       cc          0.00000  0.00000   A     3.39967e-01   3.59824e-01 ; 1.91  0.0860
 ca       ca          0.00000  0.00000   A     3.39967e-01   3.59824e-01 ; 1.91  0.0860
 na       na          0.00000  0.00000   A     3.25000e-01   7.11280e-01 ; 1.82  0.1700
 cd       cd          0.00000  0.00000   A     3.39967e-01   3.59824e-01 ; 1.91  0.0860
 hc       hc          0.00000  0.00000   A     2.64953e-01   6.56888e-02 ; 1.49  0.0157
 h1       h1          0.00000  0.00000   A     2.47135e-01   6.56888e-02 ; 1.39  0.0157
 ha       ha          0.00000  0.00000   A     2.59964e-01   6.27600e-02 ; 1.46  0.0150
 h4       h4          0.00000  0.00000   A     2.51055e-01   6.27600e-02 ; 1.41  0.0150

[ moleculetype ]
; Name            nrexcl
Protein             3

2

然后就是添加离子了 现在有一种快速添加的方法,即平衡为中性,如下:

1
 gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral

当然还有手动平衡的方法如下:

1
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -nn 4

-nn 表示的添加阴离子数目

选择15,添加在水溶剂里面

会有如下话:

1
2
3
4
5
Back Off! I just backed up topol.top to ./#topol.top.2#
Replacing solvent molecule 1173 (atom 5282) with CL
Replacing solvent molecule 3252 (atom 11519) with CL
Replacing solvent molecule 2563 (atom 9452) with CL
Replacing solvent molecule 10648 (atom 33707) with CL

表示插入进去了

11.能量最小化

参数文件如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
; LINES STARTING WITH ';' ARE COMMENTS
title       = Minimization  ; Title of run
define          = -DFLEXIBLE

; Parameters describing what to do, when to stop and what to save
integrator  = steep     ; Algorithm (steep = steepest descent minimization)
emtol       = 500.0     ; Stop minimization when the maximum force < 5.0 kJ/mol
emstep      = 0.01      ; Energy step size
nsteps      = 50000     ; Maximum number of (minimization) steps to perform
energygrps  = Protein TCD   ; Which energy group(s) to write to disk

; 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
ns_type         = grid      ; Method to determine neighbor list (simple, grid)
rlist           = 1.0       ; Cut-off for making neighbor list (short range forces)
coulombtype     = PME       ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0       ; long range electrostatic cut-off
rvdw            = 1.0       ; long range Van der Waals cut-off
constraints         = none
pbc             = xyz       ; Periodic Boundary Conditions

**能量组为energygrps = Protein TCD,这里TCD需要根据自己的体系进行修改,下同 **

1
energygrps   = Protein TCD
1
gmx grompp -f em_real.mdp -c solv_ions.gro -p topol.top -o em.tpr

运行模拟

1
gmx mdrun -v -deffnm em

个人认为在运行模拟之前要搞懂参数的意义,所幸这里有中文参考教程,因为要去吃饭了~这里就使用默认的先

我们发现700步就收敛了,所以我决定改成emtol= 1000修改成500再跑```` 1888步收敛了~所以我尝试删掉emtol= 1000 但是通过前面的教程,其也是500来步就已经收敛了~

结果如图

PS:作图方法如下:

1
gmx energy -f em.edr -o potential.xvg

选择 10 0

potential

再用共轭梯度法,mdp文件参数如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
; LINES STARTING WITH ';' ARE COMMENTS
title       = Minimization  ; Title of run
define          = -DFLEXIBLE

; Parameters describing what to do, when to stop and what to save
integrator  = cg        ; Algorithm (steep = steepest descent minimization)
emtol       = 100.0     ; Stop minimization when the maximum force < 10.0 kJ/mol
emstep      = 0.01      ; Energy step size
nsteps      = 50000     ; Maximum number of (minimization) steps to perform
nstcgsteep      = 1000
energygrps  = Protein TCD   ; Which energy group(s) to write to disk

; 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
ns_type         = grid      ; Method to determine neighbor list (simple, grid)
rlist           = 1.0       ; Cut-off for making neighbor list (short range forces)
coulombtype     = PME       ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0       ; long range electrostatic cut-off
rvdw            = 1.0       ; long range Van der Waals cut-off
constraints         = none
pbc             = xyz       ; Periodic Boundary Conditions

 

1
gmx grompp -f em2.mdp -c em.gro -p topol.top -o em2.tpr

运行模拟

1
gmx mdrun -v -deffnm em2

做出来如这个图,图是否正确待考

potential2

 

12.NVT平衡(温度平衡)

限制配体:

1
gmx genrestr -f TCD.gro -o posre_TCD.itp -fc 1000 1000 1000

在选择组过程中个人觉得选择system或者TCD均可 3

增加配体限制:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
; Include ligand topology
#include "TCD.itp"


; Ligand position restraints 
#ifdef POSRES 
#include "posre_TCD.itp" 
#endif

; Include water topology
#include "amber99sb.ff/tip3p.itp"

[注:插入的为]

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

其实这里还可以分开限制(具体这里还是没有特别理解~待我看一下参考资料)

热浴

–温度耦合需要蛋白质和配体在一起?

1
gmx make_ndx -f em2.gro -o index.ndx

选择蛋白和配体分为一组

4

然后就可以热浴了~热浴的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
title       = Protein-ligand complex NVT equilibration 
define      = -DPOSRES  ; position restrain the protein and ligand
; Run parameters
integrator  = md        ; leap-frog integrator
nsteps      = 50000     ; 2 * 50000 = 100 ps
dt          = 0.002     ; 2 fs
; Output control
nstxout     = 200       ; save coordinates every 1.0 ps
nstvout     = 200       ; save velocities every 1.0 ps
nstenergy   = 200       ; save energies every 1.0 ps
nstlog      = 200       ; update log file every 1.0 ps
energygrps  = Protein TCD
; Bond parameters
continuation    = no            ; first dynamics run
constraint_algorithm = lincs    ; holonomic constraints 
constraints     = all-bonds     ; all bonds (even heavy atom-H bonds) constrained
lincs_iter      = 1             ; accuracy of LINCS
lincs_order     = 4             ; also related to accuracy
; Neighborsearching
cutoff-scheme   = Verlet
ns_type         = grid      ; search neighboring grid cells
nstlist         = 10        ; 20 fs, largely irrelevant with Verlet
rcoulomb        = 1.4       ; short-range electrostatic cutoff (in nm)
rvdw            = 1.4       ; 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
tcoupl      = V-rescale                     ; modified Berendsen thermostat
tc-grps     = Protein_TCD Water_and_ions    ; two coupling groups - more accurate
tau_t       = 0.1   0.1                     ; time constant, in ps
ref_t       = 300   300                     ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl      = no        ; no pressure coupling in NVT
; Periodic boundary conditions
pbc         = xyz       ; 3-D PBC
; Dispersion correction
DispCorr    = EnerPres  ; account for cut-off vdW scheme
; Velocity generation
gen_vel     = yes       ; assign velocities from Maxwell distribution
gen_temp    = 300       ; temperature for Maxwell distribution
gen_seed    = -1        ; generate a random seed

 

1
2
gmx grompp -f nvt.mdp -c em2.gro -p topol.top -n index.ndx -o nvt.tpr
gmx mdrun -deffnm nvt

完成后我的析出如图:

1
2
3
4
5
6
7
NOTE: The GPU has >20% more load than the CPU. This imbalance causes
      performance loss, consider using a shorter cut-off and a finer PME grid.

               Core t (s)   Wall t (s)        (%)
       Time:     3799.500      518.255      733.1
                 (ns/day)    (hour/ns)
Performance:       16.672        1.440

可以看出我的GPU太渣,暂时没有想到解决办法

我们来看一下NVT耦合的怎么样

1
gmx energy -f nvt.edr -o template.xvg

选择15 0

template

13.NPT平衡 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
title       = Protein-ligand complex NPT equilibration 
define      = -DPOSRES  ; position restrain the protein and ligand
; Run parameters
integrator  = md        ; leap-frog integrator
nsteps      = 50000     ; 2 * 50000 = 100 ps
dt          = 0.002     ; 2 fs
; Output control
nstxout     = 200       ; save coordinates every 1.0 ps
nstvout     = 200       ; save velocities every 1.0 ps
nstenergy   = 200       ; save energies every 1.0 ps
nstlog      = 200       ; update log file every 1.0 ps
energygrps  = Protein TCD
; Bond parameters
continuation    = yes           ; first dynamics run
constraint_algorithm = lincs    ; holonomic constraints 
constraints     = all-bonds     ; all bonds (even heavy atom-H bonds) constrained
lincs_iter      = 1             ; accuracy of LINCS
lincs_order     = 4             ; also related to accuracy
; Neighborsearching
cutoff-scheme   = Verlet
ns_type         = grid      ; search neighboring grid cells
nstlist         = 10        ; 20 fs, largely irrelevant with Verlet
rcoulomb        = 1.4       ; short-range electrostatic cutoff (in nm)
rvdw            = 1.4       ; 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
tcoupl      = V-rescale                     ; modified Berendsen thermostat
tc-grps     = Protein_TCD Water_and_ions    ; two coupling groups - more accurate
tau_t       = 0.1   0.1                     ; time constant, in ps
ref_t       = 300   300                     ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl      = Parrinello-Rahman             ; pressure coupling is on for 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
; Dispersion correction
DispCorr    = EnerPres  ; account for cut-off vdW scheme
; Velocity generation
gen_vel     = no        ; velocity generation off after NVT

注意,耦合组为Protein_TCD,自己的体系需要修改成自己的,下同

1
2
3
gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -n index.ndx -o npt.tpr

gmx mdrun -deffnm npt

先看一下压力

1
gmx energy -f npt.edr -o pressure.xvg

选择 16 0

pressure

再来看一下密度

1
gmx energy -f npt.edr -o density.xvg

选择 22 0

density

可以看出在1000左右,符合水的密度,个人觉得tip3p效果比spec水模型效果好

14.成品MD

一般MD需要100ns以上,但是我的电脑我怕受不了,就先50ns吧~大约需要3天-4天左右

 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
title       = Protein-ligand complex MD simulation 
; Run parameters
integrator  = md        ; leap-frog integrator
nsteps      = 25000000    ; 2 * 25000000 = 50000 ps (50 ns)
dt          = 0.002     ; 2 fs
; Output control
nstxout             = 0         ; suppress .trr output 
nstvout             = 0         ; suppress .trr output
nstenergy           = 5000      ; save energies every 10.0 ps
nstlog              = 5000      ; update log file every 10.0 ps
nstxout-compressed  = 5000      ; write .xtc trajectory every 10.0 ps
compressed-x-grps   = System
energygrps          = Protein TCD
; Bond parameters
continuation    = yes           ; first dynamics run
constraint_algorithm = lincs    ; holonomic constraints 
constraints     = all-bonds     ; all bonds (even heavy atom-H bonds) constrained
lincs_iter      = 1             ; accuracy of LINCS
lincs_order     = 4             ; also related to accuracy
; Neighborsearching
cutoff-scheme   = Verlet
ns_type         = grid      ; search neighboring grid cells
nstlist         = 10        ; 20 fs, largely irrelevant with Verlet
rcoulomb        = 1.4       ; short-range electrostatic cutoff (in nm)
rvdw            = 1.4       ; 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
tcoupl      = V-rescale                     ; modified Berendsen thermostat
tc-grps     = Protein_TCD Water_and_ions    ; two coupling groups - more accurate
tau_t       = 0.1   0.1                     ; time constant, in ps
ref_t       = 300   300                     ; reference temperature, one for each group, in K
; Pressure coupling 
pcoupl      = Parrinello-Rahman             ; pressure coupling is on for 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        ; assign velocities from Maxwell distribution

注意:nsteps = 25000000 ; 2 * 25000000 = 50000 ps (50 ns) 这里修改时长,同时注意修改耦合组成为自己的

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

经过2d12h的计算,终于完成了,即高兴又紧张。

5

由于这个贴实在太长,另外开一贴

参考文献:自动计算ESP和RESP电荷(AMBER and G09)

参考文献:小分子在GAFF立场GROMCA中的操作

参考文献:使用AmberTools+ACPYPE+Gaussian创建小分子GAFF力场的拓扑文件

参考文献:蛋白质配体复合物

Licensed under CC BY-NC-SA 4.0