新手自学：如何用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 12center选项将复合物放在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