Amber学习---小分子肽段的MD（第一天）

参考资料：1.科学网—AMBER基础教程B0：AMBER分子动力学模拟入门 - 李继存的博文 (sciencenet.cn)

2.Benjamin D. Madej & Ross Walker, An Introduction to Molecular Dynamics Simulations using AMBER 

1 使用wsl（windows的子系统linux），了解几个常用的操作

打开cmd，输入wsl启动ubuntu，终端如图所示

1.1 使用cd命令进入指定文件夹  cd /home/lxy

1.2 使用ls命令列出当前目录中的内容 ls

1.3 使用mkdir命令创建一个名为Tutorial的目录 mkdir tutorial

输入ls命令检查所创建的文件夹

1.4 使用cd命令切换不同的目录 

如cd tutorial，cd .. ，cd ~，这三个命令分别代表进入当前目录下的tutorial目录，返回上一级目录以及返回根目录

2 准备拓扑文件和坐标文件

2.1启动leap

        使用tleap命令启动leap，如图

tleap

        现在用source命令加载蛋白质力场FF19SB

source leaprc.protein.ff19SB

2.2 建立肽段，即氨基酸序列，以 ALA VAL SER PHE为例子创建，在leap中N端残基需要在前面加上N，C端加上C，如NALA，CPHE。有三种方式建立肽或者蛋白质序列。

        1.标准氨基酸残基表示：{ALA VAL SER PHE}

        2.N或C标识：{NALA VAL SER CPHE}

        3.保护基团ACE/NME：{ACE ALA VAL SER PHE NME}

注：1.从pdb文件载入的时候，是以第二种形式。

        2.组氨酸可以以质子化的形式，或以δ或ε位置有氢的中性形式存在。对应的组氨酸残基名为HIP、HID或HIE（不是HIS）。标准的leaprc文件将HIS命名为HIE，因此，如果载入含残基HIS的PDB文件，则将该残基命名为HIE。

        3.AMBER力场可以区分残基半胱氨酸(CYS)和参与二硫键的残基胱氨酸(CYX)。此时，必须使用bond命令定义一对半胱氨酸的二硫键，因为PDB文件不能读取该信息。此外，还需使用loadPdbUsingSeq命令加载PDB文件，在将要创建二硫键的序列中将CYX替换为CYS。

foo = sequence { NALA VAL SER CPHE } 

2.3溶剂化 

        用显式的水分子溶剂化分子，模拟中, 将添加TIP3P模型水分子到系统中。在这种类型的模拟中, 系统具有周期性边界条件。因此，周期性盒子应当足够大, 即肽段周围有足够的水, 使肽段分子不与其自身的周期性映像发生相互作用。

        使用solvatebox命令对系统进行溶剂化.

source leaprc.water.tip3p
solvatebox foo TIP3PBOX 10.0

        其中第一行加载水力场，TIP3PBOX指定溶剂化水盒子的类型. 10.0表示水分子在肽段和周期性盒子之间至少有10 埃的缓冲区。溶剂化肽段后，终端如图所示，包含density，盒子大小等信息。

2.4 保存prmtop和inpcrd输入文件

        保存prmtop和inpcrd文件到当前工作目录. 现在Unit foo包括肽段，水分子和模拟所需的周期性盒子信息. 参数根据FF19SB力场指定。

        使用saveamberparm命令保存prmtop和inpcrd文件

saveamberparm foo prmtop inpcrd

2.5 退出LEaP, 使用quit命令

3 准备AMBER MD的sander输入文件

        首先进行能量最小化，其次在正则系综NVT下，20ps升温至298K，最后在标准状态下，NPT系综下，进行60ps的MD。每2 ps保存一次轨迹并写入输出文件一次，使用Langevin恒温器控制温度, 使用随机种子初始化随机数发生器。要控制这些设置, 需要使用文本编辑器编写一个简单的sander输入文件. 我们将使用ubuntu上的gedit文本编辑器。

