GPUMD 二维材料应力与真空高度

本例使用鸢算科技集群的 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. 准备二维结构与两个真空高度#

本例使用 128 原子石墨烯、Lindsay 2010 参数化的 Tersoff 型 C 势（Graphene_Lindsay_2010.txt），x/y 周期、z 非周期。初始晶格常数 a=2.46 Å；x 方向做 0.5% 应变，并用 0.45%、0.55% 两个邻近状态检查能量导数。

head -n 5 inputs/graphene-20.0-0.005/model.xyz
cat inputs/graphene-20.0-0.005/CASE.json
cat inputs/graphene-20.0-0.005/run.in

第二行的关键字段为：

Properties=species:S:1:pos:R:3:vel:R:3
pbc="T T F"

实际完整晶胞见下载包。20 Å 与 40 Å 两套输入只改变 z 盒高和片层的整体 z 位置；面内坐标、势函数、应变和零速度完全相同。这里的“真空高度”指盒子 Lz，不是给石墨烯指定物理厚度。

2. 计算后提取压力与二维应力#

重新计算完成后：

cat work/graphene-20.0-0.005/thermo.out
cat work/graphene-40.0-0.005/thermo.out

当前 GPUMD 的数字行有 18 列，前三列是 T、KE、PE，接着六列是压缩为正的压力，最后九列是晶胞。忽略 # 注释后再读：

import numpy as np
t = np.atleast_2d(np.loadtxt('thermo.out'))
Pxx_GPa = t[-1, 3]
sigma_2D_N_m = -Pxx_GPa * Lz_A * 0.1

乘 0.1 是 GPa×Å 到 N/m 的单位换算。负号把压缩为正的压力改为拉伸为正的应力。

盒高加倍，GPa 压力减半，二维应力保持一致；两套势能都为 -1020.1870049 eV。这里结果仍是压应力：0.5% 是相对 a=2.46 Å 结构的拉伸，不代表已相对该经验势的零应力晶格拉伸。要研究杨氏模量，应先求这个势的面内零应力参考，再做小应变对照。

3. 用能量随应变变化交叉检查#

在固定 y 方向的简单单轴变形下，本例用 E(0.0055)-E(0.0045) 除以 0.001、参考面内面积，再乘 16.021766208，得到 -3.169739 N/m，与 virial 提取相符。

脚本同时保持 vel=0、检查势能与原子势能之和、核对逐原子 virial 和原生压力。若进行有限温度动力学，要另外区分动能应力、瞬时波动与充分时间平均。

4. 迁移到其他二维体系#

先选择适合材料的势，准备 x/y 周期、z 非周期的结构，再设置零应力参考和面内变形。无相互作用的真空增大不应改变面内总能量和归一化二维应力；若改变，检查周期边界、原子位置、相互作用范围和截断。

不要随意给二维体系指定厚度后直接与三维 GPa 混比；若要报告“等效三维模量”，必须说明采用的物理厚度。

参考：GPUMD 模型文件、热力学输出格式。