静电相互作用设置

coulombtype

Cut-off: 具有对列表半径 rlist 和库仑截止 rcoulomb 的平面截止，其中 rlist>=rcoulumb，计算需要更加精准可以增加数值大小

Ewald: 经典的 Ewald sum 静电学。实空间截止 rcoulomb 应等于 rlist。使用例如 rlist=0.9，rcoulomb=0.9。在 reciprocal space 中使用的波矢量的最高幅度由 fourierspacing 控制。direct/reciprocal space 的相对精度由 ewald rtol 控制。注：Ewald 扩展为 O(N 3/2)，因此对于大型系统来说速度非常慢。它主要用于参考，在大多数情况下，PME将表现得更好。

静电作用计算的公式，由于静电作用属于长程作用，就是随着距离增加衰减的较慢，这样cutoff就需要设置很多，对于周期性体系n也必须设置很大，所以需要引入一定的算法进行处理。

Ewald的处理是将长程作用与短程作用分开计算，长程用在倒空间计算，短程作用则在实空间计算

因为短程在实空间中收敛较容易，但是长程不易收敛，则用倒空间进行处理，Vdir是实空间计算，Vrec是倒空间计算。

β参数用于决定直接空间加和与倒易空间加和之间的相对权重。

由于截断会带来一定的误差，因此需要加入补余误差函数V0。

这样, 我们可以对直接空间加和使用较短的截断(数量级为1 nm), 对倒易空间加和也可以使用较短的截断(例如每个方向10个波矢)。不幸的是, 倒易空间加和的计算量以N2增加(或N3/2, 如果使用更好一点的算法), 因此用于大的体系是不现实的。

每个带正电荷的离子都被周围带负电荷的离子所包围，反过来也是一样。这此外围相反电荷形成的静电屏蔽直接削弱了中心电荷与远程其它电荷的静电相互作用，这样就可以将之前电荷之间的长程相互作用看作有限长度的短程相互作用。Ewald加和方法就是以此为出发点，人为地在每个电荷周围引入一个与该电荷相反的等值屏蔽的电荷分布，并且使每个屏蔽电荷分布都有一个相反的补偿电荷分布水维持原来的电荷分布。

该方法将点电荷之间的相互作用，转化为高斯球分布的等量的异性点电荷，目的是为了加速相互作用的收敛，因为长程作用收敛速度慢，

计算1：中心电荷与近距离截断范围内电荷的静电相互作用

计算2：中心电荷与所有周期性补偿电荷分布之间的静电相互作用

计算3：需要扣除自己的高斯静电势的影响

不要使用Ewald方法, 除非你绝对确定你需要使用. 几乎对所有的情况, 下面的PME方法表现都更好. 如果你依然坚持使用精度的Ewald加和方法, 在你的.mdp文件中输入以下内容, 如果盒子的边长3 nm左右:

coulombtype  = Ewald

rvdw            = 0.9

rlist          = 0.9

rcoulomb    = 0.9

fourierspacing = 0.6

ewald-rtol      = 1e-5

盒子尺寸的比例和fourierspacing参数决定了每个方向使用的波矢 mx my mz的最大振幅. 例子中对3 nm的立方盒子每个方向将使用11个波矢(从–5到5). 

ewald-rtol参数为在截断处静电相互作用的相对强度. 减少此值可以得到更加精确的直接空间加和, 但倒易空间加和的精度会略微降低.

PME: 用于 electrostatics（具体指静电相互作用或库仑力）的 Fast smooth Particle-Mesh Ewald (SPME) 。 Direct space 类似于Ewald sum，而 reciprocal space 使用 FFT 执行。网格尺寸由 fourierspacing 控制，插值顺序由 pme-order 控制。在网格间距为 0.1 nm 和三次插值的情况下，静电力的精度为 2-3*10 -4。由于 vdw 截止的误差大于此值，您可以尝试 0.15 nm。当并行运行时，插值的并行性优于FFT，因此尝试在增加插值的同时减少网格维数。

经典 Ewald 求和计算倒易空间部分能量的时间复杂度为 O(N2)，这是由离散 Fourier 变换（discrete Fourier transform，DFT）的基本性质决定的。对于较大体量的模拟系统，使用基本的离散 Fourier 变换处理的计算效率较低，如果可以将点电荷通过插值算法分布到空间中的均匀网格上，就取得了使用快速 Fourier 变换（fast Fourier transform，FFT）的条件，将时间复杂度降为 O(Nlog2N)，从而极大地提升计算效率。PME 求和算法正是在这种思想的指导下被开发出来。这种方法不直接对波矢进行加和, 而是使用内插方法将电荷分配到网格上。GROMACS实现的PME方法使用了基数B样条插值, 通常被称为平滑PME(SPME)。先使用3D FFT算法对格点进行傅里叶变换, 在k空间中利用对格点的单个加和就可以得到倒易空间的能量项。格点上的势能可利用逆变换进行计算, 通过使用内插因子就可以得到每个原子上的力。

coulombtype= PME

rvdw = 0.9

rlist = 0.9

rcoulomb = 0.9

fourierspacing = 0.12

pme-order = 4

ewald-rtol = 1e-5

fourierspacing参数决定了FFT格点的最大间距(即, 格点数的最小值), 

pme-order控制插值的阶数，使用四阶(立方)内插与指定的间距，静电能应该能精确到5*10-3 

由于Lennard-Jones能量的精确度更低一些, 还可以轻微地增加间距

PME可用于压力缩放, 

但当心, 在一些体系中各项异性缩放可能会导致虚假的有序结构.

P3M-AD: 带用于长程静电相互作用的解析导数（analytical derivative, AD）的 Particle-Particle Particle-Mesh（P3M）算法。

