电荷的“艺术”：分子动力学力场中电荷拟合方法

大家好，我是刘建川，14 年深耕 GROMACS 分子模拟领域，始终专注分子模拟实战。我将拆解分子模拟核心原理，帮你理清技术脉络、避开选择误区，真正做到学以致用。

在分子动力学模拟的世界里，力场是决定模拟成败的核心“宪法”。而在力场的众多参数中，原子电荷无疑是最敏感、也最具争议性的存在。如果说键长、键角决定了分子的骨架，那么原子电荷则决定了分子如何“感知”世界——它控制着静电相互作用，主导着溶质与溶剂、配体与受体之间的识别与结合。

然而，原子在分子中并非孤立的离子，电荷是一种量子化学概念，无法直接通过实验测量。那么，力场中的这些数值究竟从何而来？不同的力场为何对电荷有不同的“偏好”？本文将带你揭秘各种力场背后的电荷拟合方法及其原理。

一、为什么要拟合电荷？——静电作用的基石

在分子力场中，总势能通常表示为键合项（键长、键角、二面角）与非键合项（范德华、静电）之和。对于静电部分，经典力场通常采用库仑定律描述。

其中，qi和qj就是原子中心点电荷。这套模型假设电子云是球形对称且集中在原子核上的，这显然是近似处理。因此，如何选取的值，使其能最大程度重现分子周围的真实静电势（ESP），就成了电荷拟合的核心任务。

二、主流力场的电荷拟合方法

1、AMBER、GAFF力场

AMBER、GAFF力场使用的是RESP (Restrained ElectroStatic Potential) 电荷拟合方法。RESP区别于早期静电势拟合方法（如CHELPG、MK）的关键。如果没有约束，最小二乘法有时会为了极小化静电势误差，给非极性原子（尤其是氢原子）分配出极不合理的电荷值（例如让甲基上的氢带0.3以上的电荷）。这在物理上不符合化学直觉，也会在后续模拟中产生不真实的相互作用。

RESP方法的核心是“两步走”策略。首先，通过高精度的量子化学计算（如HF/6-31G*）得到分子在真空中的波函数，计算出分子周围三维空间格点上的静电势（ESP）。然后，通过最小二乘法调整原子中心点电荷，使这些点电荷产生的静电势与量子化学计算的静电势尽可能一致。在这个过程中，RESP引入了一个关键的“约束项”（restraint），特别是对非极性原子（如氢原子）的电荷施加限制，防止其拟合出不符合化学直觉的过大数值，从而保证了电荷的化学合理性。为了解决分子构象依赖性问题，RESP通常会拟合分子多个低能构象的平均电荷。RESP能极好地重现分子在气相中的静电性质，对氢键相互作用描述准确。

一个分子可以有多种低能构象，不同构象的电子云分布不同，拟合出的电荷也会有所差异。如果只对单一构象拟合电荷，当分子在模拟中发生构象变化时，该电荷模型可能不再准确。RESP的常用策略是多构象拟合：选取分子在气相中的多个低能构象，分别进行量子化学计算，然后将所有构象的静电势格点合并成一个大的数据集，统一进行一次拟合。最终得到的电荷是一个平均意义上的最优值，能够在多种构象下都保持较好的静电势重现能力。

2、OPLS-AA力场

OPLS (Optimized Potentials for Liquid Simulations) 力场的电荷方法则有所不同。OPLS力场的电荷来源主要有：

1） 键偶极叠加 OPLS力场的一个核心特点是，其原子部分电荷并非直接拟合而来，而是根据力场参数文件中定义的键偶极计算得到的。对于一个键，力场文件会存储一个偶极矩（D，单位Debye）和自然键长（L0，单位Å）。该键上分配的电荷Q通过以下公式计算：Q = 0.2082 × (D / L0)。一个原子上的总电荷，就是它所有成键贡献的总和，并通过形式电荷进行修正。这种方法确保了分子整体的电中性，并具有良好的可传递性。

2） 用于小分子的高级方法：1.14CM1A/CM5 对于小分子，OPLS-AA力场的更常采用一种名为 1.14CM1A或CM5的电荷模型。这是一种半经验方法，通过将Mulliken等布居分析电荷与经验校正因子结合，快速、准确地预测原子电荷，特别适合药物分子等有机小分子的高通量模拟。该电荷方法通过键偶极定义电荷，参数化过程高效，模型具有良好的可传递性，能很好地重现有机液体的凝聚相性质（如密度、汽化热）。

3） 此外，通常还用1.2*CM5电荷，源于Truhlar等人发展的CM5模型，通过将原始电荷值统一乘以1.2倍，在数值上增强了原子电荷的强度，以此等效地模拟溶剂环境对溶质产生的极化作用。

CM5 (Charge Model 5) 是Truhlar课题组于2012年提出的一种第四类（Class IV）原子电荷模型。它的目标是解决传统布居分析电荷（如Mulliken）精度不足的问题，同时避免静电势拟合电荷存在的构象依赖性和“埋藏原子”不稳定性。CM5描述的是分子在真空（气相）中的电荷分布。当分子处于溶液中时，极性溶剂分子（如水）的电场会诱导溶质分子的电子云发生重新分布，使原子电荷的幅值增大——这就是极化效应。1.2CM5正是通过整体放大原子电荷的数值，来等效地描述这种溶剂环境对溶质的极化作用。它相当于在固定电荷模型的框架内，用“增强的静电势”来隐式地包含溶剂极化对溶质的影响。缩放因子1.2并非随意选取，而是基于2017年发表在J. Phys. Chem. B.上的系统测试。1.2CM5电荷的本质是一种兼顾精度与效率的经验性修正：它用CM5保证了气相静电性质的准确性，再用1.2倍缩放以固定电荷模型的方式隐式地包含了溶剂极化，从而为OPLS-AA力场提供了一套可靠、普适的原子电荷参数。

后续想学习GROMACS哪方面的实操内容，欢迎大家在留言区告诉我！