学习amber教程A17：伞形采样，绘制丙氨酸三肽的势能面

一.内容简介

在本次练习模拟中，我们的模拟目标是采集到发生顺反异构的丙氨酸三肽的足够样本，计算其平均力势。

问题1：由于在MD模拟中，分子很难跨过高势能垒而发生很大的构象改变，一般模拟过程下，分子都稳定在一个较低的能量状态下不思进取。或者说，发生分子跨过高能垒得到一个变构的构型的概率时很低很低的，运行几ps的模拟也采不到这样的构象。本次的练习模拟中，丙氨酸三肽发生顺反异构现象也需要跨越高势能垒。

解决1：我们给发生顺反异构的酰胺键添加一个人工的简谐势能k(X-X0)^2，强迫这个键在我们要求的角度范围内旋转，得到我们想要的二面角（顺反异构等效于将中间的酰胺键扭转180度）。

这个人工的简谐势能是一个二次函数，会强迫酰胺键在X0这个角度周围小范围的旋转（角度变得离这个X0太远，势能也就越大，酰胺键受到的阻力也越大）

问题2：但是光有一个酰胺键扭转180度和没扭转两个构象是不够的，顺反异构是一个连续变化的过程（从0度到180度），这两个构象中的酰胺键都不可能会在模拟中发生从0度到180度的变化。

解决2：两个模型不够，就做多个模型。将0~180度，每隔3度做一个模型，就是对酰胺键扭转3度（间隔不能太大，太大会导致采样不精准），得到91个构象进行模拟，这样的模拟就采集到顺反异构过程中二面角的变化数据了。

问题3：我们是强行添加一个不存在的简谐势能，这样得到的平均力势(PMF)靠谱吗？

解决3：我们每隔3度构建一个模型，这样可以保证模型之间的采样（数据）会出现重叠。有重叠，就可以用加权直方图分析方法去处理结果, 从而去除限制偏置并恢复PMF。利用WHAM可以实现这点（原来如此，完全听不懂）。

二.建立和弛豫初始结构

1.建立模型：

编写如下的tleap.in文件

source leaprc.protein.ff19SB
source leaprc.water.tip3pmol=sequence {ACE ALA ALA NME}solvateoct mol TIP3PBOX 10.0#保存用于模拟的inp和top文件：
saveamberparm mol mol.parm mol.crd#保存为pdb文件，可以用来对轨迹中的平移和旋转进行修正：
savepdb mol mol.pdbquit

在终端输入：

tleap -f tleap.in

2.对模型进行弛豫：

分别编写如下的最小化，升温和弛豫的输入文件：

#最小化模拟1000步：
minimize&cntrl  imin=1, maxcyc=1000, ncyc=500,  cut=8.0, ntb=1, ntp=0,  ntc=2, ntf=2,/#升温，时间步长为2fs，模拟20ps，从0K到300K:
heat NVT 20ps&cntrl  imin=0, irest=0, ntx=1,  nstlim=10000, dt=0.002,  ntc=2, ntf=2, cut=8.0,  ntt=3, gamma_ln=1.0,  ntb=1, ntp=0,  tempi=0.0, temp0=300.0,  ntpr=100, ntwx=100, ntwr=10000,  ioutfm=1,/#弛豫，模拟100ps：
equil NPT 100ps&cntrl  imin=0, irest=1, ntx=5,  nstlim=50000, dt=0.002,  ntc=2, ntf=2, cut=8.0,  ntt=3, gamma_ln=1.0,  ntb=2, ntp=1, tautp=1.0,  temp0=300.0,  ntpr=500, ntwx=500, ntwr=10000,  ioutfm=1,/

分别命名为01min.in、02heat.in和03equil.in

编写sh脚本md.sh进行模拟：

