DFT-D4参数、净电荷与力的提取

下载本篇输入、完整脚本、数据与高清图。解压后在新目录运行，reference保存鸢算科技实际运行的参考结果。

1 选择与DFT计算一致的参数#

先按DFT-D4色散能教程准备环境和XYZ。下面分别演示改变泛函及总电荷，每项放在独立目录，避免JSON被覆盖：

dftd4 input.xyz --func pbe --charge 0 --json --grad
dftd4 input.xyz --func pbe0 --charge 0 --json --grad
dftd4 input.xyz --func b3lyp --charge 1 --json --grad

我们对相同的2.9 Å双水构型实际计算了PBE、PBE0、B3LYP，每种总电荷取−1、0、+1。中性构型的D4修正分别为−0.001018327866、−0.000857018206、−0.001580506778 Hartree。

改变电荷会影响D4的EEQ电荷模型。这组差值不能解释为电离能或电子亲和能，也不是由DFT电子密度直接得到的原子电荷。正式计算必须选择与主程序一致的总电荷和泛函参数。

2 从梯度得到力#

import json
import numpy as np
x = json.load(open('dftd4.json'))
g = np.array(x['gradient']).reshape(-1, 3)
forces = -g * 27.211386245988 / 0.529177210903
np.savetxt('forces-eV-A.txt', forces)

每个原子有x、y、z三个分量。不要漏掉负号，也不要将Hartree/Bohr直接标成eV/Å。

3 用有限差分检查梯度#

下载包d4_examples.py对6个原子的18个坐标分量逐一加、减一个小位移，重新计算能量。由[E(x+h)−E(x−h)]/(2h)得到数值导数，统一单位后与程序梯度比较：

这验证了该小体系中梯度符号、单位和文件读取方式。它不能替代DFT主程序的力收敛检查。与VASP、CP2K或ORCA联用时，还应确认主程序是否已经包含同一色散修正，避免重复添加。

参考：DFT-D4官方实现。