GROMACS 能量分组、能量矩阵与 LIE 公式
为溶质和水保存分组非键能,提取enemat矩阵并用真实温度对照演示lie,明确非结合自由能。
阅读本篇目录
下载输入、真实结果与绘图资料包,解压进入资料包目录。inputs 可修改,reference 是集群已计算的参考结果;脚本每次在 work 下建立独立目录。计算在分配好的资源内执行,可用 sbatch run.sbatch 跑整组,结果路径写入 LAST-RUN.txt。
本例的用途与范围#
这里用一个中性 LJ 溶质加 883 个水分子,分别运行 280 K 和 320 K 各 6 ps,演示分组非键能、enemat 和 lie 的输入。温度对照只是公式练习,体系没有受体,结果不是结合自由能。
bash run.sh md
inputs/interactions-280 和 inputs/interactions-320 含完整输入。下面以 320 K 目录为例。
1. 在运行前定义能量组#
索引把第一个 LJ 原子命名为 Solute,其余水为 Water;另保留完整 System 温度组。能量组名和分子类型名称可以不同。本例使用一个完整温度耦合组,避免将单个溶质原子单独作为很小的温度耦合组。
完整拓扑:
#include "amber99sb-ildn.ff/forcefield.itp"
[ atomtypes ]
LJX 18 39.948 0.0 A 0.34 0.997
[ moleculetype ]
LJ 1
[ atoms ]
1 LJX 1 LJ LJ 1 0.0 39.948
#include "amber99sb-ildn.ff/tip3p.itp"
[ system ]
Neutral LJ solute in water
[ molecules ]
LJ 1
SOL 883
完整 MD 输入:
integrator = md
dt = 0.001
nsteps = 6000
cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
rcoulomb = 1.0
rvdw = 1.0
pbc = xyz
constraints = h-bonds
continuation = no
tcoupl = v-rescale
tc-grps = System
tau-t = 0.1
ref-t = 320
pcoupl = no
gen-vel = yes
gen-temp = 320
gen-seed = 20261348
energygrps = Solute Water
nstenergy = 10
nstcalcenergy = 10
nstxout-compressed = 100
ld-seed = -663491073
energygrps=Solute Water 必须在生成 TPR 前设置。它不会从一个已经不含分组能量的旧 EDR 中补出缺失数据。本例用 CPU 非键执行路径保存分组能;普通 GPU 生产轨迹若未保存这些量,可以另规划 CPU 重算并核对设置。
2. 运行并按名称提取能量#
gmx grompp -f md.mdp -c start.gro -p topol.top -n index.ndx -o md.tpr
gmx mdrun -deffnm md -ntmpi 1 -ntomp 2 -nb cpu -bonded cpu -update cpu
printf 'LJ-SR:Solute-Water\nCoul-SR:Solute-Water\n0\n' | gmx energy -f md.edr -b 2 -o solute-water.xvg
这里丢弃前 2 ps,仅作为短流程演示。读取 XVG 必须按 legend 对应列名:原生输出可能按其内部能量项顺序排列,并不一定按终端输入顺序。本例实际第一能量列是 Coul-SR,第二列才是 LJ-SR。中性溶质对应的 Coul-SR 为零。
3. 输出能量矩阵#
groups.dat:
2
Solute
Water
gmx enemat -f md.edr -groups groups.dat -b 2 -emat matrix.xpm -etot total-groups.xvg -nocoul14 -nolj14
实际生成 LJ-SRmatrix.xpm、Coul-SRmatrix.xpm、totalmatrix.xpm 等文件。XPM 色阶会量化数据,大的溶剂自相互作用还会压缩小交叉项的颜色对比;要读准确交叉能量,优先提取 EDR 到 XVG。

4. 用 lie 检查线性组合#
280 K 参考段溶质—水平均 LJ-SR 为 −9.06670448 kJ/mol,Coul-SR 为零。在 320 K 的目录执行:
gmx lie -f md.edr -b 2 -ligand Solute -Elj -9.06670448 -Eqq 0 -Clj 0.181 -Cqq 0.5 -o linear-interaction-contrast.xvg
本例逐帧结果等于:
0.181 × (E_LJ,320(t) − <E_LJ,280>)
+ 0.5 × (E_Coul,320(t) − <E_Coul,280>)
平均线性对照约 −0.00156233 kJ/mol;与按列名独立复算的最大差小于 5.1e-6 kJ/mol。短时间平均非常接近零也不能据此声称温度效应为零或误差已收敛。
真实 LIE 结合自由能需要合适的结合态与溶液参考态、充分采样和经过研究的系数。不要把这里的两个温度替换成“结合/游离”标签后直接作科学结论。