带电体系的计算问题之浅见

aaa

关于带电体系，尤其是带电荷的周期性体系的计算，网上有很多人提问。在此，记录一下我针对这个问题查阅的资料以及一些肤浅的理解。

（1）cluster体系。

Gaussian中对带电体系的处理：高斯的input file的分子说明部分的第一行就是电荷和多重度。

我觉得关于这个部分的电荷的合理解释应该是：高斯计算需要知道体系的总电子数，由用户提供起来比较费事而且容易出错，所以分子说明部分的电荷和多重度的收集就是为了让高斯方便自动生成体系的电子数目。所以，该部分的电荷是体系的总电荷数。

比如：吸附带电团簇或离子的计算，你是没办法把电荷指给某个原子或者某个片段的。用英语来讲就是you could not assign charge to a specific atom in Gaussian input file。所以，这个电荷引起的变化只能是增加或减少体系的电子数。随后的scf迭代会不断地将电子进行重新分布直到一个合理的值（也就是收敛），随之电荷也进行了合理的分布。我们可以通过波函数分析计算原子电荷并查看每个原子的上的电荷分布，也可以以此来判断计算的合理性。

这个问题：社长在这个帖子中贴了两张ppt，是对这个问题最精准的描述。http://bbs.keinsci.com/thread-3970-1-1.html

还有就是高斯里的Charge关键词，用于指定包含背景电荷（background charge）分布的计算。

关于什么是背景电荷，先贴上社长的帖子：http://sobereva.com/579

我没有用过charge这个关键词，我理解的背景电荷就是用户在特定坐标规定了一个具有电荷特征的值，从而这个位置将会影响体系的静电势等一系列性质。

这个背景电荷的分布是由一系列点电荷组成。

总的来说，带电的cluster体系处理起来相对容易。

（2）周期性体系。

带电的周期性体系就没那么简单了，所以网上的提问基本上也是针对周期性体系的。

我先说下我自己对“为什么带电的周期性体系不能简单处理的理解”，这个段落完全是我和一个没做过第一性原理的朋友，拍脑袋想出来的，请大家多多批评指正。

假设一个由100个原子组成的二维单层平面，我们让这个体系带1e-的电量，对于1mol的周期来说，这个体系所带的电量就是6.02×10^23/100=6.02×10^21个e-的电量，我们把它换算成库伦就是：1.602×10^-19×6.02×10^21C=9.64×10^2C。

