Back

GROMACS下蛋白配体复合物模拟(GAFF力场)

1. 环境介绍

终端: Windows Terminal

系统: Windows下wls2 (Ubuntu-20.04)

配置: R5-4650U 16G 512SSD

2. Conda安装

我安装的Miniconda

下载地址: https://docs.conda.io/en/latest/miniconda.html

1
2
3
4
5
6
7
8
9
#进入的用户主目录
cd /home/kangsgo
#下载安装包
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
#赋权限
chmod u+x Miniconda3-latest-Linux-x86_64.sh
#执行
bash Miniconda3-latest-Linux-x86_64.sh
#后续全部默认的安装

设置环境变量:

1
vi .bashrc

新增导出路径: 需要根据自己用户名修改:

1
2
3
...

export PATH=/home/kangsgo/miniconda3/bin:$PATH

对conda进行激活:

1
conda init

3. Ambertools安装

1
2
3
4
5
6
7
8
9
#新增加环境变量:
conda create --name ambermd
#进入环境
conda activate ambermd
#安装
conda install -c conda-forge ambertools=21 compilers
conda update -c conda-forge ambertools
#openbabel安装
conda install openbabel 

4. 小分子准备

小分子之前是用autodock对接的,因为对接结果没有氢键,可以用UCSF Chimera进行加氢或者使用openbabel,因为我Chimera一直没有太安装成功,所以用的openbanel

首先用pymol把pdbt的格式改为pdb格式,构建如下简单的python脚本:

python 版本3.8.10

openbabel 版本 3.3.1 ?

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
#参考1: http://openbabel.org/dev-api/classOpenBabel_1_1OBMol.shtml#ad3bab31bf64ff6cb18f6ee259b3b8c11
#参考2:http://openbabel.org/docs/current/UseTheLibrary/PythonDoc.html
from openbabel import openbabel

obConversion = openbabel.OBConversion
#设置输入输出格式
obConversion.SetInAndOutFormats("pdb", "pdb")

mol = openbabel.OBMol()
obConversion.ReadFile(mol,"lig.pdb")

#加氢
mol.AddHydrogens()

#保存文件
obConversion.WriteFile(mol,"addlig.pdb")

用高斯打开addlig.pdb,另存为lig.pdb (做这一步的原因是我用高斯打开显示有问题,重新保存之后再打开就没有问题了。) 注意加氢的正确率不高,可能需要手动用高斯或者pymol进行修正

(PS:一直align没有成功,若可以用pymol的align把原分子的格式弄过去就更加棒了,目前有一个想法是可以用CONECT和上面的原子序号进行修正,之后替换坐标,不知道这种是否可行?)

小分子我用的高斯准备,clean 一下之后另存为pdb格式或者mol2格式(PS:我保存的mol2格式一直没跑通我同学的分子····)

准备的小分子为lig.pdb 输出的文件为lig.gjf

1
antechamber -i lig.pdb -fi pdb -o lig.gjf -fo gcrt -pf y -gm "%mem=10GB" -gn "%nproc=30" -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 lig.gesp -gv 1

运行高斯: 高斯我放的实验室电脑跑的

1
g09 lig.gjf

拟合resp,输入文件为lig.gesp ,输出为 mol2格式,配体名称被命名为LIG,可以自己更改

1
antechamber -i lig.gesp -fi gesp -o lig_a.mol2 -fo mol2 -pf y -c resp -rn LIG

生成确实参数文件

1
parmchk2 -i lig_a.mol2 -f mol2 -o lig.frcmod

构建leap.in文件:

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

PS: 这里我卡壳了一会儿,因为之前是不需要oldff文件夹的,提示找不到文件,原因是新版本把过时的力场放到了oldff文件夹里面了。终端提示了

1
2
3
4
5
-I: Adding /home/kangsgo/miniconda3/envs/ambermd/dat/leap/prep to search path.
-I: Adding /home/kangsgo/miniconda3/envs/ambermd/dat/leap/lib to search path.
-I: Adding /home/kangsgo/miniconda3/envs/ambermd/dat/leap/parm to search path.
-I: Adding /home/kangsgo/miniconda3/envs/ambermd/dat/leap/cmd to search path.
-f: Source leap.in.

