vasptool极化与压电实操

下载本篇完整输入、实际输出和绘图脚本。解压后执行 bash run.sh，可先从实际输出重建结果和四张图，不会启动VASP。

1 先分清两种响应#

压电计算研究应变引起的极化变化。本例对四原子纤锌矿GaN施加x方向的−1%、0、+1%应变，观察z方向极化，演示proper压电响应e31。

• 固定分数坐标：晶胞形变，原子随晶胞仿射移动，不允许额外内部位移。
• 内部弛豫：每个应变下固定晶胞，优化原子位置，然后在最终结构上重新静态计算极化，包含内部位移贡献。

两者使用同一个已优化原子位置的参考晶胞。这里没有进行晶胞优化，晶格常数a=3.19Å、c=5.189Å；参考体积45.729417ų。优化力阈值0.02eV/Å，电子收敛1E-6。数字用于学习流程，不能当作已收敛GaN材料常数。

2 查看参考结构与静态设置#

module load vasptool/1.5.2
cat inputs/piezo-gan/POSCAR
cat inputs/piezo-gan/INCAR

Ga和N按下面顺序排列，POTCAR也必须为Ga_d、N：

Ga N                                    
   1.00000000000000     
     3.1899999999999999    0.0000000000000000    0.0000000000000000
    -1.5950000000000000    2.7626210380723593    0.0000000000000000
     0.0000000000000000    0.0000000000000000    5.1890000000000001
   Ga   N 
     2     2
Direct
0.0000000000000000    0.0000000000000000    0.0000000000000000
     0.3333333333333357    0.6666666666666643    0.5000000000000000
     0.3333333333333357    0.6666666666666643    0.1229999999999976
     0.0000000000000000    0.0000000000000000    0.6229999999999976

静态INCAR：

SYSTEM = piezo-gan
ENCUT = 500
EDIFF = 1E-6
PREC = Accurate
ISMEAR = 0
SIGMA = 0.05
ISPIN = 1
NSW = 0
IBRION = -1
LREAL = .FALSE.
LWAVE = .FALSE.
LCHARG = .FALSE.
ISIF = 2
EDIFFG = -0.02
LCALCPOL = .TRUE.
ISYM = 0

LCALCPOL=.TRUE.计算Berry相位极化，ISYM=0在这组形变计算中关闭对称性。KPOINTS使用与参考相同的5×5×5 Gamma网格。极化分支分析针对绝缘体；换材料后应先确认各个应变点仍保持绝缘，不把金属体系的输出套入本流程。

3 先生成应变结构#

GaN的INPUT.in与前文单轴Si示例相同：

1
3D
3
-0.01 0.00 0.01
1 0 0 0 0 0

在解压根目录执行：

python3 prepare_native.py
ls prepared/piezo-gan
ls prepared/piezo-gan-relaxed

脚本调用的原生功能是：

vasptool --feature piezo-polarization

它生成对应的三个应变结构，不会替代VASP计算。piezo-gan使用静态INCAR；piezo-gan-relaxed使用以下优化设置：

NSW = 60
IBRION = 2
ISIF = 2
EDIFFG = -0.02

ISIF=2固定当前应变晶胞，只优化原子。若让程序自由改变晶胞，会改变原先指定的应变条件。

4 内部弛豫必须有最终静态步骤#

提交包内模板前先检查账户资源和队列：

cat submit-cpu.slurm
sbatch submit-cpu.slurm
squeue -u "$USER"

模板对内部弛豫系列自动执行：优化 → 保存optimization目录 → 复制CONTCAR为POSCAR → 改为NSW=0、IBRION=-1 → 最终Berry静态计算。关键文件操作是：

mkdir optimization
cp POSCAR INCAR CONTCAR OUTCAR OSZICAR optimization/
cp CONTCAR POSCAR

随后恢复静态INCAR并用正式作业模板运行VASP。不要仅从结构优化OUTCAR里出现的偶极值就判断它对应最终结构。本流程明确用最终静态OUTCAR提取极化，也检查最终力和电子收敛。

新材料先优化参考结构的原子位置，再生成应变系列；可参考 inputs/gan-reference-relax/INCAR。本包参考结构已经完成这一步，不能把它当作其他材料的通用优化结果。

5 提取偶极、统一分支并画图#

全部结束后：

bash analyze_results.sh
cat prepared/piezo-gan/polarization-independent.csv
cat prepared/piezo-gan-relaxed/polarization-independent.csv

脚本先把INPUT.in第一行临时改为2，调用原生piezo-polarization，提取两个表：

DIPOLE_MOMENTS_ELETRONIC.dat
DIPOLE_MOMENTS_IONIC.dat

其中ELETRONIC是当前程序的实际文件名。然后从原始OUTCAR分别读取 p[elc] 和 p[ion]，只进行一次eÅ到C/m²的换算。本机441摘要对这组数据出现重复换算，附件不直接采用该摘要斜率。

原子越过周期边界时，同一个物理结构可能改变原始离子偶极的晶格分支。本包按参考结构的最近周期像和OUTCAR记录的Ga13/N5价电荷检查离子位移，确认原始值与连续值相差整数晶格量子；原始数据保留在JSON中，不覆盖OUTCAR。电子Berry分支也必须在小应变序列中保持连续，不能用任意跳变拟合。

对于形变矩阵F，本例使用连续约化偶极 F⁻¹p/V₀ 的应变导数获得proper响应，避免直接对变化晶胞中的P=p/V求导引入几何项。这里用的是小正应变，不包含剪切系数的工程应变约定；不能直接复制为所有张量分量的通用脚本。方法见Vanderbilt关于proper压电响应的说明。

6 本例得到什么#

这里的正负号对应本包的晶胞取向与极性方向，和其他资料比较前必须统一坐标、极性与proper定义。内部弛豫后，离子项改变了总响应；只取电子项会得到另一种量。

两组非零应变最终静态最大力分别为0.007459和0.016999eV/Å，均在本例0.02阈值内。三个点、1%应变以及这个力阈值只用于操作教学；精确压电常数需要更紧的力、更多小应变点，以及k网格和截断能收敛。拟合应查看残差，不能因程序输出一个数就认为材料常数可靠。

最后得到 GAN-VALIDATION.json、GAN-RELAXED-VALIDATION.json、两份CSV和两张GaN图；JSON同时保留原始偶极、连续偶极、晶格分支修正、最终力和拟合结果。