pymatgen：应变系列的proper压电响应

本例接着Berry分支处理。讨论的是x应变引起的z极化变化，即e31；不是电场响应、Born电荷、结构优化能垒或自发极化。

下载输入、真实参考表与操作脚本。

unzip pymatgen-polarization-wave667-*.zip
cd pymatgen-polarization
bash results.sh
bash plot.sh

reference-data/保存实际计算的数值，figures/reference.png为重绘图；读取和绘图不会启动VASP。若要运行本篇全部后处理，普通账号提交sbatch job.slurm。脚本申请cpu-8173m独占整节点1任务×56核，小例内部单线程，计费仍按整节点；不要在管理节点启动计算。

1. 选择参考晶胞与正确单位#

参考体积V0=45.7294174073 ų，a=3.19 Å，c=5.189 Å。x方向应变改变体积，z方向晶格矢量保持不变。直接把每点偶极都除以各自体积再求斜率，会混入体积变化的几何项；本例采用参考体积处理z分量的proper响应。

import numpy as np
from scipy.constants import elementary_charge
# same[:,2]来自原生分支处理，单位是电子电荷e<0乘Å。
# 本例z正交且c不变；其他形变需先做完整的F^-1向量变换。
V0 = structures[1].volume
Pz = -elementary_charge * 1e20 * same[:,2] / V0
strain = np.array([-0.01, 0.0, 0.01])
delta_Pz = Pz - Pz[1]
e31 = np.polyfit(strain, delta_Pz, 1)[0]
print('e31 / C m^-2:', e31)

1e20来自Å/ų到m⁻²。这里负号对应pymatgen已经转成的电子电荷约定；如果直接读取OUTCAR的|e|Å打印值，则需要匹配那套符号，不能混用两种输入和同一个换算符号。

对一般斜晶胞或包含z伸缩/剪切的形变，先恢复笛卡尔向量，再使用形变梯度F^-1映射到参考晶胞，不能照搬“仅除以V0”的本例简式。

2. 固定分数坐标与内部弛豫#

两套结果都来自已有真实DFT数据，本次只换用pymatgen完成分支与单位处理。数值与原GaN压电实例一致。符号不同来自内部坐标响应，并非把一套数据乘常数构造另一套。

固定分数坐标时，离子随晶胞仿射移动，不包含额外内部位移；内部弛豫时，每个应变下固定晶胞、优化原子，再在最终结构上计算Berry极化。电子收敛1E-6，自由原子力0.02 eV/Å。本例只有三个±1%点，参数和应变区间尚不构成全面收敛的材料常数。

3. 看结果表里的每一列#

head reference-data/GaN-pymatgen-polarization.csv
cat reference-data/GaN-pymatgen-piezo-response.csv

第一张表保留原始读入偶极、连续分支、整数分支位移、参考体积极化、相对变化与逐点体积极化。第二张表给出两套e31。不要删除分支整数记录后只留最终斜率，否则以后难以核对周期性跳变。

Polarization还提供原生μC/cm²转换和all_in_polar映射；本篇刻意保留eÅ中间量，再显式转换到C/m²和选定参考晶胞，因此没有把原生默认终态晶胞口径与这里的零应变参考体积混为一谈。无论选择哪种接口，都要把参考晶胞、单位、坐标基底和分支约定说清楚。