近年来随着分子动力学模拟方法在化学、生物及相关学科研究中的地位逐渐凸显，其得到越来越多的关注。自1977年Karplus等人首次完成牛胰蛋白酶抑制剂的分子动力学模拟以来^[1]，Web of Science数据库中题目或摘要含有分子动力学模拟(Molecular Dynamics Simulation)的文章数目飞速增长。尤其是2000年后随着计算能力的提升，分子模拟相关研究呈直线上升趋势(图 1)。同时，我们也关注到近年来美国化学教育杂志(Journal of Chemical Education)也报道了大量在物理化学或其他基础化学课程中引入分子动力学模拟实验的文章报道^[2–8]。例如，Chelsea Sweet等人^[4]在物理化学实验中设计了单原子分子体系的分子动力学模拟实验，以氩气为例采用分子动力学方法计算其在不同条件下的行为及宏观性质，并比较与理想气体的差别。Laura J. Kinnaman等人^[5]在物理化学实验中引入分子动力学模拟方法，采用不同水模型对纯水体系进行分子动力学计算，考察了不同水模型计算出的结构性质和动力学性质(径向分布函数、扩散系数)与实验值的差别。Logan H. Eckler等人^[6]在物理化学实验中采用分子动力学方法对正己烷在不同温度条件下进行模拟，并考察其扩散系数和粘度的变化。设计这些实验的初衷一方面是给学生介绍这一有力工具，另一方面是帮助学生理解一些抽象的化学概念，例如分子间作用力、溶液的微观结构等等。

目前国内鲜有分子动力学模拟相关的实验教学案例报道，笔者了解的有武汉大学侯华教授等^[9]在“分子模拟实验”中设计的对泛素分子在真空中进行的分子动力学模拟实验，考察生物大分子的空间结构随模拟时间的变化；南开大学孙宏伟教授^[10]在“量子化学与分子力学”课程中以多肽分子为例简要介绍了分子动力学模拟的基本原理和步骤。已发表的大多数计算化学实验均是基于量子化学方法研究分子的电子结构性质、化学反应等^[11–20]。例如厦门大学袁汝明等人^[19]设计的H3反应势能面的构建、气相分子的标准摩尔生成焓的预测、CO在Pd(111)表面的氧化实验等；清华大学王溢磊等人^{[13, 19]}开设的“计算化学实验”中的诸多项目；南开大学许秀芳等人^[14–16]在“计算化学”课程中设计的系列实验等。山东大学自2015年开始给化学专业大三下学期本科生开设分子模拟实验课程^[20–23]，通过设计一系列典型的分子模拟实验帮助学生从分子水平上理解化学物质的结构-性质关系、动力学性质和反应特性等，让学生了解如何从微观角度描述宏观现象，并养成通过计算化学的思维方式来解决化学问题的能力。课程内容涵盖了量子化学方法和分子力学方法。本文我们对其中的一个最基本的单分子及团簇的分子动力学模拟实验进行介绍。

分子动力学模拟主要研究体系中所有粒子的运动状态随时间的演变，在一定的统计力学系综下通过对相空间进行系综平均(时间平均)，获得体系的物理性质和化学性质。从理论上讲，分子动力学模拟才应该算是真正的“计算机实验”。分子动力学模拟的基本流程包括如下的步骤：

1) 为所研究问题建立合适的“模型体系”(modeling)；

2) 确立模型体系的初始状态：包括粒子的坐标参数、速度分布、环境状态等(initialization)；

3) 为模型体系赋予力场参数(force field)；

4) 求解牛顿运动方程，直到体系处于平衡状态(equilibrium)；

5) 从平衡态出发，继续求解运动方程，记录粒子的运动轨迹(production)；

6) 对平衡后的轨迹进行统计分析(analyze)。

本实验我们以丙氨酸模型二肽分子为例，简单说明分子动力学模拟的基本原理和步骤，并初步探讨温度、溶剂等外部因素对分子构象的影响。