#!/usr/bin/bash
pmemd -O -i 01min.in -o 01min.out -p mol.parm -c mol.crd -r 01_min.rst
echo -e "\033[1m minimize is finish! \033[0m"
pmemd -O -i 02heat.in -o 02heat.out -p mol.parm -c 01_min.rst -r 02heat.rst -x 02heat.nc
echo -e "\033[1m heat is finish! \033[0m"
pmemd -O -i 03equil.in -o 03equil.out -p mol.parm -c 02heat.rst -r 03equil.rst -x 03equil.nc
echo -e "\033[1m equilibrium is finish! \033[0m"

运行脚本：

sh md.sh

三.添加简谐势，模拟人为约束酰胺键扭转120度的构象

1.编写如下的局部能量最小化、弛豫和产品采样in文件：

#能量最小化的in文件：
2000 step minimization for 120 deg&cntrlimin = 1,maxcyc=2000, ncyc = 500,ntpr = 100, ntwr = 1000,ntf = 1, ntc = 1, cut = 8.0,ntb = 1, ntp = 0,nmropt = 1,&end&wttype='END',&end
DISANG=disang.120#弛豫的in文件：
50 ps NPT equilibration for 120 deg&cntrlimin = 0, ntx = 1, irest = 0,ntpr = 5000, ntwr = 50000, ntwx = 0,ntf = 2, ntc = 2, cut = 8.0,ntb = 2, nstlim = 50000, dt = 0.001,tempi=0.0, temp0 = 300.0, ntt = 3,gamma_ln = 1.0,ntp = 1, pres0 = 1.0, taup = 5.0,nmropt = 1, ioutfm=1,&end&wttype='END',&end
DISANG=disang.120#成品模拟的in文件：
100 ps NPT production for 120 deg&cntrlimin = 0, ntx = 5, irest = 1,ntpr = 10000, ntwr = 0, ntwx = 10000,ntf = 2, ntc = 2, cut = 8.0,ntb = 2, nstlim = 100000, dt = 0.001,temp0 = 300.0, ntt = 3,gamma_ln = 1,ntp = 1, pres0 = 1.0, taup = 5.0,nmropt = 1, ioutfm = 1,&end&wttype='DUMPFREQ', istep1=50,&end&wttype='END',&end
DISANG=disang.120
DUMPAVE=dihedral_120.dat

DISANG=disang.120：disang.120是下面描述的文件，内容是我们添加的简谐势信息

disang.120务必放在sh脚本的同级目录下。

DUMPAVE=dihearal_120.dat：dihearal_120.dat是一个输出文件，记录模拟过程中角度变化值

nmropt=1：读取之后&wt区域的内容，在type=“DUMPFREQ”中，写下的istep1=50表示每隔50步记录一次二面角度数值

2.编写disang.120文件，这里告知pmemd添加简谐势能的信息：

Harmonic restraints for 120 deg&rstiat=9,15,17,19,r1=-60.0, r2=120.0, r3=120.0, r4=300.0,rk2=200, rk3=200,&end

3.编写md脚本md_120.sh，进行模拟：

pmemd -O -i mdin_min.120 -o ala_tri_min_120.out -p ala_tri.prmtop -c 03_equil.rst -r ala_tri_min_120.rst
pmemd -O -i mdin_equi.120 -o ala_tri_equil_120.out -p ala_tri.prmtop -c ala_tri_min_120.rst -r ala_tri_equil_120.rst
pmemd -O -i mdin_prod.120 -o ala_tri_prod_120.out -p ala_tri.prmtop -c ala_tri_equil_120.rst -r ala_tri_prod_120.rst -x ala_tri_prod_120.nc

4.运行模拟：

sh md_120.sh

四.将120度的构象进行限制，进行扭转60度的模拟

根据这个模板，我们编写由120度扭转到60度的in文件，进行酰胺键扭转60度的模拟和采样

1.编写对应的四个in文件（01min.in、02heat.in、03equil.in和disang60.in）：

