GPUMD 铜的 EOS 与剪切模量

本例使用鸢算科技集群的 GPUMD 5.9.1 P100 构建（GCC 13.3、CUDA 12.9），用 ASE 3.29、NumPy 和 Matplotlib 建模、提取和画图。每个原生任务只用一张 P100。所用数据是经典势实际计算结果。

下载本篇配套输入、真实结果与脚本，解压进入目录。先查看已有结果：

bash run.sh

analysis/RESULTS.json 为结果汇总，CSV 可继续作图，figures/ 为高清图片。这个步骤不需要 GPU。要重新计算，另解压一份，进入新目录：

sbatch compute.slurm

脚本按小体系静态计算、短热接触的顺序运行，申请一张 P100、一个 CPU 核心。inputs/ 是输入，results/ 是配套参考输出，work/ 是你新运行的结果。不要把新输出追加到原参考目录。完成后自动提取 CSV 和图片。

1. 先优化参考晶胞#

本例是 4×4×4 的 FCC Cu 常规胞超胞，共 256 个原子，使用 GPUMD 官方 Zhou 2004 EAM 势。完整结构和势文件已在包内，无需手写 256 行坐标。

cat inputs/cu-reference/run.in
head -n 6 inputs/cu-reference/model.xyz
cat inputs/cu-reference/potential.txt

实际的 run.in：

potential potential.txt
minimize fire 0.00001 5000 1 1
ensemble nve
time_step 0
dump_thermo 1
dump_xyz 1 snapshot.xyz precision double mass velocity force potential virial
run 1
compute_elastic 0.005

minimize 最后两个 1 分别开启晶胞变化和静水应变；原子位置同时优化。若只优化原子、保持晶胞不动，用 minimize fire 0.00001 5000 0。compute_elastic 0.005 用 0.5% 的有限应变计算矩阵，输出 elastic.out，单位 GPa。

结构文件显式包含 vel:R:3，每个原子初速度均为 0 0 0。这里 time_step 0 只用于静态取样；不要拿它进行真实动力学。此版本 velocity 0 不接受；如果不写速度，程序会默认初始化 300 K，静态应力就会混入动能项。

2. 固定五个不同体积分别算能量#

先完成参考晶胞优化，再以它为基准，将体积取为 0.94、0.97、1.00、1.03、1.06 倍。三条晶格矢量和原子坐标按体积比的立方根同步缩放。

cat inputs/cu-volume-0.94/run.in
head -n 3 inputs/cu-volume-0.94/model.xyz

这些输入保留零速度和静态输出，不再开启晶胞优化，否则不同体积会全部回到同一个平衡体积。当前扫描使用 FCC 单原子基元的均匀缩放；有内部自由度的材料，可在每个固定体积下继续优化原子位置，但不要放开体积。

analyze.py 读取每个目录 thermo.out 的势能列及快照晶胞，用 ASE EquationOfState(..., eos='birchmurnaghan') 拟合：

from ase.eos import EquationOfState
eos = EquationOfState(volume_per_atom, energy_per_atom,
                      eos='birchmurnaghan')
v0, e0, bulk_modulus = eos.fit()
bulk_modulus_GPa = bulk_modulus * 160.21766208

实际结果：拟合平衡体积 11.818998 ų/原子，体模量 136.413 GPa。原生弹性矩阵得到的体模量约 138.910 GPa。这是五个有限体积点的教学对照；提高精度应加密体积点、比较拟合窗口和应变幅度。

3. 做正负简单剪切#

本例定义变形梯度 Fxy=γ：x'=x+γy，y'=y，z'=z，同时变换晶胞和原子位置。工程剪切 γ 为 ±0.01；它是小应变张量 εxy 的两倍。

cat inputs/cu-shear-0.01/CASE.json
head -n 3 inputs/cu-shear-0.01/model.xyz

用正、负剪切势能与未变形势能计算：

Cxyxy = [E(+γ)+E(-γ)-2E(0)] / (V0 γ²)

上式得到 eV/ų，乘 160.21766208 变成 GPa。体积不变，使用参考总体积 V0；不要混用总能量和每原子体积。

GPUMD 的热力学压力以压缩为正，六分量顺序为 xx、yy、zz、yz、xz、xy。因此拉伸为正的剪应力需要取相反号。源程序压力换算常数为舍入后的 160.2177；脚本按这个常数核对原生输出，物理单位独立换算则用上面的 SI 数值。

4. 查看导出文件#

head analysis/cu-eos.csv
cat analysis/RESULTS.json
ls figures

这些数值属于本例 Cu 经验势和有限窗口，不是 DFT 结果，也不是所有 Cu 势都必须得到的材料常数。继续看原子与晶胞优化、弹性应变幅度、拉伸及二维应力。

依据：GPUMD 最小化参数、弹性计算。