一.内容简介
在本次练习模拟中，我们的模拟目标是采集到发生顺反异构的丙氨酸三肽的足够样本，计算其平均力势。
问题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.tip3p
mol=sequence {ACE ALA ALA NME}
solvateoct mol TIP3PBOX 10.0
#保存用于模拟的inp和top文件：
saveamberparm mol mol.parm mol.crd
#保存为pdb文件，可以用来对轨迹中的平移和旋转进行修正：
savepdb mol mol.pdb
quit
在终端输入：
tleap -f tleap.in2.对模型进行弛豫：
分别编写如下的最小化，升温和弛豫的输入文件：
#最小化模拟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
&cntrl
imin = 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.120
#弛豫的in文件：
50 ps NPT equilibration for 120 deg
&cntrl
imin = 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
&wt
type='END',
&end
DISANG=disang.120
#成品模拟的in文件：
100 ps NPT production for 120 deg
&cntrl
imin = 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
&wt
type='DUMPFREQ', istep1=50,
&end
&wt
type='END',
&end
DISANG=disang.120
DUMPAVE=dihedral_120.datDISANG=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
&rst
iat=9,15,17,19,
r1=-60.0, r2=120.0, r3=120.0, r4=300.0,
rk2=200, rk3=200,
&end3.编写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.nc4.运行模拟：
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
&cntrl
imin = 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.60
#equilibration的in文件：
50 ps NPT equilibration for 60 deg
&cntrl
imin = 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
&wt
type='END',
&end
DISANG=disang.60
#produce的in文件：
100 ps NPT production for 60 deg
&cntrl
imin = 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
&wt
type='DUMPFREQ', istep1=50,
&end
&wt
type='END',
&end
DISANG=disang.60
DUMPAVE=dihedral_60.dat
#restrain的in文件：
Harmonic restraints for 60 deg
&rst
iat=9,15,17,19,
r1=-120.0, r2=60.0, r3=60.0, r4=240.0,
rk2=200, rk3=200,
&end2.同样是用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 it
print "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
&cntrl
imin = 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
EOF
close MDINFILE;
}
sub write_mdin1 {
open MDINFILE,'>', "mdin_equi.$dihed";
print MDINFILE <<EOF;
50 ps NPT equilibration for $dihed deg
&cntrl
imin = 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
EOF
close MDINFILE;
}
sub write_mdin2 {
open MDINFILE,'>', "mdin_prod.$dihed";
print MDINFILE <<EOF;
100 ps NPT production for $dihed deg
&cntrl
imin = 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
&wt
type='DUMPFREQ', istep1=50,
&end
&wt
type='END',
&end
DISANG=disang.$dihed
DUMPAVE=dihedral_${dihed}.dat
EOF
close 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
&rst
iat=9,15,17,19,
r1=$left, r2=$min, r3=$min, r4=$right,
rk2=${force}, rk3=${force},
&end
EOF
close 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_$dihed
cd $wd
pmemd -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
EOF
print 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
EOF
close 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.dat
xmgrace pmf.dat查看势能面图像：
可见，在二面角在180度时，势能最低
参考文章：
3121
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
