分子模拟/长程力处理

在分子模拟中处理长程相互作用在计算需求方面提出了严峻的挑战。通过使用周期性边界条件，可以使用有限单元内的有限数量的分子有效地计算样品的整体性质。单元中的每个分子都会与其他所有分子发生相互作用，包括同一单元内的分子和无限周期性单元图像内的分子。这些相互作用包括色散、排斥以及任何其他静电相互作用。本质上，这意味着在模拟过程中必须计算无数的分子间相互作用，这对实际计算来说是一个问题。

与分子模拟中常见的做法一样，Lennard-Jones 势用于描述色散和排斥相互作用。相对于其他静电相互作用，这些相互作用在一定距离 $r_{}$ 下衰减得相当快，因此在计算中可以相当容易地处理。静电相互作用（例如电荷-偶极、偶极-偶极、四极-四极等）在给定距离 $r_{}$ 下通常要强得多，并且随着 $r_{}$ 的增加以更慢的速度衰减，并且衰减速度随着 $r_{}$ 指数的增加而增加。这使得它们的处理变得更加复杂。

截断

未修改的 Lennard-Jones 势（蓝色）和使用平滑势截断的 Lennard-Jones 势，使用切换函数（红色）。切换函数在灰色区域内处于活动状态。切换势从蓝色区域开始为零。

非键截断用于节省计算资源，同时不牺牲结果的有效性。这些截断用于模拟模拟中原子或分子之间的长程相互作用，并且本质上在特定距离处截断相互作用。因为分子间相互作用势随着 $r\to \infty$ 收敛到特定值（即零），因此相互作用势可以在一定距离后有效地调整到该值。对于色散和排斥相互作用，这些相互作用由 Lennard-Jones 势建模，在模拟中使用切换函数 $S(r)$。此切换函数使相互作用的衰减在所需的截断距离处达到零，此时参与分子间相互作用的两个分子或原子根本不通过 Lennard-Jones 势发生相互作用。

$\nu (r)={\begin{cases}\nu _{LJ}(r)&{\text{ for }}r\leq r_{switch}\\\nu _{LJ}(r)S(r)&{\text{ for }}r_{switch}\leq r\leq r_{cutoff}\\0&{\text{ for }}r>r_{cutoff}\end{cases}}$

切换函数开启的分子间距离是切换距离，通常比截断距离短约 2 Å。在切换距离和截断距离之间，切换函数逐渐将相互作用势缩放到零。这种相互作用势的逐渐变化保持了模拟的物理基础，而不使用相互作用的突然、非物理截断（如硬球模型中的排斥相互作用）。这种方法避免了由于突然截断相互作用势函数而导致的相互作用势的不连续性。

这种对长程相互作用的便捷处理可以轻松应用于由 Lennard-Jones 势建模的色散和排斥相互作用，但并不适用于所有情况。对于静电相互作用，必须采用其他方法（例如 Ewald 求和，见下文）。这是因为静电相互作用通常比色散和排斥相互作用的长程得多，并且通过使用简单的切换函数进行截断会导致结果准确性严重下降。这些相互作用比色散或排斥更强，必须更仔细地处理。

Ewald 分解

虽然使用切换函数和非键截断使 Lennard-Jones 势的色散和排斥相互作用的截断变得相对简单，但处理长程静电相互作用并不那么直接，需要采用不同的方法。这些相互作用的评估通常使用 Ewald 方法进行，该方法将静电相互作用分为短程和长程成分。短程成分在真实空间中迅速衰减到零，因此可以使用切换函数进行截断，就像 Lennard-Jones 势一样。长程成分更难处理，但可以在倒空间中有效地计算。在 Ewald 求和方法中，使用高斯误差函数 (erf) 和互补误差函数 (erfc)。 这些函数可以看作是

Ewald 分解通过使用 erfc(r) 和 erf(r) 恒等式将库仑相互作用分解为短程真实空间成分和长程倒空间成分。真实空间 erfc(r)/r 成分在短距离处收敛到零，避免了截断误差。