1) 了解分子动力学的基本原理；

2) 能够区分比较量子化学计算与分子动力学模拟的区别、优势和适用范围；

3) 掌握分子动力学计算的一般步骤，并能够进行单分子及团簇的分子动力学模拟计算；

4) 掌握并运用简单的分子动力学模拟结果分析方法。

计算机，Materials Studio软件(Visualizer、Forcite模块)。

乙酰基和甲胺基封端的丙氨酸二肽(N-acetyl-N’-methyl-L-alanylamide，NANMA)是一个常用来研究蛋白质中氨基酸的构象行为的典型模型分子^{[7, 24]}，其分子结构如图 2所示。首先采用Visualizer模块构建其三维结构。点击菜单栏File-New...在弹出的对话框中选择3D Atomistic。首先需要画出分子中所有骨架原子(即非氢原子)，在工具栏点击按钮，此时光标会变成铅笔状(右侧下拉条可选O原子)，在*.xsd文件上单击并移动即可画出多个原子，按ESC键终止。非氢原子完成后单击工具栏完成加氢，再单击其附近的Clean功能调整优化分子结构。完成后可将文件名重命名为NANMA.xsd。

通常在运行分子动力学计算之前需要对初始模型进行能量极小化(Minimization)，即几何优化(Geometry Optimization)，以防止人为搭建的原子模型间有不合理的接触。在上一步构建分子模型时已采用Visualizer模块的Clean功能进行了初步的构型优化，故可以直接进行分子动力学的模拟。

分子动力学计算的流程可以概括为：①根据麦克斯韦分布给初始构型随机分配速度；②根据力场参数和原子坐标计算各原子的受力；③根据牛顿运动方程更新下一步原子的位置和速度；④应用控温控压算法调节体系的温度压力；⑤记录轨迹等信息；⑥循环②–⑤直至达到预定的模拟时间。

具体操作如下：打开Forcite计算模块，在Setup标签页将任务设置为Dynamics即动力学计算，再点击右侧More...按钮设置具体模拟参数，如图 3所示，在NVT系综下计算(粒子数、温度、体积恒定)，随机分配初速度，温度300 K，动力学计算积分步长1 fs (通常积分步长应小于体系中振动频率最快的运动模式的十分之一，氢原子由于质量最小，所以通常和氢原子相连的键振动频率最快，周期约为10 fs，故积分步长不应大于1 fs)，总共模拟10 ns，每间隔10000步记录一次轨迹信息，其中控温算法选择Berendsen方法。

在Energy标签页将力场设置为COMPASS (该力场由上海交通大学孙淮等人于90年代开发，普适性较强)，对于非键相互作用(包括静电相互作用和范德华相互作用)的加和方式均采用Atom based加和法。Job Control标签页可以设置并行计算的核数和任务描述，应采用尽可能多的CPU核数并行以加速计算，如图 4所示。将构型文件NANMA*xsd窗口激活，点击Run即可提交计算，计算进程会在下方任务栏显示。

1) 动力学轨迹的可视化分析。

动力学计算完成后会生成一条记录原子每一步坐标、速度和受力信息的轨迹文件NANMA.xtd，通过Animation工具栏可以对体系在动力学计算过程中的变化进行可视化分析，调整好视角之后亦可以导出为视频文件(File-Export-*.avi)。

2) 骨架二面角随时间的变化和分布。

在获得轨迹的基础上，我们可以对体系中感兴趣的量进行定量统计分析。以二肽分子中的骨架二面角phi和psi (图 2)随时间的变化和分布为例：选择工具栏中角度和距离测量工具，标记出感兴趣的二面角phi和psi，此时该二面角的度数会以蓝色显示，单击选中该二面角后呈黄色高亮显示，打开工具栏Edit-Edit Sets-New…对话框，将二面角分别命名为phi和psi，如图 5所示。

