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
2
3
|
...
export PATH=/home/kangsgo/miniconda3/bin:$PATH
|
对conda进行激活:
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
|
运行高斯: 高斯我放的实验室电脑跑的
拟合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等更为新的力场。
构建:
这样我们获得了amber的GAFF力场
输出可能如下:
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
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,并将前面的数字增加(配体加受体总原子数)
底部如下:
将前面得到的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
|
; 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作为输出
|