使用Amber对偶氮苯二面角进行伞状采样(Umbrella Sampling)

Amber Umbrella Sampling(伞状采样)教程

Umbrella Sampling 伞状采样

参考链接：
AMBER高级教程A17: 伞形采样计算平均力势PMF
MoBioChem

前言

  通常我们可能想要知道自由能随特定反应坐标的分布情况. 这种分布被称为平均力势, 它对识别过渡态, 中间产物和终点的相对稳定性是非常有用的, 一开始人们可能认为能够通过运行一次MD模拟来产生自由能随特定反应坐标的变化, 然后查看样品采样的概率. 然而, 通常感兴趣的能量位垒是kbT的数倍, 因此MD模拟将仍然处在它一开始的局部最小值或者跨越不同的最小值, 但极少会采样到过渡态. 采样的缺乏使得我们不能产生一个精确的PMF. 事实上, 甚至对于很小的系统我们也需要超过一毫秒的MD去获得可用的采样. 这显然很难达到, 因此需要使用更好的方法.
  其中一种方法是伞形采样. 这项工作可以通过分解反应坐标成为一系列窗口, 然后采用限制迫使反应坐标维持在靠近窗口的中心. 显然对于理想采样, 最好的方法是用自由能形貌自身作为伞形限制. 这将使能量表面足够平坦以至统一采样. 然而这不可行, 因为它需要我们在计算它之前知道能量表面形状. 相反, 我们通常所做的是用一系列简谐势去限制反应坐标到某个值. 这一系列表示采样的反应坐标值的直方图能够被产生, 假设窗口重叠, 这就能够使用加权直方图分析方法去处理结果, 从而去除限制偏置并恢复PMF.
  本教程利用amber程序采用伞状采样方法研究偶氮苯在水中的热异构化(研究不同二面角情况下偶氮苯的平均力势)。
本文所需文件

具体步骤

系统准备

1. 偶氮苯结构文件获取 进入ZINC数据库下载偶氮苯(Azobenzene)的 mol2文件

2. 在本机创建工作路径并处理偶氮苯

   # 将准备下载的偶氮苯结构文件命名为Azobenzene.mol2
   # 使用antechamber将mol2文件转换为gaussian文件
   antechamber -i Azobenzene.mol2 -fi mol2 -o lig.com -fo gcrt
   # 将gaussian文件转换为gaff原子格式，bcc电荷的mol2文件
   antechamber -i lig.com -fi gcrt -o lig.mol2 -fo mol2 -c bcc -at gaff
   # 生成偶氮苯的topol文件
   parmchk2 -i lig.mol2 -f mol2 -o lig.frcmod

   ! 使用两次antechamber是为了确保文件内格式的正确，以顺利的模拟。

   antechamber命令解释：

   • -i -o: 输入输出文件
   • -fi -fo: 指定输入输出文件格式
   • -c bcc: 指定电荷类型为AM1-BCC
   • -at gaff: 指定原子类型

初始MD

  在进行伞状采样之前，我们需要对偶氮苯系统进行平衡，使偶氮苯处于平衡状态之后，在进行伞状采样。

# 将偶氮苯的mol2文件和topol文件移到当前工作路径
cp ../1_init/lig.mol2 .
cp ../1_init/lig.frcmod .
# 生成偶氮苯溶剂化系统
tleap -f tleap.in

# tleap.in文件
# 导入力场
source leaprc.gaff2
source leaprc.water.tip3p
# 导入偶氮苯文件
MOL=loadmol2 lig.mol2
# 导入偶氮苯topol
loadamberparams lig.frcmod
# 生成偶氮苯lib
saveoff MOL lig.lib
# 生成含有偶氮苯的水盒子，八面体盒子，盒子边长为10A
solvateoct MOL WAT 10
# 保存系统的topol和坐标文件
saveamberparm MOL system.parm7 system.rst7
# 保存系统的pdb文件
savepdb MOL system.pdb
quit

  上述步骤，我们将偶氮苯放在了一个边长为10A的八面体盒子中，并对其溶剂化。生成的下一步使用的文件为system.parm7和system.rst7   生成以上文件之后，可以使用如下命令查看结构

