在本文中，我将再次深入探讨化学领域。 许多软件包，包括商业的和开源的，都可以用于进行量子水平的化学计算。 我在这里介绍的是 gromacs (http://www.gromacs.org)。 它应该可以通过您的发行版的软件包管理器获得。

gromacs 软件包实际上是一整套小型应用程序，用于处理从创建初始化文件到执行计算运行，再到分析和可视化结果的所有步骤。 您可能还会注意到有多个版本可用。 例如，Ubuntu 既有单处理器版本，也有 MPI 版本，用于在并行计算机上进行计算。 如果您拥有多核机器或机器集群，这可以大大加快您的计算速度。

在开始之前，先简单介绍一下可以进行的计算类型以及 gromacs 使用的方法。 理想情况下，在计算化学中，您应该能够获取所研究系统的薛定谔方程并完全求解它们。 但是，除了只有几个原子的非常简单的系统之外，这样做很快就会变得不可能。 这意味着需要某种近似方法才能获得结果。 在 gromacs 中，近似方法是分子动力学 (MD)。 MD 模拟求解一组相互作用粒子的牛顿方程，以获得它们的轨迹。 计算来自所有其他粒子的粒子上的总力，然后除以粒子的质量以获得其加速度。 这是针对所有粒子计算的，每个粒子都根据其加速度移动。 然后，时间步进一个单位，并再次计算整个过程。 您可以提高空间和时间分辨率，但会以更多的计算机时间为代价，直到获得足够准确的结果以满足您的目的。

MD 模拟有几个局限性。 首先，它是一个经典模拟。 在大多数情况下，这很好，但是在某些情况下，您需要了解系统中原子的量子效应。 第二个局限性是电子处于基态，并且被视为瞬间移动以保持在其原子周围的轨道中。 这意味着您无法模拟化学反应或任何其他电子相互作用。 下一个局限性是长程相互作用被截断。 当处理带电粒子时，这成为一个严重的问题。 我在这里看到的最后一个局限性是，周期性边界条件被用于尝试模拟本体系统。 这与上面提到的截断相结合，意味着您可能会得到一些不物理的结果，尤其是在长时间运行时。

既然您有了一些背景知识，那么您实际上如何使用 gromacs 呢？ 您需要从输入开始：所有原子的初始位置、初始速度以及描述所有原子之间力的相互作用势。 整个系统的质心通常被定义为具有零速度，这意味着系统上没有外力。 输入此初始数据后，gromacs 需要计算力，将这些力应用于每个粒子并计算它们的新位置。 一旦它们移动了，gromacs 需要重新计算力。 这是使用蛙跳算法完成的。 为了更好地了解它的外观，让我们考虑一个具体的例子：水中的蛋白质。

第一步是提出进行计算所需的初始化文件。 这可以完全从头开始完成，但在许多情况下，没有必要。 对于蛋白质，有蛋白质数据库 (PDB)，位于 http://www.pdb.org。 确保蛋白质结构具有您在模拟中寻找的所需细节。 您也可以将其加载到 PyMOL 中并查看其外观。 当您对您的选择感到满意时，您可以使用它通过 pdb2gmx 为 gromacs 生成所需的初始化文件

pdb2gmx -f my_protein.pdb -water tip3p

其中 tip3p 是您可以选择的水模型之一。 此命令生成多个输出文件，其中最重要的是 conf.gro、topol.top 和 posre.itp。 此时，您仍然没有将水添加到初始化文件中。 为此，您首先需要定义一个盒子来容纳水和蛋白质。 为此，您可以使用以下命令编辑配置文件

editconf -f conf.gro -bt dodecahedron -d 0.5 -o box.gro

这定义了一个形状为十二面体且直径为 0.5nm 的盒子。 您也可以使用立方体或八面体形状的盒子。 现在您有了一个盒子，您可以使用以下命令添加水

genbox -cp box.gro -cs spc216.gro -p topol.top -o solvated.gro

此命令获取盒子 (box.gro) 并用文件 spc216.gro 中定义的水分子填充它。 -p topol.top 选项将此水添加到系统的拓扑中。 最后，所有这些都写入到文件 solvated.gro 中。

如果您现在尝试运行它，您可能会遇到由于引入的水分子引起的巨大力而导致的问题。 为了解决这个问题，您可以在实际进行任何计算之前最小化系统的能量。 您可以通过创建一个参数文件来定义如何进行最小化。 例如

------em.mdp------
integrator   = steep
nsteps       = 200
nstlist      = 10
rlist        = 1.0
coulombtype  = pme
rcoulomb     = 1.0
vdw-type     = cut-off
rvdw         = 1.0
nstenergy    = 10
------------------

在此示例中，通过最速下降法进行最小化，超过 200 步。 您可以在 gromacs 文档中查找所有其他选项的详细信息。 完成此操作后，您可以使用 grompp 命令完成所有必要的预处理

grompp -f em.mdp -p topol.top -c solvated.gro -o em.tpr

实际的最小化是通过以下方式处理的

mdrun -v -deffnm em

前缀 em 用于所有相关文件名，具有不同的扩展名。 这使得在您的桌面上完成所有预处理并完成初始化步骤变得更容易，然后在超级计算机上进行实际运行。

当您准备好进行最终运行时，您需要设置一个参数文件来描述详细信息。 在此示例中，您可以使用类似这样的内容

------run.mdp------
integrator   = md
nsteps       = 5000
dt           = 0.002
nstlist      = 10
rlist        = 1.0
coulombtype  = pme
rcoulomb     = 1.0
vdw-type     = cut-off
rvdw         = 1.0
tcoupl       = Berendsen
tc-grps      = protein non-protein
tau-t        = 0.1 0.1
ref-t        = 298 298
nstxout      = 1000
nstvout      = 1000
nstxtcout    = 100
nstenergy    = 100
------------------

使用此参数文件，您可以使用以下命令进行预处理

grompp -f run.mdp -p topol.top -c pr.gro -o run.tpr

实际的 MD 计算使用以下命令完成

mdrun -v -deffnm run

现在您可以去喝几天咖啡了。 实际运行时间将取决于您有多少处理器可以投入到这个问题中。

一旦运行完成，您仍然需要分析它以查看您可以从中学习到什么。 您可以将结果与来自 X 射线结构测量的实验结果进行比较。 您可以使用 g_rms 程序测量重原子与 X 射线结构的位移。 您可以使用 g_dist 和 g_hbond 程序分析距离和氢键。 您甚至可以使用 trjconv 程序制作电影

trjconv -s run.gro -f run.xtc -e 2500.0 -o movie.pdb

这将导出 2.5 纳秒的轨迹电影。 然后您可以使用 PyMOL 查看它。

这篇短文仅提供了 gromacs 功能的冰山一角。 我希望它能激发您阅读相关资料并进行一些新的和有趣的工作的兴趣。

化学图像 通过 Shutterstock.com。