本文面向Tanglab及类似的希望针对生物体系做分子动力学模拟的科研人士
本文在一定程度上参考了Amber Tutorials ，本系列教程旨在提供一些作者的使用经验。若想更为全面地掌握Amber分子动力学模拟软件，应当阅读官方教程。
本系列教程所需的预备知识（prerequisite）包括掌握Linux的基本命令，以及对分子动力学模拟理论的大致了解。
起这个标题名字的时候，感觉就像是在写实验报告捏
本系列教程从第一篇起，开始真正地上手计算。在本节中，我们使用Amber自带的建模程序tleap构建一条短肽并溶剂化。通过一段时间的分子动力学模拟观察模拟过程中的轨迹，并对模拟结果做简要的分析。通过本节的模拟，笔者带大家Amber进行Amber的初体验，力图以实操的形式阐明分子动力学模拟的基本流程。
搭建模拟体系永远是分子动力学模拟中最为重要的步骤。
正如本系列教程第0章所说的，在模拟之肇始，Amber依靠两个文件，prmtop与inpcrd[1]，来识别一个分子。prmtop文件包含了原子名称，原子类型，残基名称，力场参数等信息，inpcrd包含了原子的所有坐标，以及模拟盒子的形状。
倘若已在系统中成功安装了AmberTools，那么直接在命令行窗口输入tleap[2]
1 | COPY $ tleap |
那么你应当看到如下字样：
1 | COPY -I: Adding /smb/open-apps/amber22/amber22/dat/leap/prep to search path. |
且命令提示符由 $ 变为 > ，这意味着tleap已被成功打开。path因Amber的安装地址而异。
分子动力学模拟当然少不了力场。本节中我们要模拟一条肽链，因此选择使用专用于蛋白质的力场，ff系列。具体为ff19SB。关于力场本身，包括其加载方式，还有许多有趣的subtle。此处暂且按下不表。
1 | COPY > source leaprc.protein.ff19SB |
你应当见到如下所示的输出：
1 | COPY ----- Source: /smb/open-apps/amber22/amber22/dat/leap/cmd/leaprc.protein.ff19SB |
source即加载力场命令，leaprc.protein.ff19SB为力场名称。如果想查看source的详细用法，可以参考amber reference manual，也可以在tleap中执行。适用于其他命令。
1 | COPY > help source |
tips: tleap是大小写不敏感的（case-insensitive）。故“help source” 也可写成 “HELP source”，“help SOURCE”，“HelP SoURce”等形式。
接下来我们使用tleap创建一条肽链。命名该段肽链为polyppd，并将其保存为pdb文件。
1 | COPY > polyppd = sequence { ACE ALA MET LYS LEU GLU NME } |
tips： Amber中分子模型分为三级，自下而上分别是原子（atom）、残基（res/residue）、单元（unit）。在本例中polyppd即一个unit，其包含了七个res，每个res又含有若干个原子。
打开pdb文件，看一下我们的多肽链长啥样
创建的多肽链
相关的文件会在每节的最后给出：
接下来我们给多肽添加一个厚度为10 $\AA$ 的立方体水溶剂盒子，使用TIP3P水模型。
1 | COPY > source leaprc.water.tip3p |
最后一行说明该溶剂化命令添加了1404个水分子。
我们先检查一下搭建的模型有无参数缺失
1 | COPY > check polyppd |
众所周知，warning是向来不用管的 $~$ Warning提示我们有两对原子距离太近了。这是很trivial的。tleap只是帮我们把残基连起来，它并不会考虑结构的合理性。后面做分子动力学模拟时我们会优化结构。
保存文件
1 | COPY > saveAmberParm polyppd prmtop inpcrd |
命令提示符变为熟悉的$符号，这代表我们回到了Linux界面。看一下我们现在得到了什么。
1 | COPY $ ls |
leap.log记录了你在tleap中执行的操作。其余三个文件为我们在tleap中生成的。
以下是得到的溶剂化后的Amber输入文件：
建模部分顺利结束。接下来是正式的模拟。一般说来，一个完整的模拟流程包括最小化minimization、升温heating、成品production。根据模拟目的不同，单个流程中可以包含多段模拟。
我们需要准备参数文件，来告诉Amber我们希望做怎样的模拟。这里使用vim文本编辑器。
1 | COPY $ vim min.in |
tips: 强烈建议使用Linux端的文本编辑器。windows端的文本编辑软件（尤其是word）默认使用的换行符，及其他相关的格式有可能导致程序报错。如使用windows端文本编辑软件，推荐使用notepad++或其他文本编辑器。
但是如果没有vim使用经验的话，请先查阅相关教程，掌握基本的vim命令。
写入以下内容
1 | COPY Minimize ! title line. type anything you want |
tips: Amebr是基于Fortran语言编写的，在输入文件中，半角叹号 ! 引导注释内容。其不会被软件读取。
在参数输入文件中，第一行是标题，尽情发挥。模拟参数由namelist引导。在学习初期，我们只涉及 &cntrl 这一个namelist。每个namelist由 / 作为结束。行首空格只是排版好看。可顶格写，也可一行写入多个参数。
imin=1 $~~~~$ 单点能计算
ntx=1 $~~~~$ 从输入文件中读取坐标信息。不读取速度信息。
irest=0 $~~~~$ 不要重启模拟 (重启模拟不适用极小化)
maxcyc=2000 $~~~~$ 做2000步的极小化模拟
ntmin=1 $~~~~$ 在一开始的ncyc步使用最速下降法，然后转换到共轭梯度法
ncyc=1000 $~~~~$ 使用最速下降法做1000步的极小化
ntpr=100 $~~~~$ 每100步输出一次能量信息到输出文件
ntwx=0 $~~~~$ 不输出坐标信息
cut=8.0 $~~~~$ 非键相互作用的截断距离
极小化结束后体系还是0K的。贸然做成品MD，有可能导致体系的崩溃。故我们还需要升温。
将以下内容写入heat.in
1 | COPY Heat |
imin=0 $~~~~$ 运行分子动力学模拟，而不是最小化
nstlim=10000 $~~~~$ MD运行步数 (nstlim * dt = 模拟时长（ps）)
dt=0.002 $~~~~$ 每步的演化时长，亦即相邻两步的时间间隔。单位为皮秒，ps
ntf=2 $~~~~$ 不计算SHAKE算法所约束的化学键的受力
ntc=2 $~~~~$ 使用SHAKE约束所有含氢原子的化学键
tempi=0.0 $~~~~$ 初始的恒温器温度。单位：K
temp0=300.0 $~~~~$ 最终的恒温器温度。单位：K
ntwx=1000 $~~~~$ 每1000步将体系的坐标信息输出到Amber轨迹文件mdcrd
ntb=1 $~~~~$ 固定模拟盒子体积的周期性边界条件
ntp=0 $~~~~$ 不做压强约束
ntt=3 $~~~~$ 使用郎之万（Langevin）恒温器做温度控制
gamma_ln=2.0 $~~~~$ 郎之万恒温器碰撞频率
nmropt=1 $~~~~$ 使用核磁限制信息与读取重量变化
ig=-1 $~~~~$ 郎之万恒温器中产生随机数的随机种子。除非调试程序，强烈建议设置ig=-1
heat.in最后三行代表体系将用9000步从0K升温至300K，并在从9001步到10000步这一段维持300K的温度。
做最终的成品MD，我们选用NPT系综来模拟真实体系。对于不同的模拟目的，应谨慎选取所使用的系综。
将以下内容写入prod.in
1 | COPY Production |
ntx=5 $~~~~$ 从rst7坐标输入文件中读取坐标与速度信息。
irest=1 $~~~~$ 重启先前的MD模拟。这意味着rst7坐标输入文件中应当有速度信息，并被程序读取作为原子的初始速度
temp0=300.0 $~~~~$ 恒温器初始温度。
ntwv=500 $~~~~$ 每500步将速度信息输出到Amber速度输出文件mdvel
ntb=2 $~~~~$ 使用恒压周期性边界条件
ntp=1 $~~~~$ 使用贝伦德森（Berendsen）恒压器调节体系压强
贵组服务器资源充足，可爽用GPU/CPU。
以CPU单核运算为例。通过以下命令调用sander进行模拟计算
1 | COPY $ sander -O -i min.in -o min.out -p prmtop -c inpcrd -r min.rst |
对上命令的各参数做出说明
-O $~~~~$ Overwrite the output files if they already exist
-i min.in $~~~~$ Choose input file min.in (default mdin)
-o min.out $~~~~$ Write output file min.out (default mdout)
-p prmtop $~~~~$ Choose input parameter and topology file prmtop
-c inpcrd $~~~~$ Choose input coordinate file inpcrd
-r min.rst $~~~~$ Write output restart file with coordinates and velocities (default restrt)
-inf 0min.info $~~~~$ Write MD info file with simulation status (default mdinfo)
-x prod.crd $~~~~$ Write output trajectory file prod.crd (default mdcrd)
-v prod.vel $~~~~$ Write output velocity file prod.vel (default mdvel)
简而言之，一个模拟（无论是最小化，升温还是成品）需要三个输入文件，分别是MD参数文件mdin、分子参数与拓扑文件prmtop、分子坐标文件inpcrd/restrt。而其输出包括输出文件mdout、重启文件restrt[3]、轨迹坐标文件mdcrd、轨迹速度文件mdvel。
在运行sander（pmemd同理）的命令中，若不指定输入文件名，则程序会根据默认名称搜索当前目录下有无对应文件。若不指定输出文件名，则程序以默认文件名进行输出。
restrt、mdcrd、mdvel文件默认是二进制的。
做多核并行运算则使用sander.MPI程序。例如使用96个CPU做并行运算：
1 | COPY $ mpirun -np 96 sander.MPI -O -i min.in -o min.out -p prmtop -c inpcrd -r min.rst |
tips: 若做多核运算，使得模拟计算并行效率最大的做法是为模拟任务分配的CPU数目为2的n次幂个，如2，4，8，16，32，等。当然在计算资源充足的情况下总是核数越多，计算越快。例如总的而言96核总是快于64核。此外有极个别体系[4]，倘若不分配 $2^n$ 个CPU，将导致程序报错。此规则在模拟时不需要刻意遵守，遇到所谓的“极个别”体系时另行修正即可。
使用pmemd.cuda程序做GPU加速运算。
1 | COPY $ nvidia-smi # find out if any GPU available |
不愧是GPU，一颗4090足以薄纱96个CPU。
然而遗憾的是有很多体系不支持GPU加速。在这类体系的mdout报错信息中会有相应提示。
让我们看看得到了哪些输出文件（此处没有列出全部输出文件）吧！
min.out $~~~~$
min.rst $~~~~$
heat.out $~~~~$
heat.rst $~~~~$
prod.out $~~~~$
prod.rst $~~~~$
prod.crd
以prod.out文件的片段为例，说明输出文件中含有哪些信息。
1 | COPY NSTEP = 16500 TIME(PS) = 53.000 TEMP(K) = 299.55 PRESS = -1.6 |
目前先介绍几个基本信息
NSTEP $~~~~$ 步数，即体系在当前过程中运算到了哪一步。
TIME $~~~~$ 模拟所经历的时间。注意倘若输入文件中irest=1，则该时间将从上一段模拟中开始累加。
TEMP $~~~~$ 温度。
PRESS $~~~~$ 压强。注意其输出的是相对波动值，单位0.000001 bar[5]。
Etot、EKtot、EPtot $~~~~$ 体系的总能量、动能、势能。单位kcal/mol。
Density $~~~~$ 密度。单位$g/cm^3$。
此外在mdout文件（prod.out）末尾也会包含体系性质的平均值与涨落。
1 | COPY A V E R A G E S O V E R 100000 S T E P S |
使用amber自带的分析脚本简单看一下成品MD跑得怎么样。
1 | COPY $ mkdir analysis |
每个性质都输出了summary_avg.xxx, summary.xxx, summary_rms.xxx，代表了对应性质的平均值、列表汇总、标准差。既然成品MD是在NPT系综下进行的，我们关注体系的体积、能量变化曲线，以判断体系是否达到平衡。
模拟盒子体积曲线
对体积而言，其在21ps[6]至约50ps不断下降，此后保持稳定。这是由于系统自动的溶剂化时添加的水分子往往排列太过“稀疏”，使得模拟盒子体积过大。那么在MD模拟时就会有一个收缩的过程。倘若我们需要获知体系的平衡性质，那么显然应当只取50ps之后的轨迹。
|
|
|
对于动能曲线，可见其变化大致在200kcal/mol的范围内波动，即恒温器较好地控制了体系的温度。而势能曲线与总能量曲线在约30ps处出现了一个明显的能量下降过程，对应了模拟盒子的收缩。而在50ps后，势能曲线与总能量曲线没有再发生大的变化，也是在大致在200kcal/mol的范围内波动，由此可以认为体系最终确实到了比较平衡、稳定的状态。当然可能有人认为这很难讲，尤其是看总能量变化的时候，波动似乎还是挺大的。固然200ps的模拟时长确实略显不足。不过这只是一个演示实验，大体已经可以说明原理了。说白了，只是从性质变化曲线上判断体系平衡与否，并没有什么很统一的标准。你能说服自己即可。
VMD软件是一款非常强大的可视化软件，支持多种文件格式，可用于分子动力学模拟轨迹的可视化及其他一些后处理工作。遗憾的是，尽管VMD支持Windows与Linux两个系统，但Windows版本暂时无法识别amber输出文件。一个较为曲折的解决方案是安装WSL（windows subsystem for Linux），然后安装子系统的虚拟桌面，远程连接虚拟桌面后使用VMD。好在贵组的新服务器是有图形化界面的，并且也已部署好VMD，也就省去了诸多麻烦。
关于如何使用VMD，请参考 Using VMD with AMBER 。这篇教程写得非常详尽，笔者不再班门弄斧。这里直接列出模拟结果，不再说明VMD如何操作。
值得一提的是在VMD中，prmtop格式选择Amber 7 Parm，rst文件（包括mdrst、mdcrd）格式为NetCDF，rst7文件对应格式Amber Restart。
成品MD轨迹
使用VMD制作动图的步骤在上面列出的参考教程中亦有提及，此处不再赘述。这里直观展示一下得到的模拟轨迹。可以看到模拟初期确实有盒子（蓝色边框）缩小的过程。
多肽链最终构象
成品MD轨迹最后一帧（Frame）对应的分子构象。与刚建模时的构象相比，变化还挺大的。我们这里只是目测一下，不做进一步的定量讨论。顺带一提将rst文件转换为pdb文件的方法：
1 | COPY $ cpptraj -p prmtop -y prod.rst -x prod.rst7 |
这样的话得到的最终构象也可以通过GaussianView等软件打开了。
简而言之本教程就图一乐，展示了MD模拟的基本流程。最终我们得到了一段轨迹，以及多肽在溶液中的稳定构象。此外还通过分析体系的体积、能量变化曲线，判断其是否达到平衡状态。由于没有明确的探究问题，因此最终也没有什么robust的结论。希望读者能够掌握如何使用tleap建模，通过调节参数控制模拟条件，并最终跑一段MD轨迹。
[1] 这两个默认文件名的含义：prmtop 即 parameter & topology， inpcrd 即 input coordinate。也可自定义文件名。
[2] 本系列教程约定：除非额外说明，命令提示符为 $ 时，代表命令直接在Linux终端中输入。命令提示符为 > 时，代表其为在打开的软件中执行的命令。没有命令提示符则代表为命令的输出内容或文件内容（总之不是命令就对啦）。
[3] 重启文件restrt实时记录、更新模拟过程中轨迹的最后一帧(frame)的空间坐标。其有两个作用。第一是模拟因为意外中断时，可以以restrt文件作为坐标输入文件恢复模拟，从而节省时间；第二是上一步模拟的重启文件是下一步模拟的坐标输入文件。
[4] 在笔者近一年的使用经历中，只遇到过一个这种体系，即对非键相互作用做出修正的12-6-4模型。
[5] 存疑。待考据。
[6] 回顾成品MD的模拟中，我们设置了参数irest=1，即重启MD，因此计时也是承接升温过程的。升温模拟时长20ps，因此成品MD轨迹对应的时间段为21ps-220ps。