2024.11.26 O原子的能量组成

序

这是一个不够化学的桃子对反应力场模型的惊叹。（并不保证以下内容的正确性，仅为个人工作记录）

从一个O原子开始的模拟

如图所示，把1个O原子放进模拟盒子里，如果用反应力场（ReaxFF）模拟，会奇怪的发现它的势能不为零。

或者，如果让一个O原子逐渐远离其他原子，即使非常非常远，依然显示有势能。

这一个原子的势能是哪里来的呢？

一个O原子在盒子里的秘密

我们首先关注到了上面一同输出的还有nlp=3。众所周知，反应力场里面确实有这么一项，$E_{\mathrm{lp}}$：

这一项代表孤对电子能。显然，单独的一个O原子，周围有3对孤电子对，$E_{\mathrm{lp}}$必然做出了贡献。

不过我们实际用compute pair命令输出一下，会发现这其实不是主要的影响因素，小到完全可以忽略不计。

写个备注，lammps里面的反应力场比较独立，它的输出并没有完全和lammps的命令配套。因此像compute pe/atom 这样的命令是没法细致的区分“pair or bond or angle or dihedral or improper or kspace or fix”的。反应力场里面分分合合的分子也是没法通过mol-id来描述的。

compute pair命令自然也做不出epair or evdwl or ecoul这些区分。取而代之的是特别的reaxff设置。输出以下14项能量（全局的，thermo命令输出）：

看看结果。和ea相比，elp真是小得微不足道……

所以这里的 atom energy 又是什么呢？

无意中在文献 [Bertels 2020 (doi:10.1021/acs.jpca.0c02734)] 看到一个说法：

对于过配位和前配位的原子，有两个能量校正项：$E_{\mathrm{over}}$和$E_{\mathrm{under}}$。

$E_{\mathrm{over}}$这一项的公式中出现了BO（bond order），由于我们的系统里只有一个原子，BO=0，因此$E_{\mathrm{over}}=0$。计算可得$\Delta_i^{\mathrm{lpcor}}\approx-2.098$，进而可以得到$E_{\mathrm{under}}\approx-2.51 \mathrm{kcal/mol}$，接近lammps软件输出的$-2.46 \mathrm{kcal/mol}$。

因此文献中说

$$ ea\ (\mathrm{atom\ energy}) = E_{\mathrm{over}}+E_{\mathrm{under}} $$

没毛病。

回过头来仔细看一下$E_{\mathrm{lp}}$的描述，“…（当BO很大导致孤电子对逐渐被拆散）…This is accompanied by an energy penalty…”也就是说$E_{\mathrm{lp}}$是孤对电子被破坏的能量损失，没有受到影响的时候（比如只有一个O原子）当然$E_{\mathrm{lp}}=0$啦！

如果盒子里有一个O$_2$分子，$E_{\mathrm{lp}}=0$就稍微明显一点了，但也只有一点点点点。依然可以忽略不计，毕竟成键数量完美，没有孤电子对受到破坏。

这个时候pe里面占比最高的是eb（bond energy），其次是ew（van del waals energy）和ea（atom energy）。且ew>0

千里之外的影响从何而来

还存有疑惑的地方是，把O原子无限远离其他原子后，pe大约-38kcal/mol,和单独一个O原子的-2.46kcal/mol不同。

暂时认为这可能是与$E_{\mathrm{coulomb}}$有关。如下图，后面的阶梯状变化是每次删除一部分其他原子，直至最后只剩我们统计的那一个原子：

这不科学……查文献！！！

参考一篇专门研究ReaxFF的论文，QEq电荷平衡方法是EEM (Electronegativity Equalization Method) 的变体，确实有被称为“unphysical long-range charge transfer”问题，对此前人已有些研究。

在lammps中，有两个命令fix acks2/reaxff command和fix qtpie/reaxff command，ACKS2原理上基于DFT，更准确更物理，但需要更复杂的求解器，更难算；QTPIE好算一些，但是基于QEq加入了基于直觉而非严格推导的修正。前者于2021年加入lammps，后者是最近刚刚被加到lammps的新命令：

手册中明确提到，“ACKS2 impedes unphysical long-range charge transfer sometimes seen with QEq”，以及 QTPIE “penalize long-range charge transfer seen with the QEq charge equilibration scheme”。也就是说，用ACKS2和QTPIE代替QEq都可能解决这个问题。

对比三种命令：

*ACKS2的报错：

WARNING: Fix acks2/reaxff BiCGStab convergence > failed after 200 iterations at step 15 (../fix_acks2_reaxff.cpp:606)
ERROR on proc 0: step 14: bondchk failed: i=7448 end(i)=203608 str(i+1)=203607 (../reaxff_forces.cpp:97)

参考网址中的讨论以及acks2相关论文和补充材料，ACKS2的问题可能是与ReaxFF的参数不匹配。相比之下，QTPIE需要的高斯指数可以在Chen的论文中找到部分。

使用fix qtpie/reaxff command之后，c_PE的变化还是与$E_{\mathrm{coulomb}}$有关的：

现在势能和电荷量都与0的差距小了一个数量级，更符合直观认识（一个原子远离其他原子时，与其他原子的相互作用越来越小，势能趋于单个原子在盒子中的结果）

其他

1. compute pe/atom和compute pair里面的nsub是使用了hybrid之后的第几个势函数模型

2. ReaxFF中，

   • bond order的公式在reaxff_bond_order.cpp
   • bond energy在reaxff_bonds.cpp
   • lone pair energy的公式(7-8)在reaxff_bond_order.cpp，其余在reaxff_multi_body.cpp
   • overcoordination和undercoordination在reaxff_multi_body.cpp
   • lone pair energy的公式(8)里面int是抹掉小数（向0取整），且第一个int前面有-（负号）（2008年Kimberly ChenowethAdri和C. T. van Duin等人文章的附件说明里漏了这个负号）
   • lone pair energy的公式(9)里面$n_{\text{lp,opt}} = \frac{1}{2} \cdot \left( Val_i^\text{e} - Val_i \right)$

参考资料：

1. https://pubs.acs.org/doi/10.1021/jp709896w
2. https://www.scm.com/doc/ReaxFF/ffield_descrp.html
3. https://pubmed.ncbi.nlm.nih.gov/32501686/
4. https://github.com/lammps/lammps/tree/develop/src
5. Chen Jiahao, Martínez Todd J. QTPIE: Charge transfer with polarization current equalization. A fluctuating charge model with correct asymptotics. Chemical Physics Letters 2007;438(4–6):315–20. Doi: 10.1016/j.cplett.2007.02.065.
6. 其他前面可以点击的链接……不一一罗列了