GROMACS 软核自由能与BAR
软核自由能与BAR，包含终端操作、输入文件、实际数据和适用范围。
阅读本篇目录
本文以鸢算科技集群 GROMACS 2026.4 实际结果为例。下载输入、参考轨迹、分析与绘图资料包，资料包约 529 MB，包含参考轨迹。解压后进入 gromacs-examples。参考输出在 cases；新计算由 new-case.py 建立独立目录。
本篇任务与操作要点#
本例在水中放置一个自定义中性 Lennard–Jones 测试粒子，σ=0.34 nm、ε=0.997 kJ/mol、质量 39.948。它用于演示 λ 文件和 BAR，不是经过校准的某种真实溶质模型。
五个 λ 状态从完全相互作用逐步变到关闭相互作用，使用软核避免端点粒子重叠导致发散。calc-lambda-neighbors = -1 输出所有状态的能差；每窗口 50 ps，BAR 舍弃前 10 ps。若输入改成带电溶质，还要分别设计电荷/LJ路径、处理电荷有限尺寸和标准态问题。
gmx bar 读取各窗口 DHDL 文件计算相邻自由能差与总和。正反向分布需要充分重叠，软件给出的误差条也依赖采样质量；50 ps 只用于演示，不能把总和当作收敛水合自由能。去耦与耦合方向相反，引用数值前先核对符号。
1 准备与提交#
python new-case.py lj-solute
cd work/lj-solute
sbatch run.sh
run.sh 申请资源后执行 workflow.sh。不要在登录节点直接执行完整模拟。下列每个算例可单独新建副本并提交。
| 算例目录 | 本次内容 |
|---|---|
lj-solute |
完整输入、前置文件与真实计算输出 |
fep-0 |
完整输入、前置文件与真实计算输出 |
fep-1 |
完整输入、前置文件与真实计算输出 |
fep-2 |
完整输入、前置文件与真实计算输出 |
fep-3 |
完整输入、前置文件与真实计算输出 |
fep-4 |
完整输入、前置文件与真实计算输出 |
bar |
实际完成的分析输入与结果 |
最后的 wham/bar 依赖前面全部窗口，必须等所有窗口完成后再提交。清单中的相对路径读取 work 下的新结果；不能只创建分析目录就直接运行。
2 实际终端命令#
以下是在对应算例目录中实际执行并生成结果的命令。分析命令在该算例的 analysis 子目录执行；交互回答已用 printf 写入，组名需随自己的体系核对。
lj-solute
gmx solvate -cp solute.gro -cs spc216.gro -o solv.gro -p topol.top
gmx grompp -f em.mdp -c solv.gro -p topol.top -o em.tpr
gmx mdrun -deffnm em -ntmpi 1 -ntomp 4 -nb gpu -pme cpu
fep-0
gmx grompp -f run.mdp -c start.gro -p topol.top -o run.tpr
gmx mdrun -ntmpi 1 -ntomp 4 -nb gpu -pme gpu -update cpu -deffnm run
bar
gmx bar -f ../fep-0/run.xvg ../fep-1/run.xvg ../fep-2/run.xvg ../fep-3/run.xvg ../fep-4/run.xvg -b 10 -o bar.xvg -oi bar-integral.xvg -oh overlap.xvg
3 输入文件与关键设置#
lj-solute/topol.top 完整输入（点击展开）
#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
fep-0/topol.top 完整输入（点击展开）
#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
fep-0/run.mdp 完整输入（点击展开）
cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
rcoulomb = 1.0
rvdw = 1.0
pbc = xyz
constraints = h-bonds
integrator = md
dt = 0.002
nsteps = 25000
nstenergy = 50
nstcalcenergy = 50
nstlog = 1000
nstxout-compressed = 50
tcoupl = v-rescale
tc-grps = System
tau-t = 0.5
ref-t = 300
pcoupl = no
gen-vel = yes
gen-temp = 300
gen-seed = 88100
continuation = no
pcoupltype = isotropic
tau-p = 5
ref-p = 1
compressibility = 4.5e-05
nstxout = 0
nstvout = 0
nstfout = 0
free-energy = yes
init-lambda-state = 0
vdw-lambdas = 0 0.25 0.5 0.75 1
couple-moltype = LJ
couple-lambda0 = vdw-q
couple-lambda1 = none
couple-intramol = no
sc-alpha = 0.5
sc-power = 1
sc-sigma = 0.3
calc-lambda-neighbors = -1
nstdhdl = 50
fep-1/topol.top 完整输入（点击展开）
#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
fep-1/run.mdp 完整输入（点击展开）
cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
rcoulomb = 1.0
rvdw = 1.0
pbc = xyz
constraints = h-bonds
integrator = md
dt = 0.002
nsteps = 25000
nstenergy = 50
nstcalcenergy = 50
nstlog = 1000
nstxout-compressed = 50
tcoupl = v-rescale
tc-grps = System
tau-t = 0.5
ref-t = 300
pcoupl = no
gen-vel = yes
gen-temp = 300
gen-seed = 88101
continuation = no
pcoupltype = isotropic
tau-p = 5
ref-p = 1
compressibility = 4.5e-05
nstxout = 0
nstvout = 0
nstfout = 0
free-energy = yes
init-lambda-state = 1
vdw-lambdas = 0 0.25 0.5 0.75 1
couple-moltype = LJ
couple-lambda0 = vdw-q
couple-lambda1 = none
couple-intramol = no
sc-alpha = 0.5
sc-power = 1
sc-sigma = 0.3
calc-lambda-neighbors = -1
nstdhdl = 50
fep-2/topol.top 完整输入（点击展开）
#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
fep-2/run.mdp 完整输入（点击展开）
cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
rcoulomb = 1.0
rvdw = 1.0
pbc = xyz
constraints = h-bonds
integrator = md
dt = 0.002
nsteps = 25000
nstenergy = 50
nstcalcenergy = 50
nstlog = 1000
nstxout-compressed = 50
tcoupl = v-rescale
tc-grps = System
tau-t = 0.5
ref-t = 300
pcoupl = no
gen-vel = yes
gen-temp = 300
gen-seed = 88102
continuation = no
pcoupltype = isotropic
tau-p = 5
ref-p = 1
compressibility = 4.5e-05
nstxout = 0
nstvout = 0
nstfout = 0
free-energy = yes
init-lambda-state = 2
vdw-lambdas = 0 0.25 0.5 0.75 1
couple-moltype = LJ
couple-lambda0 = vdw-q
couple-lambda1 = none
couple-intramol = no
sc-alpha = 0.5
sc-power = 1
sc-sigma = 0.3
calc-lambda-neighbors = -1
nstdhdl = 50
fep-3/topol.top 完整输入（点击展开）
#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
fep-3/run.mdp 完整输入（点击展开）
cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
rcoulomb = 1.0
rvdw = 1.0
pbc = xyz
constraints = h-bonds
integrator = md
dt = 0.002
nsteps = 25000
nstenergy = 50
nstcalcenergy = 50
nstlog = 1000
nstxout-compressed = 50
tcoupl = v-rescale
tc-grps = System
tau-t = 0.5
ref-t = 300
pcoupl = no
gen-vel = yes
gen-temp = 300
gen-seed = 88103
continuation = no
pcoupltype = isotropic
tau-p = 5
ref-p = 1
compressibility = 4.5e-05
nstxout = 0
nstvout = 0
nstfout = 0
free-energy = yes
init-lambda-state = 3
vdw-lambdas = 0 0.25 0.5 0.75 1
couple-moltype = LJ
couple-lambda0 = vdw-q
couple-lambda1 = none
couple-intramol = no
sc-alpha = 0.5
sc-power = 1
sc-sigma = 0.3
calc-lambda-neighbors = -1
nstdhdl = 50
fep-4/topol.top 完整输入（点击展开）
#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
fep-4/run.mdp 完整输入（点击展开）
cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
rcoulomb = 1.0
rvdw = 1.0
pbc = xyz
constraints = h-bonds
integrator = md
dt = 0.002
nsteps = 25000
nstenergy = 50
nstcalcenergy = 50
nstlog = 1000
nstxout-compressed = 50
tcoupl = v-rescale
tc-grps = System
tau-t = 0.5
ref-t = 300
pcoupl = no
gen-vel = yes
gen-temp = 300
gen-seed = 88104
continuation = no
pcoupltype = isotropic
tau-p = 5
ref-p = 1
compressibility = 4.5e-05
nstxout = 0
nstvout = 0
nstfout = 0
free-energy = yes
init-lambda-state = 4
vdw-lambdas = 0 0.25 0.5 0.75 1
couple-moltype = LJ
couple-lambda0 = vdw-q
couple-lambda1 = none
couple-intramol = no
sc-alpha = 0.5
sc-power = 1
sc-sigma = 0.3
calc-lambda-neighbors = -1
nstdhdl = 50
4 实际结果与画图#
本次测试粒子去耦方向的 BAR 总和为 −8.23 ± 0.84 kJ/mol（程序对这段短数据的估计），不是经采样收敛的真实溶质水合自由能。
用附件中的 Python 脚本重画所有参考图：
# 回到 gromacs-examples 根目录
module load ase/3.29.0
python plot.py --root .
plot.py 使用 NumPy 和 Matplotlib，直接读取 XVG/COLVAR/坐标数值，图片保存在 figures。对单项新结果，可照着脚本对应段替换文件路径；不要把不同模型或单位的数据直接拼在一起。
参数查阅#
以本集群 gmx 命令名 -h 为准；GROMACS 2026.4 命令手册和MDP 参数手册可查完整选项。
需要根据你的体系选择资源？联系讨论计算需求。实际可用软件与权限，以当前账户和运行环境为准。