ASE 能量数值求导、力与应力
畸变Cu4比较三个差分步长，独立核对解析受力、应力、剪切及单位。
阅读本篇目录
本篇使用鸢算科技集群实际运行的ASE3.29.0小例子，全部为单线程CPU。Ar3用两套Lennard–Jones参数说明接口；Cu4用EMT经验势，不是DFT或GPU性能测试。
下载完整Python脚本、实际结构、数据与绘图包，在新目录：
unzip ase-calculator-wrappers-practical*.zip -d ase-example
cd ase-example
module load ase/3.29.0
bash run.sh
默认读取随附实算数值生成4图和analysis/RESULTS.json，不重算势能。要自己重算，在已分配的CPU资源中执行：
bash run.sh compute
cat work/RESULTS.json
输出包括自己的work/、27项数值核对和4图。脚本不会覆盖已有work；第二次重算请另解压一份。其他电脑已安装ASE3.29、NumPy与Matplotlib时，可用PYTHON=python3 bash run.sh。已有work需要重画时用bash run.sh analyze。
1. 从一份有非零受力和应力的Cu4开始#
先建立FCC铜常规胞，再轻微移动一个原子并施加晶胞畸变。这样各量不全是对称性给出的数值零，适合检查求导。
import numpy as np
from ase.build import bulk
from ase.calculators.emt import EMT
atoms = bulk('Cu', 'fcc', a=3.62, cubic=True)
atoms.positions[0] += [.04,-.03,.02]
atoms.set_cell(atoms.cell @ np.array([
[1.01,.006,0],[0,.99,.003],[0,0,1.005]]), scale_atoms=True)
atoms.calc = EMT()
print(atoms.get_potential_energy())
print(atoms.get_forces())
print(atoms.get_stress())
完整初始结构在reference/cu4-distorted.extxyz，自己的重算生成work/cu4-distorted.extxyz。EMT只是这个金属小例的经验势，不用于据此判断其他材料的DFT精度。
2. 用能量差分获得力和应力#
from ase.calculators.fd import FiniteDifferenceCalculator
numerical = atoms.copy()
numerical.calc = FiniteDifferenceCalculator(
EMT(), eps_disp=1e-4, eps_strain=1e-4,
force_consistent=False)
print(numerical.get_forces())
print(numerical.get_stress())
eps_disp是Å，eps_strain是无量纲应变。力是能量对位移的负导数；应力还涉及晶胞变形和体积归一化，不能把两个步长的单位混用。force_consistent=False在这里选择通常势能路径；其他含电子温度的计算器应先确认能量、自由能和解析受力的定义是否一致。
3. 三个步长的实际数值#
| 位移Å / 应变无量纲 | 最大力差 / eV·Å⁻¹ | 最大应力差 / eV·Å⁻³ |
|---|---|---|
| 1e-03 | 3.77005e-07 | 7.82991e-07 |
| 1e-04 | 3.77980e-09 | 7.83028e-09 |
| 1e-05 | 8.97399e-10 | 8.28473e-11 |
这个光滑小体系中减小步长后差异减小；并非步长越小在所有程序中都更准确。DFT的SCF噪声、积分网格、势函数截断或近邻变化都可能让过小步长放大数值误差。
cat work/VALIDATION.json
python - <<'PYDATA'
import json
x=json.load(open('work/RESULTS.json'))['derivatives']
print(x['analytic_stress_eV_A3'])
for row in x['steps']:
print(row['step_A_and_strain'],row['force_max_error_eV_A'])
PYDATA
4. 独立核对剪切与应力符号#
ASE Voigt顺序为xx, yy, zz, yz, xz, xy，本例单位eV/ų。脚本另建六种对称形变；剪切矩阵的两个对称非对角元各取1/2，再用(Eplus-Eminus)/(2hV)与解析应力比较，不直接复用数值包装器内部的求导过程。
如果要以GPa表示，1 eV/ų ≈ 160.217662 GPa。不要把软件中“正应力”和“正压力”的符号约定混为一谈。需要用应力做晶胞优化时，继续看几何与应力格式和状态方程。
来源：ASE数值求导计算器。