学习Amber T3.3：隐式溶剂模型(GB)的MD

一、简介

隐式溶剂模型，这个模型里只有溶质分子没有溶剂分子（水，离子都没有），这样的模型好处是减少MD的计算量（毕竟原子数量大大减少），而且可以合理考虑溶剂分子对溶质的势能变化影响。

Amber中常用的隐式溶剂模型是广义的Born模型（GB模型），这个模型的算法在力场ff99SB和ff14Bonlysc中有不错的测试结果。

二、过程

1、准备模拟的蛋白质：1L2Y

下载链接：https://www.rcsb.org/structure/1L2Y

对下载的pdb文件进行处理：

pdb4amber 1l2y.pdb > xx.pdb

tleap

source leaprc.protein.ff19SB

set default PBradii mbondi3   #保证添加的H原子间距符合一个GB模型

protein=loadpdb xx.pdb

saveamberparm protein xx.parm7 xx.rst7

savepdb protein xx.pdb

得到进行模拟的top和crd文件，还有加了H的pdb文件。 

2、准备模拟的输入文件： 

进行能量最小化：

我们用tleap或是charmm-GUI进行体系加水，蛋白质加H时，原子和原子之间可能出现不好的接触。若是不处理这种不好的接触，模拟过程中就会出现局部能量过高最终破坏体系。能量最小化的作用就是找到一个能量最低的构象，消除这些不好的接触。

编写01min.in文件：

energy minimization
&cntrl
  imin = 1,                        #1进行最小化，0不进行最小化
  maxcyc=100,                 #最小化最多进行100次
  ntx = 1,                          #从没有速度的rst文件中读入坐标
  ntwr = 100, ntpr = 10,    #前者设置输出到rst的步数，后者设置输出到out的步数
  ioutfm=0, ntxo=1,          #前者指出nc文件的格式二进制，netcdf，后者指定rst文件格式ascii
  cut = 1000.0,                 #非键截断距离（单位埃）
  ntb=0,                            #没有周期性边界（毕竟没有水嘛）
  igb = 8,                          #指定这是一个隐式溶剂模型 GBneck2

  gbsa=3, surften=0.007, #指定溶剂可及表面积和表面张力
  saltcon = 0.0,
&end

进行升温：

我们不能直接在想要的温度下进行模拟，而是由一个低温向高温慢慢爬上去，Amber一般用到的恒温器是Langevin thermostat，设置的碰撞常量γ(gamma=1.0ps-1)，一般目标温度是300K。这个过程中，我们会对主链原子(CA,C,O,N)加上限制力进行约束。

heating with backbone restraints
&cntrl
   ntx = 1,                            
   ntwx = 5000,ntwr = 500, ntpr = 5000,

   ntwe = 0,                  #ntwe步输出能量和温度，设为0就没有          
   ioutfm=0, ntxo=1,     #二进制格式的nc和ascii格式的rst                
   imin = 0,
   nstlim = 500000, dt = 0.002, #dt：每步的时长，单位是ps
   ntt = 3,gamma_ln=1.,temp0 = 100.0,   #ntt=3，用Langevin thermostat
   nscm = 1000,
   ig = -1,         #随机数种子
   ntc= 2, ntf = 2, #对H原子键振动的约束（shake算法）
   cut = 1000,    #对原子之间非键作用力的限制，距离大于1000A就忽略这两原子间的作用力

   igb=8,
   gbsa=3, surften=0.007,
   ntb = 0,
   saltcon=0.,
   ntr=1,restraintmask="@CA,N,C,O",  restraint_wt=10.0, #ntr=1：将对指定原子进行约束
   nmropt=1,    #1：读取下面的&wt区域参数                         
&end
&wt
   TYPE="TEMP0", istep1=0, istep2=500000,
   value1=100., value2=300.,        #0步时，温度是100K，50000步(1ns)时温度是300K    
&end
&wt
   TYPE="END",
&end 

进行平衡：

这个模拟中，体系在我们设置的最终温度下进行模拟，模拟时间长，分两个阶段，依次减小对主链原子的限制力。

编写03equil01.in，跑250ps：

&cntrl
   ntx = 5,                              

   ntwx = 5000, ntwe = 0, ntwr = 500, ntpr = 5000,                 

   ioutfm=0, ntxo=1, 

   imin = 0, nstlim = 125000, dt = 0.002, 

   ntt = 3, gamma_ln=1.0, temp0 = 300.0,  
   nscm = 1000, ig = -1, irest = 1,        #Iirest=1：重启模拟，从rst文件中读取速度和坐标  

   ntc= 2, ntf = 2, cut = 1000,           

   igb=8, gbsa=3, surften=0.007, ntb = 0,
   saltcon=0., 

   ntr=1, restraintmask="@CA,N,C,O",     
   restraint_wt=1.0, 
/

编写03equil02.in：

&cntrl
   ntx = 5,                             
   ntwx = 5000, ntwe = 0, ntwr = 500,    
   ntpr = 5000,                 

   ioutfm=0, ntxo=1,                     

   imin = 0, nstlim = 125000, dt = 0.002, 

   ntt = 3, gamma_ln=1., temp0 = 300.0,  
   nscm = 1000, ig = -1, irest = 1,                 

   ntc= 2, ntf = 2, cut = 1000,           

   igb=8, gbsa=3, surften=0.007, ntb = 0,
   saltcon=0., 

   ntr=1, restraintmask="@CA,N,C,O",     
   restraint_wt=0.1, 
/

 两个平衡模拟的区别在于restraint_wt参数值的减小，从10到1到0.1，表示对主链原子的限制力越来越小，运动会更加自然。restraint_wt的单位是 kcal/mol/angstrom2，越大，原子的运动受限越厉害。

3、进行模拟：

编写sh脚本：

#！/usr/bin/bash

export CUDA_VISIBLE_DEVICES=3  #指定跑模拟的GPU型号

export path=/home/packages/amber20/bin  #说明路径为全局变量
$path/pmemd.cuda -O -i 01min.in -p xx.parm7 -c xx.rst7 -o min.out\
-x min.crd -inf min.info -r min.rst7

echo  -e “\033[1;5;31;40m minimize is finished \033[0m"

$path/pmemd.cuda -O -i 02heat.in -p xx.parm7 -c min.rst7 -ref min.rst7\
-o heat.out -x heat.crd -inf heat.info -r heat.rst7

echo -e “\033[1;5;31;40m heating is finished \033[0m"

$path/pmemd.cuda -O -i 03equil01.in -p xx.parm7 -c heat.rst7 -ref heat.rst7\
-o eq1.out -x eq1.crd -inf eq1.info -r eq1.rst7

echo -e “\033[1;5;31;40m equilibrium 01 is finished \033[0m"

$path/pmemd.cuda -O -i 03equil02.in -p xx.parm7 -c eq1.rst7 -ref eq1.rst7\
-o eq2.out -x eq2.crd -inf eq2.info -r eq2.rst7

echo -e “\033[1;5;31;40m equilibrium 02 is finished \033[0m"

三、参考链接：

Equilibration of Implicit Solvent