本手册已过时, 不再更新. 如果需要最新手册, 请加入下方QQ群.
1.1 计算化学与分子建模
GROMACS是一个用于分子动力学模拟和能量最小化的计算引擎. 分子动力学模拟和能量最小化是计算化学和分子建模领域众多技术中的两种. 计算化学 是计算技术在化学中应用, 涉及范围从分子的量子力学到复杂大分子聚集体的动力学. 分子建模 是用实际的原子模型描述复杂化学体系的一般方法, 其目的是以原子尺度的详细知识为基础, 理解和预测物质的宏观性质. 通常, 分子建模被用于设计新的材料, 因为在这过程中需要对实际体系的物理性质进行准确的预测.
宏观的物理性质可以分为 (a) 静态平衡性质, 如一种酶与其一种抑制剂的结合常数, 系统的平均势能, 或液体的径向分布函数；(b) 动态或非平衡性质, 如液体的粘度, 膜中的扩散过程, 相变的动力学, 反应动力学, 或在晶体中缺陷的动力学. 计算方法的选择取决于要回答的问题以及该方法在当前技术发展水平下得到可靠结果的可行性. 理想情况下, (相对论)含时薛定谔方程能以很高的精度描述分子体系的性质, 但这种 从头算 的方法水平只能处理几个原子的平衡态计算. 因此, 采用近似是必要的. 体系的复杂度越高, 待研究过程的时间跨度越长, 计算所需的近似度就越高. 到某个阶段(远早于人们所希望的), 从头算方法必须借助于 经验 的参数化模型, 或被其完全取代. 当体系的复杂度导致基于原子间相互作用的物理原理的模拟仍然失败时, 分子建模就完全建立在对已知结构和化学数据进行相似性分析的基础之上. QSAR方法(定量结构-活性关系)和许多基于同源蛋白的蛋白质结构预测就属于这类方法.
宏观性质始终是分子体系某种代表性统计系综(平衡或非平衡)的系综平均值. 对分子建模, 这有两个重要的推论:
仅有单一结构的信息是远远不够的, 即使这个结构是全局能量最小点. 为计算宏观性质, 必须产生指定温度下体系的代表性系综. 但对于计算那些与自由能相关的热力学平衡性质, 如相平衡, 结合常数, 溶解度, 分子构象的相对稳定性等等, 这仍然不够. 自由能和热力学势的计算需要对分子模拟技术进行特殊地推广.
原则上, 分子模拟提供了体系结构和运动的原子尺度的细节, 然而, 这些细节常常与那些感兴趣的宏观性质无关. 这为简化对相互作用的描述并对无关细节进行统计平均提供了可能. 统计力学为这些简化提供了理论框架, 并可以使用一系列不同等级的方法, 涉及到的方法从将一组原子视为一个单元, 以少量的集合坐标描述运动, 对溶剂分子利用平均力势和随机动力学进行统计平均[7], 一直到 介观动力学. 介观动力学直接描述密度而不是原子, 并以描述流(作为对热力学梯度的响应)取代描述速度或加速度(作为对力的响应)[8].
有两种方法可用于产生一个具有代表性的平衡系综: (a) 蒙特卡洛(Monte Carlo, MC)模拟 和(b) 分子动力学(Molecular Dynamics, MD)模拟. 对于非平衡系综的产生和动态事件的分析, 只有第二种方法是适用的. 与MD相比, 蒙特卡洛模拟更简单(不需要计算力), 然而, 在一定的计算时间内, MC并不能得到比MD明显更好的统计结果. 因此, MD是更通用的技术方法. 如果初始构型离平衡态非常远, 体系的受力可能过大, 以致MD模拟失败. 在这种情况下, 需要先对体系进行稳健的 能量最小化. 另一个进行能量最小化的原因是去除体系的所有动能: 如果要对动力学模拟过程中几个“快照”进行比较, 能量最小化能够减少结构和势能中的热噪声, 因此比较结果可能更好.
1.2 分子动力学模拟
对含 \(N\) 个相互作用原子的系统, MD模拟求解其牛顿运动方程:
\[m_i \oppar {\boldsymbol r_i} t = \boldsymbol F_i, i=1, \cdots N \tag{1.1}\]
原子所受的力是一个势能函数导数的负值:
$$ \bi F_i = - \opar V {\bi r_i} \tag{1.2} $$
模拟时, 会以很小的时间步长同时求解这些方程, 使系统演化一段时间, 同时小心地使温度和压力维持在设定值, 并以一定的时间间隔将原子坐标写入输出文件中.
| 键的类型 | 振动类型 | 波数(cm-1) | |
|---|---|---|---|
| C-H, O-H, N-H | 伸缩 | 3000-3500 | |
| C＝C, C＝O | 伸缩 | 1700-2000 | |
| HOH | 弯曲 | 1600 | |
| C-C | 伸缩 | 1400-1600 | |
| H2CX | 剪式, 摇摆 | 1000-1500 | |
| CCC | 弯曲 | 800-1000 | |
| O-H...O | 摆动 | 400-700 | |
| O-H...O | 伸缩 | 50-200 | |
随着时间变化的坐标代表了体系的一个 轨迹. 经过初始的一些变化之后, 体系通常会达到 平衡态. 通过对平衡后的轨迹进行平均, 就可以从输出文件中提取出许多宏观性质.
此时, 了解MD模拟的局限性是有益的. 用户应该清楚这些局限性, 并始终通过对已知实验性质进行检查来评估模拟的准确度. 下面列出了MD所采用的近似.
模拟是经典的
使用牛顿运动方程自然意味着我们使用 经典力学 描述原子的运动. 对处于正常温度下的大多数原子来说, 这是合适的, 但也有例外: 氢原子质量非常小, 质子的运动有时具有显著的量子力学特性. 例如, 在沿氢键传递过程中, 质子可能通过 隧穿效应 越过势垒. 经典力学不可能正确地处理这样的过程！低温下的液氦是经典力学失效的另一个例子. 或许我们不是很关心氦, 但共价键的高频振动却应该让我们担心！当共振频率 \(\n\) 接近或超过 \(k_BT/h\) 时, 经典谐振子与真正的量子谐振子的统计力学明显不同. 室温下当 \(h\n=k_B T\) 时, 波数 \(\s=1/\l=\n/c\), 大约为200 cm-1. 因此, 所有高于, 比如说100cm-1的频率, 在经典模拟中的行为都可能不正确. 这意味着所有键和键角的振动实际上都是可疑的, 甚至氢键键合基团的运动, 比如平动或氢键的摆动都超出了经典的极限(见表1.1). 我们该怎么办呢？
除了进行真正的量子动力学模拟, 我们可以采取下面两种方法中的一种:
(a) 如果进行MD模拟时对键采用了谐振子近似, 我们应该对总内能 \(U = E_{kin} + E_{pot}\) 和比热 \(C_V\) (以及熵 \(S\) 自由能 \(A\) 或 \(G\), 如果要计算它们) 进行校正. 对频率为 \(\n\) 的一维振子, 能量和热容的校正分别为[9]:
\[U^{QM} = U^{cl} + kT \left( {1\over2} x - 1 + {x \over e^x-1} \right) \tag{1.3}\]
\[C_V^{QM} = C_V^{cl} + k \left( {x^2 e^x \over (e^x-1)^2} - 1 \right) \tag{1.4}\]
其中 \(x =h\n/kT\). 当高频量子振子处于基态时, 具有 \({1\over2}h\n\) 的零点能, 经典谐振子会吸收过多的能量(\(kT\)).
(b) 我们可以将键(和键角)视为运动方程的 约束. 其理由是, 与经典振子相比, 处于基态的量子振子与被约束的键更相似. 作这种选择有一个很实用的理由: 当最高频率被移除后, 可以在算法中使用更大的时间步长. 实际中, 相比把键当作谐振子, 当键被约束时的时间步长可以取为原先的4倍[10]. GROMACS对键和键角提供了这种方案. 灵活的键角约束是相当重要的, 它使得运动更加真实, 并能覆盖构型空间[11].
电子处于基态
在MD中, 我们使用 保守力场, 它只是原子位置的函数. 这意味电子的运动被忽略了: 我们假定电子能够瞬间调整自己的运动状态以适应原子位置发生的变化(玻恩-奥本海默近似), 并始终处于它们的基态. 这个假定几乎总是正确的, 但显然不能用于电子转移过程和电子激发态, 也不能正确处理化学反应. 当然, 暂时回避化学反应还有着其他原因.
力场是近似的
力场提供了力. 它们实际上不是模拟方法的一部分. 当需要或对体系有了更多了解时, 用户可以修改力场中的参数. 但是, 在一个特定的程序中所能使用的力的形式是有限制的. GROMACS中可用的力场将在第4章进行说明. 在目前的版本中, 力场是对势累加(pair-additive)的(除了长程库仑力), 所以不能包含极化, 也不包含可以微调的键相互作用. 由此还导致了下面列出的一些限制. 除此以外, 力场对水溶液中的生物相关的大分子是相当有用, 相当可靠的！
力场是对势累加的
这意味着, 所有的 非键 力来源于非键对势相互作用的加和. 非对势累加相互作用, 其中最重要的是通过原子极化产生的相互作用, 要通过 有效对势 进行描述. 这里仅考虑了平均的非对势累加的贡献. 这也意味着, 对相互作用(pair interactions)并不纯粹, 也就是说, 它们对孤立的原子对, 或是与参数化时基于的模型有明显差距的体系, 是无效的. 实际上, 有效对势在实际使用中并没有那么糟糕. 但忽略极化也意味着原子中的电子并不能提供正常的介电常数. 举例来说, 实际的液体烷烃具有略大于2的介电常数, 这会降低(部分)电荷之间的长程静电相互作用. 因此, 模拟将高估远程库仑相互作用. 幸运的是, 下一近似稍稍弥补了这一缺陷.
长程相互作用被截断
在本版本中, GROMACS总是对Lennard-Jones相互作用使用截断半径, 对库伦作用有时也使用截断半径. GROMACS使用的“最小映像约定”要求: 对一个对相互作用(pair interaction)中的每个粒子, 只考虑其在周期性边界条件下的一个映像, 所以截断半径不能超过盒子大小的一半. 对大的体系来说, 这仍然相当大, 但仅仅对于含有带电粒子的体系会有麻烦. 这时可能会发生许多糟糕的事, 如电荷会在截断边界处累积, 或计算的能量完全错误！对于这样的系统, 你应该考虑包含了那些长程静电相互作用的算法, 如粒子网格Ewald (particle mesh Ewald, PME)[12, 13].
边界条件是不自然的
由于体系尺寸很小(即便含10000个粒子的体系也还很小), 粒子组成的团簇与其环境(真空)之间具有许多不需要的边界. 如果想模拟块状体系, 我们必须避免这种情况. 对此, 我们采用周期性边界条件来避免实际的相边界. 由于液体不是晶体, 一些不自然情况仍然存在. 我们最后提及这一方面是因为它的影响最小. 对大体系其误差很小, 但对于具有很强内部空间相关性的小体系, 周期性的边界可能会增强内部的相关性. 在这种情况下, 要多加小心并注意测试系统尺寸的影响. 当使用格点加和方法计算远程静电相互作用时, 这一点尤其重要, 因为我们知道这样有时会引入额外的有序.
1.3 能量最小化和搜索方法
在1.1节已经提到, 许多情况下需要进行能量最小化. GROMACS提供了许多方法进行局部能量最小化, 我们将在3.10节进行详细讨论.
(大)分子体系的势能函数是多维空间中一个超曲面, 具有非常复杂的形貌. 它有一个最低点, 全局极小点 和大量的 局部极小点. 在这些点上, 势能函数对所有坐标的导数为零, 对所有坐标的二阶导数都是非负的. 二阶导数组成的矩阵, 即所谓的 Hessian矩阵, 具有非负的本征值；只有对应于孤立分子的平动和转动的集约坐标具有零本征值. 在局部极小值之间有 鞍点, 其Hessian矩阵仅有一个负的本征值. 通过这些点, 体系可以从一个局部极小值到转移到另一个局部极小值.
拥有所有局部极小点, 包括全局最小点和所有鞍点的信息, 我们就可以描述有关的结构和构象, 它们的自由能, 以及结构转化的动力学. 不幸的是, 构型空间的维数和局部极小点的数目非常多, 我们不可能进行足够的采样来获得全部的信息. 特别是, 没有一种最小化方法能够保证在实际可承受的时间范围内找到全局极小点. 存在一些不很实用的方法, 速度有快有慢[14]. 然而, 给定一个初始构型, 可以找到它的 最近极小点. “最近”在这里并不总是意味着结构意义上的“最近”(即二者坐标差距的平方和最小), 而是指可以通过系统地沿着局部梯度最陡方向往下移动能够到达的极小点. 很不幸, GROMACS所能为你做的只是找到这种最近极小！如果你想找到其他极小点, 并希望在此过程中找到全局极小点, 最好尝试下温度耦合的MD: 先在高温下模拟体系一段时间, 然后慢慢下降到所需的温度, 并不断重复这一过程！如果存在熔点或玻璃化转变温度, 可以先在稍低于该温度的条件下模拟一段时间, 再使用某些聪明的方案慢慢降温, 这样的过程被称为 模拟退火. 由于这一过程不需要与物理真实过程相对应, 你尽可以发挥你的想象力以加速这一进程. 一个经常可行的小技巧是增加氢原子的质量(可增加到10左右): 虽然这将减缓氢原子的快速运动, 却几乎不会影响到体系的慢运动, 同时能使时间步长增加为原来的3倍或4倍. 在搜索过程中你也可以修改势能函数, 例如, 通过移除势垒(去除二面角函数或以 软核 势取代排斥势[15]), 但当慢慢恢复势能函数的时候要小心. 允许结构剧烈变化的最好搜索方法是在四维空间中进行漫游[16], 但这超出了GROMACS的标准功能, 需要一些额外的编程工作.
有三类可能的能量最小化方法:
只需要计算函数值的方法. 如单纯形法及其变种. 前进的每一步都是在前面计算结果的基础上做出的. 如果可以利用导数信息, 这类方法要差于那些使用导数信息的方法.
使用导数信息的方法. 由于在MD程序中, 势能相对于所有的坐标的偏导数(等于力的负值)都是已知的, 由MD程序稍加修改得到的这类方法非常合适.
还使用二阶导数信息的方法. 这些方法在接近极小点时的收敛性非常好: 一个二次势能函数的最小化只需要一步！问题是, 对于N粒子的体系, 必须对 \(3N \times 3N\) 的矩阵进行计算, 存储, 求逆. 除非用额外的编码获得二阶导数, 对大多数系统, 这类方法都不可行. 此外, 还有一些实时建立Hessian矩阵的中间方法, 但它们也需要大量的额外存储空间. 因此, GROMACS没有使用这类方法.
GROMACS使用的 最速下降法 属于第二类方法. 它只是简单地沿负梯度方向(因此, 在力的方向)前进, 而不考虑前面的已有的任何步骤. 搜索过程中, 步长可以进行调整, 以便使搜索快速进行, 但搜索永远沿着能量减小的方向. 这是一个简单稳健但有些愚蠢的方法: 它的收敛速度可能相当慢, 特别是处于局部极小值附近时！收敛更快的 共轭梯度法 (见[17])使用了前面步骤的梯度信息. 一般情况下, 最陡下降法能够很快地接近最近极小点, 而共轭梯度法能够 非常 接近局部极小点, 但当远离极小点时效果较差. GROMACS也支持L-BFGS最小化, 这是几乎与共轭梯度法相当的方法, 但在某些情况下, 其收敛速度更快.