采用Forcite模块的Analysis功能，选择Torsion evolution，Sets分别选择上一步中定义的phi和psi，激活轨迹文件NANMA.xtd后点击Analyze进行计算，输出结果为* Torsion Evolution.xtd和* Torsion Evolution.xcd，分别为数据文件和绘图文件。类似的方法选择Torsion distribution功能可计算两个骨架二面角的分布，结果如图 6所示。

参考4.2中的方法可进一步对丙氨酸二肽模型分子在500 K温度下的构象变化进行模拟，统计骨架二面角phi和psi随时间的变化，考察其与300 K温度下的异同。显然可以发现在500 K温度下，分子的构象变化得更为剧烈，两个二面角的分布也更加宽泛。这也说明升高模拟温度可以帮助分子跨越一些常温条件下难以跨越的势垒，更加充分地对势能面进行采样。

注：分子力学方法采用谐振式描述原子间的成键相互作用，不必担心高温引起原子间化学键的断裂。高温情况下，分子运动较为剧烈，为了确保分子动力学的稳定，步长建议不超过1 fs。

通过对两个温度下进行模拟可以发现在低温时分子倾向于在极小点附近振动，在高温时分子则能轻易地跨越势能面上的势垒。因此我们可以通过高温分子动力学结合能量极小化的办法找到势能面上的一批极小点结构(local minimum)，在分子动力学中这种方法称为淬火。

具体操作如下(图 7)：打开Forcite计算模块，在Setup标签页将任务设置为Quench即淬火，再点击右侧More...按钮设置具体淬火参数，每10000步淬火一次，在NVT系综下计算，随机分配初速度，模拟温度500 K，动力学计算积分步长1 fs，共模拟1 ns，控温算法选择Berendsen方法。

结果文件中*.xtd文件为动力学轨迹文件，* Quench.xtd文件为经能量极小化后的轨迹文件，*Quench.std文件为淬火动力学得到的分子极小点结构和其他信息汇总表。在此基础上可进行一系列的性质分析。

通过淬火动力学方法我们能够找到一批极小点的结构，但是我们更关心的是分子的全局最小点结构(global minimum)，显然在高温动力学的基础上，我们缓慢地降温，则分子极有可能收敛到该结构，如果周期性地循环升温-降温则可以找到势能面上的一批极小点结构，这批结构中则有可能存在分子的全局最小点结构。这种周期性升温-降温的方法称作模拟退火动力学方法。

具体操作如下(图 8)：打开Forcite计算模块，在Setup标签页将任务设置为Anneal即周期性模拟退火，再点击右侧More...按钮设置具体退火参数，共进行5轮退火，初始温度1 K，最高温度500 K，升温或降温阶梯数500 (为了确保能收敛到极小点，降温过程应尽可能慢，升温过程可以较快)，每个温度阶梯模拟的动力学步数1000步(相当于每1 ps升温/降温1 K)，在NVT系综下计算，随机分配初速度，动力学计算积分步长1 fs，控温算法选择Berendsen方法。

其中* Anneal.xtd文件为5轮退火得到的分子结构组成的轨迹，* Anneal.std为5轮退火得到的分子结构和其他信息汇总表。由于分子力学相对较为粗糙，为了得到精确的极小点结构，可以以模拟退火得到的结构为初始构型，采用量子化学方法进行进一步的几何优化，并计算能量。

采用类似的方法可对丙氨酸二肽模型分子形成的二聚体或与水分子形成的团簇结构进行分子动力学模拟，并考察其组装结构。

注：原子数增多后需要计算的原子对间的相互作用的数目也会增加，导致计算量增加，达到同样的模拟时间的目标下，为了节约耗时，可以在Setup标签页将模拟步长设置为2 fs，同时勾选上Fix bonds，这样在动力学计算过程中程序会约束住与H原子相连的键的振动，可以使用较大的时间步长。

