第三章　　算法3.1　简介

在本章中$,$我们首先介绍GROMACS用到的一些基本概念： 周期性边界条件 (3.2节)和 组 (3.3节)。MD算法将在3.4节介绍：首先给出算法的整体形式，然后在后面的小节中进行详细说明。(简单的)EM(Energy minimization, 能量最小化)算法在3.10节说明。用于特殊目的动力学的一些其他算法在此后进行说明。

有几个问题是人们普遍关心的。在所有情况下，必须定义由分子组成的 体系 。而分子又是由定义了相互作用形式的粒子组成。关于分子 拓扑结构 和 力场 的详细描述以及力的计算将在第四章讨论。在本章中，我们只讨论算法的其他方面，如配对列表(pair list)的生成，速度和位置的更新，与外部温度和压力的耦合，约束的保持。对MD模拟数据的 分析 将在第八章讨论。

3.2　周期性边界条件

对有限的体系，减小边缘效应的经典方法是应用 周期性边界条件 。将待模拟体系的原子放到一个空间填充(space-filling)的盒子中，周围是体系自身平移后的副本(图 3.1)。因此，体系不存在边界，由孤立团簇的多余边界导致的问题现在被由周期性条件引起的问题所取代。如果体系具有晶体结构，这样的边界条件符合预期(尽管运动自然被限制于周期性运动，且其波长与盒子适应)。如果希望模拟非周期性体系，如液体或溶液，周期性就会带来误差。此误差可以通过比较不同尺寸体系的模拟结果进行评估。比起由不自然的真空边界引起的误差，周期性边界引起的误差更小一些。

用于空间填充的元胞的有几种可能的形状。其中的一些，如 菱形十二面体 与 截角八面体 [18]比立方体更接近于球形，并因此更适合用于研究溶液中近似球形的大分子，这是由于对给定的大分子映像的最小距离，填充这两种形状的盒子需要的溶剂分子更少。同时，菱形十二面体和截角八面体是 三斜晶胞 的特殊情况，而三斜晶胞是最一般的空间填充元胞，组成了所有可能的空间填充形状[19]。出于这个原因，GROMACS是基于三斜晶胞的。

GROMACS使用周期性边界条件，并结合了 最小映像约定 ：对短程非键相互作用，每个粒子只考虑它的一个最近的映像。对长程静电相互作用，这有时不够准确，因此GROMACS还采用了点阵加和的方法，如Ewald加和，PME和PPPM。

GROMACS支持任何形状的三斜盒子。模拟盒(元胞)由3个盒矢量 a,b,c 定义。盒矢量必须满足以下条件：

通过旋转盒子，总可以使得等式3.1成立，不等式(3.2)和(3.3)也总可以通过增减盒矢量使其成立。

即便使用三斜盒子进行模拟，为了提高计算效率，GROMAC始终将粒子放在长方体空间内，对二维体系的图示见图3.1。因此，从输出的轨迹来看，模拟好像是在长方体盒子中进行的。trjconv程序可用于转换轨迹，以便在不同的元胞中显示。

模拟时也可不使用周期性边界条件，但当模拟孤立的分子团簇时，使用一个很大的周期性盒子通常效率更高，这是因为快速的格点搜索只能在周期性体系中使用。

3.2.1 一些有用的盒子类型

表3.1列出了用于模拟溶液体系的三种最常用的盒子类型。菱形十二面体(图3.2)是最小，最规则的空间填充元胞。12个映像元胞的每一个具有相同的距离，其体积是具有相同映像距离的立方体体积的71%。当模拟溶液中球形或柔性的分子时，这样可以节省约29%的CPU时间。有两种不同取向的菱形十二面体满足方程3.1，3.2和3.3。程序editconf生成与xy平面有正方截面的菱形十二面。之所以选择这个取向，是因为该取向的菱形十二面的前两个盒矢量与x轴和y轴重合，因此容易理解。另一取向可用于膜蛋白的模拟。这时它在xy平面的截面是一个六边形，其面积比具有相同的映像距离的正方形的面积小14%。为获得最佳间距，可以改变盒子的高度(cz)。这种盒子形状不仅可以节省CPU时间，也可以使蛋白质的分布更加均匀。