#minimize的in文件：
2000 step minimization for 60 deg&cntrlimin = 1,maxcyc=2000, ncyc = 500,ntpr = 100, ntwr = 1000,ntf = 1, ntc = 1, cut = 8.0,ntb = 1, ntp = 0,nmropt = 1,&end&wttype='END',&end
DISANG=disang.60#equilibration的in文件：
50 ps NPT equilibration for 60 deg&cntrlimin = 0, ntx = 1, irest = 0,ntpr = 5000, ntwr = 50000, ntwx = 0,ntf = 2, ntc = 2, cut = 8.0,ntb = 2, nstlim = 50000, dt = 0.001,tempi=0.0, temp0 = 300, ntt = 3,gamma_ln = 1.0,ntp = 1, pres0 = 1.0, taup = 5.0,nmropt = 1,&end&wttype='END',&end
DISANG=disang.60#produce的in文件：
100 ps NPT production for 60 deg&cntrlimin = 0, ntx = 5, irest = 1,ntpr = 10000, ntwr = 0, ntwx = 10000,ntf = 2, ntc = 2, cut = 8.0,ntb = 2, nstlim = 100000, dt = 0.001,temp0 = 300.0, ntt = 3,gamma_ln = 1,ntp = 1, pres0 = 1.0, taup = 5.0,nmropt = 1, ioutfm = 1,&end&wttype='DUMPFREQ', istep1=50,&end&wttype='END',&end
DISANG=disang.60
DUMPAVE=dihedral_60.dat#restrain的in文件：
Harmonic restraints for 60 deg&rstiat=9,15,17,19,r1=-120.0, r2=60.0, r3=60.0, r4=240.0,rk2=200, rk3=200,&end

2.同样是用pmemd进行模拟和采样。

五.批量进行0~90度，90~150度，150~180度的模拟和采样

这里直接修改教程里提供的perl脚本就可以了，命名为md.perl

#!/usr/bin/perl -w
use Cwd;
$wd=cwd;print "Preparing input files\n";$name="ala_tri";   #这里修改inp和top文件的路径
$incr=3;
$force=200.0;&prepare_input();exit(0);sub prepare_input() {$dihed=0;                  #这里修改模拟的角度最小值while ($dihed <= 90) {     #这里修改模拟的角度最大值if ($dihed != 60) {      #这里修改要跳过的模拟，分别是60,120和180#skip 60 degrees since we already ran itprint "Processing dihedral: $dihed\n";&write_mdin0();&write_mdin1();&write_mdin2();&write_disang();&write_batchfile();system("sh run.pbs.$dihed");     #这里改为用sh运行}$dihed += $incr;}
}sub write_mdin0 {open MDINFILE,'>', "mdin_min.$dihed";print MDINFILE <<EOF;
2000 step minimization for $dihed deg&cntrlimin = 1, maxcyc=2000, ncyc = 500,ntpr = 100, ntwr = 1000,ntf = 1, ntc = 1, cut = 8.0,ntb = 1, ntp = 0,nmropt = 1,&end&wt type='END',&end
DISANG=disang.$dihed
EOFclose MDINFILE;
}
sub write_mdin1 {open MDINFILE,'>', "mdin_equi.$dihed";print MDINFILE <<EOF;
50 ps NPT equilibration for $dihed deg&cntrlimin = 0, ntx = 1, irest = 0,ntpr = 5000, ntwr = 50000, ntwx = 0,ntf = 2, ntc = 2, cut = 8.0,ntb = 2, nstlim = 50000, dt = 0.001,tempi=0.0, temp0 = 300, ntt = 3,gamma_ln = 1.0,ntp = 1, pres0 = 1.0, taup = 5.0,nmropt = 1, ioutfm=1,&end&wt type='END',&end
DISANG=disang.$dihed
EOFclose MDINFILE;
}
sub write_mdin2 {open MDINFILE,'>', "mdin_prod.$dihed";print MDINFILE <<EOF;
100 ps NPT production for $dihed deg&cntrlimin = 0, ntx = 5, irest = 1,ntpr = 10000, ntwr = 0, ntwx = 10000,ntf = 2, ntc = 2, cut = 8.0,ntb = 2, nstlim = 100000, dt = 0.001,temp0 = 300, ntt = 3,gamma_ln = 1.0,ntp = 1, pres0 = 1.0, taup = 5.0,nmropt = 1, ioutfm=1,&end&wttype='DUMPFREQ', istep1=50,&end&wt type='END',&end
DISANG=disang.$dihed
DUMPAVE=dihedral_${dihed}.dat
EOFclose MDINFILE;
}
sub write_disang {$left  = sprintf "%4.1f", $dihed - 180;$min   = sprintf "%4.1f", $dihed;$right = sprintf "%4.1f", $dihed + 180;open DISANG,'>', "disang.$dihed";print DISANG <<EOF;
Harmonic restraints for $dihed deg&rstiat=9,15,17,19,r1=$left, r2=$min, r3=$min, r4=$right,rk2=${force}, rk3=${force},&end
EOFclose DISANG;
}sub write_batchfile {open BATCHFILE, '>', "run.pbs.$dihed";print BATCHFILE <<EOF;
#PBS -l nodes=1:ppn=2
#PBS -l walltime=1:00:00
#PBS -N run_$dihedcd $wdpmemd -O -i mdin_min.$dihed -p ${name}.prmtop -c ala_tri_prod_60.rst -r ${name}_min_${dihed}.rst  -o ${name}_min_${dihed}.out
pmemd -O -i mdin_equi.$dihed -p ${name}.prmtop -c ${name}_min_${dihed}.rst  -r ${name}_equil_${dihed}.rst -o ${name}_equil_${dihed}.out
pmemd -O -i mdin_prod.$dihed -p ${name}.prmtop -c ${name}_equil_${dihed}.rst -r ${name}_prod_${dihed}.rst -o ${name}_prod_${dihed}.out -x ${name}_prod_${dihed}.mdcrd
EOFprint BATCHFILE "\necho \"Execution finished\"\n";close BATCHFILE;
}

