pymatgen：压电张量与响应分解

压电响应往往包含电子和离子两部分，两者可能明显抵消。本例复用鸢算科技集群已经完成的ZnS DFPT结果，用pymatgen读取和转换张量，再重建离子分量。

前置数据使用VASP6.6.1 CPU版本、oneAPI2024.2生成，EDIFF=1E-6；本页后处理不会重新算DFPT。它是小体系操作参考，没有宣称精密材料常数或完整收敛研究。

下载实际响应输入、数据、分析与绘图脚本。

unzip pymatgen-piezo-thermodata-wave722-*.zip
cd pymatgen-piezo-thermodata
module load pymatgen/2026.9.24
bash results.sh
bash plot.sh

后处理环境为pymatgen2026.9.24、Python3.11、NumPy/SciPy与Matplotlib。run.sh在新建的analysis重新读取已有数据，不提交DFT；原数据在reference-data，绘图输出到figures。

1. 先明确实际输入#

inputs/INCAR#

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

inputs/POSCAR#

S  Zn
 1.0000000000000000
     0.0000000000000000    2.7250000000000001    2.7250000000000001
     2.7250000000000001    0.0000000000000000    2.7250000000000001
     2.7250000000000001    2.7250000000000001    0.0000000000000000
 S   Zn 
   1   1
Direct
  0.2500000000000000  0.2500000000000000  0.2499999999999999
  0.0000000000000000  0.0000000000000000 -0.0000000000000000

inputs/KPOINTS#

Gamma mesh
0
Gamma
7 7 7
0 0 0

IBRION=8和LEPSILON开启响应计算，NSW=1在这里不是做普通结构优化。结构元素顺序为S、Zn；下载包不含POTCAR。

若需要重算，在已有授权势函数环境准备并在有可用CPU资源时提交：

bash prepare-potcar.sh
cd inputs
sbatch job.slurm

模板使用vasp/6.6.1-8173m CPU模块。CPU分区未开放时，先使用包内数据学习后处理，等待适用CPU资源再提交。两种CPU/GPU构建的功能范围不能仅按同一版本号推断。

2. VASP与标准Voigt顺序不同#

VASP的6列为XX、YY、ZZ、XY、YZ、ZX；pymatgen张量的标准Voigt顺序为xx、yy、zz、yz、xz、xy。用专门转换接口，避免手动把列号当成物理分量：

python - <<'PY'
import json
from pymatgen.analysis.piezo import PiezoTensor
with open('reference-data/ZnS-response-inputs.json') as f:
    d = json.load(f)
tensor = PiezoTensor.from_vasp_voigt(d['electronic_vasp_voigt_C_m2'])
print('3x3x3完整张量:', tensor.shape)
print('e(x,yz):', tensor[0,1,2])
print('标准Voigt:', tensor.voigt)
PY

张量完整形式的后两个指标对称。本例e(x,yz)的VASP位置是x行、YZ列，电子分量为−0.69510C/m²。

当前VASP输出的最终电子压电表含including local field effects字样；本例读取这个最终对称化C/m²表。不能把较早版本读器没有匹配到表头时留下的全零数组当作实际电子响应。

3. Born、内部应变和力导数各是什么#

下载包的实际OUTCAR保存：

• Born有效电荷Z*，单位|e|；Zn约+2.06282、S约−2.06282。
• 内部应变响应，打印为每个原子的3×6表，单位eV/Å。
• 原子位移的力导数矩阵，单位eV/Ų。

Born有效电荷不是Bader/Mulliken静态电荷。内部应变表也需要按VASP剪切列转换成完整张量。

本例打印的未对称化力导数存在约0.043562eV/Ų的非对称残差。后处理保留原矩阵，再取对称部分并投影掉三个刚体平移方向；对这个两原子胞，得到三个声学零模和三个正的物理刚度模。它是已有响应矩阵的对称/声学和处理，不是新DFPT结果。

4. 原生重建与体积换算#

python - <<'PY'
import json, numpy as np
from pymatgen.core import Structure
from pymatgen.analysis.piezo_sensitivity import get_piezo
with open('reference-data/ZnS-response-inputs.json') as f:
    d = json.load(f)
s = Structure.from_dict(d['structure'])
Z = np.array(d['Born_charge_e'])
L = np.array(d['internal_strain_eV_A'])
F = np.array(d['force_derivative_symmetric_ASR_eV_A2'])
blocks = F.reshape(2,3,2,3).swapaxes(1,2)
ionic = get_piezo(Z,L,blocks)/s.volume
print('体积/A^3:', s.volume)
print('离子e(x,yz)/C m^-2:', ionic[0,1,2])
PY

这里F沿用VASP力导数符号，物理刚度为−F。原生函数对其做适用的广义逆，去掉声学方向。按本例输入单位，还须显式除以晶胞体积40.46965625ų，才能与VASP打印的C/m²分量比较。

不要直接把Vasprun.force_constants不加单位检查地传进来：该属性可来自XML中带质量/频率单位的dynamical-matrix表，本例明确使用OUTCAR的eV/Ų力导数。

两个较大数值抵消后总量很小，不能从最后几位差别推断新的材料效应。这里的重建与打印值相差约8.4×10⁻⁷C/m²。

5. 对称约束和敏感性操作#

python - <<'PY'
import json, numpy as np
from pymatgen.core import Structure
from pymatgen.symmetry.site_symmetries import get_site_symmetries
from pymatgen.analysis.piezo_sensitivity import BornEffectiveCharge
with open('reference-data/ZnS-response-inputs.json') as f:
    d = json.load(f)
s = Structure.from_dict(d['structure'])
pointops = get_site_symmetries(s,precision=0.01)
bec = BornEffectiveCharge(s,np.array(d['Born_charge_e']),pointops)
print([(a,b,len(operations)) for a,b,operations in bec.get_BEC_operations()])
PY

完整脚本还实际调用内部应变和力导数的对称关系接口，并检查对已经满足声学和的矩阵再次投影时结果不变。

敏感性示例把Z*、内部应变或力刚度单独乘0.9/1/1.1，再调用同一个重建函数。前两种使离子响应按比例变化，刚度增大则使响应按反比变化；这只是算术控制，没有假定真实材料能独立改变这些张量而其他性质保持不变。

cat reference-data/ZnS-piezo-components.csv
cat reference-data/ZnS-piezo-algebraic-sensitivity.csv

若使用模块的随机张量工具，生成的是受给定约束的模型张量，不是新DFPT预测。真实材料研究仍应核对结构、张量来源、对称性、声学和、刚度稳定性及参数收敛。

接口见pymatgen压电响应模块，前置计算说明见VASP极化、介电与压电。