DP-GEN弹性常数自动测试

本篇接自动优化与EOS，使用同一个实例包。先完成基础晶胞优化，再对优化结构施加应变、计算应力。应变后的形状和体积保持固定，只允许内部原子位置调整。

1 输入与命令#

cat templates/elastic.json
bash elastic.sh

{
  "structures": [
    "confs/Ni-fcc"
  ],
  "interaction": {
    "type": "eam_alloy",
    "model": "__ROOT__/NiAlH_jea.eam.alloy",
    "type_map": {
      "Ni": 0
    }
  },
  "properties": [
    {
      "type": "elastic",
      "norm_deform": 0.005,
      "shear_deform": 0.005
    }
  ]
}

本例正应变、剪切应变幅度参数均为0.005，生成24个任务。relax_pos=true、relax_shape=false、relax_vol=false维持所施加的应变；如果在每个任务中重新自由优化晶胞，应力应变关系就会被改变。

ls confs/Ni-fcc/elastic_00/task.000000
cat confs/Ni-fcc/elastic_00/task.000000/strain.json
cat confs/Ni-fcc/elastic_00/result.out

每个任务的strain.json是实际应变张量，in.lammps、conf.lmp是实际模拟输入，result_task.json是提取结果。最终result.json里的elastic_tensor是36个按行排列的元素，不能当成36个互相独立的立方弹性常数。

2 本例实际数值#

泊松比为0.276609，无单位。对于立方结构可用B=(C11+2C12)/3作直观核对。C11、C12、C44是单晶张量分量，BV、GV等是Voigt平均得到的量，不要混写成同一个物性。

3 读取矩阵和画图数据#

source env.sh
python - <<'PY'
import json
import numpy as np
p = 'confs/Ni-fcc/elastic_00/result.json'
r = json.load(open(p))
C = np.array(r['elastic_tensor']).reshape(6, 6)
print('弹性矩阵(GPa):')
print(C)
print('BV(GPa):', r['BV'])
PY

只查看附件参考值时用bash results.sh或读取弹性矩阵CSV，无需再次执行24个计算。正式研究还应改变应变幅度、优化精度和超胞检查线性区，选用适合目标材料的势函数。

当前原生程序把LAMMPS应力整理成结果张量并输出GPa弹性常数；不要再对result.json里的最终弹性常数重复乘1000或10000。参见DP-GEN官方弹性教程。