Gromacs结合自由能计算教程（三）
以下文章来源于网络，如有侵权,请联系管理员删除!
一、准备拓扑文件
对于一个有一定经验的分子模拟工作者而言，应该知道如何准备拓扑文件。如果你恰好要处理的是GAFF力场下小分子的模拟，那么你可以参考这篇文章来学习如何准备拓扑文件。如果你在mdp文件中定义了耦合基团为“ligand”（couple-moltype = ligand，如本例），那么在top文件中你就必须单独建立一个ligand分子描述模块。在本例中我在复合物的top文件中用了#include "ligand.itp"来包含了单独描述ligand的itp文件。top和itp文件的格式可以参照本例的complex.top ligand.itp。在这里还定义了用于施加position restrain的posre_ligand.itp文件，其内容是这样的
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
2 1 1000 1000 1000
3 1 1000 1000 1000
4 1 1000 1000 1000
5 1 1000 1000 1000
6 1 1000 1000 1000
7 1 1000 1000 1000
8 1 1000 1000 1000
9 1 1000 1000 1000
10 1 1000 1000 1000
11 1 1000 1000 1000
12 1 1000 1000 1000
13 1 1000 1000 1000
14 1 1000 1000 1000
15 1 1000 1000 1000
16 1 1000 1000 1000
17 1 1000 1000 1000
18 1 1000 1000 1000
19 1 1000 1000 1000
20 1 1000 1000 1000
21 1 1000 1000 1000
22 1 1000 1000 1000
23 1 1000 1000 1000
24 1 1000 1000 1000
25 1 1000 1000 1000
26 1 1000 1000 1000
27 1 1000 1000 1000
其中第一个数字指的是施加限制势的原子编号，1-27就是将整个分子（本例-樟脑）的所有原子都施加了限制，第二个数字是限制势函数类型，之后3个数字是限制势的力常数，1000已经够用。
需要注意的是限制势的信息一定要添加在施加限制的分子的随后，不一定非得用include一个分开的文件，也可以直接加在分子描述的后面。
二、准备MDP文件
和其它的分子模拟操作类似，自由能计算也包含了能量最小化、NVT平衡、NPT平衡、生产相模拟。本例中所有的mdp文件都可以在压缩包的MDP文件夹内找到，注释都很详尽了，在此不作赘述，其中举一个生产相模拟的MDP文件为例：
;====================================================
; Production simulation
;====================================================
;----------------------------------------------------
; RUN CONTROL
;----------------------------------------------------
integrator = sd ; stochastic leap-frog integrator
nsteps = 2000000 ; 2 * 2,000,000 fs = 4000 ps = 4ns
dt = 0.002 ; 2 fs
comm-mode = Linear ; remove center of mass translation
nstcomm = 100 ; frequency for center of mass motion removal
;----------------------------------------------------
; OUTPUT CONTROL
;----------------------------------------------------
nstxout = 0 ; don't save coordinates to .trr
nstvout = 0 ; don't save velocities to .trr
nstfout = 0 ; don't save forces to .trr
nstxout-compressed = 1000 ; xtc compressed trajectory output every 1000 steps (2 ps)
compressed-x-precision = 1000 ; precision with which to write to the compressed trajectory file
nstlog = 1000 ; update log file every 2 ps
nstenergy = 1000 ; save energies every 2 ps
nstcalcenergy = 100 ; calculate energies every 100 steps
;----------------------------------------------------
; BONDS
;----------------------------------------------------
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; hydrogens only are constrained
lincs_iter = 1 ; accuracy of LINCS (1 is default)
lincs_order = 4 ; also related to accuracy (4 is default)
lincs-warnangle = 30 ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
continuation = yes
;----------------------------------------------------
; NEIGHBOR SEARCHING
;----------------------------------------------------
cutoff-scheme = Verlet
ns-type = grid ; search neighboring grid cells
nstlist = 20 ; 20 fs (default is 10)
rlist = 1.0 ; short-range neighborlist cutoff (in nm)
pbc = xyz ; 3D PBC
;----------------------------------------------------
; ELECTROSTATICS
;----------------------------------------------------
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
ewald_geometry = 3d ; Ewald sum is performed in all three dimensions
pme-order = 6 ; interpolation order for PME (default is 4)
fourierspacing = 0.10 ; grid spacing for FFT
ewald-rtol = 1e-6 ; relative strength of the Ewald-shifted direct potential at rcoulomb
;----------------------------------------------------
; VDW
;----------------------------------------------------
vdw-type = PME
rvdw = 1.0
vdw-modifier = Potential-Shift
ewald-rtol-lj = 1e-3
lj-pme-comb-rule = Geometric
DispCorr = EnerPres
;----------------------------------------------------
; TEMPERATURE & PRESSURE COUPL
;----------------------------------------------------
tc_grps = System
tau_t = 1.0
ref_t = 300
pcoupl = Parrinello-Rahman
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2 ; time constant (ps)
ref_p = 1.0 ; reference pressure (bar)
compressibility = 4.5e-05 ; isothermal compressibility of water (bar^-1)
;----------------------------------------------------
; VELOCITY GENERATION
;----------------------------------------------------
gen_vel = no ; Velocity generation is off (if gen_vel is 'yes', continuation should be 'no')
gen_seed = -1 ; Use random seed
gen_temp = 300
;----------------------------------------------------
; FREE ENERGY CALCULATIONS
;----------------------------------------------------
free-energy = yes
couple-moltype = ligand
couple-lambda0 = none
couple-lambda1 = vdw-q
couple-intramol = no
separate-dhdl-file = yes
sc-alpha = 0.5
sc-power = 1
sc-sigma = 0.3
init-lambda-state = 0
coul-lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.00 0.00 0.00 0.0 0.2 0.4 0.5 0.6 0.7 0.8 0.9 1.0
vdw-lambdas = 0.0 0.2 0.4 0.5 0.6 0.7 0.8 0.85 0.9 0.95 0.97 0.99 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
nstdhdl = 100
calc-lambda-neighbors = -1
关于自由能计算的参数信息含义在上一节中已经讲得很清楚了。这里模拟时间设置为4ns，对于自由能计算，一般来说时间越长，采样越充分，所得模拟结果越可靠，当然这也需要根据你的体系大小进行权衡，但2ns大概是最低标准了。事实上往往还要进行好几次独立的模拟，数据才有说服力。在静电计算上采用了较高精度的PME格点，这一块的计算是模拟中最耗时的一部分，但为了自由能计算的准确度还是不能省的。
你可能需要根据自己的体系调整mdp文件，你可以在mdp文件夹内直接修改那些mdp文件，其中保留住"init-lambda-state = $LAMBDA$"关键字供脚本识别修改，然后运行脚本“bash mk_mdp.sh coul”这里指的是按照coul的lambda向量来生成对应mdp文件。如此便会自动地在ENMIN，NVT, NPT, PROD文件夹下产生对应的mdp文件
for mdp in *.mdp
do
Nlambdas=`cat $mdp | grep $1-lambdas | wc -w`
Nlambdas=`expr $Nlambdas - 2`
filename=$(basename "$mdp")
filename="${filename%.*}"
newdirname=`printf $filename | tr '[:lower:]' '[:upper:]'`
newdir=`printf "%s" $newdirname`
echo "Making directory $newdirname"
mkdir -p ${PWD}/$newdirname
i=0
while [ $i -lt $Nlambdas ]
do
new_filename=`printf "%s.%02d.mdp" $filename $i`
echo "Writing file $new_filename"
sed "s/\\\$LAMBDA\\\$/${i}/" $mdp > $newdirname/$new_filename
i=`expr $i + 1`
done
done
三、运行模拟
在根目录下执行run.sh脚本
for d in lambda.*/; do
d1=$(basename $d)
lam="${d1##*.}"
cd $d
mkdir ENMIN
cd ENMIN
gmx grompp -f ../../MDP/ENMIN/enmin.$lam.mdp -c ../../ligand.gro -p ../../ligand.top -o enmin.tpr
gmx mdrun -v -stepout 1000 -s enmin.tpr -deffnm enmin
cd ../
mkdir NVT
cd NVT
gmx grompp -f ../../MDP/NVT/nvt.$lam.mdp -c ../ENMIN/enmin.gro -p ../../ligand.top -o nvt.tpr -r ../../ligand.gro
gmx mdrun -stepout 1000 -s nvt.tpr -deffnm nvt
cd ../
mkdir NPT
cd NPT
gmx grompp -f ../../MDP/NPT/npt.$lam.mdp -c ../NVT/nvt.gro -t ../NVT/nvt.cpt -p ../../ligand.top -o npt.tpr -r ../../ligand.gro
gmx mdrun -stepout 1000 -s npt.tpr -deffnm npt
cd ../
mkdir PROD
cd PROD
gmx grompp -f ../../MDP/PROD/prod.$lam.mdp -c ../NPT/npt.gro -t ../NPT/npt.cpt -p ../../ligand.top -o prod.tpr
gmx mdrun -stepout 1000 -s prod.tpr -deffnm prod -dhdl dhdl
cd ../../
done
这个脚本会打开MDP文件夹内的mdp文件，然后根据根目录内lambda.XX文件夹，自动运行从lambda0到lambda20（本例）的副本并保存在lambda.XX文件夹中。如果要中断模拟然后续算，那么就可以将已经算完的模拟文件夹置于其它地方，这样脚本就会从当前目录的最小序号模拟开始计算。当然，你需要根据自己的模拟调整参数；例如本例中的gro文件为ligand.gro，你需要按照自己的实际情况调整。
需要注意的是这里有一项“-r ligand.gro”，这是在Gromacs 2018中新加入的项目，用于给施加position restrain时作为参考坐标，参考坐标文件可以是原文件。
取决于你的计算体系大小和算力，一套模拟下来可能需要一天以上的时间。待全部模拟完成后，在根目录下执行cp_xvg.sh脚本
#!/bin/bash
if [ ! -d "dHdl_files" ]; then
mkdir dHdl_files
fi
for i in $(ls | grep "lambda.*"); do
lam="${i##*.}"
cp ./$i/PROD/dhdl.xvg ./dHdl_files/dhdl.$lam.xvg
done
这个脚本会将目录下所有lambda下的dhdl.xvg文件收集到创建的dhdl文件夹中并且重命名，之后便可进行下一步的数据分析。
