ABACUS Berry极化与连续路径实操

下载完整输入、赝势、原生结果与绘图脚本，解压进入附件根目录。先查看已有数据：

bash analyze.sh

图和表在reference/figures与reference/analysis。从头运行则在计算资源环境中提交：

sbatch submit.slurm

附件run.sh会先建立独立work目录，按前置SCF→后续计算顺序运行，再生成work/analysis与work/figures。已有work时请新解压一份，避免覆盖原结果。模板申请1张P100，SCF使用GPU，Berry/Wannier接口指定CPU；每步1个MPI进程、1个线程。接口阶段仍占所申请的资源，长批次可拆分提交。你也可以只选一个结构手动计算。

ABACUS为3.10.1，模块abacus/3.10.1-p100，GNU13.3、HPC-X2.50/OpenMPI5.0.10rc2环境；本模块CUDA12.8，接口CPU路径不代表GPU加速。独立Wannier90为4.0.3串行版，Intel2024.2构建，单独子环境运行，避免与ABACUS的MPI环境混装；后处理使用ASE3.29.0及NumPy。

1 先确定你要比较的两个结构#

本例使用5原子的PbTiO3固定晶胞：a=b约3.837Å、c约4.787Å。先建立中心对称参考，再沿同一个结构路径逐步移动Ti和O，λ=0、0.25、0.50、0.75、1。这个拉长晶胞来自教学输入，没有在本篇做结构优化；所得数值不作为平衡晶体的实验自发极化。

λ=1的实际STRU如下，原子位置为分数坐标：

ATOMIC_SPECIES
Pb 207.2 Pb_ONCV_PBE-1.0.upf
Ti 47.867 Ti_ONCV_PBE-1.0.upf
O 15.9994 O_ONCV_PBE-1.0.upf

LATTICE_CONSTANT
1.889725988579

LATTICE_VECTORS
3.8369998932         0.0000000000         0.0000000000
0.0000000000         3.8369998932         0.0000000000
0.0000000000         0.0000000000         4.7870001793
	
ATOMIC_POSITIONS
Direct

Pb
0.0
1
 0.000000000         0.000000000         0.000000000 0 0 0


Ti
0.0
1
 0.500000000         0.500000000         0.555999975 0 0 0

 
O
0.0
3
 0.500000000         0.500000000         0.194999994 0 0 0
 0.000000000         0.500000000         0.671000004 0 0 0
 0.500000000         0.000000000         0.671000004 0 0 0

中心对称参考中Ti的z=0.5、轴向O的z=0、另外两个O的z=0.5；其余坐标相同。附件保留全部五个结构、镜像结构及各自SCF，不能把一个结构的密度直接交给另一个结构。

2 每个结构先自洽#

本例完整SCF输入：

INPUT_PARAMETERS
pseudo_dir ../../PP_ORB
ntype 3
ecutwfc 50
scf_nmax 200
symmetry 1
scf_thr 1e-6
calculation scf
basis_type pw
out_chg 1
orbital_dir ../../PP_ORB
device gpu
ks_solver cg
suffix tutorial
precision double
cal_force 0
cal_stress 0
out_bandgap 1

K_POINTS
0
Gamma
4 4 4 0 0 0

scf_thr 1e-6是本例电子自洽阈值；本篇固定结构，没有离子优化步骤。先看SCF原生输出中的收敛提示和E_bandgap。五个路径点的采样KS能隙分别约1.2243、1.4904、1.5967、1.7738、2.0468eV，均为绝缘态；这些是当前k网格结果，不是基本带隙的精密计算。

3 读取相同结构的密度，计算Berry相位#

λ=1的完整Berry输入与k网格：

INPUT_PARAMETERS
suffix tutorial
pseudo_dir ../../PP_ORB
orbital_dir ../../PP_ORB
ntype 3
basis_type pw
calculation nscf
ecutwfc 50
symmetry -1
device cpu
precision double
init_chg file
pw_diag_thr 1e-10
nbands 32
smearing_method fixed
berry_phase 1
gdir 3
out_chg 0
cal_force 0
cal_stress 0

K_POINTS
0
Gamma
2 2 8 0 0 0

关键设置：

• calculation nscf、init_chg file：读取同一结构SCF密度。
• smearing_method fixed：本例为绝缘体固定占据，不沿用金属高斯展宽。
• berry_phase 1、gdir 3：沿第三个实空间晶格矢计算。本例为正交晶胞，对应z方向。
• gdir取1/2/3；电场教程的efield_dir取0/1/2且指倒格矢方向，两者不能按同一编号复制。
• 第三方向每串8点，横向2×2；完整均匀网格与能带高对称路径不是一回事。

已有前置SCF后，可在附件根目录执行下面这段实际操作：

cd work/cases/pbt-berry-1.00
mkdir -p OUT.tutorial
cp ../pbt-polar-plus-scf/OUT.tutorial/SPIN1_CHG.cube OUT.tutorial/
module purge
module load abacus/3.10.1-p100
mpirun -np 1 abacus > berry.log

运行后查看OUT.tutorial/running_nscf.log内的VALUES OF POLARIZATION，包括离子相位、电子相位和三种单位。C/m^2一行的括号里是笛卡尔分量。

4 按连续路径处理分支#

不能只减去两个结构屏幕上打印的极化数。当前3.10.1无自旋、成对占据路径打印的周期为2.1748536C/m²；通常定义的单电子极化量子eR/Ω在本例为约1.0874268C/m²。原生实现的自旋成对相位周期与常见单电子量子要区分。下表按原生周期逐步连续展开；每个相邻变化约0.3C/m²，也小于常见单电子量子的一半，所以本路径不会因这两个周期约定而选出不同连续变化。

因此这条固定结构路径的相对极化变化为−1.2731022C/m²；λ=1原生打印的+0.9017514C/m²属于另一个等价分支。若路径出现闭合能隙、相邻变化接近半个量子，应增加结构路径点并重新判断，不能让脚本猜出唯一物理解。

读取附件已经展开的数据：

python - <<'PY'
import csv
with open('reference/analysis/polarization-path.csv') as f:
    for r in csv.DictReader(f):
        print(r['lambda_path'],r['P_raw_Cm2'],r['delta_P_Cm2'])
PY

5 做方向和k点控制#

镜像极性结构与正向结构总能相差约1.76×10⁻⁸eV，极化相加在当前分支约−3.5×10⁻⁶C/m²；横向两个分量均约10⁻⁷C/m²。将z串点从8改为16，原生端点由0.9017514变为0.9001721C/m²，变化−0.0015793C/m²。横向仍是2×2，还不能称为完整布里渊区收敛。

ABACUS使用的旧版单位常数与今天某些SI常数略有差别，且原生输出仅保留有限小数；核对原生相位和单位时应允许输出精度差异，不要人为平移或改写数据。

参阅官方Berry说明及本例版本原生实现。