在两个平衡阶段都已完成后，体系已经在一定的温度和压力下平衡得很好了。我们现在可以释放位置限定，执行MD来收集数据。我们会在grompp中用到一个检查点文件（在此例中包含恒压耦合信息）。我们会运行一个10ns的MD，用到的.mdp文件在此处（www.mdtutorials.com/gmx/complex/Files/md.mdp）。

正如任何在周期性边界条件下进行的模拟一样，分子可能会被盒子截断或者反复穿过盒子。为了将蛋白重新置于盒子中央并且将晶胞恢复成菱形十二面体，使用trjconv模块：

(资料图)

根据提示centering选择“Protein”，输出选择“System”。注意，将复合物（蛋白-配体，蛋白-蛋白）置于中心对于更长的模拟可能是困难的，会有很多在周期性边界反复横跳的情况。对于那样的情况（尤其是在蛋白-蛋白复合物中），根据蛋白的活性位点或者复合物相互作用界面的残基，创建一个单独的组用于将分子置于中心是必要的。

若想提取轨迹的第一帧（t = 0 ns），对置于中心的轨迹使用trjconv -dump：

对于更加平滑的可视化，最好要进行旋转和平移的调整。执行trjconv如下：

选择“Backbone”对蛋白骨架进行最小二乘拟合，选择“System”作为输出。注意，同时进行周期性边界条件恢复和坐标拟合在数学上是不相容的。如果想要进行拟合的话（对于可视化很有用，但对于大多数的分析是不必要的），分别执行这些坐标调整。

本教程不可能涵盖你想要做的所有分析方法，一些基础的分析在此进行说明。

2-丙基苯酚配体能够与Gln102侧链产生氢键。GROMACS的hbond模块可以容易地计算出在任意两组原子之间形成的氢键数量，但在我们这个例子中，这个值只有1或0。为了能够更加仔细地观察配体与Gln102间的相互作用，我们可以计算JZ4的羟基与Gln102羰基氧原子之间的距离。对于判断是否形成氢键，通常的标准是供体与受体原子之间的距离小于等于3.5A（0.35nm）。用distance模块可以计算整个轨迹中的距离，gmx中的selection语法选择原子。

平均距离是0.31±0.05 nm，符合氢键形成的条件。

第二个通常用于判断氢键是否存在的标准是供体、氢原子和受体原子之间形成的键角。计算键角有不同的惯例。在GROMACS的hbond模块中，键角被定义为氢原子-供体-受体的夹角，这个角应当小于等于30°。若要进行这一分析，首先要为供体原子(必须一同包括供体重原子和连接的氢原子)和受体原子创建编号组。

用angle模块分析这三个原子之间的键角：

注意到angle模块没有-s声明，这和其他大多GROMACS分析模块不同。键角计算不需要质量或周期性信息，原子名称等等。所以轨迹可以在没有.tpr文件或坐标文件的情况下被处理。这个命令返回的结果大约为23±9°。这个结果有些意外，因为我们是按照OAB,H12,OE1的顺序创建编号组的，这对应着供体-氢原子-受体的顺序，我们期待的结果是接近180度的。打开index文件你就会发现：

make_ndx自动地将原子序数从低到高排列，因此计算的结果是受体-供体-氢原子之间的键角，与hbond模块本应计算的角是一致的。所以结果符合氢键形成的条件，因为确实小于30度。为了获得供体-氢原子-受体之间的键角，我们就要人工地用文本编辑器打开编号组，重新对原子序号进行排序（2624 2636 1616）。用这个编号组计算键角得到的结果是147±11°。

最后，我们可能会感兴趣如何定量地表达配体结合模式在整个模拟过程中的变化。创建一个新的编号组，计算JZ4重原子的RMSD值：

执行rms模块，选择“Backbone”进行最小二乘拟合，选择"JZ4_Heavy"进行RMSD计算。这样的话，蛋白整体的旋转和平移通过拟合被去除了，报告的RMSD值是JZ4相对于蛋白结构的变化，这可以很好地说明结合模式在整个模拟过程中是否能保持。

计算的RMSD值应该小于0.1nm（1A），说明配体位置只有很小的变化。

为了定量地描述JZ4和T4溶菌酶之间相互作用的强度，计算二者之间的非键相互作用作用能。GROMACS能够将任意组之间的短程非键能分解。值得注意的是这个这个值并不是自由能或结合能。事实上，大多数力场的这个值都是没有实际物理意义的。CHARMM能够描述和水之间的量子力学相互作用能，所以这个得到的结合能是有意义的。

结合能的计算通过.mdp文件中的energygrps关键词计算。尽管是.mdp关键词，相互作用能的计算不应被视作常规模拟的一部分。短程能的分解与GPU上的MD运算是不相容的，同时会使运算速度发生不必要的减慢。mdrun模块不需要做这份额外的工作来完成一次有效的模拟。因此，只将计算相互作用能作为你分析中的一部分，而不是你的动力学实验。从包含energygrps = ProteinJZ4定义的.mdp文件出发创建一个新的.tpr文件，例如这个：

接着，根据已有的模拟轨迹加入-rerun选项启动mdrun：

注意用-deffnm读取ie.tpr并且将所有输出文件命名为ie.*。-rerun选项选择你用于重新计算能量的轨迹，-nb cpu告诉mdrun你只想在cpu上跑。正如之前所说，这种计算不能在GPU上完成。重新计算的过程应该很快，只需几分钟。

通过energy模块提取感兴趣的能量项。我们感兴趣的是Coul-SR：Protein-JZ4和LJ-SR：Protein-JZ4。

平均短程Coulombic相互作用能是-20.5 ± 7.4 kJ mol-1，而短程Lennard-Jones相互作用能是 -99.1 ± 7.2 kJ mol-1。我们很想从这两个数据中得出一些结论，但即使是CHARMM也很难给出有意义的结论。没有实验的办法能够验证这些值，所以没法知道它们是否有意义。而总的结合能在此例中是有意义的。值是 -119.6 ± 10.3 kJ mol-1（误差部分根据标准公式计算）