3.1能量最小化

        使用gedit打开文本编辑器，输入以下内容，并将文件命名为01_min.in。

 Minimize
 &cntrl
  imin=1,
  ntx=1,
  irest=0,
  maxcyc=2000,
  ncyc=1000,
  ntpr=100,
  ntwx=0,
  cut=8.0,
 /
 

  设置含义如下：

• imin=1: 选择运行能量最小化

• ntx=1: 从ASCII格式的inpcrd坐标文件读取坐标, 但不读取速度

• irest=0: 不重启模拟

• maxcyc=2000: 最小化的最大循环数

• ncyc=1000: 最初的0到ncyc循环使用最速下降算法, 此后的ncyc到maxcyc循环切换到共轭梯度算法

• ntpr=100: 每ntpr次循环写入Amber mdout输出文件一次

• ntwx=0: 不输出Amber mdcrd轨迹文件，能量最小化不需要

• cut=8.0: 以埃为单位的非键截断距离( 不要使用低于8.0的值. 较高的数字略微提高精度, 但是大大增加计算成本)

3.2 NVT升温

          使用gedit打开文本编辑器，输入以下内容，并将文件命名为02_heat.in。

Heat
 &cntrl
  imin=0,
  ntx=1,
  irest=0,
  nstlim=10000,
  dt=0.002,
  ntf=2,
  ntc=2,
  tempi=0.0,
  temp0=298.0,
  ntpr=100,
  ntwx=1000,
  cut=8.0,
  ntb=1,
  ntp=0,
  ntt=3,
  gamma_ln=2.0,
  nmropt=1,
  ig=-1,
 /
 &wt type='TEMP0', istep1=0, istep2=9000, value1=0.0, value2=298.0,/
 &wt type='TEMP0', istep1=9001, istep2=10000, value1=298.0, value2=298.0,/
 &wt type='END',/

设置含义如下:

• imin=0: 选择运行分子动力学(MD)，不运行能量最小化

• nstlim=10000: 要运行的MD步数(运行时间长度为nstlim*dt, 单位ps)

• dt=0.002: 以皮秒(ps)为单位的时间步长. 每一MD步骤的时间长度

• ntf=2: 不计算受SHAKE约束的键所受的力

• ntc=2: 启用SHAKE来约束所有包含氢的键

• tempi=0.0: 初始恒温器的温度, 单位K

• temp0=298.0: 最终恒温器的温度 单位K

• ntwx=1000: 每ntwx步输出Amber轨迹文件mdcrd一次，也就是每2ps输出一次

• ntb=1: 等容的周期性边界

• ntp=0: 无压力控制

• ntt=3: 使用Langevin恒温器控制温度

• gamma_ln=2.0: Langevin恒温器的碰撞频率

• nmropt=1: 读入NMR限制和权重变化

• ig=-1: 随机化伪随机数发生器的种子

• 通过NMROPT控制的恒温器温度：上面输入文件的最后三行允许恒温器在整个模拟过程中改变其目标温度. 对于前9000步, 温度将从0 K增加到290 K. 对于9001至10000步, 温度将保持在298 K.

        注意：本MD的NTPR和NTWX设置得非常低, 这样才可以对这个很短的模拟进行分析. 使用这样的设置进行更长时间的MD模拟会产生非常大的输出文件和轨迹文件, 并且比常规MD设置更慢. 对于真正的成品MD, 需要增加NTPR和NTWX的值。 

3.3 NPT生产

         使用gedit打开文本编辑器，输入以下内容，并将文件命名为03_prod.in。

Production
 &cntrl
  imin=0,
  ntx=5,
  irest=1,
  nstlim=30000,
  dt=0.002,
  ntf=2,
  ntc=2,
  temp0=290.0,
  ntpr=100,
  ntwx=100,
  cut=8.0,
  ntb=2,
  ntp=1,
  ntt=3,
  gamma_ln=2.0,
  ig=-1,
 /