接下来我们看看库伦这个单位有多大：A coulomb is an enormous charge - two 1 C charges that are 1 m apart exert a force of 9 x 10^9 newtons (see Coulomb's law). That's over two million tonnes, ~720x as much as the thrust of a space shuttle solid rocket booster during liftoff （摘自https://energyeducation.ca/encyclopedia/Coulomb）.用中文来讲：两个相距1m的1C电荷会产生9 x 10^9牛顿的力，大约是200万吨的航天飞机固体火箭助推器在升空期间推力的720倍。（备注：后半句的翻译不知道对不对。这个后半句，我没太理解。）

我们的slab平面带了这么多电，必然导致费米能级剧烈上升，接着导致电子逃逸。也就是不管在理论中还是现实中，这个带电的slab都根本不可能存在。

以上纯属个人瞎猜。

插播一句网友们对这个原因的解释：这是由于带电的slab的周期性复制造成的。

我是没理解这句话。希望理解了这句话的老师、同学帮忙解答。

接着看看来自VASP的官方说法：

我没找到官方关于不能简单粗暴计算带电slab的原因，但是找到了这样做会导致的后果的解释及一些解决方案。

VASP的官方说法：https://www.vasp.at/wiki/index.php/FAQ

Reason: If you set the number of electrons in the INCAR file for a slab calculation you end up with a charged slab. The electrostatic energy of such a slab is however only conditionally convergent and worse, in practice, even infinite (basic electrostatics). Therefore, no method whatsoever exists to correct the error in the electrostatic energy. E.g. the energy converges towards infinity, when the vacuum width is increased. You can try to validate this, by simply increasing the vacuum width in VASP for a charged slab. You will find that the energy increases or decreases linearly with the vacuum width.

我的理解：如果我们像在高斯里那样，直接通过设置NELECT去模拟一个带电体系是不现实的，大部分情况下无法收敛。VASP也提供了一种可能会产生收敛值的解决方案：在固定真空层厚度的同时，不断增加slab平面的大小，直到平面无限大的时候得到一个外推能量值。我想这种做法的原理是不是就是使得slab体系的电荷不那么密了。不过VASP官方明显不推荐这种做法。

解决方案：

文献：Mauricio Chagas da Silva, Michael Lorke, Bálint Aradi, Meisam Farzalipour Tabriz, Thomas Frauenheim, Angel Rubio, Dario Rocca, Peter Deák, Self-consistent potential correction for charged periodic systems, Phys. Rev. Lett. 126, 076401 (2021). https://doi.org/10.1103/PhysRevLett.126.076401

文章中说SCPC方法已经被VASP 6.2.0 implemented。

github上这个方法的使用：https://github.com/isobestico/SCPC-Method/tree/master/VASP

再分享一个文献：Jian Liu, Shuming Wen, Xiumin Chen, Shaojun Bai, Dan Liu, Qinbo Cao, DFT computaion of Cu adsorption on S atoms of sphalerite (110), Minerals Engineering, 2013, 46-47, 1-5.

这个文献较早了，他的作者应该是在小*虫针对这个问题进行过提问。

审稿人提出了周期性带电体系无法给出合理的能量值的问题。

http://mu****g.com/t-5347282-1

他的解决方案应该是通过与审稿人协商，作者指出了计算的能量只是给出一个定性描述。

The above adsorption and substitution energies calculated from DFT cannot represent the real energies and have little physical meaning, due to the fact the calculations are performed in vacuum and the energies are different for different calculation systems. The calculated energies only provide a qualitative criterion for the adsorption and substitution of Cu2+ ion.

在那个年代也算解决问题了，但是不知道这样的回答是不是准确。

卡开发发

granvia 发表于 2022-3-25 18:43
PBC模型未必就是三维无限延展的。实际计算采取的是Born—von Karman近似，相当于用环来模拟周期性平移体 ...

确实，Born-von Karman的要求是使用一个首位相接的环来讨论这个问题。那么这几个问题我觉得可以思考，你想一下BZ积分为啥是积分？为啥要做收敛性测试？

granvia

再补充一个带电周期性体系是的问题： 当周期性模型中是带电粒子在金属表面上方时，就会产生一个镜像电荷，这就引入了一种（不真实的）原来电荷与其镜像电荷的作用。

granvia

卡开发发 发表于 2022-3-25 10:55
随意做点讨论，我以前和别人说也有说错的时候，下面的说法都只是我自己看到的资料做的一些简单总结，一些地 ...

PBC模型未必就是三维无限延展的。实际计算采取的是Born—von Karman近似，相当于用环来模拟周期性平移体系（只有环无限大才严格等价）。由于实际计采用算有限个k点，所以这个环的大小是有限的

卡开发发

额外的注释：
1、修正方法其实不只是SCPC，只是其他方法没有在VASP中实现，比如
10.1103/PhysRevB.96.075448
10.1103/PhysRevB.73.115407
10.1103/PhysRevB.90.245101
其中有些方法在QE中实现。

2、

我想这种做法的原理是不是就是使得slab体系的电荷不那么密了。

实际上电荷也不该那么密，比如在z方向。

3、表面带电吸附的情况，我印象中vasp的论坛有人也专门提到，如果不做校正或外推，算出的吸附能没有意义，最多同一个体系在c相同的情况定性看看吸附能来判断吸附强弱的相对趋势。

卡开发发

随意做点讨论，我以前和别人说也有说错的时候，下面的说法都只是我自己看到的资料做的一些简单总结，一些地方可能不大严谨，也没有展开说，就仅供参考和讨论，希望是真的能用在讨论上吧。

1、怪异的PBC
可曾想过，真正古怪的不是那个带电的slab，而是PBC（周期性边界）模型本身？如果你假想PBC的体系不是模型而是真的存在，一个晶体就会三维无限延伸充满全世界，这显然是荒谬的，因此要理解这个模型如何近似以及这种近似的好坏才比较重要。

2、能量大爆炸？
PBC计算中指定的电荷是每个晶格中的电荷（一般显示出来是其中一个，即单胞，也可能是多个，即超胞），因为这种无限重复，实际上过剩电荷之间的相互作用是无穷多个，这个Coulomb作用是发散的，但是实际计算处理会引入凝胶近似给背景补上均匀的正电荷来克服能量发散问题，即vasp手册：

If the number of electrons is not compatible with the number derived from the valence and the number of atoms a homogeneous background-charge is assumed.

。

3、不发散就行了？
上述处理后，晶格的能量（每个晶格平均的）保证不发散了，但并不意味着结果有意义。由于使用平面波方法，无论是晶体中低浓度的带电缺陷还是表面上低浓度的带电吸附，这种三维PBC是强制使用的，我们希望使用一个足够大但有限的晶胞或真空层分别取模拟上述两种情况，而此时体系总能量的变化会随着每一个方向上晶格尺寸L呈现1/L的关系：

A final note has to be made: besides the electrostatic interactions there exist also spurious elastic interactions between neighboring cells which (according to a simple “elastic dipole lattice model”) should scale like 1/L^3 (leading order). Therefore, the corrected values may still show a certain variation with respect to the supercell size.

，以至于收敛过于缓慢。

4、外推和校正
这种1/L的渐进收敛缓慢，我们根据这种特性其实可以做不同L的结果然后对1/L做外推。实际操作上这样的方法还是存在困难，对于大晶胞的计算本身比较缓慢，在外推过程中甚至体系的稳定结构也可能发生变化。究其根源，还是在人为引入的周期性边界上，故才有后续的这些校正方法。

sobereva

带电体系的周期性计算怎么考虑也看具体算什么问题。很多问题都依赖的是整体能量在不同状态间求差进行计算，虽然带电会造成与周期镜像的虚假的静电互斥作用，但如果前后状态的这种虚假的作用抵消得很充分，那么也是可以照常算的。为稳妥起见，做个晶胞尺寸的收敛性测试或外推总是有益的。