3.2.2 截断限制

最小映像约定意味着，用于截断非键相互作用的截断半径不能超过最短盒矢量的一半：

否则将会有多个映像出现在力的截断距离内。当研究溶液中的大分子——如蛋白质——时，仅考虑这个限制是不够的：原则上，单个的溶剂分子应该无法同时“看到”大分子的两侧。这意味着，每个盒矢量的长度都必须超过大分子在该方向的长度 再加上 两倍的截断半径 Rc 。但人们经常不严格遵守这个限制，以便稍微减少溶剂层从而节约计算成本。出于效率的考虑，三斜盒子的截断更受限制。对格点搜索的额外限制弱一些：

对简单搜索的额外限制更强：

每个元胞(立方，长方或三斜)都被26个平移的映像所包围。因此一个特定的映像始终可以用指向27个 平移矢量 之一的索引来指认，并通过对索引矢量进行平移来构建(见3.4.3)。约束式(3.5)保证了只需要考虑26个映像。

3.3　组的概念

GROMACS的MD和分析程序可对用户自定义的原子 组 进行一些操作。组的最大数目为256，但每个原子最多只能属于六类不同的组。这六类组如下：

温度耦合组 对每个温度耦合组可单独定义温度耦合参数(参考温度，时间常数，自由度数目，参见3.4.4)。例如，在高分子溶液中，相比高分子，溶剂分子(力和积分的误差使其容易产生更多热效应)可以使用更短的时间常数与热浴耦合，又如，可以将表面的温度维持得低于比吸附分子的温度。 可以定义许多不同的温度耦合组。也请参看下面质心组。

冻结组 属于冻结组的原子在模拟过程中始终保持静止。在平衡体系的过程中这是很有用的，例如，可避免不当放置的溶剂分子对蛋白质的原子产生不合理的碰撞，尽管对受保护原子添加约束势可以得到同样的效果。如果需要，冻结选项可仅仅用于原子的一个或两个坐标，从而将原子冻结在一个平面或一条直线上。当一个原子被部分冻结时，它所受约束仍然可以使它移动，即便是在冻结方向上。一个被完全冻结的原子不能被它所受的约束移动。可定义许多冻结组。冻结坐标不受压力缩放影响；在某些情况下，这可能会产生无用的结果，尤其是与约束同用时(在这种情况下，你会得到非常大的压力)。因此，我们建议避免将冻结组与约束和压力耦合混合使用。为了平衡体系，可以先在使用冻结进行等体积模拟时，然后再使用位置约束与恒压模拟。

加速组 加速组的每个原子上会被加上一个加速度 ag ，这等同于受到一个外力。利用这个特性可驱使体系进入非平衡态，并进行非平衡MD模拟以计算输运性质。

能量监测组 在模拟中，所有能量监测组之间的交叉相互作用都会被考虑，且对Lennard-Jones项和库仑项的计算是分开进行的。原则上，最多可定义256组，但这样将会计算256×256项相互作用！最好保守地使用.

能量监测组之间的所有非键相互作用都可以被排除在外(参见7.3节)。在能量监测组中被排除的粒子对不会被放入配对列表，当不需要计算体系中的某些相互作用时，这样做可以显著提高模拟速度。

质心组 GROMACS可以移除质心(center of mass, COM) 的运动，无论是整个体系的质心还是或原子组的质心。后者(移除原子组的质心)是有用的，比如，对于存在有限的阻止质心运动的摩擦的体系(如气体体系)。对温度耦合与质心运动移除使用相同的组是合理的。

压缩位置输出组 为了进一步减小压缩轨迹文件(.xtc或.tng)的大小，可以仅仅存储一部分粒子的轨迹。所有标记为压缩组的会被保存，其他则不会。如果没有指定这样的输出组，所有的原子坐标都将被保存到压缩轨迹文件。

GROMACS工具中组的使用将在8.1节详述。

◆图片表格公式代码完整版请参看GROMACS中文手册第3章：算法 1-2节◆