GROMACS 小溶质自由能实操

下载完整输入、真实轨迹、dhdl、BAR结果与绘图脚本。解压进入资料包目录，先处理已有结果：

bash analyze.sh

reference/figures 保存数据图，reference/analysis 保存CSV、完整JSON与NPZ。想重新计算，在新解压目录只需：

bash submit.sh

脚本提交16个λ状态的数组任务，每个任务用1张P100、2个CPU核心顺序跑两组种子；全部结束后自动提交一个小约束对照与分析任务。不会让你手工串接一长串作业号。新结果在 work。这里是实际32个λ/种子条件加1个约束对照，不是33种独立性质。

本例使用 gromacs/2026.4-p100-threadmpi，GNU13.3、CUDA12.8、thread-MPI/OpenMP环境，生产阶段非扰动短程相互作用和PME使用GPU，扰动相关部分按原生程序的CPU路径执行，更新/约束使用CPU；不是所有计算都在GPU上。双精度或MPI/PLUMED构建与这里不同。数据处理使用ASE3.29.0环境中的Python、NumPy、SciPy与Matplotlib。不要在登录节点运行模拟；已有结果的提取与绘图无需重新跑MD。

1 这个算例具体算什么#

在2.6 nm立方周期盒中放一个中性LJ测试粒子与577个SPC/E水分子。粒子标记为AR，质量39.948、σ=0.340 nm、ε=0.997 kJ/mol；它用于学习流程，并非经实验校准的某种溶质水合模型。采用OPLS几何混合规则，0.9 nm截断并平移LJ势，不加长程色散尾修正。

λ=0时粒子与水完整相互作用；λ=1时完全去耦。这里求的是固定盒体积下，从相互作用态到不相互作用态的自由能差。其符号与把粒子耦合进水中相反；在本模型下，去耦ΔG对应过量化学势的负值。它不是配体结合自由能，也不是默认1 bar下的实验水合自由能。

2 从拓扑和坐标开始#

每个任务先用 gmx solvate 从相同的单粒子坐标建立水盒，实际加入577个水。拓扑核心如下：

#include "oplsaa.ff/forcefield.itp"
[ atomtypes ]
ARLJ 18 39.948 0.0 A 0.340 0.997
[ moleculetype ]
AR 1
[ atoms ]
1 ARLJ 1 AR AR 1 0.0 39.948
#include "oplsaa.ff/spce.itp"
[ system ]
One neutral Lennard-Jones Ar in SPC/E water
[ molecules ]
AR 1
SOL 577

inputs里的初始拓扑只有AR行，gmx solvate -p topol.top 自动补水数；上面展示的是补齐后的版本。不要重复添加SOL。AR是分子类型，ARLJ是原子类型，二者不能混用到 couple-moltype。

在已分配资源、从某个输入副本目录执行的实际命令为：

module purge
module load gromacs/2026.4-p100-threadmpi
gmx solvate -cp argon.gro -cs spc216.gro -box 2.6 2.6 2.6 -o solvated.gro -p topol.top
gmx grompp -f em.mdp -c solvated.gro -p topol.top -o em.tpr
gmx mdrun -deffnm em -ntmpi 1 -ntomp 2 -nb gpu -pme cpu -update cpu

先能量最小化到最大力小于500 kJ mol⁻¹ nm⁻¹，再进入各自λ状态的平衡。能量最小化与动力学阶段的PME设备设置不同，不要把上面的EM命令机械套成生产命令。

3 λ、软核和dhdl怎样设置#

以下是 lambda00-seed1/md.mdp 的完整实际配置：

integrator = sd
dt = 0.002
tcoupl = no
tc-grps = System
tau-t = 1.0
ref-t = 300
pcoupl = no
nstcalcenergy = 100
nstenergy = 1000
nstlog = 10000
nstxout = 0
nstvout = 0
nstfout = 0
cutoff-scheme = Verlet
nstlist = 20
rlist = 0.9
coulombtype = PME
rcoulomb = 0.9
fourierspacing = 0.12
pme-order = 4
vdwtype = Cut-off
vdw-modifier = Potential-shift
rvdw = 0.9
DispCorr = no
pbc = xyz
constraints = h-bonds
constraint-algorithm = lincs
lincs-iter = 1
lincs-order = 4
free-energy = yes
init-lambda-state = 0
fep-lambdas = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.94 0.97 1
couple-moltype = AR
couple-lambda0 = vdw
couple-lambda1 = none
couple-intramol = no
sc-alpha = 0.5
sc-power = 1
sc-sigma = 0.3
sc-coul = no
calc-lambda-neighbors = -1
nstdhdl = 100
dhdl-derivatives = yes
dhdl-print-energy = potential
separate-dhdl-file = yes
nsteps = 500000
ld-seed = 2921600
continuation = yes
gen-vel = no
nstxout-compressed = 2500

每个状态改变 init-lambda-state，索引为0–15；具体λ数组是0、0.1、0.2、0.3、0.4、0.5、0.6、0.65、0.7、0.75、0.8、0.85、0.9、0.94、0.97、1。实际路径不是均匀分布，TI积分时必须用这些真实间距。

couple-lambda0 = vdw 到 couple-lambda1 = none 关闭这个中性粒子的LJ相互作用。sc-alpha = 0.5 等参数启用软核，缓和接近去耦端点时的重叠奇异性。软核使Hamiltonian随λ非线性变化，不能只用 dH/dλ × Δλ 代替真实的相邻状态能差来做BAR。

calc-lambda-neighbors = -1 输出全部16个外部λ状态的ΔH；nstdhdl = 100 对应每0.2 ps记录一次。dhdl-print-energy = potential 多写一列势能，它不是dH/dλ列，读数据时先看XVG图例。

4 平衡、生产、后处理依次做什么#

每个λ先做20 ps随机动力学平衡，再从检查点接续1 ns生产；300 K、2 fs、固定盒体积。integrator = sd 给完全去耦粒子也提供随机热浴，本例用它采样平衡自由能，不能将这类轨迹用于解释真实动力学输运系数。

gmx grompp -f eq.mdp -c em.gro -p topol.top -o eq.tpr
gmx mdrun -deffnm eq -ntmpi 1 -ntomp 2 -nb gpu -pme gpu -bonded cpu -update cpu -pin off
gmx grompp -f md.mdp -c eq.gro -t eq.cpt -p topol.top -o md.tpr
gmx mdrun -deffnm md -ntmpi 1 -ntomp 2 -nb gpu -pme gpu -bonded cpu -update cpu -pin off

上面是工作流内部实际执行的命令；平时直接 bash submit.sh 即可。MD输出 md.xvg 为dhdl，md.xtc为结构轨迹，md.edr为能量，md.tpr保存运行拓扑，md.log记录运行与参数。

两次重复分别改变初速度与随机动力学种子，各状态也使用独立种子。本例各条件采样平均温度299.54–300.70 K；固定体积下各自平均压力约−111到30 bar，并未控制为1 bar。若要其他密度、压力或真实溶质，需要重新设计平衡与生产条件。

5 最后得到什么#

BAR、TI、重叠与误差分析给出实际自由能和判断采样是否充分的方法；约束与标准态说明为什么去耦自由能不能直接写成绝对结合自由能。完整原始dhdl、真实结构、能量和短命令都随包提供。

λ与dhdl参数定义见 GROMACS官方MDP选项。这里是小模型流程参考，正式性质需要对时间、λ密度、盒大小、相互作用模型和统计误差继续检查。