VASP 弹性张量的方向模量与SI导出

下载实际DFT数据、导出脚本和参考图，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd vasp-multiscale-reference-examples
bash results.sh
bash plot.sh

本包复用本站已经完成的VASP数据，只做参照、单位和张量方向的整理，没有新增DFT或有限元计算。分析使用NumPy和Matplotlib，集群 ase/3.29.0 环境已提供所需依赖。

1. 准备材料张量#

输入 inputs/Si-source-material.json 来自Si有限应变弹性实例。原始DFT使用PBE PAW64、550 eV、EDIFF=1E-6，应变点固定晶胞、内部原子弛豫到0.02 eV/Å；VASP6.6.1的P100构建采用NVHPC25.3、CUDA12.8及配套MPI。

本次采用弛豫离子张量，C11=153.14954474 GPa、C12=56.73880305 GPa、C44=72.78602494 GPa，密度2280.906 kg/m³。输入JSON中的矩阵已经乘10⁹转为Pa，不能再乘一次。

cat inputs/Si-source-material.json
bash run.sh
head -n 12 results/Si-directional-compliance.csv

results已有时脚本拒绝覆盖。新输出还包括 Si-material-SI.json，可传给你后续采用同一约定的连续介质模型。

2. 工程剪切约定要成对使用#

本例Voigt顺序为 xx yy zz yz xz xy，应力向量的剪切分量是σyz等；应变向量则是 εxx εyy εzz 2εyz 2εxz 2εxy。计算柔度 S=C⁻¹后，要把工程剪切应变除以2才能拼回3×3应变张量。

S = np.linalg.inv(C)  # C in Pa
nx,ny,nz = n
stress = np.array([nx*nx,ny*ny,nz*nz,ny*nz,nx*nz,nx*ny])
strain = S @ stress
epsilon = np.array([
    [strain[0],strain[5]/2,strain[4]/2],
    [strain[5]/2,strain[1],strain[3]/2],
    [strain[4]/2,strain[3]/2,strain[2]]])
E = 1.0 / (n @ epsilon @ n)

这里对单位应力 n⊗n求轴向应变，所以E的单位为Pa。n必须是单位向量。换成其他软件的Voigt顺序或剪切定义时，不能只把矩阵原样复制。

3. 实际方向结果#

本次在(001)平面取 n=(cosθ,sinθ,0)，θ从0到90°、每5°一个点，共19个方向。这些是由同一真实DFT张量计算出的方向响应，没有为每个方向重新运行VASP。

结果同时与立方晶体的独立解析公式比较。立方对称使[100]与[010]等价，但并不要求[110]的模量相同。把一个多晶平均E填入单晶所有方向，会丢掉这类各向异性。

4. 接到自己的更大尺度模型#

保留晶轴与器件坐标的关系，选择单晶张量或与你模型相符的平均模量，再传递密度、单位和温度条件。本次只是静态零温附近的DFT弹性输入，没有补入缺陷、晶界、温度或表面效应。

cat results/Si-material-SI.json
bash plot.sh --data-dir results

想接着计算具体结构响应，可参考Si悬臂梁的实际参数传递例子。那篇采用解析梁模型，也不等同于已经完成任意几何的有限元求解。