模拟计划准备

进行分子动力学模拟的步骤：

1.准备模拟需要的文件：

autodock对接完成得分高的复合物构象。

文本编辑器内对构象中分子分离，得到:

1): 带有受体和配体的pdb文件

2): 不含有配体的受体pdb文件

注意：在进行tleap操作时，为了避免报错，需要删掉蛋白和配体的H原子，使用pdb4amber可以完成这一点：

pdb4amber -i complex.pdb -o complex_notH.pdb --dry --nohyd
 
pdb4amber -i protein.pdb -o protein_notH.pdb --dry --nohyd
 
#--dry：去掉pdb里的WAT信息
#--nohyd：去掉体系里的H原子信息

进行tleap进行配置力场，加水盒子，计算电荷，加抗衡离子，将复合物保存为inp和top文件：

#导入leaprc配置文件：
source leaprc.protein.ff19SB
source leaprc.water.tip3p
 
#设置原子半径：
#set default PBRadii mbondi3   #不是隐式溶剂模型
 
#导入复合物pdb文件，tleap会补H和缺失的重原子
mol = loadpdb complex_notH.pdb 
 
#添加水盒子
solvateBox mol TIP3PBOX 10 
 
#在添加抗衡离子前，用charge命令查看当前体系的电荷总和：
 charge mol
 
#若是显示整体带正电，就添加CL-，相反，若是体系是带负电，就添加Na+
#这里添加Cl-，使体系整体电荷为0，这里注意命令的大小写
addIons2 mol Cl- 0
 
#保存为parm7和crd文件：
saveAmberParm mol complex.parm7 complex.rst7
 
#保存为pdb：
savepdb mol complex.pdb
 
------------------------------------------
#实际操作上，发现168的SER和小配体的NTYR多出些原子，
#我们用tleap的命令查看对应残基具体的原子名：
desc SER.SER
desc NTRY.NTYR
#对应原子名，在pdb文件里删除或修改

同样的操作，保存一个不含配体的蛋白溶剂体系。

建议使用ambpdb进行对parm7和rst7文件进行整合保存为pdb，检查pdb，确保拓扑和坐标文件正确：

#使用ambpdb合成pdb：
ambpdb -p complex.parm7 < complex.rst7 > complex_solvated.pdb

2.对体系进行最小化和控温，恒容，恒压模拟：

目前体系里的水分布还是很不真实的（一个立方体盒子），需要经过最小化、升温、恒容-恒温模拟、恒压-恒温模拟调整水的密度达到真实水平（1g/cm3）

这里的能量最小化分两步走：①约束溶质所有非H原子；②约束溶质所有重原子。

在参考文献里：

使用最陡下降和共轭梯度法进行能量最小化（常规操作）

采用 2 fs积分步骤并使用蛙跳算法更新时间步长（也是常规操作）

静电相互作用使用PME(Particle-mesh Ewald)求和方法处理（是默认吧）。

01.编写第一次最小化的01min.in文件：

 
Minimization with Cartesian restraints for the solute
 &cntrl
  imin=1, 
  maxcyc=1000,ncyc=500,dt=0.002,
  ntpr=5, 
  ntr=1,  #该关键字告知需要进行约束
  ntr=1,restraintmask=":1-292&!@H=",  restraint_wt=10.0,
 &end

这里进行1000步最小化，前500步用最陡下降法，后500步用共轭梯度法(maxcyc=1000,ncyc=500)

常用2fs的积分步长(dt=0.002ps)。

这里对蛋白和配体进行约束，只是让水分子放开跑。

02.看看能不能跑：

第一次最小化，对蛋白和配体进行约束，试试看能不能跑：

#export CUDA_VISIBLE_DEVICES=3  #指定跑模拟的GPU型号为3号（具体看情况）
 
#pmemd -O -i min.in -p complex.parm7 -c complex.rst7 -r complex_first_min.rst \
#-o complex_first_min.out -ref complex.rst7
 
#模拟用的是sander：
 
sander -O -i min.in -p complex.parm7 -c complex.rst7 -r complex_first_min.rst \
-o complex_first_min.out -ref complex.rst7

注意：因为对原子进行约束，所以，我们需要提供约束位置的参考信息，做法就是将坐标信息(*.rst7)文件接在关键词-ref后

03.编写第二次最小化的02min.in文件：

Minimization of the entire molecular system
 &cntrl
  imin=1, maxcyc=1000,dt=0.002，ncyc=500,
  ntpr=5,
  ntr=1,  #该关键字告知需要进行约束
  ntr=1,restraintmask=":1-292&@CA,N,C,O",  restraint_wt=10.0,
 &end

第二次约束复合物的主链原子。跑完后查看out文件里能量是否有下降。

04.看看能不能继续跑：

export CUDA_VISIBLE_DEVICES=3
pmemd -O -i min.in -p complex.parm7 -c complex_first_min.rst -r \
complex_second_min.rst -o complex_second_min.out

用第一次跑的文件作为启动坐标进行模拟，跑完后查看out文件的能量信息。

