3.4.2 差分电荷密度

差分电荷密度（charge density difference，电荷密度差分，一个意思），可以简单理解为同一个结构在成键（自洽计算）前后的电荷密度差值，用来直观定性地判断体系成键时不同区域的电子流向（A区域的电子流向B区域，使得成键）。

以CO分子为例：

（1）手动写一个POSCAR文件：先构建一个大胞（晶格大于10Å，尽可能减弱周期性带来的分子之间相互作用），在合适的位置放置C和O原子各一个，并且保证它们的距离合适。网上查找CO实验键长1.1283Å，所以建模型时键长也取1.128Å（取Cartesian坐标比较方便）。POSCAR内容如下

CO
1
10 0 0
0 10 0
0 0 10
C O
1 1
Cartesian
5 5 5
5 5 6.128

（2）（这一步结构优化可以省略）再做ISIF=2的结构优化（固定晶格，只优化原子位置）来得到理论键长（PBE泛函）。其中，KPOINTS只取111即可，ENCUT=520eV（C的ENMAX=400，这里取1.3倍比较安全）。POTCAR中C考虑了4个价电子，O考虑了6个价电子。基于先验知识，CO分子没有孤立电子，所以应该没有磁性，ISPIN=1即可。为了检验磁性，也可以在结构优化完，自己设置ISPIN=2（为简单起见，不额外设置初始磁矩，让程序自动设置）做一次scf计算，看OSZICAR中是否有磁矩来初步简单地判断是否有磁性。

确认结构优化成功后，从CONTCAR中读取得到的CO键长是1.143Å（比实验偏大，有可能是计算精度不够，不过PBE泛函的老毛病就是通常会高估键长和晶格）

（3）在新的目录下，用结构优化好的CONTCAR作为新的POSCAR，做一次静态自洽计算，并确保INCAR中设置了LCHARG=.TRUE.以保证输出自洽电荷密度文件。注：这一步的CHGCAR也可以直接用relax输出的CHGCAR。这一步我们获得了CO分子自洽的电荷密度文件，记为CHARGE-scf

（4）在新的目录下，做非自洽计算(NSCF)计算。这里的INCAR和上一步静态自洽计算的INCAR基本一样，唯一区别的参数就是ICHARG=12 （12的含义是生成初始电荷密度后，在后续迭代中不更新电荷密度）。这一步计算结束，得到了nscf的CHGCAR。这里的CHGCAR是VASP程序初猜的体系波函数，可以理解为未成键时的CHGCAR。

（5）将scf和nscf计算得到的CHGCAR文件都用vesta可视化，直观感受一下差别（文件有点大，打开需要点时间）

左边是scf的CHGCAR，右边是nscf的，可以看到左边经过自洽迭代得到的电荷密度上头大，下头小（红色是氧O原子），与右边不同。

（6）做scf和nscf的电荷密度差分：

方法一：可以在vesta里面做差，请自行网上搜索“vesta做电荷差分”关键词

https://blog.csdn.net/xk6891/article/details/104779701

方法二：借助辅助程序，这里使用vaspkit 的31 → 314 功能来处理：

                   Structural Utilities
1)  VASP Input-Files Kit      2) Mechanical Properties
3)  K-Pathfor Band-Structure  4) Structure Editor
5)  Catalysis-ElectroChem Kit 6) Symmetry Analysis
8)  Advanced Structure Models
               === Electronic Utilities
11) Density-of-States         21) DFT Band-Structure
23) 3D Band-Structure         25） Hybrid-DFT Band-Structure
26) Fermi-Surface             28) Band-Structure Unfolding
31）（ Charge & Spin Density    42) Potential Analysis
51) )Wave-FunctionAnalysis    62) Magnetic Properties
65） Spin-Texture              68） Transport Properties
                      Misc Utilities
71）（ Optical Properties       72) Molecular-Dynamics Kit
73) VASP2other Interface      74) User Interface
91) Semiconductor Kit         92) 2D-Material Kit
(0 Quit
        -->>

                    Charge Options
311) Charge Density
312） Spin Density
313) Spin-Up & Spin-down Densitv
314)( Charge-Density Difference
315) 2D Linear-Average Charge-Density
316) 1D Planar-Average Charge-Density
317) 1D Macroscopic-Average Charge-Density
318) Export CHGCAR/PARCHG to XcrySDen .xsf format
319) Export CHGCAR/PARCHG to Gaussian .cube format
320) Redefine Cell Size and Shape of CHGCAR
                    Advance Options ===
325)STM Simulation
(0  Quit
9)   Back
        --->>

                     File Options
Input the Names of Charge/Potential Fileswith Space:
(e.g., to get AB-A-B type:~/AB/CHGCAR ./A/CHGCAR ../B/CHGCAR)
./CHGCAR R./CHGCAR

注意是用scf的CHGCAR减去nscf的CHGCAR，所以两个CHGCAR的路径要注意前后顺序（有可能这一步会卡住，一直没有输出，当遇到这种情况时建议更换服务器或者用方法一）。将得到的CHGDIFF.vasp直接用vesta打开即可（下方最右侧）：

黄色区域对应的是电荷增多，蓝色对应减少（vesta默认设置）。这意味着CO分子中C与O在成键时，C和O（红色）的电子会往中间区域和上下两侧转移（相比于未成键时的孤立原子来说），中间黄色区域增加的电子可以理解为CO分子中的共用电子对。

有时候需要在vesta里面调整Objects-Properties-Isosurfaces来调整显示效果，请自行研究。

另外一个例子是MgB2（层状材料，常压下超导转变温度39K的超导体）。可以作为差分电荷密度的练习。MgB2的晶体结构数据如下：

mgensium diboride
1.0
3.068802008 -0.000000000 0.000000000
-1.534401004 2.657660498 -0.000000000
-0.000000000 -0.000000000 3.525151116
Mg B
1 2
Direct
0.000000000 0.000000000 0.000000000
0.333333343 0.666666687 0.500000000
0.666666627 0.333333313 0.500000000

这里直接给出差分电荷结果（橙色小球为Mg原子）：

黄色区域对应的是此处电荷增多，蓝色对应减少。很明显，定性上看，Mg的电子向B层转移了。有时候需要在vesta里面调整Objects-Properties-Isosurfaces来调整显示效果，请自行研究。

注意，这里画的差分电荷密度都是原子成键之后自洽的电荷密度，与每个原子之间未成键时的孤立原子电荷密度叠加，之间的差。实际上，差分电荷密度有很多种定义。取决于研究的问题，电荷密度差分的其它情况，请参考刘锦程博士的blog：https://blog.shishiruqi.com/2019/07/12/chgdiff/