运行这个脚本，就会编写好一个，跑一个，直到跑完规定的角度：

perl md.perl

跑完所有的模拟，编写如下脚本，将得到的一系列dihedral_${dihed}.dat存到一个文件里：

#!/usr/bin/python3
import os
for i in (3,180,3):os.system("cat dihedral_"+i+".dat > all_dihedral.dat")
del i

用xmgrace进行绘图：

这里的x轴标签改为step ，由于我是隔5度一个模型进行模拟的，反应坐标会少一些，图像就会显得度数的覆盖程度不高，看起来较稀疏。

六.运行WHAM

1.WHAM的tgz包的下载：WHAM的下载链接

2.使用脚本建立meta.dat文件：

#!/usr/bin/perl -w
use Cwd;
$wd=cwd;print "Preparing meta file\n";$name="meta.dat";
$incr=3;
$force=0.12184;&prepare_input();exit(0);sub prepare_input() {$dihed=0;while ($dihed <= 180) { print "Processing dihedral: $dihed\n";&write_meta();$dihed += $incr;}
}sub write_meta {open METAFILE,'>>', "$name";print METAFILE <<EOF;
dihedral_$dihed.dat $dihed.0 $force
EOFclose MDINFILE;
}

这里的meta.dat里指定了要调用的dihedral文件，所以生成的91个dihedral_${name}.dat要放在当前目录下，运行wham：

wham P -1 181 61 0.01 300 0 meta.dat result.dat#P: 反应坐标是周期性变化的，即0度=360度
#-1 181: 从-1度到180度
#61 0.01: 采样61个模型，能量的最单位是0.01kcal
#300: 温度是300K
#0: 不加数据点

提取result的前两列，xmgrace绘制PMF图：

cat result.dat | awk '{print$1,$2}' > pmf.datxmgrace pmf.dat

查看势能面图像：

可见，在二面角在180度时，势能最低

参考文章：

amber教程A17