分子模拟中静电力计算方法的研究.pdfVIP

分子模拟中静电力计算方法的研究

摘要：

分子模拟是现今科学家们探索物质特性的微观机制与内在联系的重要手段，承担着沟

通模型与理论、与实验的关键作用。通过分子力场建模，人们可以设定分子内、分子间的

化学键作用和非键作用。其中被赋予了部分电荷的原子之间的非键且长程的静电作用因收

敛缓慢、有效距离长，而常常占用着大量的计算时间、制约着模拟尺度的扩大。

在静电算法的发展过程中，Ewald3D求和方法首先使精确计算成为可能，而后PME的

提出更使原子个数在百万量级的模拟成为可能，也揭开了分子模拟广泛应用于生物分子体

系的序幕。随着研究内容向多样性发展，静电算法也向着专门化、特异化发展，这里面的

主要原因是体系的维度决定了静电作用的具体形式。因此，相对于百年前的

Ewald3D-tinfoil，晚近才诞生的Ewald2D拥有一套独立的数学表达，但也不该忽视两者之

间物理图像上或数学原理上的相似性，这正是各种近似方法（Ewald3DC、Ewald3DLC等）

具有等价性的基础。

这篇论文的第一个原创性研究工作正是进一步探索Ewald2D的数学基础，并以此严格

地证明了几种近似方法的物理图像的完备。我们先把Ewald2D虚部能量写成傅里叶变换的

形式，利用简单的梯形法则做近似，再借助留数定理将原式与近似式化作复平面积分的形

式，最后通过估计两复平面积分的差值给出一个误差界限，并将此与实际计算误差进行对

比、发现二者始终相当，故从数值的角度证明了此种方法之可行。这个过程其实是将Mori

等人的数值估计与误差处理的经验应用到Ewald2D的分析中。尤其复平面路径积分中的奇

点贡献恰好对应了Ewald3DLC中的多余镜像层的静电作用，因而从精确数值的角度肯定了

后者物理含义的正确。

既然Ewald2D公式天然地分成若干项，我们优化算法也就从对这些分量的相对大小和

计算时间的消耗对比入手。论文中简单的举例说明了一个算法优化的准则，即找到计算中

不怎么影响精度又占用大量计算资源的那些成分、直接约去以显著提高计算效率。这其实

对参数的设置提出了一定的要求，所以操作之前要大致了解算法的原理。

将同样的复平面路径积分方法运用于Ewald1D首先面对的问题却是标准值的计算，即

公式中特殊积分的精确计算。以往依赖于指数积分及其递推关系的数值方法通常产生较大

的误差，这是因为这种递推关系会随阶数的增加而累积误差，所以为准确起见，我们借助

Harris等人对渗漏含水土层函数研究的经验，进而得到可靠的标准值。在Ewald1D虚部的

复平面积分形式中，由于二次积分复杂的误差形式，我们放弃了理论误差的解析表达，但

留数修正项的发现却弥补了前人近似算法的缺失。

通过将特殊积分被积函数泰勒展开，我们探索了一种新的拆分方式，即将无限积分中

因缓慢收敛而难于做梯形法则近似的部分抽出来、求解其解析结果，而余下的部分则采取

简单的近似策略。几个数值的例子说明了这种分解是行之有效的。

可惜的是我们暂时没能完成对Ewald1D算法的进一步数值分析，因而也没能完成

Ewald1D算法的优化问题，但这一题目是可以仿照我们优化Ewald2D的过程进行的，比如

分析Ewald1D各分量在不同体系参数设置下的相对大小和计算耗时对比。

通过对Ewald算法的研究，我们发现在Ewald1D、Ewald2D中起重要作用的虚部项，其

在Ewald3D中对应项反而往往被轻易地省略了、并应用于模拟中。这促使我们把研究目光

重新聚焦在Ewald3D上，反复思考其物理与数学的含义。由于其他非Ewald算法（如MMM、

Leknersum等）在三维周期性条件下的公式往往暗含一种其他的特殊边界条件，这样不利

于我们寻找Ewald3D的真正意义，所以我们采取了另一种原始的方法——以晶胞为单位的

截断法作为对照方法。这其实是Evjen在1932年求解晶体静电能相关问题时采取的策略，

I

依据的是体系静电能的多极展开（即数学上的多维泰勒展开）高阶项贡献衰减迅速、低阶

项可能因体系对称性而消失。

将Ewald3D中发散项泰勒展开，取其条件收敛项，去其高阶项及常数项，我们可以得

到Ewald3D最完整且最简约的形式，只是条件收敛的物理意义并不明确。再将该项写作傅

里叶变