如何对量化软件中输出的原子实R_i对r处的势U(R_i,r)去除非库伦贡献

赵小壮

我在阅读该篇文献（https://doi.org/10.1021/jz900268s）时，遇见这段表述：

"Here, \rho(r) is the valence electron density,  R_i is the position of pseudonuclei i, and Z_i is the pseudonuclear charge in the pseudopotentials (PP), with ZO = +6|e| and ZH = +1|e|. By default, VASP instead outputs the negative of eq(2) after replacing the last term with the entire local PP U(R_i, r). The short-range, nonCoulombic contribution to U(Ri, r) should be removed when computing φ. This ambiguity arises from the use of PPs and is absent in all-electron calculations."

第一，我想问一下我的理解对不对。R_i位置处的原子实对位置r处的作用势在理论上以及全电子计算中如HF里全部都是静电势，但是在DFT里由于使用泛函所以还有非静电势的项。所以在VASP里使用DFT/PBE输出的就包含非静电项，如果我只想得到纯静电势，就需要减去VASP输出的原子实部分，自己按eq(2)最后一个求和项coding算纯静电势。


第二，如果上边理解无误的话，我想问软件输出的r处的势V(r)应该来源既有原子实也有电子，然我如何分辩出原子实的贡献呢？
第三，如果不能分辨，那我可能需要全部自己coding算，原子实的位置是知道，但如何得到体系中电子的电荷网格呢，在CP2K中可以直接得到吗，谢谢大家

万里云

静电势包括Hartree和原子核两部分。Hartree就是电荷密度积分那一项，原子核是带Z那一项。采用赝势后，实际计算用的不是图里的公式，而是赝势U(r)。作者表达的就是这个意思。

赵小壮

万里云 发表于 2023-10-26 11:21
静电势包括Hartree和原子核两部分。Hartree就是电荷密度积分那一项，原子核是带Z那一项。采用赝势后，实际 ...

谢谢，请问实际计算中Hartree贡献的部分和原子核的部分是统一用赝势算的吗，还是是分开算？分开算的话软件里算的Hartree部分就是图片里电荷密度积分那一项吧，只是第二项原子核的部分有差异，请问我理解的对吗，谢谢

卡开发发

赵小壮 发表于 2023-10-26 12:07
谢谢，请问实际计算中Hartree贡献的部分和原子核的部分是统一用赝势算的吗，还是是分开算？分开算的话软 ...

如果你要单独计算电子Hartree贡献，那么只计算你所谓的电子密度积分那一项。需要注意的是，在使用赝势的情况下电子密度会被替换为赝电子密度（仅包含价电子且在原子芯附近更光滑）。电子Hartree部分可以通过求解Poisson方程来实现，办法比较多，当然其中有一些细节问题这里没办法讨论到。VASP的LOCPOT没办法单独打印电子Hartree的部分。

赵小壮

卡开发发 发表于 2023-10-26 13:57
如果你要单独计算电子Hartree贡献，那么只计算你所谓的电子密度积分那一项。需要注意的是，在使用赝势的 ...

非常感谢您的解惑，“在使用赝势的情况下电子密度会被替换为赝电子密度（仅包含价电子且在原子芯附近更光滑）”，那请问如果我不使用赝势的话，计算出来的结果会包括价电子之外的东西的贡献吗？这个之外的东西会一般有什么呀？我想得到一个纯粹的价电子的贡献的话是必须使用赝势是吗？谢谢大佬，祝大佬长命百岁

卡开发发

赵小壮 发表于 2023-10-26 15:33
非常感谢您的解惑，“在使用赝势的情况下电子密度会被替换为赝电子密度（仅包含价电子且在原子芯附近更光 ...

这个价电子密度实际上也不是真实的密度，前面提到在原子芯内其实因为使用赝势的原因，波函数和密度相对于全电子情况（单电子近似下）在原子芯附近的密度要更为光滑，只有远离原子芯的位置才能保证与全电子情况价电子密度的贡献一致。当然，引入赝势使赝密度光滑在数值处理上是很有必要的，否则求解Poisson方程也要非常小心。

如果你真的需要全电子（冻芯）情况下价电子密度，可以借助PAW框架，用关键字LAECHG=T，这样会输出AECCAR2，这对应于此时的价电子密度，但这个泊松方程在数值技巧上应当小心处理。

卡开发发

赵小壮 发表于 2023-10-26 15:33
非常感谢您的解惑，“在使用赝势的情况下电子密度会被替换为赝电子密度（仅包含价电子且在原子芯附近更光 ...

如果真要处理那个全电子的价电子密度，你也许可以参考类似于CP2K中GAPW的手段，即便仅处理赝价电子密度的Poisson方程也存在一些trick。你方便的话可以透露一下你要这方面信息的意图，也许并不需要折腾这些。

赵小壮

卡开发发 发表于 2023-10-26 23:07
如果真要处理那个全电子的价电子密度，你也许可以参考类似于CP2K中GAPW的手段，即便仅处理赝价电子密度的 ...

谢谢大佬的解惑，我是想计算一个充满水的box里边每个位置上的感受到的库伦势。并且想看看两种方法的差异，
方法一  ：
就是先使软件直接输出，然而按照文献所说直接输出的有非库伦项，需要减去公式中原子实的贡献，替换为自己计算的原子实的贡献，
而自己计算原子实贡献时候，原子实是完全的O+6|e|, H +1|e|, 这就要求第一部分必须完全是价电子的贡献。然而经过上边的讨论按我的理解
软件直接输出的时候第一部分还不完全是价电子，那我进行替换后计算的东西的物理图像就不是我想要的 O+6|e|, H +|e|, 四个价电子-1|e|了。
方法二：
我想让软件输出box价电子密度分布，软件内部应该储存的是离散的电荷网格，此时原子实的位置知道，那么我就可以完全自己coding算，我想完全复现
方法一的结果，来证明软件中确实是这么算的。
谢谢大佬

卡开发发

赵小壮 发表于 2023-10-27 08:20
谢谢大佬的解惑，我是想计算一个充满水的box里边每个位置上的感受到的库伦势。并且想看看两种方法的差异 ...

我不确定你算这个东西要分析或者作为哪方面的用途，其实这是不太常用的信息，所以大部分程序输出的都是包含离子势的那个部分。并且实际计算当中涉及到一些trick，不见得是你计算想讨论的问题。

软件自身用的是赝密度处理的，全电子的密度是借助PAW框架变换过来的，所以一般计算过程应该不包含全电子（包括此时价电子）静电势计算的信息。

赵小壮

卡开发发 发表于 2023-10-27 10:50
我不确定你算这个东西要分析或者作为哪方面的用途，其实这是不太常用的信息，所以大部分程序输出的都是包 ...

我们其实是做Ewald summation的组，对于这种周期性体系，本身我们有完全exact的计算公式，但是我阅读文献发现这个作者用的是软件直接输出的，然后他说软件里就是按这个公式这么算的，所以我觉的不然就是软件错了，不然就是软件里不是按他说这个公式算的。因此想看看软件这这个公式对不对的上