vmd system.parm7 system.rst7

Minimization

sander -O -i min.in -p system.parm7 -c system.rst7 -o min.out -r min.rst7 -inf min.info

# min.in
Minimization
&cntrl
	imin=1  # 进行最小化
	maxcyc=4000 # 最大最小化步数 
	ncyc=1000 # 做完1000步最速下降法后转为共轭梯度法
	ntpr=2   # 每两步打印一次能量
/

Heating

sander -O -i heat.in -p system.parm7 -c min.rst7 -o heat.out -x heat.crd -r heat.rst7 -inf heat.info

# heat.in
Heat
&cntrl
	imin=0, # 关闭能量最小化
	ntx=1,  # 初始输入文件不带速度
	irest=0, # 不重启MD
	nstlim=10000, # 运行步数为10000
	dt=0.002,  # 步长为0.002ps
	ntf=2, # 不为SHAKE算法计算力
	ntc=2, # 使用SHAKE限制所有氢键
	tempi=0.0, # 初始温度
	temp0=300.0, # 终温
	ntpr=1, # 每一步都打印温度
	ntwx=10, # 每十步输出一次轨迹
	ioutfm=0, # 以二进制文件输出轨迹
	iwrap=1, 
	cut=8.0, # 非键截断 A
	ntb=1, # 固定体积周期性边界
	ntp=0, # 不控压
	ntt=3, # Langevin控温
	gamma_ln=2.0, # Langevin控温更新频率
	nmropt=1, # NMR限制
	ig=-1, # 随机种子数
/
&wt
	type='TEMP0',   # 升温程序 从0到8000步升温到300k
	istep1=0,istep2=8000,
	value1=0.0,value2=300.0,
/
&wt
	type='TEMP0', # 恒温8000步之后稳定300k
	istep1=8000,
	istep2=10000,
	value1=300.0, value2=300.0,
/
&wt
type='END'
/

Production

sander -O -i prod.in -p system.parm7 -c heat.rst7 -o prod.out -r prod.rst7 -x prod.crd -inf prod.info &
# 使用Amber内置脚本处理out文件得到运行结果
process_mdout.perl prod.out

# prod.in
Prod
&cntrl
    imin=0
    irest=0
    ntx=7   # 初始文件含有速度
    ntb=2   # 周期性边界条件
    ntp=1
    barostat=1 # Berendsen热浴
    pres0=1.0  # 初始压力
    taup=2.0  # 弛豫时间 ps
    cut=8.0
    ntc=2
    ntf=2
    tempi=300.0
    ntt=3
    gamma_ln=1.0
    nstlim=40000
    dt=0.002
    ntpr=200 # 每200步输出一次能量
    ntwx=200 # 每200步输出一次轨迹
    ntwr=200 # 每200步输出一次重启文件
    ntxo=1   # 最终输出格式为二进制
    ioutfm=0 
    iwrap=1
/

Umbrella Sampling

​ 限制二面角的MD

sander -O -i prod.in -p system.parm7 -c prod.rst7 -r prod_180.rst7 -x prod_180.crd -inf prod_180.info &

# disang
Harmonic restraints for 180 deg
&rst
    iat=4,5,6,7,  # 二面角涉及原子
    # 当二面角处于(r1,r2)偏置势为200，具体说明请见amber 25.1 章节
    r1=60.0,r2=180.0,r3=180.0,r4=300.0, 
    rk2=200,rk3=200,
/

# prod.in
Prod
&cntrl
    imin=0
    irest=0
    ntx=7
    ntb=2
    ntp=1
    barostat=1
    pres0=1.0
    taup=2.0
    cut=8.0
    ntc=2
    ntf=2
    tempi=300.0
    ntt=3
    gamma_ln=1.0
    nstlim=2000
    dt=0.002
    ntpr=10
    ntwx=10
    ntwr=20
    ntxo=1
    ioutfm=0
    nmropt=1
/
&wt
    type='DUMPFREQ',  # 每50步输出一次二面角值
    istep1=50,
&end
&wt
    type='END',
