pymatgen弹性张量与应力应变
从硅的实际应变数据拟合三项立方弹性常数，计算柔度和多晶平均，辨别剪切因子与单位。
阅读本篇目录
本篇接着VASP弹性实算教程，用同一批20个硅计算结果练习 pymatgen 2026.9.24 的弹性分析。固定离子和弛豫离子两组分别处理，不能混入同一张量。
1. 下载与运行#
下载实际输入数据、脚本和绘图包，解压进入 pymatgen-elastic，在已分配的CPU作业内执行：
bash run.sh
也可用包内 sbatch job.slurm 提交；CPU分区按56核整节点分配和计费，分析只用一个线程。脚本读取已有应力和结构，不会启动VASP；Python路径已指向集群现有pymatgen环境。
输入为 inputs/actual-stress-strain.json 和实际未变形的8原子 inputs/Si-POSCAR。20条数据是两种内部原子处理方式×两种应变模式×五个应变点；每条保留结构标签、应变、六分量应力和能量。
head -n 28 inputs/actual-stress-strain.json
cat inputs/Si-POSCAR
2. 先统一应力、剪切和晶体近似#
本包应力已经由ASE转换为拉伸为正的GPa，顺序为 xx,yy,zz,yz,xz,xy，不要再按VASP原始OUTCAR反号。直接读OUTCAR时，压缩为正的kbar需要另外处理，10 kbar=1 GPa。
工程剪切γ对应对称小应变εxy=εyx=γ/2。pymatgen的 Strain.voigt 对剪切使用工程分量，Stress.voigt 的剪切应力不乘2。
这里沿原教程采用立方近似，从xx的σxx、横向σyy/σzz和xy的剪切σxy分别拟合C11、C12、C44。实际晶胞有微小非立方残差；这不是对任意晶系21个独立弹性常数的完整拟合，不能把未采样分量当作独立算过的结果。拟合保留截距，以免把残余应力直接除以应变。
3. 建立ElasticTensor并计算平均模量#
核心形式如下，完整脚本会从输入实际拟合，而不是硬编码这三个近似值：
import numpy as np
from pymatgen.analysis.elasticity.elastic import ElasticTensor
C11, C12, C44 = 153.14954474, 56.73880305, 72.78602494
C = np.zeros((6, 6))
C[:3, :3] = C12
np.fill_diagonal(C[:3, :3], C11)
C[3:, 3:] = np.eye(3) * C44
elastic = ElasticTensor.from_voigt(C)
print(elastic.k_vrh, elastic.g_vrh, elastic.y_mod / 1e9)
compliance_tensor.voigt 为C的逆，单位GPa⁻¹。B、G和Cij以GPa表示；y_mod 返回Pa，输出GPa必须除以10⁹。
cat results/moduli.csv
| 量 / GPa | 固定离子 | 弛豫离子 |
|---|---|---|
| C11 | 153.150 | 153.150 |
| C12 | 56.739 | 56.739 |
| C44 | 99.163 | 72.786 |
| Hill体积模量 | 88.876 | 88.876 |
| Hill剪切模量 | 74.237 | 61.704 |
| Hill杨氏模量 | 174.207 | 150.324 |
左图是原始DFT点，右图比较三项重建常数。C44的变化展示内部原子弛豫的作用；本例不能推导成“所有材料的弛豫修正都相同”。
4. 应变能、有限变形与重新提取#
calculate_stress 按线性Hooke关系计算应力，energy_density 默认给eV/ų；选择 convert_GPa_to_eV=False 时给等价GPa能量密度。小应变使用对称ε，有限变形的Green–Lagrange应变则为 (F.T @ F - I)/2，两者不能在大变形时直接混用。
/public/toolkit/miniconda3/envs/pymatgen-2026.9.24/bin/python -B analyze.py
完整张量、柔度和拟合数组保存在 results/full-fields.npz，新图在 plots/。更换材料时准备与对称性匹配的完整应变集合，重新检查k点、应变幅度和内部力收敛，不能只替换材料名。下一篇继续算方向模量与声速；接口见官方弹性模块。