ABACUS+Bader charge 分析教程

作者：黄一珂，邮箱：huangyk@aisi.ac.cn

最后更新时间：2024 年 7 月 13 日

前言

我们在软件开发时不可避免遇到化学相关应用场景的需求，其中一个可能重点关注的信息便是“价态”。尽管从原理上考虑，“价态”这个量的定义并不明确，但化学场景的使用者总希望找到各式各样的方法来契合他们的“化学直觉”，便诞生了各种布居分析（Population Analysis）方法：如 Mulliken、Hirshfeld、Hirshfeld-Iterative (Hirshfeld-I)、ADCH、DDEC6、CM5、Mayer、Bader 等，以及各种电子的定域化分析方法，如 Wannier、ELF 等。尽管他们各自存在或致命或物理上定义不明的问题，如 Mulliken 电荷具有非常显著的基组相关效应，Hirshfeld 电荷经常对过渡元素计算出携带负电，Bader 对于电荷分布极其不均匀的体系会得出异常结果等，但仍然有足够普遍的受众群体。

Bader charge 是将电子密度进行类似于 Voronoi tesellation 的剖分，在原子间以电荷密度的梯度为 0 处创建截面，然后将原子周围所有截面包络出的 cell 进行实空间积分，即为 Bader charge。考虑实际计算，我们总能对价电子的密度进行切分，但对近核电子的处理则有所差别：CP2K 支持构造出核电荷分布，VASP 支持从 PAW 势中重建核电荷分布，而 QE 和 ABACUS 则不支持以上处理方式。

Bader 分析软件部署

网站：https://theory.cm.utexas.edu/henkelman/code/bader/

我们可以直接运行 linux 命令以下载 Bader 电荷分析软件的可执行文件版本，也可以选择下载源码然后手动编译：

wget https://theory.cm.utexas.edu/henkelman/code/bader/download/bader_lnx_64.tar.gz
tar -zxvf bader_lnx_64.tar.gz
chmod +x bader

下载源码：

wget https://theory.cm.utexas.edu/henkelman/code/bader/download/bader.tar.gz
tar -zxvf bader.tar.gz
cd bader
make
chmod +x bader

使用示例

在了解 bader charge 的原理后，我们首先从 nspin 1（only one spin channel）开始讲解具体使用。

nspin 1 case

生成所需 cube 文件

bader 程序需要读取以 cube 方式，或者 CHGCAR 格式存储的实空间电荷密度，因此运行 ABACUS 时，只需要添加

out_chg 1

即可输出电荷密度。

运行 bader 电荷分析

只需要以

/path/to/bader Si2-ELECTRON_DENSITY.cube

即可，屏幕会有如下输出：

GRID BASED BADER ANALYSIS  (Version 1.05 08/19/23)

  OPEN ... Si2-ELECTRO
  GAUSSIAN-STYLE INPUT FILE
  DENSITY-GRID:   60 x  60 x  60
  CLOSE ... Si2-Si2.cube-ELECTRO
  RUN TIME:    0.02 SECONDS

  CALCULATING BADER CHARGE DISTRIBUTION
                 0  10  25  50  75  100
  PERCENT DONE:  **********************

  REFINING AUTOMATICALLY
  ITERATION: 1
  EDGE POINTS:         97653
  REASSIGNED POINTS:    7908

  RUN TIME:       0.23 SECONDS

  CALCULATING MINIMUM DISTANCES TO ATOMS
                 0  10  25  50  75  100
  PERCENT DONE:  **********************
  RUN TIME:    0.02 SECONDS

  WRITING BADER ATOMIC CHARGES TO ACF.dat
  WRITING BADER VOLUME CHARGES TO BCF.dat

  NUMBER OF BADER MAXIMA FOUND:              8
      SIGNIFICANT MAXIMA FOUND:              8
                 VACUUM CHARGE:         0.0000
           NUMBER OF ELECTRONS:        8.00001

