摘要：
分子模拟是现今科学家们探索物质特性的微观机制与内在联系的重要手段，承担着沟通模型与理论、与实验的关键作用。通过分子力场建模，人们可以设定分子内、分子间的化学键作用和非键作用。其中被赋予了部分电荷的原子之间的非键且长程的静电作用因收敛缓慢、有效距离长，而常常占用着大量的计算时间、制约着模拟尺度的扩大。
在静电算法的发展过程中，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、Lekner sum等）在三维周期性条件下的公式往往暗含一种其他的特殊边界条件，这样不利于我们寻Ewald3D的真正意义，所以我们采取了另一种原始的方法——以晶胞为单位的截断法作为对照方法。这其实是Evjen在1932年求解晶体静电能相关问题时采取的策略，
I
济南丝足依据的是体系静电能的多极展开（即数学上的多维泰勒展开）高阶项贡献衰减迅速、低阶项可能因体系对称性而消失。
将Ewald3D中发散项泰勒展开，取其条件收敛项，去其高阶项及常数项，我们可以得到Ewald3D最完整且最简约的形式，只是条件收敛的物理意义并不明确。再将该项写作傅里叶变换的形式，并利用高斯定理，我们重新得到该项的实空间表述模样。通过几组数值对比各种条件下的Ewald3D和截断法，
我们更直观地解释了Ewald3D的物理图像，并在结论中对“无穷大”的概念做了简单的探讨。
关键词：
分子模拟，静电算法，Ewald求和
太原郝建秀>白城师范学院学报
II
A study on electrostatics algorithms in molecular simulations
Abstract：
Nowadays, more and more scientists choose molecular simulations to investigate the microscopic physics and intrinsic mechanism of materials, because this technology can both check the validity of the model and the corresponding theory(or hypothesis). By molecular force field modeling, people can set up the inter- and intra-molecular chemical bonded or non-bond interactions, in which the long-ranged Coulombic interactions between atoms that are partially charged are always considered as the computational barrier in one simulation for its slowly converging effect.
During the development of algorithms for electrostatics, Ewald3D was firstly established to accompli
sh the accurate calculation, and then PME enabled the simulations that are composed by millions of atoms, which also opened the gate for this tool to be extensively used in biological systems. As the researches become specified and specialized, the electrostatic algorithms also change their expressions as a result of different formula determined by the dimension and periodicity of the system. Thus, compared to Ewald-tinfoil method that existed for a century or so, the much younger Ewald2D method has a totally different representation. However, it should not be neglected that these two methods have some similarities in both physical picture and mathematical formula, which is also the fundamental of the equivalence between different approximations(Ewald3DC, Ewald3DLC etc).
This dissertation illustrates its first innovation on the mathematics of Ewald2D and rigorously demonstrates the validity of several approximations. To achieve these goals, we firstly rewrite the the reciprocal part of Ewald2D energy into a form that is similar to Fourier transform, and approximate the integral by trapezoidal rule. Then, we transform the integral and its quadrature into path integrals on the complex plane by Cautchy theorem. At last, we estimate the error bounds of this approximation by subtracting these two complex integral, and find that the theoretical errors always agree with the realistic computational error which can be viewed as a numerical proof of this approac
h. Especially, the contribution by the singularities in the path integral just corresponds to the electrostatic interactions by excess image layers in Ewald3DLC method, which states that the physics in Ewald3DLC is at least numerically reasonable.
Since Ewald2D formula naturally separate into several terms, we optimize the algorithm by analyzing the relative magnitude and computational complexity of each ingredient. We exemplify to explain one rule in optimizing algorithms, that is to find the parts which nearly harmless to the accuracy but time consuming badly and ignore it completely to strikingly increase the calculation efficiency. Obviously, this strategy depends on how to set up the parameters. Thus it is better to understand the algorithm before optimizing it.
Exerting the same tactic on Ewald1D, we immediately confront the problem of how to calculate the special integral in it. The old formula using exponential
包场中学III
integral and its recurrence relation usually generate large errors. So for accurately evaluating the special integral, we reform it into a sum of leaky aquifer function that has been discussed by Harris and other mathematicians. By the methods suggested by them, we can get the precise Ewald1D ene
rgy. However, the path integral form of the reciprocal energy is a twice integral which is too complex to handle with. Though we cannot provide the analytical expression of the error bounds, the newly found residue corrections still remedy the old method to a better accuracy. After Taylor expansion exerted on the integrand, we split the integral into two terms. One of them has the slowly converging integrand but can be integrated, and the other has the fast decaying integrand and thus can be approximated simply by trapezoidal rule. We also present several numerical cases to justify this new split pattern. But, we do not optimize Ewald1D temporarily which we believe can be dealt with by the same way as we optimize Ewald2D.
As we study the mathematics and physics in Ewald sums, we find that the reciprocal terms in Ewald1D and Ewald2D are very important but the same term in Ewald3D is always left out without convincing reasons. More seriously, this simplified Ewald3D (or Ewald3D-tinfoil) is just the standard method in the simulation community. The ambiguity needs to be clarified, which is just the last research in my dissertation. As far as we are concerned, other methods, like MMM type or Lekner sum type, always imply some extra boundary condition in its summing order, so it is better to choose another very clumsy method as the reference, that is the truncating the crystal by shells of unit cells. This strategy has been used in 1932 by Evjen when he studied on the electrostatic energy of finite cr
老家的树ystal by expanding the energy into multipoles in which the higher ordered terms decay fast and the lower ordered terms vanish because of charge neutrality of the cell.
Using Taylor expansion on the diverging term of Ewald3D, we keep the conditioned converging term and abandon the higher ordered terms and the constant term. As a result, it gives rise to a simple formula without any loss. But we do not know the full information the formula shows. Thus we rewrite this conditioned term in Fourier form, and by Gauss theorem re-express it into a real space form. By some numerical comparisons between Ewald3D under various conditions and corresponding truncation methods, the physical picture can be directly explained. At last, we simply discuss the definition of “infinite”.
Keywords:
molecular simulations, electrostatic algorithms, Ewald sums
IV
目录
第1章绪论 (1)
1.1分子模拟的简介 (1)
1.1.1计算机模拟的由来、早期发展及意义 (1)
1.1.2分子模拟的模型 (2)
1.1.3周期性边界条件 (2)
1.1.4粒子间相互作用 (4)
1.1.5MC原理简述 (8)
1.1.6MD原理简述 (9)
1.1.7这篇论文与分子模拟的联系 (10)
1.2静电作用的精确计算 (11)
1.2.1Ewald型算法 (12)
1.2.2MMM型方法 (23)
联通宝视通
1.3基于Ewald求和的近似算法 (28)
1.3.1体相静电算法：Particle Mesh Ewald (28)
1.3.2界面静电算法：Ewald3DC (31)
1.3.3界面静电算法：Ewald3DLC (39)
1.4第1章小结 (41)
第2章优化Ewald2D方法 (43)
2.1Ewald2D数值近似方法概要 (43)
2.2Mori的误差分析方法原理 (45)
2.3特殊积分的数值计算与误差分析（其一） (49)
2.3.1积分路径与误差限 (51)
2.4特殊积分的数值计算与误差分析（其二） (62)
2.5Ewald2D与Ewald3DLC (72)
2.6计算精度举例 (73)
2.7优化Ewald2D (80)
2.8第2章小结 (88)
第3章Ewald1D方法的再发展 (91)
V
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议