会在这些文件夹里面搜索,所以找了一下

1
find ./miniconda3/ -name "leaprc.*"

若有兴趣可以使用ff14SB等更为新的力场。

构建:

1
tleap -f leap.in

这样我们获得了amber的GAFF力场

1
ls

输出可能如下:

image-20211116163719009

5. 转换为gromacs可读力场文件

安装acpype:

1
conda install -c conda-forge acpype

转换:

1
acpype -p lig.prmtop -x lig.inpcrd -d

我将前面得到的所有文件,都移动到了LIG的子文件内。方便下一步。

6. 安装GROMACS

因为我最后跑是在服务器上跑的,所以在自己电脑上安装一个cpu版本即可

1
wget https://ftp.gromacs.org/gromacs/gromacs-2021.4.tar.gz

执行:

1
2
3
4
5
6
7
8
9
tar xfz gromacs-2021.4.tar.gz
cd gromacs-2021.4
mkdir build
cd build
cmake .. -DGMX_BUILD_OWN_FFTW=ON -DREGRESSIONTEST_DOWNLOAD=ON
make
make check
sudo make install
source /usr/local/bin/GMXRC

7. 处理复合物

1
cat complex.pdb

我的文件结果如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
ATOM      1  C1  RES d   1       1.569   4.295  -1.522  0.00  0.00           C
ATOM      2  C2  RES d   1       0.511   3.793  -0.688  0.00  0.00           C
ATOM      3  C3  RES d   1      -0.842   4.328  -0.878  0.00  0.00           C
ATOM      4  C4  RES d   1      -1.045   5.270  -1.813  0.00  0.00           C
ATOM      5  C6  RES d   1       2.877   3.761  -1.326  0.00  0.00           C
...
ATOM    733  H8   DA A  25      -0.046   5.444   2.135  1.00  0.59           H
ATOM    734  H1'  DA A  25      -3.003   3.634   2.817  1.00  0.76           H
ATOM    735  H2'  DA A  25      -2.883   4.867   5.348  1.00  0.73           H
ATOM    736  H3'  DA A  25      -2.556   7.093   4.556  1.00  0.75           H
ATOM    737  H4'  DA A  25      -4.418   6.552   2.265  1.00  0.70           H
ATOM    738  H5'  DA A  25      -1.811   8.069   2.413  1.00  0.89           H
ATOM    739  H61  DA A  25       3.164   1.975   6.139  1.00  0.70           H
ATOM    740  H62  DA A  25       3.280   3.119   4.820  1.00  0.65           H
ATOM    741 H2''  DA A  25      -4.366   4.653   4.430  1.00  0.98           H
ATOM    742 H5''  DA A  25      -2.785   7.842   0.945  1.00  0.79           H
TER
END

可以看到小分子配体为RES标签

获取受体的结构

1
2
grep -v RES complex.pdb > repoter.pdb
cat repoter.pdb

结果如下:

1
2
3
4
5
6
7
ATOM     27  C5'  DT A   4       7.122  -9.923 -11.819  1.00  0.65           C
ATOM     28  O5'  DT A   4       6.723  -8.905 -12.741  1.00  0.75           O
ATOM     29  C4'  DT A   4       8.597  -9.841 -11.505  1.00  0.59           C
ATOM     30  O4'  DT A   4       8.854 -10.536 -10.257  1.00  0.66           O
ATOM     31  C3'  DT A   4       9.514 -10.491 -12.546  1.00  0.69           C
ATOM     32  O3'  DT A   4      10.717  -9.728 -12.693  1.00  0.69           O
....

7.5 增加力场

核酸力场amber现在建议的ff-nucleic-OL15, 但是找的参考文件使用的ff12SB力场,故我也是用的ff12SB力场,可以在此下载:

1
2
https://www.gromacs.org/Downloads/User_contributions/Force_fields
#wget https://www.gromacs.org/@api/deki/files/253/=FF12SB_FF14SB_GROMACS_TEST.tar.bz2

解压拷贝ff12SB文件到您实验目录即可。

