九、限制配体
gmx genrestr -f Ligand.gro -o posre_Ligand.itp -fc 1000 1000 1000
出现以下信息：
Reading structure fileSelect group to position restrainGroup0 ( System) has 68 elementsGroup1 (Other) has 68 elementsGroup2 ( MOL) has 68 elementsSelect a group:
选哪个都行，其实是一样的。我个人选择了2
得到posre_Ligand.itp
输出配体施加位置限制的itp文件
找到topol.top中的以下部分：
; Include Position restraint file#ifdef POSRES#include "posre.itp"#endif; Include ligand topology#include "Ligand.itp"; Include water topology#include "amber99sb-ildn.ff/spc.itp"
在配体拓扑的下面、水拓扑的前面插入配体限制：
; Ligand position restraints#ifdef POSRES#include "posre_19.itp"#endif
千万别放错地方了
十、体系整合
gmx make_ndx -f solv_ions.gro -o index.ndx
出现以下信息：
0 System : 125394 atoms1 Protein : 3265 atoms2 Protein-H : 1644 atoms3 C-alpha : 202 atoms4 Backbone : 606 atoms5 MainChain : 809 atoms6 MainChain+Cb : 993 atoms7 MainChain+H : 1000 atoms8 SideChain : 2265 atoms9 SideChain-H : 835 atoms10 Prot-Masses : 3265 atoms11 non-Protein : 122129 atoms12 Other: 68 atoms13 MOL : 68 atoms14 NA :3 atoms15 Water: 122058 atoms16 SOL : 122058 atoms17 non-Water : 3336 atoms18 Ion :3 atoms19 MOL : 68 atoms20 NA :3 atoms21 Water_and_ions: 122061 atomsnr : group'!': not 'name' nr name 'splitch' nr Enter: list groups'a': atom '&': and 'del' nr 'splitres' nr 'l': list residues't': atom type '|': or 'keep' nr 'splitat' nr 'h': help'r': residue 'res' nr 'chain' char"name": group 'case': case sensitive 'q': save and quit'ri': residue index>
选择你要整合的部分。在本实验中，我将蛋白和小分子绑定在一起。
> 1 | 13> q
使用的mdp文件内容如下：
title = Protein-ligand complex NVT equilibration define= -DPOSRES ; position restrain the protein and ligand; Run parametersintegrator = md ; leap-frog integratornsteps= 50000; 2 * 50000 = 100 psdt= 0.002; 2 fs; Output controlnstxout= 500 ; save coordinates every 1.0 psnstvout= 500 ; save velocities every 1.0 psnstenergy = 500 ; save energies every 1.0 psnstlog= 500 ; update log file every 1.0 psenergygrps = Protein MOL; Bond parameterscontinuation = no ; first dynamics runconstraint_algorithm = lincs ; holonomic constraints constraints= h-bonds; all bonds (even heavy atom-H bonds) constrainedlincs_iter= 1 ; accuracy of LINCSlincs_order= 4 ; also related to accuracy; Neighborsearchingcutoff-scheme = Verletns_type = grid; search neighboring grid cellsnstlist = 10 ; 20 fs, largely irrelevant with Verletrcoulomb = 1.4 ; short-range electrostatic cutoff (in nm)rvdw = 1.4 ; short-range van der Waals cutoff (in nm); Electrostaticscoulombtype= PME ; Particle Mesh Ewald for long-range electrostaticspme_order = 4 ; cubic interpolationfourierspacing = 0.16; grid spacing for FFT; Temperature couplingtcoupl= V-rescale ; modified Berendsen thermostattc-grps= Protein_MOL Water_and_ions ; two coupling groups - more accuratetau_t = 0.1 0.1 ; time constant, in psref_t = 300 300 ; reference temperature, one for each group, in K; Pressure couplingpcoupl= no ; no pressure coupling in NVT; Periodic boundary conditionspbc = xyz ; 3-D PBC; Dispersion correctionDispCorr = EnerPres ; account for cut-off vdW scheme; Velocity generationgen_vel= yes ; assign velocities from Maxwell distributiongen_temp = 300 ; temperature for Maxwell distributiongen_seed = -1 ; generate a random seed
如果有报错，可以根据提示在constraints、energygrps或tc-grps查找问题。一般报错的话就出在这几个地方。
然后我们就即将开始NVT、NPT平衡和正式的分子动力学模拟了。在此之前我想特别提醒一下，强烈建议各位使用北京超算！！！NVT和NPT平衡还好，一两个小时能做出来，但是成品MD模拟总共50万步，用自己的电脑跑真的非常非常非常慢。北京超算注册条件又简单又有200元免费额度，实际上买起来价格也不太贵，RTX3090区是4.5多一点/小时，Tesla性能高的区稍微贵一些，这个看个人需求，总之租3090和Tesla V100都是很好的选择。按照手册登录平台，选Linux系统，把GROMACS安装包和你已经做好的文件上传上去，打开web ssh就可以正常在超算平台上使用GROMACS了，贼给力！
十一、NVT平衡
gmx grompp -f GMXtut-5_nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr
如果出现
Fatal error:Cannot find position restraint file restraint.gro (option -r).From GROMACS-, you need to specify the position restraint coordinate filesexplicitly to avoid mistakes, although you can still use the same file as youspecify for the -c option.
则说明你没有用上面那行命令中的-r来指定.gro文件。在早期版本的教程没有关于这个问题的描述，官网最新版教程可以看看。
然后就可以开始运行NVT平衡了！！
gmx mdrun -deffnm nvt -v
敲下回车，开始坐等
完成时出现以下信息：
Writing final coordinates.Core t (s) Wall t (s) (%)Time: 25564.0006391.000400.01h46:31(ns/day) (hour/ns)Performance: 1.352 17.752
十二、NPT平衡
gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -r nvt.gro -p topol.top -n index.ndx -o npt.tpr
使用的mdp文件如下：
title = Protein-ligand complex NPT equilibration define = -DPOSRES ; position restrain the protein and ligand; Run parametersintegrator = md ; leap-frog integratornsteps = 50000; 2 * 50000 = 100 psdt = 0.002; 2 fs; Output controlnstenergy= 500 ; save energies every 1.0 psnstlog = 500 ; update log file every 1.0 psnstxout-compressed= 500 ; save coordinates every 1.0 ps; Bond parameterscontinuation = yes ; continuing from NVT constraint_algorithm = lincs; holonomic constraints constraints = h-bonds ; bonds to H are constrained lincs_iter = 1 ; accuracy of LINCSlincs_order = 4 ; also related to accuracy; Neighbor searching and vdWcutoff-scheme = Verletns_type = grid; search neighboring grid cellsnstlist = 20 ; largely irrelevant with Verletrlist = 1.2vdwtype = cutoffvdw-modifier = force-switchrvdw-switch = 1.0rvdw= 1.2 ; short-range van der Waals cutoff (in nm); Electrostaticscoulombtype = PME ; Particle Mesh Ewald for long-range electrostaticsrcoulomb= 1.2pme_order= 4 ; cubic interpolationfourierspacing= 0.16; grid spacing for FFT; Temperature couplingtcoupl = V-rescale ; modified Berendsen thermostattc-grps = Protein_MOL Water_and_ions ; two coupling groups - more accuratetau_t = 0.1 0.1 ; time constant, in psref_t = 300 300 ; reference temperature, one for each group, in K; Pressure couplingpcoupl = Berendsen ; pressure coupling is on for NPTpcoupltype = isotropic ; uniform scaling of box vectorstau_p = 2.0 ; time constant, in psref_p = 1.0 ; reference pressure, in barcompressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1refcoord_scaling = com; Periodic boundary conditionspbc = xyz ; 3-D PBC; Dispersion correction is not used for proteins with the C36 additive FFDispCorr= no ; Velocity generationgen_vel = no ; velocity generation off after NVT
然后就可以运行
gmx mdrun -deffnm npt -v
来进行最后的NPT平衡了！
（坐等again）
Writing final coordinates.step 50000, remaining wall clock time:0 sCore t (s) Wall t (s) (%)Time: 10568.0002642.000400.044:02(ns/day) (hour/ns)Performance: 3.270 7.339
得到一堆npt开头的文件
至此，两步平衡结束，准备开始最终模拟！
十三、成品MD
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md_0_1.tpr
得到md_0_1.tpr
gmx mdrun -deffnm md_0_1 -v
-v可以让你看到现在运行到哪一步
然后就慢慢等呗
