OPLS-AA使用记录

Background

近年来，全原子分子动力学模拟在聚合物及凝聚态材料研究中得到广泛应用。相比粗粒化模型，全原子模型保留了具体的化学结构，因此能够描述局部构象、官能团相互作用及链段运动。

为了比较不同经典力场的适用性，我们首先考察 OPLS-AA 的参数来源、势能函数和自动参数化方法。

Introduction

OPLS-AA（Optimized Potentials for Liquid Simulations - All Atom）是其中一种面向有机分子和凝聚态体系的经典全原子力场， 最早源于由 William L. Jorgensen 及其合作者在 20 世纪 80 年代发展的 OPLS-UA，，最初用于模拟有机液体的分子间相互作用。早期 OPLS-UA（United Atom）模型将部分氢原子合并到相应的重原子相互作用中心。

1996 年，Jorgensen、Maxwell 和 Tirado-Rives 系统建立了 OPLS-AA 力场。该模型显式描述氢原子，并利用量子化学构象能与液体实验数据确定二面角及非键相互作用参数。此后，OPLS-AA 不断扩展和修订了不同官能团及化学环境的参数。

2015 年，OPLS-AA/M 针对肽和蛋白质的构象能进行了优化，重点修订主链与侧链二面角参数。2019 年，该系列进一步扩展并改善了 RNA 的描述。

2017 年，LigParGen 提供了自动生成 OPLS-AA 参数及分子拓扑的方法，并支持 CM1A-LBCC 原子电荷。LigParGen 的 Python 封装已有公开实现，但核心参数分配依赖 BOSS 非开源程序，拥有诸多限制。

2019 年，PolyParGen 采用分子片段化方法，为较大的分子和聚合物生成力场参数。

此外，Schrödinger 于 2021 年发表了商业力场 OPLS4，针对更广泛的化学空间改进了参数和预测精度。

OPLS-AA 形式

Jorgensen、Maxwell 和 Tirado-Rives 于 1996 年在 Journal of the American Chemical Society 发表了 Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids。

该工作主要借鉴 AMBER 全原子力场的键伸缩和键角参数，使用 RHF/6-31G*//RHF/6-31G* 计算的旋转能曲线拟合二面角参数，并结合实验液体密度和汽化焓等数据优化非键相互作用。

OPLS-AA 的总势能可以写成：

$$U=U_{\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{d}}+U_{\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}}+U_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{r}}+U_{\mathrm{i}\mathrm{m}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{r}}+U_{\mathrm{L}\mathrm{J}}+U_{\mathrm{C}\mathrm{o}\mathrm{u}\mathrm{l}}.$$

其中前四项描述成键相互作用，后两项描述非键相互作用。

Bond stretching

键伸缩采用谐振子近似：

$$U_{\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{d}}=\sum _{(ij)\in B}\frac{k_{ij}}{2}(r_{ij}-r_{ij}^{(0)}{)}^2.$$

其中，$r_{ij}^{(0)}$ 为平衡键长，$k_{ij}$ 为力常数。只考虑了势能在平衡键长附近的二阶展开。属于非常简单的谐振近似。

这里的 $1/2$ 系数属于参数定义约定。不同力场文件可能将该系数吸收到力常数中。

Angle bending

键角采用类似的谐近似：

$$U_{\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}}=\sum _{(ijk)\in A}\frac{k_{ijk}}{2}({\theta }_{ijk}-{\theta }_{ijk}^{(0)}{)}^2.$$

其中，${\theta }_{ijk}^{(0)}$ 为平衡键角。其忽略了较高阶非谐项及键长与键角之间的显式耦合。

Proper torsion

经典 OPLS 二面角势使用截断至第四谐波的傅里叶级数：

$$\begin{array}{rl}{U}_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{r}}(\varphi )=& \frac{V_1}{2}(1+\cos \varphi )\\ & +\frac{V_2}{2}(1-\cos 2\varphi )\\ & +\frac{V_3}{2}(1+\cos 3\varphi )\\ & +\frac{V_4}{2}(1-\cos 4\varphi ).\end{array}$$

其中，$V_1,\dots ,V_4$ 控制不同二面角构象的相对能量和旋转势垒。 与键伸缩和键角弯曲相比，许多单键绕键旋转所需的能量较低，因此能够产生多种热力学可访问的构象。因此对于柔性聚合物，主链二面角的能量分布会显著影响链的局部柔性、构象统计和构象转变速率。

Improper torsion

Improper 二面角主要用于限制局部几何结构，例如维持芳香环或其他平面官能团的平面性。 根据具体力场参数和程序实现，可以采用谐势或周期势。例如谐势形式为：

$$U_{\mathrm{i}\mathrm{m}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{r}}=\sum _{\xi }\frac{k_{\xi }}{2}(\xi -{\xi }_0{)}^2.$$

其中，$\xi$ 是 improper 二面角，${\xi }_0$ 为参考角度。

Nonbonded interactions

OPLS-AA 的非键相互作用由 Lennard-Jones 势和库仑势组成。

Lennard-Jones 势为：

$$U_{\mathrm{L}\mathrm{J}}=\sum _{i<j}{f}_{ij}^{\mathrm{L}\mathrm{J}}4{ϵ}_{ij}[{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^{12}-{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^6].$$

其中，${ϵ}_{ij}$ 为势阱深度，${\sigma }_{ij}$ 为 LJ 势能过零位置对应的距离参数。

传统 OPLS-AA 对分子内相互作用只描述 1-4（即间隔三个共价键）及以上的 LJ ，非相互作用参数使用几何混合规则：

$$\boxed{{\displaystyle {\sigma }_{ij}=\sqrt{{\sigma }_i{\sigma }_j},{ϵ}_{ij}=\sqrt{{ϵ}_i{ϵ}_j}.}}$$

静电相互作用采用固定点电荷模型：

$$U_{\mathrm{C}\mathrm{o}\mathrm{u}\mathrm{l}}=\sum _{i<j}{f}_{ij}^{\mathrm{C}\mathrm{o}\mathrm{u}\mathrm{l}}\frac{q_i{q}_j}{4\pi {ϵ}_0{r}_{ij}}.$$

标准 OPLS-AA 对 1–4 (指通过三个共价键相连) 非键相互作用采用特殊的$1/2$缩放：

$$U_{14}=0.5U_{14}^{\mathrm{L}\mathrm{J}}+0.5U_{14}^{\mathrm{C}\mathrm{o}\mathrm{u}\mathrm{l}}.$$

也就是 1-4 非键相互作用 LJ 和 Coulomb 均需乘 0.5，其与 proper torsion 项共同定义了二面角。

因此，对某个二面角进行参数拟合时，应比较完整的势能变化，而不能只考虑显式的 Fourier 二面角项。特别要注意 1–4 相互作用的缩放系数属于力场定义的一部分。修改这些系数后，通常需要重新检查或拟合相关二面角参数。