1) 请分别描述300 K和500 K温度下丙氨酸二肽模型分子轨迹的特点。(提示：300 K温度下分子基本维持在特定结构，500 K温度下分子构型波动较为剧烈。)

2) 请根据分子动力学模拟的结果，以phi和psi为xy坐标，分别绘制出300 K和500 K温度下丙氨酸二肽模型分子的势能面，并对比二者的区别。(注：能量可由Forcite分析模块Potential energy components获得。)

3) 请基于分子力学，利用系统格点搜索法(Conformers模块)对丙氨酸二肽模型分子的两个骨架二面角同时做构象搜索，并绘制出势能面，与上一题中的势能面相比有何不同，并推测原因。(提示：系统格点搜索法能够得到完整的势能面，而分子动力学方法无法跨越所有的能垒，难以在有限的时间内遍历势能面。)

4) 请编写Perl脚本采用量子化学半经验方法(VAMP模块)对两个骨架二面角做势能面扫描，并与分子力学方法得到的势能面比较。

5) 请采用量子化学方法对周期性退火得到的极小点结构做几何优化，并计算能量，并分析分子为何倾向于采用该构象。(提示：形成分子内氢键，能量降低。)

6) G. N. Ramachandran发现氨基酸的骨架二面角只能取某些有限的范围，后人对大量蛋白质分析并绘制了标准的拉氏图(Ramachandran map) ^[25]，请将丙氨酸二肽模型分子与标准拉氏图对比，分析其可能的二级结构。(提示：beta折叠。)

7) 将中间C原子上相连的甲基替换为H，结果会有何不同。(提示：没有侧链psi和phi的取值范围更广。)

8) 请构建丙氨酸模型二肽分子与水分子的团簇结构，并进行分子动力学模拟，考察在有水分子存在的条件下势能面有何不同。(提示：水分子存在的情况下会阻止分子内氢键的形成。)

9) 请比较并总结量子化学计算与分子动力学模拟的区别、优势、及适用范围。

本实验由课上两个学时完成，实验前主讲教师需对分子动力学方法的基本原理和流程进行简要说明，并进行主要实验流程的演示。实验内容4.1–4.3为基本分子动力学模拟流程，由学生在课上完成；实验内容4.4–4.7可由学生在课后选择性完成，亦可根据课时安排或计算资源在课上选择性完成。

根据山东大学分子模拟实验的课程安排，本实验是学生学习完量子化学实验后(分子模型创建与优化、分子性质计算、势能面构建、热力学量计算、过渡态搜索、分子光谱计算)向分子动力学实验过渡的第一个实验，因此在设计时我们采用了最简单的单分子模型，不涉及周期性边界条件、非键相互作用处理方式等特殊处理方法，而与量子化学计算实验中的孤立体系的计算较为类似，仅计算方法由量子化学方法过渡到了分子力学方法，易于学生类比和迁移。

通过实践我们发现学生能基本熟悉分子动力学方法的基本原理和流程，为后续开展均相体系(溶液)、非均相体系(界面)等复杂体系分子动力学模拟打下了基础。同时由于分子力学方法是基于量子化学方法构造的势能面来描述原子间的相互作用，因此通过本实验学生对于势能面的物理意义以及通过量子化学方法构造出的势能面的应用也有了更深刻的认识。本实验与量子计算化学方法对孤立体系的静态性质研究不同，通过本实验学生第一次直观地看到分子结构的动态变化过程，也算是真正意义上的计算机实验，故均表现出了浓厚的兴趣，也是学生感觉分子模拟实验有趣有用之处。

通过本实验学生能够初步了解分子动力学模拟的基本原理及流程，学会处理和分析模拟得到的数据。同时通过本实验的学习与上机操作，能够将涉及到的物理化学中的势能面等抽象概念理解更为深入，并学以致用。即使对于没有专门开设计算化学实验课程的高校，在物理化学实验中增设此实验亦能开阔学生的视野，并激发学生的科研兴趣。