py4vasp 压电与内部应变

压电张量、内部应变耦合和晶胞的应变不是同一个量。本篇用Si、NaCl与ZnS三个实际DFPT例子读取前两者；其中Si、NaCl的压电为零，ZnS提供非零实例。

1. 获取数据与选择操作#

下载py4vasp数据、输入与绘图实例包：

unzip py4vasp-quantities-wave151-*.zip -d py4vasp-examples
cd py4vasp-examples
bash run.sh list

后处理模块为 py4vasp/0.11.3，使用独立Python3.11环境；结构和静态图使用 ase/3.29.0 的Python3.11.16、NumPy2.3.5、Matplotlib3.11.2。只需已分配资源内1个CPU，不需要GPU。

bash run.sh 1    # 39个数据入口、版本和文件前提
bash run.sh 2    # 介电、Born、压电和内部应变
bash run.sh 3    # 电子收敛、带隙及CSV

也可以在新解压目录运行 bash run-all.sh，得到 analysis/ 与 figures/。另一份新目录运行 bash analyze.sh，结果写入 analysis-reference/、figures-reference/。这两条命令只处理已经提供的数据，不重新运行VASP。输出目录保护会拒绝覆盖已有结果。

8份输入包含两类实际计算：5份已有VASP6.6.1/P100 HDF5结果（NVHPC25.3、CUDA12.8），以及3份Si/NaCl/ZnS的CPU DFPT结果（vasp/6.6.1-8173m-hdf5，InteloneAPI2024.2、MPI/OpenMP、HDF51.14.6）。新CPU计算采用1个MPI进程×4个OpenMP线程；读取HDF5与计算平台分开，py4vasp本身不替代VASP求解。

2. 为什么再准备一个ZnS原胞#

ZnS使用闪锌矿两原子原胞，常规晶格参数5.45Å。它在这里不是纤锌矿结构。完整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

完整INCAR：

SYSTEM = zns-response
ENCUT = 550
PREC = Accurate
EDIFF = 1E-6
EDIFFG = -0.02
IBRION = 8
NSW = 1
ISIF = 2
ISMEAR = 0
SIGMA = 0.05
LREAL = .FALSE.
LASPH = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
NELM = 150
ALGO = Normal
ISYM = 2
LORBIT = 11
NBANDS = 24
LH5 = .TRUE.
LEPSILON = .TRUE.
LRPA = .FALSE.
NPAR = 1

Small response tutorial mesh; not a convergence study
0
Gamma
4 4 4
0 0 0

用静态响应教程中的赝势准备与 bash run-native.sh 可实际重新求解。只做提取则直接 bash run.sh 2。本篇是固定结构、4×4×4采样的操作例，未给出收敛的压电材料常数。

3. 先读取完整数组，不急着挑一个数字#

module load py4vasp/0.11.3
python - <<'PYCODE'
import py4vasp
c = py4vasp.Calculation.from_path("inputs/zns-response")
p = c.piezoelectric_tensor.read()
i = c.internal_strain.read()
print("piezo electronic shape:", p["clamped_ion"].shape)
print("piezo total shape:", p["relaxed_ion"].shape)
print("internal strain shape:", i["internal_strain"].shape)
print("atom order:", i["structure"]["elements"])
print(c.piezoelectric_tensor)
print(c.internal_strain)
PYCODE

压电的两个数组都是 (3,3,3)；前两个轴是应变指标，最后一个轴是电场/极化方向。clamped_ion为电子项，relaxed_ion为电子加离子总项。单位是C/m²。

内部应变数组是 (原子数,3,3,3)：依次为原子、位移方向、应变的两个指标；它表示能量对原子位移与应变的混合导数，单位eV/Å。它不是一个“实际应变大小”数组，也不能直接当成压电系数。

4. 按OUTCAR顺序导出3×6表#

VASP表格按 XX YY ZZ XY YZ ZX 排列，这与一些软件采用的 XX YY ZZ YZ ZX XY 顺序不同。用明确的索引转换：

python - <<'PYCODE'
import numpy as np
import py4vasp
c = py4vasp.Calculation.from_path("inputs/zns-response")
pairs = [(0,0),(1,1),(2,2),(0,1),(1,2),(2,0)]
p = c.piezoelectric_tensor.read()
for name in ["clamped_ion", "relaxed_ion"]:
    t = p[name]
    table = np.array([
        [0.5*(t[i,j,k]+t[j,i,k]) for i,j in pairs]
        for k in range(3)
    ])
    np.savetxt(f"piezo-{name}.csv", table, delimiter=",")
t = c.internal_strain.read()["internal_strain"]
for atom, tensor in enumerate(t, 1):
    table = np.array([
        [0.5*(tensor[d,i,j]+tensor[d,j,i]) for i,j in pairs]
        for d in range(3)
    ])
    np.savetxt(f"internal-strain-ion{atom}.csv", table, delimiter=",")
PYCODE

压电CSV的行是x/y/z电场方向；内部应变CSV的行是该原子的x/y/z位移方向。不要仅因为它们都能写成3×6表，就把行和单位当成相同。

5. 这三个实际例子说明什么#

Si和NaCl的压电张量为零；这是有实际张量数据并符合对称性的零，和入口索引里的“空返回/没有计算数据”不同。

Si第一原子的内部应变剪切耦合约8.423eV/Å，虽然它的压电为零。这正说明内部位移–应变耦合不能直接替代电极化–应变耦合。

ZnS的x电场–YZ应变分量，本例电子项为 -0.803108 C/m²，电子加离子总项为 0.027621 C/m²。两项会明显抵消，不能把 clamped_ion 单独当成总压电，也不要把 relaxed_ion 误当成纯离子项。坐标轴、原子顺序与结构取向改变后，分量及符号也需按相同约定比较。

6. 最后检查哪些文件#

ls response
cat response/zns-response.json
# 若执行的是全包后处理，改看：
ls analysis/response

配套JSON保留3×6压电电子/总项和每个原子的内部应变表；原始OUTCAR也在各数据目录中。用相应单位的最终累计表对照，不取stdout中单次扰动的中间数字。

这些具体结构的对称性检查只适用于本例。处理自己的低对称材料时，直接用上面的通用读取/索引代码，并按实际结构重新选择应检查的对称关系。

接口参见py4vasp calculation目录；相关标签参见VASP LEPSILON。