GPUMDkit局域位移与极化实操

本篇把BaTiO₃结构依次变成邻居表 → 局域位移 → 八面体偏角 → 近似局域极化 → 空间矢量图。输入来自鸢算科技VASP实际计算，Born电荷由同一参考结构的DFPT获得。9帧文件是静态软模路径，不能据此计算有限温度涨落或介电常数。

下载输入、结果与绘图脚本。解压后在附件根目录操作；不需要安装额外依赖：集群gpumdkit/1.5.8已配套ferrodispcalc0.1.0。

1. 准备工作目录和结构#

module load gpumdkit/1.5.8
gpumdkit.sh -doctor
python new-case.py ferro
cd work/ferro
head -2 model.xyz

model.xyz为3×3×3复制的135原子参考结构；movie.xyz含9帧，每帧27个Ti中心。原始5原子晶胞如下，晶格常数4.02Å，顺序Ba、Ti、O：

BaTiO3 cubic reference
1.0
4.02 0 0
0 4.02 0
0 0 4.02
Ba Ti O
1 1 3
Direct
0.0 0.0 0.0
0.5 0.5 0.5
0.5 0.5 0.0
0.5 0.0 0.5
0.0 0.5 0.5

这条路径按软模方向共同移动原子，u=0.06Å是路径振幅，不是单独将Ti平移0.06Å。其Ti相对氧笼中心的实际位移为0.07566059Å。各原子精确坐标见附件polar-path/frame-0至frame-8的POSCAR。

2. 建立两份邻居表#

gpumdkit.sh -calc nlist -i model.xyz -c 2.7 -n 6 -C Ti -E O
gpumdkit.sh -calc nlist -i model.xyz -c 3.7 -n 8 -C Ti -E Ba
head -3 nl-Ti-O.dat
head -3 nl-Ti-Ba.dat

-C指定中心Ti，-E指定邻居元素，-c为Å单位搜索半径，-n为每个中心的邻居数。本例第一份27行×7列，第二份27行×9列；首列为Ti编号，其余是邻居编号，编号从1开始。

这组半径只用于本例。自己的结构应检查每个Ti确实连接6个O和8个Ba；不能用加大半径强凑配位数。后续所有帧必须保持相同原子顺序，重排原子后应重建邻居表。

3. 提取位移和八面体偏角#

gpumdkit.sh -calc disp -i movie.xyz -n nl-Ti-O.dat -o displacements.dat
gpumdkit.sh -calc oct-tilt -i movie.xyz -n nl-Ti-O.dat -o tilt.dat
head -3 displacements.dat
head -3 tilt.dat

位移输出单位Å，每行x、y、z三个分量。243行按“9帧×每帧27个Ti”排列。它表示Ti坐标减去六个O的周期最小镜像平均坐标，不能直接当极化。

这条软模路径只发生极性位移，O–O轴不倾斜，因此tilt.dat为0°。该工具输出O–O轴相对于笛卡尔轴的无符号夹角，不等同于带正负号的Glazer转角或Euler角。作为几何对照，绕z整体转10°的刚性八面体得到(10,10,0)°，而非(0,0,10)°。跨周期边界的原子按最小镜像处理。

4. Born电荷从哪里来#

本次使用VASP6.6.1普通CPU构建与oneAPI2024.2，在5原子参考晶胞执行DFPT。附件born-dfpt保留输入和OUTCAR。主要INCAR如下：

ENCUT = 550
PREC = Accurate
EDIFF = 1E-6
EDIFFG = -0.02
IBRION = 8
NSW = 1
ISIF = 3
ISMEAR = 0
SIGMA = 0.05
LREAL = .FALSE.
LASPH = .TRUE.
NELM = 150
ALGO = Normal
ISYM = 2
LCALCEPS = .FALSE.
LEPSILON = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.

Gamma mesh
0
Gamma
7 7 7
0 0 0

这里IBRION=8是响应计算，不是结构优化；EDIFFG不是本次停止准则。自己的体系先完成结构优化，再进行响应计算。使用个人有授权的POTCAR；程序不会由几何文件自动生成Born电荷。

查看实际张量块：

cd ../../born-dfpt
sed -n '/BORN EFFECTIVE CHARGES/,+24p' OUTCAR
cd ../work/ferro

本例Ba对角元2.73777，Ti对角元7.52747；沿Ti–O键的O对角元−5.98674，垂直方向−2.13925。pol-abo3当前每元素只接收一个标量，因此本例用O张量迹的平均−3.421746666666667，并满足Ba+Ti+3O=0。标量平均丢失O电荷的方向差异；形式价态Ba²⁺、Ti⁴⁺、O²⁻不能替代这些动态有效电荷。

5. 算近似局域极化并作对照#

gpumdkit.sh -calc pol-abo3 -i movie.xyz \
  --nl-ba nl-Ti-Ba.dat --nl-bo nl-Ti-O.dat \
  --bec Ba=2.73777 Ti=7.52747 O=-3.421746666666667 \
  -o polarization.dat
head -3 polarization.dat

每行x、y、z，单位C/m²；每个Ti局域单元分配1个Ti、八个Ba各1/8、六个O各1/2。本例各单元等价，以下列出z方向平均值：

所有结果以中间非极性参考帧为零。完整张量线性估计采用ΔPα=e/Ω Σκβ Z*καβ Δuκβ；附件polar_compare.py可从OUTCAR和各POSCAR重算。Berry结果来自这些静态结构的实际VASP计算，并沿连续路径处理极化分支。振幅增大时，参考构型的线性Born响应与非线性Berry结果逐渐偏离；不能挑选其中更接近预期的一列作为“正确极化”。研究畴结构、非等价氧或大畸变时，应保留完整张量并核对参考态与分支。

6. 输出空间矢量图#

gpumdkit.sh -get_frame movie.xyz 8
gpumdkit.sh -plt plane-grid -i frame_8.xyz -d displacements.dat \
  -e Ti -g 3 3 3 -f 7 \
  --select-xy 0 --select-xz 0 --select-yz 0 -o plot
ls plot

先用从1计数的get_frame提取第8帧，再用从0计数的-f 7选同一帧位移（本例u=+0.18Å）。当前绘图脚本的-f只选择位移数据，直接传多帧结构会默认读取末帧；显式提取单帧可避免坐标与位移错配。选择该帧和各方向第0层，输出XY_0.png、XZ_0.png、YZ_0.png。这是一组沿z的位移，所以XY平面内箭头接近零，XZ/YZ更直观。-g必须匹配Ti网格；任意无序结构不能强行设成3×3×3。

回到附件根目录，重绘本篇两幅比较图及离子电导图：

cd ../..
python plot_results.py

进一步阅读：GPUMDkit极性分析入口、ferrodispcalc上游。