05.编写升温文件01heat.in，进行升温-恒容模拟：

Heating up the system equilibration stage 1
 &cntrl
  nstlim=5000, dt=0.002,  "进行10ps模拟"
  ntx=1, irest=0,  #irst=0 不使用先前的速度；ntx=1 使用坐标
  ntpr=500, ntwr=5000, ntwx=5000, #每500步输出out信息，5000步输出能量信息和轨迹坐标
 
  tempi =0.0, temp0=300.0, 
 
#ntt=1或ntt=3都是设置升温算法：
 " ntt=1, #使用Berendsen coupling algorithm \
  tautp=2.0, #算法里的时间常量（tautp）设置为2ps"
 
  ntt=3, gamma_ln=2.0, #用Langevin thermostats和对应参数
 
  ig=209858, #随机种子，用于分配初始速度，要是用同一个数字，会得到相同的轨迹（官推ig=-1）
  ntb=1, ntp=0,  #周期性边界，无压力控制（恒容）
  ntc=2, ntf=2, cut=8.0, #SHAKE algorithm限制和H原子形成的共价键,常量设置为2
  nmropt=1,   #1：读取下面的&wt区域参数
 /
 &wt TYPE='TEMP0', istep1=0, istep2=5000,
  value1=0.1, value2=300.0, 
 /
 &wt TYPE='END' 
 /

注意：这个升温很快（就10ps），out文件里只有10step和轨迹里只有1帧 

 使用Langevin thermostats控制温度，进行升温、恒容模拟，让体系收敛。 

pmemd -O -i 01heat.in -p complex.parm7 -c complex_second_min.rst  \
-r complex_first_v.rst -x complex_first_v.crd -o complex_first_v.out

从out文件中提取温度变化信息：

grep TEMP complex_first_v.out | awk '{print $6, $9}' > temp.dat

在参考文献里，研究进行了100ps的恒温-恒容模拟，编写02heat.in文件，进行恒温-恒容(NVT)模拟：

06.编写02heat.in，进行恒温-恒容（NVT）模拟：

Constant volume constant temperature equilibration stage 2
 &cntrl
  imin=0, irest=1, ntx=5,  #irest=1，从rst文件中读取初始坐标和速度
  nstlim=50000, dt=0.002,  #进行100ps模拟
  ntc=2, ntf=2,
  ntt=3, gamma_ln=2.0,
  tempi=300.0, temp0=300.0,
  ntpr=500, ntwx=50000, #得到一个1帧的轨迹和100步的out信息
  ntb=1, ntp=0,  #周期性边界，无压力控制（恒容）
  ntc=2, ntf=2,  #SHAKE algorithm限制和H原子形成的共价键,常量设置为2
  cut=8.0,
 &end

注意：这里使用rst文件里的速度，就不需要用ig设置随机种子。

pmemd -O -i 02heat.in -p complex.parm7 -c complex_first_v.rst \
-r complex_second_tv.rst -x complex_second_tv.crd -o complex_second_tv.out

同样，利用grep命令从out文件里提取TEMP的信息，输入dat文件。

下面和文献一样，进行100ps的恒温-恒压(NPT)模拟，保证水分子的密度达到实验值。

07.编写03equil.in文件，进行恒温-恒压(NPT)模拟：

Constant pressure constant temperature equilibration stage 3
 &cntrl
  imin=0,ntx=5, irest=1,
  nstlim=50000, dt=0.002, #100ps
  ntpr=500, ntwr=500, ntwx=50000, #得到记录100次的out信息
  tempi=300.0,temp0=300.0,  
  ntt=3, gamma_ln=2.0, #控温
  ntb=2, #周期性边界
  ntp=1,   "ntp=1使用Berendsen恒压器进行恒压模拟"
  taup=2.0, "压力松弛时间是2.0ps"
  ntc=2, ntf=2,cut=8.0,
&end

稳妥起见，还是先试试跑个10ps或是5ps再考虑跑100ps。 在vmd里查看配体在受体里的情况。

pmemd -O -i 03equil.in -p complex.parm7 -c complex_second_tv.rst \
-r complex_third_tp.rst -x complex_third_tp.crd -o complex_third_tp.out

 最后用grep提取out里的势能和密度信息，绘制曲线。

在文献里，完成100ps的NVT和NPT模拟后，进行一个长度为60ns的模拟（...），实际只完成了20ns。（怎么跑完1000ns的。。。）

3.进行最后的成品模拟：

编写一个prod.in的文件，内容和03equil.in类似，只是步长要换成如下：

nstlim=2500000  #积分步长dt是2fs,这里250w步，就是5ns
ntpr=5000, ntwx=5000, #一帧的长度是10ps（教程如此设置），5ns共500帧，out信息中输出500项

编写跑produce的脚本：

#!/usr/bash
sander -O -i prod01.in -o prod01.out -p xx.parm7 -c xx.rst7 -r xx2.rst7 \
-x xx.mdcrd
...
#gzip -9 prod*.mdcrd  #对所有的文件压缩成*.tar.gz