&end
DISANG=disang   # 二面角设置
DUMPAVE=dihedral_180.dat

伞状采样

  进行采样需要多个窗口，创建脚本，生成伞状采样的窗口。
  此脚本使用二面角为180的结构作为初始文件。依次运行二面角逐渐减少3度的模型。意味着一共会产生60个窗口。

bash umbrella.sh
# 脚本运行完成之后运行如下命令
for i in { 1..9 }
do
mv md$i.crd md0$i.crd
done
# 将不同窗口帧合成一个轨迹文件
parm system.parm7
trajin *.crd 1 last 20
autoimage :MOL
trajout final.crd
go
exit
vmd system.parm7 final.crd
# 查看不同窗口之间是否有重合
xmgrace dihedral{1..180}.out

# umbrella.sh
#!/bin/bash -f

delta=$(echo "-3" | bc)
rini=$(echo "180" | bc)
limit=61

n=1

while [ $n -le $limit ]
do

r=$(echo "$rini+$delta*($n-1)" | bc)
rmin=$(echo "$r-30.0" | bc)
rmax=$(echo "$r+30.0" | bc)

[ -e md$r.in ] && rm md$r.in
[ -e dihedral$r.dat ] && rm dihedral$r.dat

# Print MD file
echo umbrellatest                    >> md$r.in
echo \&cntrl                         >> md$r.in
echo imin=0                          >> md$r.in
echo irest=1                         >> md$r.in
echo ntx=7                           >> md$r.in
echo ntb=2                           >> md$r.in
echo ntp=1                           >> md$r.in
echo barostat=1                      >> md$r.in  
echo pres0=1.0                       >> md$r.in 
echo taup=2.0                        >> md$r.in 
echo cut=8.0                         >> md$r.in 
echo ntc=2                           >> md$r.in 
echo ntf=2                           >> md$r.in 
echo tempi=300.0                     >> md$r.in 
echo temp0=300.0                     >> md$r.in 
echo ntt=3                           >> md$r.in 
echo gamma_ln=1.0                    >> md$r.in 
echo nstlim=4000                     >> md$r.in 
echo dt=0.002                        >> md$r.in 
echo ntpr=10                         >> md$r.in 
echo ntwx=5                          >> md$r.in 
echo ntwr=5                          >> md$r.in 
echo ntxo=1                          >> md$r.in 
echo ioutfm=0                        >> md$r.in 
echo nmropt=1                        >> md$r.in 
echo /                               >> md$r.in 
echo \&wt                            >> md$r.in 
echo type=\'DUMPFREQ\',              >> md$r.in 
echo istep1=1,                       >> md$r.in 
echo \&end                           >> md$r.in 
echo \&wt                            >> md$r.in 
echo type=\'END\',                   >> md$r.in 
echo \&end                            >> md$r.in 
echo DISANG=dihedral$r.dat           >> md$r.in 
echo DUMPAVE=dihedral$r.out          >> md$r.in 

echo Harmonic restraints for $r dihedral  >> dihedral$r.dat
echo \&rst                                >> dihedral$r.dat
echo iat=4,5,6,7                          >> dihedral$r.dat
echo r1=$rmin,r2=$r,r3=$r,r4=$rmax        >> dihedral$r.dat
echo rk2=200,rk3=200,                     >> dihedral$r.dat
echo /                                    >> dihedral$r.dat

rnext=$(echo "$r + $delta" | bc)

sander -O -i md$r.in -o md$r.out -p system.parm7 -c ini$r.rst7 -r md$r.rst7 -ref ini$r.rst7 -x md$n.crd -inf md$r.info

cp md$r.rst7 ini$rnext.rst7

((n++))
done

使用WHAM生成PMF

  WHAM软件可以去WHAM – Grossfield Lab下载，安装很简单。
  使用WHAM之前，我们需要使用脚本生成meta.dat。此文件有三列，由二面角文件名，二面角值和力常数组成。

Please remember that Amber uses k(x − x₀)² with k in kcal/mol/rad² whereas the WHAM program uses 0.5k(x-x₀)² with k in kcal/mol/deg² therefore we need to multiply our force constant by 2(π/180)².