写此文章的最新的力场gromacs格式可以去如下网站下载:

1
https://fch.upol.cz/ff_ol/downloads.php

当您很多年后查看的时候,请前往:

1
http://ambermd.org/AmberModels.php

查看amber现在建议的力场。(2021/11/17)

1
2
3
(base) kangsgo@FQI5F177Y1UC7IR:~/zjh/MD$ cp -r /mnt/d/shiyan/zjh/amber12sb.ff/ amber12sb.ff
(base) kangsgo@FQI5F177Y1UC7IR:~/zjh/MD$ ls
 amber12sb.ff  bak  complex.pdb   repoter.pdb

8. 构建复合物拓扑结构

首先构建受体的拓扑结构,我这里是核酸:

1
gmx pdb2gmx -f repoter.pdb -o repoter_processed.gro -water tip3p -noignh

其中-noignh 是忽略氢原子,并进行重建,因为这个结构是autodock得到的,一般结构都不是标准结构,所以进行了忽略,让GROMACS进行重构标准的H。

将前面的LIG_GMX.gro文件的坐标拷贝入repoter_processed.gro,并将前面的数字增加(配体加受体总原子数)

image-20211117110802715

底部如下:

image-20211117110823453

将前面得到的LIG_GMX.tio文件最底下的[system][molecules]删除,修改名字为LIG.itp 放入开始生成的top.top文件内:位置位于核酸和水的位置限制之间

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

; Include ligand topology
#include "LIG.itp"

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

在top的底部再增加一个分子:

1
2
3
4
[ molecules ]
; Compound        #mols
DNA_chain_A         1
LIG					    1

将前面整合的repoter_processed.gro文件重命名为complex.gro方便理解:

1
mv repoter_processed.gro complex.gro

将顶部的[ defaults ]删除,将[ atomtypes ]内容剪切,放入topol.top文件中:

位于:

1
2
; Include forcefield parameters
#include "./amber12sb.ff/forcefield.itp"

1
2
3
4
5
6
7
[ moleculetype ]
; Name            nrexcl
DNA_chain_A         3