同时生成了 ACF.dat，AVF.dat 和 BCF.dat 文件，ACF.dat 中即是 bader 切分该电荷密度后在“Voronoi cell”切分的结果：

#         X           Y           Z       CHARGE      MIN DIST   ATOMIC VOL
 --------------------------------------------------------------------------------
    1    3.653762    2.109500    1.491642    3.999891     1.952476   137.965434
    2    0.000000    0.000000    0.000000    4.000124     1.952474   137.962879
 --------------------------------------------------------------------------------
    VACUUM CHARGE:               0.0000
    VACUUM VOLUME:               0.0000
    NUMBER OF ELECTRONS:         8.0000

我们将“CHARGE”列数据和赝势中 PP_HEADER 部分的 z_valence 值相减，即可得到该原子所携带电荷量。

nspin 2 case

ABACUS 目前分 spin channel 进行电荷密度的输出，若 nspin 2，则会在目录下输出 SPIN1_CHG.cube 和 SPIN2_CHG.cube。然而，直接对每一个 spin channel 进行体积剖分，然后按照原子序号相加的结果可能是错误的，因为两个 spin channel 中不一定每个原子均具有位置相同的“0 梯度面”。因此对同一原子，可能在两个 spin channel 中具有不同大小的 cell。为了避免这种情况，可以首先将两 cube 文件进行加和，之后再调用 bader.x。

Cube manipulator

Cube manipulator 是处理 Gaussian cube 格式 3D 实空间格点数据的小工具，可以在 ABACUS 代码包 tools 文件夹下找到（https://github.com/deepmodeling/abacus-develop/blob/develop/tools/plot-tools/cube_manipulator.py）。

通过 --help 可以查看可用选项：

myaccount@mycomputer:~/abacus-develop/examples/spin_polarized/AFM/OUT.ABACUS# python3 cube_manipulator.py --help
usage: cube_manipulator.py [-h] [-i INP] [-o OUT] [-s SCALE] [--p1d P1D] [--s2d S2D] [-p PLUS] [-m MINUS]

manipulate the Gaussian cube format volumetric data.

options:
  -h, --help            show this help message and exit
  -i INP, --inp INP     the input Gaussian cube file.
  -o OUT, --out OUT     the output file.
  -s SCALE, --scale SCALE
                        scale the Gaussian cube file by a factor.
  --p1d P1D             integrate the Gaussian cube file in 2D, followed by the axis: 'x', ...
  --s2d S2D             slice the Gaussian cube file along one axis, followed by string like 'x=0.0', 'y=0.0', 'z=0.0'. Note: should be fractional coodinate.
  -p PLUS, --plus PLUS  plus the two Gaussian cube files.
  -m MINUS, --minus MINUS
                        minus the two Gaussian cube files.

Once meet any problem, please submit an issue at: https://github.com/deepmodeling/abacus-develop/issues

生成所需 cube 文件

我们选择 abacus-develop/examples/spin_polarized/AFM 算例，将 INPUT 文件中 out_chg flag 修改为 1，得到两个 cube 文件分别对应于两个 spin channel。此时我们只需要将两个 SPIN*_CHG.cube 相加即可：

python3 cube_manipulator.py -i SPIN1_CHG.cube -p SPIN2_CHG.cube -o ELECTRONIC_DENSITY.cube

运行 bader 电荷分析

对文件 ELECTRONIC_DENSITY.cube 使用 bader.x，得到：