该方法及其代码与 SPME 相同，只是影响函数（influence function）针对网格进行了优化，使精度略有提高。

在GROMACS中, 也可使用Hockney和Eastwood发展的粒子-粒子粒子-网格方法来处理长程静电相互作用.

 动力学计算中, 尽管P3M方法是第一个用于分子模拟的高效长程静电方法, 但在原子模拟中, P3M方法很大程度上已经被平滑PME方法(SPME)所取代。

原始P3M方法计算时有一个缺点: 它需要3次3D-FFT后变换来获得粒子上的力. 但P3M并不需要这样, 可以通过对势能进行数值微分得到所需的力, 像在PME中一样. 这样的方法被称为P3M-AD。

P3M-AD和PME唯一的不同在于: PME对晶格Green影响函数进行了优化以尽量减小误差. 然而, 在2012年有研究表明, 可以修改SPME影响函数得到P3M. 

这意味着P3M-AD误差最小化的优点在PME中可以同样的计算代价和同样的代码获得, 只需增加几行代码修改影响函数。

然而, 对最佳参数设置, P3M-AD误差最小化的效果小于10%. 对于隔行(也称为交错)的格点, P3M-AD确实展现了很大的精度提升, 但GROMACS(目前)并不支持这种处理（暂时不知道新版本是否支持）

Reaction-Field: 具有库仑截止 rcoulomb 的反应场处理方法，且 rlist>=rvdw。超过截止值的介电常数为 epsilon-rf。介电常数可以通过设置 epsilon-rf＝0 而设置为无穷大。

这个原子在截断半径内的都进行计算

半径外的全部用恒定的介电常数来处理

半径内的也是用介电常数来处理

半经内的数目会经常发生变化

导致计算后的数据不连续

介电常数需要提前知道

所以很麻烦

但是效果就是计算的速度更快

User: 注意可以使用的版本

gmx mdrun 现在希望找到一个文件表table.xvg，其具有用户定义的势能函数, 包括排斥, 色散和库仑相互作用。

当存在对相互作用时，gmx mdrun 还希望找到一个文件 tablep.xvg 用于对相互作用。当相同的相互作用应用于 non-bonded 及 pair interactions时，用户可以为这两个选项指定相同的文件。

这些文件应包含7列：x值, f(x), -f'(x), g(x), -g'(x), h(x), -h'(x)，其中 f(x) 是库仑函数，g(x) 色散函数，h(x) 排斥函数。

如果未将 vdwtype 设置为 User，则忽略 g、-g'，h 和 -h' 的值。

对于非键相互作用，x值应从0到最大截止距离 + table-extension，且应间隔均匀。

对于配对相互作用，将使用文件中的表长度。

用于非用户表格的最佳间距在混合精度下运行时为 0.002 nm，在双精度下运行后为 0.0005 nm。

x=0 时的函数值不重要。更多信息见手册。

PME-Switch: Direct space 部分的切换函数（switch function）（见上文）同 PME 的组合，rcoulomb 可以小于 rlist。

PME-User: PME 同用户表的组合（见上文），rcoulomb 可以小于 rlist。通过gmx mdrun 从用户表中减去 PME 栅格贡献，这使得用户表应该包含大约10位小数。

PME-User-Switch: PME-User 和 Switch 功能的组合（见上文）。切换函数应用于终末粒子-粒子相互作用，例如同时用于用户提供的函数和 PME 栅格校正部分。

Generalized-Reaction-Field: 2022.2中已弃用。库仑截止为 rcoulomb 的广义反应场，其中 rcoulumb≥rlist。超过截止值的介电常数为 epsilon-rf。离子强度由带电（即非零电荷）电荷基团的数量计算。GRF 电势的温度用 ref-t[K] 设定。

Reaction-Field-zero: 2022.2中已弃用。在GROMACS中，cutoff-scheme=group 的正常反应场静电导致较差的能量守恒。Reaction-Field-zero 通过使电势零点超过截止值来解决这个问题。它只能用于无限介电常数（epsilon-rf=0），因为只有该值的力在截止处消失。rlist 应比 rcoulomb 大 0.1 至 0.3 nm，以适应电荷组的大小和相邻列表更新间的扩散。这一现象以及使用表查找代替分析函数的事实，使得 Reaction-Field-zero 在计算上比正常反应场更耗费算力。

Shift: 2022.2已弃用。类似于 vdwtype 中的 Shift。你可能想用 Reaction-Field-zero 来代替，它具有相似的势形状，且有物理解释，并且由于排除校正项，它具有更好的能量。

Encad-Shift: 2022.2不再支持。使用 Encad 模拟软件包中的定义，库仑电势在整个范围内被降低。

Switch: 2022.2已弃用。类似于 vdwtype 中的 Switch。切换库仑电势可能会导致严重的伪影，建议使用 Reaction-Field-zero。

coulomb-modifier

Potential-shift: 将库仑电势增减一个常数，使其在截止点处为零。这使得势成为力的积分。请注意，这不会影响力或采样过程。

None: 使用未修改的库仑势。这在将能量与用其他软件计算的能量进行比较时非常有用。

Potential-shift-Verlet: 2022.2中已弃用。cutoff-scheme=Verlet 时启用 Potential-shift，因为它（几乎）是无限制的；cutoff-scheme=group 时启用 None。

rcoulomb-switch (0) [nm]

设置从何处开始切换库仑电势，仅在使用力或电势切换时有意义。

rcoulomb (1) [nm]

库仑截止距离。请注意，对于PME，此值可以通过 gmx mdrun 中的 PME 调优以及 PME 栅格间距来增加。

epsilon-r (1)

相对介电常数，0表示无穷大。

epsilon-rf (0)

反应场的相对介电常数。这仅用于 reaction-field，0表示无穷大。