分子动力学模拟的主要步骤

本文来源：“分子动力学”公众号

原文链接：https://mp.weixin.qq.com/s/htX-p_gTb6TkdyoGuenjTw

1、确定起始构型

一个能量较低、结构合理的起始三维构型是进行MD模拟的第一步。生物大分子的起始构型主要来自实验数据，如X射线晶体衍射法或核磁共振波谱法测定的分子结构，或根据已知分子结构通过同源建模得到的结构等。体系中有机小分子的构型也可通过量子化学计算得到。

2、选用适当力场和模拟软件

选择适当的力场是进行MD模拟的基础。不同的力场具有不同的适用范围和局限性，需要根据所研究的体系和问题适当选取。同时使用多种力场时，应当注意所用力场间的兼用性。力场的选择与最终模拟结果的准确性息息相关。软件的选择往往与所用的力场有关，应着重考虑所需的算法、软件的运行速度和并行计算能力。

3、构建体系和能量最小化

已经有了研究对象分子的起始构型，接着要根据研究对象所处的环境（如气相、水溶液或跨膜环境等）构建模拟体系，在生物大分子周围加上足够的溶剂分子。体系的大小和模拟元胞的形状依具体体系而定，需要兼顾合理性和可行性。

初步建立的体系中常常存在局部不合理性（如相邻原子间隔太近），不能马上进行动力学模拟，需要先进行体系能量最小化（energy minimization）。比较常用的能量最小化方法有最徒下降法和共轭梯度法。

最徒下降法是一次求导法，利用当前位置的导数作为直线搜索的方向，进行直线搜索求势能面的极小值。这种方法的优点是优化速度比较快，能迅速调整扭曲分子的起始结构；缺点是在极小值附近收敛性较差。最徒下降法常常与共轭梯度法结合使用。

共轭梯度法也是一种一次求导方法，与最徒下降法不同的是，它的直线搜索运用的是前一步的结果，所以能较好地预测下一步搜索的方向，收敛性比最徒下降法好。但共轭梯度法对于不够合理的初始结构，则搜索速度较慢。因此，一般进行能量最小化时，先用最徒下降法迅速优化分子结构，再继续用共轭梯度法优化，直至收敛为止。

4、平衡过程

体系构建好之后要赋予各个原子初始速度，这一速度是根据一定温度下的玻尔兹曼分布随机生成的，然后对各个原子的运动速度进行调整，使得体系总体在各个方向上的动能之和为零，即保证体系没有平动位移。通常在低温下生成初始速度（避免初始温度过大，引起原子碰撞、体系不稳定等），然后在NVT条件下约束住溶质进行逐渐升温，以防止升温过程中损坏溶质的合理构型。升至所需的温度后，接着在NPT条件下进行模拟，调整体系的压强和密度。这一过程中需要对体系的能量、温度、压强、密度等进行监控，看是否收敛，直至体系达到平衡。

5、数据采集过程

体系达到平衡之后则可以进行长时间的模拟，从这个过程中采集样本进行分析。采样要在系统平衡后进行，可记三条轨迹，即体系中粒子的坐标、速度和能量随时间的变化。能量中又包括了不同能量项，以及体系的温度和压力等随时间的变化，如静电能、范德华能，体系的动能、势能和总能量。记轨迹的频率可以根据具体的模拟体系、模拟步长及所研究的现象和性质来选择，并兼顾模拟时间的长短和输出文件的存储需求。如果记轨迹不够频繁，容易缺少信息，不能很好地观测一些相关时间较短的现象；记轨迹过于频繁则会生成巨大的轨迹文件，给后期处理和数据备份造成不必要的麻烦，并且可能超出可用的存储空间，导致模拟中断。模拟时间要尽量长一些，以确定所研究的现象或性质能够被观测到，并且需要确保此现象出现的可重复性。

6、MD结果分析

结果分析主要是通过系综平均得到可与实验结果相比较的宏观物理量。数据分析可用的工具有多种，常用MD软件包内包含一些结果分析程序，也可以根据需要自编一些结果分析软件。以下介绍几种常用的结果分析方法。

1）平均能量及其涨落随模拟时间的变化

2）体系平均原子坐标及其涨落随模拟时间的变化

3）计算体系的结合自由能

4）由几何判据分析分子内及分子间的氢键结构及其氢键占有率

5）蛋白质二级结构及其随时间变化

6）计算体系的动力学性质

文章摘录自王存新等《蛋白质模拟：原理、发展和应用》