[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  typeB    chargeB      massB
; residue   4 DT  rtp DT5  q -0.3

之间,不可放错,最后得到的结果如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
; Include forcefield parameters
#include "./amber12sb.ff/forcefield.itp"

[ atomtypes ]
;name   bond_type     mass     charge   ptype   sigma         epsilon       Amb
 ca       ca          0.00000  0.00000   A     3.39967e-01   3.59824e-01 ; 1.91  0.0860
....

[ moleculetype ]
; Name            nrexcl
DNA_chain_A         3

9.构建盒子并填充水

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

创建一个最小的mdp文件用于平衡离子ions.mdp

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

; 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

; 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	    = cutoff	; Treatment of long range electrostatic interactions
rcoulomb	    = 1.0		; long range electrostatic cut-off
rvdw		    = 1.0		; long range Van der Waals cut-off
pbc             = xyz 		; Periodic Boundary Conditions

生成tpr文件

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

用KCl进行中和离子

1
2
 gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname K -nname CL -neutral
 #我选择的是替换SQL原子

输出如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
Replacing solvent molecule 685 (atom 2818) with K
Replacing solvent molecule 2010 (atom 6793) with K
Replacing solvent molecule 5260 (atom 16543) with K
Replacing solvent molecule 5683 (atom 17812) with K
Replacing solvent molecule 1930 (atom 6553) with K
Replacing solvent molecule 1454 (atom 5125) with K
Replacing solvent molecule 644 (atom 2695) with K
Replacing solvent molecule 151 (atom 1216) with K
Replacing solvent molecule 1416 (atom 5011) with K
Replacing solvent molecule 580 (atom 2503) with K
Replacing solvent molecule 6002 (atom 18769) with K
Replacing solvent molecule 6244 (atom 19495) with K
Replacing solvent molecule 6250 (atom 19513) with K
Replacing solvent molecule 5553 (atom 17422) with K
Replacing solvent molecule 2379 (atom 7900) with K
Replacing solvent molecule 5380 (atom 16903) with K
Replacing solvent molecule 5221 (atom 16426) with K
Replacing solvent molecule 5674 (atom 17785) with K
Replacing solvent molecule 1641 (atom 5686) with K
Replacing solvent molecule 4133 (atom 13162) with K
Replacing solvent molecule 1839 (atom 6280) with K

10. 能量最小化

我们用梯度最小化,请注意energygrps组需要改成自己的,因为我做的DNA和LIG配体分子,所以我这里帮他们作为一组,在新的教程里面也有不做限制组的。

构建mini.mdp文件,内容如下:

 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  = DNA LIG   ; 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
2
3
gmx grompp -f mini.mdp -c solv_ions.gro -p topol.top -o em.tpr
#运行模拟
gmx mdrun -v -deffnm em

我这里399步达成了,若体系比较小,其实可以能量最小化Stop minimization更小的设置

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
Step=  396, Dmax= 5.2e-03 nm, Epot= -3.39672e+05 Fmax= 4.34592e+03, atom= 504
Step=  398, Dmax= 3.1e-03 nm, Epot= -3.39769e+05 Fmax= 3.77077e+02, atom= 509

writing lowest energy coordinates.

Steepest Descents converged to Fmax < 500 in 399 steps
Potential Energy  = -3.3976897e+05
Maximum force     =  3.7707678e+02 on atom 509
Norm of force     =  3.3953958e+01

GROMACS reminds you: "It's Not Your Fault" (Pulp Fiction)

11. NVT与NPT

即热平衡和压力平衡,为了防止在NVT和NPT的时候蛋白复合物爆开,所以要加一个力的限制,首先构建配体的力限制

将前面得到的LIG_GMX.gro重命名为LIG.gro保存在MD文件夹内:

1
2
3
4
#创建索引文件
gmx make_ndx -f LIG.gro -o index_lig.ndx
> 0 & ! a H*
> q

可以 tail index_lig.ndx 查看一下结果:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
  31   32   33   34   35   36   37   38   39   40   41   42   43   44   45
  46   47
[ LIG ]
   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
[ System_&_!H* ]
   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

给增加1000的三个方向的力

1
2
3
gmx genrestr -f LIG.gro -n index_lig.ndx -o posre_LIG.itp -fc 1000 1000 1000
#选择刚刚创建的ndx组
> 3 

topol.top文件的

1
2
3
4
5
6
7
; Include ligand topology
#include "LIG.itp"

;>这里插入

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

之间插入位置限制

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

; Include ligand Position restraint file
#ifdef POSRES_LIG
#include "posre_LIG.itp"
#endif

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

为了防止体系爆炸,我们需要把蛋白和配体设置为一组,再次需要用到ndx:

1
2
3
4
gmx make_ndx -f em.gro -o index.ndx
#由于我是核酸,没有蛋白那么多组,若是蛋白配体复合物的时候需要注意
> 1 | 2
> q

实际操作如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
  0 System              : 20812 atoms
  1 DNA                 :   716 atoms
  2 LIG                 :    47 atoms
  3 K                   :    21 atoms
  4 Other               :    47 atoms
  5 LIG                 :    47 atoms
  6 K                   :    21 atoms
  7 Water               : 20028 atoms
  8 SOL                 : 20028 atoms
  9 non-Water           :   784 atoms
 10 Ion                 :    21 atoms
 11 LIG                 :    47 atoms
 12 K                   :    21 atoms
 13 Water_and_ions      : 20049 atoms

> 1 | 2

Copied index group 1 'DNA'
Copied index group 2 'LIG'
Merged two groups with OR: 716 47 -> 763

 14 DNA_LIG             :   763 atoms

> q

构建nvt的mdp : nvt.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
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
nstenergy               = 500   ; save energies every 1.0 ps
nstlog                  = 500   ; update log file every 1.0 ps
nstxout-compressed      = 500   ; save coordinates every 1.0 ps
; Bond parameters
continuation            = no        ; first dynamics run
constraint_algorithm    = lincs     ; holonomic constraints 
constraints             = h-bonds   ; bonds to H are constrained 
lincs_iter              = 1         ; accuracy of LINCS
lincs_order             = 4         ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme           = Verlet
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 20        ; largely irrelevant with Verlet
rlist                   = 1.2
vdwtype                 = cutoff
vdw-modifier            = force-switch
rvdw-switch             = 1.0
rvdw                    = 1.2       ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
rcoulomb                = 1.2       ; short-range electrostatic cutoff (in nm)
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.16      ; grid spacing for FFT
; Temperature coupling
tcoupl                  = V-rescale                     ; modified Berendsen thermostat
tc-grps                 = DNA_LIG 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 is not used for proteins with the C36 additive FF
DispCorr                = no 
; Velocity generation
gen_vel                 = yes       ; assign velocities from Maxwell distribution
gen_temp                = 300       ; temperature for Maxwell distribution
gen_seed                = -1        ; generate a random seed

注意tc-grps = DNA_LIG Water_and_ions的热浴组合温度需要改成自己的,这个体系用的常温300K,跑的时间为100ps

执行:

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

NPT文件如下[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
47
48
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
nstenergy               = 500       ; save energies every 1.0 ps
nstlog                  = 500       ; update log file every 1.0 ps
nstxout-compressed      = 500       ; save coordinates every 1.0 ps
; Bond parameters
continuation            = yes       ; continuing from NVT 
constraint_algorithm    = lincs     ; holonomic constraints 
constraints             = h-bonds   ; bonds to H are constrained 
lincs_iter              = 1         ; accuracy of LINCS
lincs_order             = 4         ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme           = Verlet
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 20        ; largely irrelevant with Verlet
rlist                   = 1.2
vdwtype                 = cutoff
vdw-modifier            = force-switch
rvdw-switch             = 1.0
rvdw                    = 1.2       ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
rcoulomb                = 1.2
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.16      ; grid spacing for FFT
; Temperature coupling
tcoupl                  = V-rescale                     ; modified Berendsen thermostat
tc-grps                 = DNA_LIG 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                  = Berendsen                     ; 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 is not used for proteins with the C36 additive FF
DispCorr                = no 
; Velocity generation
gen_vel                 = no        ; velocity generation off after NVT 

执行:

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

gmx mdrun -v -deffnm npt

12. 成品MD

md.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 MD simulation 
; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 400000000   ; 2 * 400000000 = 800000 ps (800 ns)
dt                      = 0.002     ; 2 fs
; Output control
nstenergy               = 50000      ; save energies every 100.0 ps
nstlog                  = 50000      ; update log file every 100.0 ps
nstxout-compressed      = 50000      ; save coordinates every 100.0 ps
; Bond parameters
continuation            = yes       ; continuing from NPT 
constraint_algorithm    = lincs     ; holonomic constraints 
constraints             = h-bonds   ; bonds to H are constrained
lincs_iter              = 1         ; accuracy of LINCS
lincs_order             = 4         ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme           = Verlet
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 20        ; largely irrelevant with Verlet
rlist                   = 1.2
vdwtype                 = cutoff
vdw-modifier            = force-switch
rvdw-switch             = 1.0
rvdw                    = 1.2       ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
rcoulomb                = 1.2
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.16      ; grid spacing for FFT
; Temperature coupling
tcoupl                  = V-rescale                     ; modified Berendsen thermostat
tc-grps                 = DNA_LIG 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 is not used for proteins with the C36 additive FF
DispCorr                = no 
; Velocity generation
gen_vel                 = no        ; continuing from NPT equilibration 

成品模拟

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

gmx mdrun -v -deffnm md_0_800

13. 轨迹处理

轨迹处理参考:

http://www.mdtutorials.com/gmx/complex/09_analysis.html

1
2
3
4
5
6
7
8
#中心并压缩
gmx trjconv -s md_0_10.tpr -f md_0_10.xtc -o md_0_10_center.xtc -center -pbc mol -ur compact
#选择protein/nucleic 作为center,system作为输出
#释放第一帧结果
gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o start.pdb -dump 0
#释放fit结果
gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o md_0_10_fit.xtc -fit rot+trans
#选择Backbone作为比对,System作为输出
Licensed under CC BY-NC-SA 4.0