GRID BASED BADER ANALYSIS  (Version 1.05 08/19/23)

  OPEN ... ELECTRONIC_DENSITY.c
  GAUSSIAN-STYLE INPUT FILE
  DENSITY-GRID:   36 x  36 x  36
  CLOSE ... ELECTRONIC_DENSITY.c
  RUN TIME:    0.00 SECONDS

  CALCULATING BADER CHARGE DISTRIBUTION
                 0  10  25  50  75  100
  PERCENT DONE:  ********************** 

  REFINING AUTOMATICALLY
  ITERATION: 1
  EDGE POINTS:         39985
  REASSIGNED POINTS:    6753

  RUN TIME:       0.07 SECONDS

  CALCULATING MINIMUM DISTANCES TO ATOMS
                 0  10  25  50  75  100
  PERCENT DONE:  **********************
  RUN TIME:    0.00 SECONDS

  WRITING BADER ATOMIC CHARGES TO ACF.dat
  WRITING BADER VOLUME CHARGES TO BCF.dat

  NUMBER OF BADER MAXIMA FOUND:            108
      SIGNIFICANT MAXIMA FOUND:            108
                 VACUUM CHARGE:         0.0000
           NUMBER OF ELECTRONS:       31.99909

该例为，呈反铁磁态，每个 Fe 原子有 16 个价电子，最终积分得到电子数量为 31.99909，因此符合预期。同样查看 ACF.dat 文件：

#         X           Y           Z       CHARGE      MIN DIST   ATOMIC VOL
 --------------------------------------------------------------------------------
    1    0.000000    0.000000    0.000000   15.999544     2.060461    76.704205
    2    2.676622    2.676622    2.676622   15.999544     2.060459    76.704205
 --------------------------------------------------------------------------------
    VACUUM CHARGE:               0.0000
    VACUUM VOLUME:               0.0000
    NUMBER OF ELECTRONS:        31.9991

发现两 Fe 原子均为电中性，同样符合预期。

附录：cube manipulator 更多的使用方法

Cube manipulator 还有除了加法之外的更多功能，例如两 cube 文件进行减法、一个 cube 文件进行乘法等，以及还有面向更实际的应用场景，2D 平面积分与 cube 3D 内容切片。

AXPY 类操作

AXPY 是 BLAS 中函数，代表了诸如：

$y=\alpha x+\beta y$

的操作。因此-p flag 即为$\alpha=1, \beta=1$，-m flag 即为$\alpha=1,\beta=-1$，-s 为$\alpha=?, \beta=0$。

例如将 SPIN1_CHG.cube 和 SPIN2_CHG.cube 相减可以获得自旋密度：

python3 cube_manipulator.py -i SPIN1_CHG.cube -m SPIN2_CHG.cube -o SPIN_DENSITY.cube

Profile1d 操作

在进行表面模型的功函数计算时，需要对 Hartree 势进行面积分，得到沿某个轴的 1 维曲线。我们分别使用 ELECTRON_DENSITY.cube 和 SPIN_DENSITY.cube 对该功能进行演示：

python3 cube_manipulator.py -i ELECTRON_DENSITY.cube --p1d x -o charge1d.dat
python3 cube_manipulator.py -i SPIN_DENSITY.cube --p1d x -o spin1d.dat

即对 yz 平面进行积分，进行简单的绘图：

图 4. charge1d.dat，0和0.5位置各有一个Fe原子，Bader电荷切分的地方接近0.25和0.75

图 5. spin1d.dat，spin density由正到负，符合反铁磁Fe2的预期

Slice2d 操作

有时我们会考察某个截面的电子密度，对原子周围电子的局域化程度进行分析。我们仍然分别使用 ELECTRON_DENSITY.cube 和 SPIN_DENSITY.cube 对该功能进行演示：

python3 cube_manipulator.py -i ELECTRON_DENSITY.cube --s2d x=0 -o charge2d_000.dat
python3 cube_manipulator.py -i ELECTRON_DENSITY.cube --s2d x=0.25 -o charge2d_025.dat
python3 cube_manipulator.py -i ELECTRON_DENSITY.cube --s2d x=0.5 -o charge2d_050.dat

python3 cube_manipulator.py -i SPIN_DENSITY.cube --s2d x=0 -o spin2d_000.dat
python3 cube_manipulator.py -i SPIN_DENSITY.cube --s2d x=0.25 -o spin2d_025.dat
python3 cube_manipulator.py -i SPIN_DENSITY.cube --s2d x=0.5 -o spin2d_050.dat