${\text{erf}}(x)={\frac {2}{\sqrt {\pi }}}\int \limits _{0}^{x}e^{-t^{2}}dt$

以及

${\text{erfc}}(x)=1-{\frac {2}{\sqrt {\pi }}}\int \limits _{0}^{x}e^{-t^{2}}dt$

其中 ${\text{erf}}(0)=0$ 以及 ${\text{erfc}}(0)=1$。因此，这两个函数互为补函数。此外，请注意

$\lim _{x\to \infty }{\text{erf}}(x)=1$

以及

$\lim _{x\to \infty }{\text{erfc}}(x)=0$

因此，根据定义，

${\textrm {erfc}}(x)+{\textrm {erf}}(x)=1$.

在使用这些高斯函数时，库仑相互作用势的 ${\frac {1}{r}}$ 项被分解为上述短程项和长程项，

${\frac {1}{r}}={\frac {1}{r}}\left[{\textrm {erfc}}(\alpha r)+{\textrm {erf}}(\alpha r)\right]={\frac {{\textrm {erfc}}(\alpha r)}{r}}+{\frac {{\textrm {erf}}(\alpha r)}{r}}$

其中 $\alpha$ 是一个可调常数，用于确保静电相互作用在所需的 $r$ 值处衰减到零。在这种情况下，${\frac {{\textrm {erfc}}(\alpha r)}{r}}$ 项对应于短程相互作用势，它通过在实空间中的评估迅速衰减到零，而 ${\frac {{\textrm {erf}}(\alpha r)}{r}}$ 项对应于长程相互作用势，它更容易在倒空间中进行评估。 该势被评估为模拟盒子的所有可能的周期图像中的所有电荷-电荷相互作用的总和。在倒空间中，这些相互作用迅速收敛，通过快速傅立叶变换的应用允许有效计算。 在模拟系统的所有周期图像中，所有长程电荷-电荷相互作用的总和可以表示为

$\sum _{\mathbf {S} }{\frac {q_{1}q_{2}{\text{erf}}(\alpha |\mathbf {r_{ij}} +\mathbf {S} |)}{|\mathbf {r_{ij}} +\mathbf {S} |}}$

其中 $\mathbf {r_{ij}}$ 是相互作用的电荷 $q_{1}$ 和 $q_{2}$ 之间的距离向量，而 $\mathbf {S}$ 是不同周期性晶胞图像之间的平移向量。等价地，这个和可以写成傅里叶级数，

${\frac {1}{L^{3}}}\sum _{\mathbf {g} }C_{\mathbf {g} }e^{i\mathbf {g} \cdot \mathbf {r} }$

其中 $C_{\mathbf {g} }$ 是倒空间展开系数，而 $\mathbf {g}$ 是倒空间向量。注意，虽然这个和仍然是无限的，但是 $C_{\mathbf {g} }\to 0$ 随着 $\mathbf {g} \to \infty$ 快速趋近于 0，允许通过应用快速傅里叶变换来有效地计算长程静电相互作用。

粒子网格 Ewald 方法

粒子网格 Ewald 方法是对 Ewald 求和的数值近似。 它涉及将电荷“平滑地”分布在一个有限的晶格点集上，这些晶格点位于放置在模拟单元上的网格（或“网格”）上。通过将电荷密度定位在这个网格上，跨越模拟单元及其无限周期图像的整体电荷分布被简化为更容易处理的形式。这使得长程静电相互作用的计算变得明显更容易和更便宜，特别是在模拟单元长度很大，因此 $\mathbf {g}$ 向量数量很多的情况下。从“平滑”电荷密度的晶格，可以使用快速傅里叶变换来计算长程相互作用势。虽然这种近似显著降低了计算成本，但收集数据的准确性仍然很高。粒子网格 Ewald 方法也相对容易集成到模拟中。 这些因素结合在一起，使其成为计算长程静电相互作用势的极其有吸引力的方法。

另见

维基百科关于 Ewald 求和的文章