新手自学：如何用gromacs对简单分子复合物进行伞形采样

1、建立体系：

1、将蛋白的pdb文件转化为gmx：

gmx pdb2gmx -f 2BEG_model1_capped.pdb -ignh -ter -o complex.gro

这个网页可以实现将多肽序列转化为pdb：

ProBuilder On-line

这个教程的蛋白2BFG包含两条链（chain A和B）

在生成的topol文件中，增加如下的内容，效果就是chain B将被添加限制力而不能移动：

#ifdef POSRES_B
#include "posre_Protein_chain_B.itp"
#endif

这里chain A做为拉动链，而chain B作为固定链，用上面的参数进行位置限制。 

这里的ifdef是一个函数定义，如果在mdp文件中，有相关的字段(define=-DPOSRES_B)，模拟就是读取这个itp拓扑文件里的数据。

1.1、将小分子-小分子复合物建立gro和top文件。

（这是教程外的，我个人的实践）多肽类小分子可以通过gmx自身完成gro和top文件的生成，文件如下：一个posre、topol和gro文件。

若是pdb中有多个分子，还会额外生成一个itp文件，并在topol中用incluede命令进行链接导入。这里pdb中只有一个小分子短肽，故其拓扑数据直接记录在topol文件（框1）中。

对于非蛋白，非核酸的分子，可以用sob老师的软件sobtop完成 gro和top生成。

这里，用sob老师生成的小分子itp数据如下图，我分为框内和框外（剩余）部分。这个文件不能直接用# include “xx.itp"的方式直接写入到topol文件中。

 在topol文件中，其逻辑是先导入atomtype、bond等信息，证据就是其在开始的注释后，是先导入了gmx的share文件夹里top的forcefield.itp文件（下图框）。打开这个forcefield.itp文件，可以看到调用的itp文件，会先读取atomtype区块。

因此，若直接将小分子的itp以# include “xx.itp"命令写入，会再读一遍atomtype（由sobtop生成的），就会出错。

 我们可以把小分子（例如CH3F）的atomtype区块剪切出来，独立作为一个itp文件，例如命名为CH3Fatomtype.itp，放在当前工作路径下，在topol读取moleculetype区块之前，写入#include ”CH3Fatomtype.itp",就不会出错。

1.2、给体系填加盒子：

进行拉动模拟（pull simulation）不像其他模拟，

拉力距离必须始终小于箱向量长度的一半，箱向量沿着箱向量进行拉力。原因是：GROMACS 计算距离，同时考虑周期性。这样，如果你有一个10纳米的盒子，你拉过一个大于5.0纳米的距离，周期性的距离变成了拉的参考距离，这个距离实际上小于5.0纳米！

这个教程里，打算在一个12nm长的箱子里对分子进行5nm的拉动模拟。盒子设置命令如下：

gmx editconf -f complex.gro -o newbox.gro -center 3.280 2.181 2.4775 -box 6.560 4.362 12

center选项将复合物放在3.28*2.18*2.48这个坐标，box选项设置盒子的三维是6.56*4.36*12.0（教程里的情况）。

实际情况还是根据体系的大小进行考虑，盒子太大对计算也是一种负荷。

在VMD里导入这个建立的newbox.gro文件可视化盒子的情况，打开tk consle窗口输入如下命令：

pbc box #这里是pbc 表示周期性元胞边界the periodic unit cell boundaries

1.2.1、借用VMD软件进行更便捷的操作： 

其实，以上的gmx editconf命令实际效果只是在gro文件下增加了如下的周期性边界信息（下图框内），这个命令不会对topol文件造成修改（命令完成后，这个文件的修改日期没有改变）。 

利用VMD的效果也是相同的，我在gro文件中对盒子的大小进行调整（对于小分子，教程里设置的盒子太大了),调整后的效果如下图所示。另外，盒子改变后，会出现复合物可能跑到盒子外去的情况，可以用VMD的快捷键8（移动分子）和R（旋转）进行调整复合物的位置：

 1.3、添加水分子：

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

以上的命令会根据盒子的大小填充水分子，效果用VMD打开如下：

 1.4、建立ions.mdp文件，添加离子：

ions.mdp文件如下：

; ions.mdp - used as input into grompp to generate ions.tpr
; 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 < 1000.0 kJ/mol/nm
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    ; Buffered neighbor searching
ns_type         = grid      ; Method to determine neighbor list (simple, grid)
rlist           = 1.4       ; Cut-off for making neighbor list (short range forces)
coulombtype     = cutoff    ; Treatment of long range electrostatic interactions
rcoulomb        = 1.4       ; Short-range electrostatic cut-off
rvdw            = 1.4       ; Short-range Van der Waals cut-off
pbc             = xyz       ; Periodic Boundary Conditions

 利用grompp生成这个mdp文件的二进制文件tpr：

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

添加NaCl：

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

在跳出的选项中选择SOL对应的序号，就会将水分子换为离子。背后的逻辑就是：离子不是直接加上去的，而是会随机夺取体系中水分子的位置。

可以查看添加离子前后在topol里水分子的数量和离子的数量。

1.5、之后是常规的能量最小化，NPT平衡：

教程的minim.mdp文件和上面的ions.mdp是一样的，生成tpr然后用mdrun跑一遍能量最小化：

gmx grompp -f minim.mdp -c solv_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em

教程里的npt模拟的npt.mdp文件：

title       = NPT Equilibration 
define      = -DPOSRES          ; position restrain the protein
; Run parameters
integrator  = md                ; leap-frog integrator
nsteps      = 50000             ; 2 * 50000 = 100 ps
dt          = 0.002             ; 2 fs
; Output control
nstxout     = 1000              ; save coordinates every 2 ps
nstvout     = 1000              ; save velocities every 2 ps
nstenergy   = 1000              ; save energies every 2 ps
nstlog      = 1000              ; update log file every 2 ps
; Bond parameters
continuation         = no        ; Initial simulation 
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
ns_type     = grid              ; search neighboring grid cels
nstlist     = 5                 ; 10 fs
rlist       = 1.4               ; short-range neighborlist cutoff (in nm)
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 is on
tcoupl      = Berendsen             ; Weak coupling for initial equilibration 
tc-grps     = Protein   Non-Protein ; two coupling groups - more accurate
tau_t       = 0.1       0.1         ; time constant, in ps
ref_t       = 310       310         ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl              = Berendsen     ; Pressure coupling on in NPT, also weak coupling
pcoupltype          = isotropic     ; uniform scaling of x-y-z box vectors
tau_p               = 2.0