可以使用脚本生成neta.dat文件也可以手动完成。

# create meta-perl
#!/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=180;
	while ($dihed>=0) {
	print "Processing dihedral:  $dihed \n"
	&write_meta();
	$dihed +=$incr;
	}
}

sub write_meta {
	open METAFILE,">>","$name";
	print METAFILE<<EOF;
dihedral$dihed.out $dihed.0 $force
EOF
	close MDINFILE;
}

# meta.dat
dihedral180.out	180.0 	0.12184
dihedral177.out	177.0 	0.12184
dihedral174.out	174.0 	0.12184
dihedral171.out	171.0 	0.12184
dihedral168.out	168.0 	0.12184
dihedral165.out	165.0 	0.12184
dihedral162.out	162.0 	0.12184
dihedral159.out	159.0 	0.12184
dihedral156.out	156.0 	0.12184
dihedral153.out	153.0 	0.12184
dihedral150.out	150.0 	0.12184
dihedral147.out	147.0 	0.12184
dihedral144.out	144.0 	0.12184
dihedral141.out	141.0 	0.12184
dihedral138.out	138.0 	0.12184
dihedral135.out	135.0 	0.12184
dihedral132.out	132.0 	0.12184
dihedral129.out	129.0 	0.12184
dihedral126.out	126.0 	0.12184
dihedral123.out	123.0 	0.12184
dihedral120.out	120.0 	0.12184
dihedral117.out	117.0 	0.12184
dihedral114.out	114.0 	0.12184
dihedral111.out	111.0 	0.12184
dihedral108.out	108.0 	0.12184
dihedral105.out	105.0 	0.12184
dihedral102.out	102.0 	0.12184
dihedral99.out	99.0 	0.12184
dihedral96.out	96.0 	0.12184
dihedral93.out	93.0 	0.12184
dihedral90.out	90.0 	0.12184
dihedral87.out	87.0 	0.12184
dihedral84.out	84.0 	0.12184
dihedral81.out	81.0 	0.12184
dihedral78.out	78.0 	0.12184
dihedral75.out	75.0 	0.12184
dihedral72.out	72.0 	0.12184
dihedral69.out	69.0 	0.12184
dihedral66.out	66.0 	0.12184
dihedral63.out	63.0 	0.12184
dihedral60.out	60.0 	0.12184
dihedral57.out	57.0 	0.12184
dihedral54.out	54.0 	0.12184
dihedral51.out	51.0 	0.12184
dihedral48.out	48.0 	0.12184
dihedral45.out	45.0 	0.12184
dihedral42.out	42.0 	0.12184
dihedral39.out	39.0 	0.12184
dihedral36.out	36.0 	0.12184
dihedral33.out	33.0 	0.12184
dihedral30.out	30.0 	0.12184
dihedral27.out	27.0 	0.12184
dihedral24.out	24.0 	0.12184
dihedral21.out	21.0 	0.12184
dihedral18.out	18.0 	0.12184
dihedral15.out	15.0 	0.12184
dihedral12.out	12.0 	0.12184
dihedral9.out	9.0 	0.12184
dihedral6.out	6.0 	0.12184
dihedral3.out	3.0 	0.12184
dihedral0.out	0.0 	0.12184

​ 运行wham软件

./wham P 0 180 60 0.01 300 0  meta.dat result.dat

P：反应坐标是周期性的
0 180 60 ：从0-180 使用60个bins生成PMF
0.01：PMF的容忍度
300：温度
0：不填充数据点
meta.dat: meta file
result.dat: output file

# 将数据提取出来
cat result.dat | awk "{print $1,$2}" > pmf.dat
# 绘制pmf随反应坐标变化图，也可以手动提取数据到熟悉的软件中作图
xmgrace pmf.dat

后话

  因为自身过于懒惰，本来想把该教程写的更加详实以及直白，但是偷懒了，最后也就这样了。在此根据该教程的使用软件情况给自己挖几个博客的坑.什么时候能够填上？有谁会知道呢？

• perl/bash语言使用教程
• wham、xmgrace/qtgace安装教程--->极致简单
• TCL语言开发VMD插件