ABACUS：用极化差分计算Born有效电荷
以SiC八原子晶胞演示位移建模、GPU自洽、CPU Berry极化和Zxx提取。
阅读本篇目录
1. 这个例子算什么#
Born有效电荷描述绝缘体中原子位移引起的宏观极化变化，不等于Bader电荷、形式价或赝势价电子数。这里固定晶胞，对SiC施加沿x方向的微小位移，通过Berry极化的中心差分得到 Zxx，没有计算完整的3×3张量。
实际模型是8原子的立方SiC常规胞，晶格常数4.36 Å，4个Si和4个C，32价电子。采用PBE、ONCV规范守恒赝势、60 Ry截断、SCF 4×4×4网格，Berry阶段8×4×4网格。它是操作参考，未做截断能、网格和晶格常数的系统收敛研究。
下载完整输入、实际SCF与Berry输出、CSV及绘图脚本，将下载文件保存为 example.zip：
python3 -m zipfile -e example.zip .
cd abacus-sic-born-fd
module load ase/3.29.0
bash read.sh
bash plot.sh
读取已有结果不需要提交计算；analysis-read 保存提取的数据，figures 保存图片。需要自己重算时，再按下面的步骤准备与提交。
2. 准备成对位移#
本例有八个结构：Si子晶格整体±0.002 Å、C子晶格整体±0.002 Å、Si子晶格整体±0.004 Å，以及所有原子整体±0.002 Å。前四个结构给出Si、C的Zxx，另外四个用于观察步长与整体平移的数值误差。
四个Si在这个晶体中等价，因此同时移动四个Si后，差分结果必须除以4才是每个Si的有效电荷。C同理。对于任意低对称结构，通常应每次移动一个待研究原子，不能把多个不等价原子的响应除以数量就当作各自的Born电荷。
下面是 Si-plus/scf/STRU 的完整内容。Si的分数x坐标增加0.002/4.36，C保持不变；末尾0表示此单点输入不做原子位置优化：
ATOMIC_SPECIES
Si 28.0855 Si_ONCV_PBE-1.0.upf
C 12.011 C_ONCV_PBE-1.0.upf
LATTICE_CONSTANT
1.8897261246257702
LATTICE_VECTORS
4.36 0 0
0 4.36 0
0 0 4.36
ATOMIC_POSITIONS
Direct
Si
0
4
0.00045871559633 0.00000000000000 0.00000000000000 0 0 0
0.50045871559633 0.50000000000000 0.00000000000000 0 0 0
0.50045871559633 0.00000000000000 0.50000000000000 0 0 0
0.00045871559633 0.50000000000000 0.50000000000000 0 0 0
C
0
4
0.25000000000000 0.25000000000000 0.25000000000000 0 0 0
0.75000000000000 0.75000000000000 0.25000000000000 0 0 0
0.75000000000000 0.25000000000000 0.75000000000000 0 0 0
0.25000000000000 0.75000000000000 0.75000000000000 0 0 0
运行准备命令会下载并核对固定版本官方示例赝势，将八组输入复制到自己的 work-new：
bash prepare.sh
# 已有相同赝势时，可改为 bash prepare.sh /你的赝势目录
ls work-new/cases
cat work-new/cases/Si-plus/scf/STRU
资料包不附带二进制程序和赝势。若外网下载不通，可按 MODEL.json 中的官方文件名和SHA256准备同一套赝势，不应悄悄替换成另一套赝势再与本页数字比较。
3. 每个结构先做SCF#
scf/INPUT：
INPUT_PARAMETERS
suffix SiCBorn
pseudo_dir ../../../pseudo
ntype 2
nspin 1
nelec 32
nbands 24
basis_type pw
ecutwfc 60
dft_functional pbe
precision double
cal_force 0
cal_stress 0
ks_solver cg
smearing_method fixed
calculation scf
device gpu
symmetry 0
scf_thr 1e-8
scf_nmax 180
mixing_type broyden
mixing_beta 0.3
out_chg 1
out_bandgap 1
scf/KPT：
K_POINTS
0
Gamma
4 4 4 0 0 0
out_chg 1 保存电荷密度，供同一个位移结构的Berry计算读取。固定占据适用于本例绝缘体；若自己的结构变成金属，不要继续把这套绝缘体Berry公式用于金属。
4. 再从同一密度计算Berry极化#
berry/INPUT：
INPUT_PARAMETERS
suffix SiCBorn
pseudo_dir ../../../pseudo
ntype 2
nspin 1
nelec 32
nbands 24
basis_type pw
ecutwfc 60
dft_functional pbe
precision double
cal_force 0
cal_stress 0
ks_solver cg
smearing_method fixed
calculation nscf
device cpu
symmetry -1
init_chg file
pw_diag_thr 1e-10
berry_phase 1
gdir 1
out_chg 0
berry/KPT：
K_POINTS
0
Gamma
8 4 4 0 0 0
这里 gdir 1 是x方向；Berry阶段关闭对称性、读取SCF电荷密度。程序要求此阶段采用CPU执行路径，因此每个已分配GPU作业内先运行GPU SCF，再运行CPU Berry，并非两个不同结构共用密度。
一个结构的核心终端操作如下，资料包的 run.py 会自动执行这些步骤：
cd work-new/cases/Si-plus/scf
mpirun -np 1 abacus > native.out
grep "charge density convergence is achieved" OUT.SiCBorn/running_scf.log
cd ../berry
mkdir -p OUT.SiCBorn
cp ../scf/OUT.SiCBorn/SPIN1_CHG.cube OUT.SiCBorn/
mpirun -np 1 abacus > native.out
grep "P =" OUT.SiCBorn/running_nscf.log
5. 一次提交八组结构#
回到解压后的资料包根目录，普通账号执行：
bash submit.sh
squeue -u "$USER"
提交器为每组结构申请V100 16G的一张卡、1个MPI进程、3个CPU核心和8 GiB内存；线程数为1。它明确指定工作目录，避免从错误目录读取输入。实际编译环境是ABACUS 3.10.1、GCC 13.3、HPC-X 2.50、OpenBLAS 0.3.30与CUDA 12.8，使用 abacus/3.10.1-v100 模块。
八个作业全部结束后，在资料包根目录提取自己的新结果：
bash read.sh work-new analysis-new
bash plot.sh analysis-new
cat analysis-new/Born-Zxx.csv
原始已完成数据和新计算目录分开保存；脚本遇到已经存在的新计算结果会停止，避免覆盖。若只需理解方法，直接读取资料包已有结果即可。
继续阅读：极化分支、差分公式与实算结果。