• ntx=5: 从无格式的inpcrd坐标文件中读取坐标和速度

• irest=1: 重新启动以前的MD运行(这意味着inpcrd文件中存在速度, 并将使用它们来提供初始原子速度)

• temp0=298.0: 恒温器温度. 在298 K运行

• ntb=2: 在恒定压力下使用周期性边界条件

• ntp=1: 使用Berendsen恒压器进行恒压模

        至此，sander输入文件准备完成，包括参数和拓扑文件prmtop, 坐标文件inpcrd和输入文件01_min.in, 02_heat.in, 03_prod.in, 接下来我们准备运行实际的最小化, 升温和成品MD. 

4 运行Amber MD模拟程序sander

4.1 运行肽段的能量最小化

        使用程序sander, Amber的通用MD引擎(也有一个高性能的版本, 称为pmemd,, 是MD引擎的最佳选择，但是是商业版)。在命令行上, 我们可以指定更多的选项, 并选择使用哪个文件用于输入。

$AMBERHOME/bin/sander -O -i 01_min.in -o 01_min.out -p prmtop -c inpcrd -r 01_min.rst -inf 01_min.mdinfo

sander对MD模拟的每一步都使用一致的语法. 以下是sander命令行选项的总结:

• -O: 覆盖输出文件, 如果它们已经存在

• -i 01_min.in: 选择输入文件

• -o 01_min.out: 输出文件

• -p prmtop: 选择参数和拓扑文件prmtop

• -c inpcrd: 选择坐标文件inpcrd

• -r 01_min.rst: 输出包含坐标和速度的重启文件

• -inf 01_min.mdinfo: 输出包含模拟状态的MD信息文件

sander运行完成后, 应该有一个输出文件01_min.out, 一个重启文件01_min.rst和一个MD信息文件01_min.mdinfo. 将使用重启文件01_min.rst来升温系统.

       使用gedit打开输出文件01_min.out，可以看到系统能量逐步降低

4.2 运行升温MD

        使用sander运行升温MD

$AMBERHOME/bin/sander -O -i 02_heat.in -o 02_heat.out -p prmtop -c 01_min.rst -r 02_heat.rst -x 02_heat.mdcrd -inf 02_heat.mdinfo

以下是sander命令行选项的总结:

• -c 01_min.rst: 从最小化的重启文件获取输入坐标

• -x 02_heat.mdcrd: MD模拟的输出轨迹文件

        使用gegit打开输出文件02_heat.out查看。在02_heat.out文件中, 能找到升温MD的输出，能够看到系统的信息, 包括每步的能量和温度. 例如在1000步的时候:

   

一些重要的数值包括:

• NSTEP: MD模拟的时间步

• TIME: 模拟的总时间(包括重新启动)

• TEMP: 系统温度

• PRESS: 系统压力

• Etot: 系统的总能量

• EKtot: 系统的总动能

• EPtot: 系统的总势能

此时发现并未升温到298步，故将时间改为80ps，并去除nmropt限制条件，再次查看最后一步的能量，温度等数据。

4.3运行生产NPT MD

        使用sander运行NPT MD，

$AMBERHOME/bin/sander -O -i 03_prod.in -o 03_prod.out -p prmtop -c 02_heat.rst -r 03_prod.rst -x 03_prod.mdcrd -inf 03_prod.info &

注意: 命令的末尾加上了&, 使sander在后台运行

        如果我们想监控NPT MD的状态，通过sander输出文件来检查运行状态。Linux程序tail可以输出文件的结尾部分。推出tail，使用ctrl+c

tail -f 03_prod.out

         当sander正在运行时, 上面的命令可以显示输出文件, 这对追踪过程很有帮助。 还可以监控mdinfo文件(cat 03_prod.info), 该文件提供了详细的性能数据以及预计的完成时间。

5 NPT MD模拟输出

        完成后, 打开输出文件检查模拟是否正常完成。本次学习到此为止，明日再学VMD软件可视化MD及MD后处理数据分析。