要说清楚多肽和LBP（脂多糖结合蛋白）之间的结合自由能怎么算，光有PDB结构才走了一小半的路。实际操作中，从体系搭建到轨迹平衡，再到采样和最终的自由能计算，每一步都藏着不少细节。下面梳理一套完整的操作流程，用的是GROMACS version 2025.2，力场选AMBER99SB-ILDN，水模型用TIP3P。这套组合在多肽-蛋白复合物体系里表现稳定，也是目前比较主流的配置之一。
先说几个核心判断：第一，起始结构质量决定了模拟的下限。拿到复合物PDB文件（比如这里用的 binder_l37_s368566.pdb）之后，千万别急着直接扔进pdb2gmx，要先看看里面有没有非标准氨基酸、残缺残基或者链标识混乱的问题。简单的排查方法就是用一行命令看看原子信息：
grep "^ATOM" binder_l37_s368566.pdb | awk '{print $4,$5,$6}' | sort | uniq -c
这一步能快速辨认出每个残基的编号和链归属，如果发现异常，比如某条链上的残基编号断档、或者有非标准残基名字出现，就需要先在原始PDB里做修正，否则后面生成的拓扑文件肯定是错的。
确认结构没问题之后，生成拓扑文件就比较流畅了：
gmx pdb2gmx -f binder_l37_s368566.pdb -o complex.gro -p topol.top -i posre.itp -ff amber99sb-ildn -water tip3p -ignh
这里的 -ignh 指令会忽略原有PDB里的氢原子，让GROMACS按照力场标准重新加氢，这一步对后续的质子化状态判断非常关键。有些多肽的末端残基在生理pH下可能是带电的，pdb2gmx 会自动处理，但是如果你对某个残基的质子化状态有特殊要求（比如组氨酸的不同质子化形式），就需要在执行前手动修改残基名称来指定。
建立模拟盒子与溶剂化
盒子形状上，十二面体相比正方体能节省大约30%的水分子量，计算效率明显更高：
gmx editconf -f complex.gro -o box.gro -c -d 1.0 -bt dodecahedron
边缘距离1.0 nm作为最小值比较稳妥，既能防止蛋白质与自己的镜像相互作用，也不会浪费太多溶剂体积。然后填充溶剂：
gmx solvate -cp box.gro -cs spc216.gro -o solv.gro -p topol.top
这里的 spc216.gro 虽然是SPC/E水模型的平衡构型，但用来做TIP3P的初始填充完全没问题，后续通过平衡步骤会迅速弛豫到正确状态。
离子添加：既要中性化，又要生理盐浓度
这一步容易被忽视，但事实证明，离子浓度不合理会导致静电差分布异常，进而影响结合自由能的计算精度。先创建一个简单的ions.mdp，只用来生成TPR，不做真正的能量最小化：
; ions.mdp - only for grompp/genion
integrator = steep
nsteps = 1000
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb = 1.2
rvdw = 1.2
pbc = xyz
执行离子化命令：
gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr -maxwarn 2
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.15 -seed 42
设置 -seed 42 是为了让随机替换过程可复现，这在后续重复实验或调试时非常有用。
能量最小化：别只盯着最大步数
很多人喜欢将 nsteps 设成很大的值，比如50000步，但实际上收敛的判断依据是 emtol。建议先跑一个短步数的测试，看看体系趋向收敛的速度。标准参数如下：
; em.mdp
integrator = steep
nsteps = 50000
emtol = 500.0
emstep = 0.001
nstlog = 100
nstenergy = 100
cutoff-scheme = Verlet
nstlist = 20
ns_type = grid
verlet-buffer-tolerance = 0.005
coulombtype = PME
pme_order = 4
fourierspacing = 0.12
rcoulomb = 1.2
vdwtype = Cut-off
rvdw = 1.2
rvdw-switch = 1.0
DispCorr = EnerPres
constraints = none
执行：
gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr -maxwarn 1
gmx mdrun -deffnm em -v -ntmpi 1 -ntomp 16
如果发现 emtol 在几千步内就达到了，那是理想情况；如果跑完50000步还没收敛，说明结构中有较大空间冲突，这时需要检查起始结构是否合理。
NVT和NPT平衡：两个阶段都不能省
NVT阶段的作用是让温度均匀分布，NPT阶段则是让密度和压力稳定下来。两个阶段的参数文件在下述基础上略有差异：
NVT（控温、固定体积）：
; nvt.mdp
integrator = md
dt = 0.002
nsteps = 500000 ; 1 ns
nstxout = 5000
nstvout = 5000
nstfout = 5000
nstlog = 1000
nstenergy = 1000
nstcalcenergy = 100
constraints = h-bonds
constraint_algorithm = LINCS
lincs_iter = 1
lincs_order = 4
cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
pme_order = 4
fourierspacing = 0.12
rcoulomb = 1.2
rvdw = 1.2
DispCorr = EnerPres
tcoupl = V-rescale
tc-grps = Protein Non-Protein
tau_t = 0.1 0.1
ref_t = 300 300
pcoupl = no
pbc = xyz
gen_vel = yes
gen_temp = 300
gen_seed = -1
comm-mode = Linear
comm-grps = Protein Non-Protein
define = -DPOSRES
执行：
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr -maxwarn 1
gmx mdrun -deffnm nvt -v -ntmpi 1 -ntomp 16
NPT阶段相比NVT的两点关键变化：打开压力耦合，以及设为 continuation = yes 来续接NVT的速度参数：
; npt.mdp
integrator = md
dt = 0.002
nsteps = 500000 ; 1 ns
nstxout = 5000
nstvout = 5000
nstfout = 5000
nstlog = 1000
nstenergy = 1000
constraints = h-bonds
constraint_algorithm = LINCS
lincs_iter = 1
lincs_order = 4
cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
pme_order = 4
fourierspacing = 0.12
rcoulomb = 1.2
rvdw = 1.2
DispCorr = EnerPres
tcoupl = V-rescale
tc-grps = Protein Non-Protein
tau_t = 0.1 0.1
ref_t = 300 300
pcoupl = C-rescale
pcoupltype = isotropic
tau_p = 2.0
ref_p = 1.0
compressibility = 4.5e-5
pbc = xyz
continuation = yes
gen_vel = no
comm-mode = Linear
comm-grps = Protein Non-Protein
define = -DPOSRES
执行命令：
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr -maxwarn 1
gmx mdrun -deffnm npt -v -ntmpi 1 -ntomp 16
生产MD模拟：采样长度直接决定自由能计算的可靠性
对于多肽和LBP这种蛋白-配体体系，结合自由能计算（无论是MM-PBSA、GBSA还是FEP）都要求配体-蛋白界面有足够多的构象采样。100 ns是一条经验下线，如果体系较大或动态较强，建议延长到200-500 ns。下面是100 ns的生产模拟参数（md.mdp），关键参数已经标出：
; md.mdp - 100 ns
integrator = md
dt = 0.002
nsteps = 50000000 ; 100 ns
nstxout = 50000 ; 全精度的坐标写少了不影响分析
nstvout = 50000
nstfout = 50000
nstlog = 5000
nstenergy = 5000
nstcalcenergy = 100
nstxout-compressed = 5000 ; XTC每10 ps写一帧
compressed-x-grps = System
constraints = h-bonds
constraint_algorithm = LINCS
lincs_iter = 1
lincs_order = 4
cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
pme_order = 4
fourierspacing = 0.12
rcoulomb = 1.2
rvdw = 1.2
DispCorr = EnerPres
tcoupl = V-rescale
tc-grps = Protein Non-Protein
tau_t = 0.1 0.1
ref_t = 300 300
pcoupl = C-rescale
pcoupltype = isotropic
tau_p = 2.0
ref_p = 1.0
compressibility = 4.5e-5
pbc = xyz
continuation = yes
gen_vel = no
comm-mode = Linear
comm-grps = Protein Non-Protein
如果计算资源有限，也可以先跑50 ns（nsteps = 25000000），但那样后续的收敛性检验就尤其重要。启动命令推荐用 screen 来避免终端意外断开：
screen -S md_prod1
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr -maxwarn 1
gmx mdrun -deffnm md -v -ntmpi 1 -ntomp 16 -gpu_id 0 2>&1
如果模拟中途意外中断，可以使用续跑命令：
gmx mdrun -deffnm md -cpi md.cpt -append -v -ntmpi 1 -ntomp 16 -gpu_id 0 -maxh 47.5 | tee -a md_run.log
轨迹后处理
生产模拟结束后，先用 gmx trjconv 对轨迹做去周期性处理和聚类分析：
gmx trjconv -f md.xtc -s md.tpr -o cluster.xtc -pbc cluster -n index.ndx
这里 -pbc cluster 会将分子重新聚拢到主盒子中，消除由于周期性边界造成的“碎片化”假象，对后续的氢键分析、RMSD计算和结合自由能提取都有直接影响。
最后再强调一句：以上流程仅仅是一个标准MD模板，实际场景中还需要根据具体的蛋白-配体变化做很多针对性调整——比如是否需要溶剂边界修正、是否要引入伞形采样或者自由能微扰来精确计算结合自由能，这些都需要对照你的研究目标来决定。
