1. 为什么伞形采样值得你花时间啃下来做分子动力学模拟的人迟早会碰到一个绕不开的需求某个分子要从一个位置移动到另一个位置中间隔着一道能量势垒常规模拟在有限时间内根本翻不过去。比如配体从蛋白口袋中解离、离子穿过膜通道、构象从一种稳态翻到另一种稳态这些过程的自由能变化才是真正有价值的信息但直接跑模拟你只能看到分子在原地抖动。伞形采样就是解决这类问题的经典手段。它的核心思路并不复杂沿着你关心的反应坐标人为施加一系列简谐势也就是“伞”把体系限制在反应坐标的不同窗口里让每个窗口都能充分采样最后再用WHAMWeighted Histogram Analysis Method或类似方法把各窗口的数据拼起来还原出整条自由能曲线PMFPotential of Mean Force。我接触Gromacs伞形采样有好几年了从最开始照着教程跑通、到后来自己设计pull code、调窗口、处理WHAM报错踩过的坑可以说能写一本小册子。这篇内容就是把我这些年积累的实战经验整理出来从原理到操作到排错尽量讲透。适合已经会跑基础Gromacs模拟、想进阶做自由能计算的读者也适合跑过一两次伞形采样但结果总是不对劲、想搞清楚问题出在哪的人。整篇内容围绕一条主线怎么用Gromacs把伞形采样做对而不是仅仅做出来。做出来和做对之间差的就是那些教程里不会写的细节。2. 伞形采样的整体设计与方案选型2.1 反应坐标怎么选决定了后面一半的成败伞形采样最重要的一步不是写mdp文件而是选反应坐标Reaction Coordinate。反应坐标选错了后面窗口拉得再漂亮、WHAM收敛得再好得到的PMF也没有物理意义。反应坐标的本质是你要用一个或少数几个几何变量来描述整个过程的进展。常见的选择有距离两个原子组之间的质心距离适合配体解离、离子穿透这类过程。二面角适合构象翻转、旋转异构化。配位数适合描述溶剂化壳层变化、离子配对。向量投影沿某个方向的位移分量。我个人的经验是能用距离就别用角度能用一维就别上二维。原因很直接一维反应坐标的窗口数量、采样时间、WHAM收敛难度都是可控的一旦上二维窗口数呈平方增长计算量直接爆炸而且WHAM的收敛判断会变得非常痛苦。选反应坐标时有一个自检标准沿着这条坐标体系从初态到末态应该是单调推进的中间不应该出现坐标值回退但物理状态也回退的歧义。举个反例如果你用配体到蛋白某个残基的距离做坐标但配体解离路径上会先靠近另一个残基再远离那这条坐标就不是单调的PMF会出现假势垒。2.2 窗口数量和分布别拍脑袋定窗口怎么分布是新手最容易拍脑袋的地方。常见错误是均匀撒20个窗口就开跑结果发现势垒附近采样不够PMF在关键区域抖得厉害。正确的做法是先跑一条短时间的探索性模拟粗略估计自由能曲线的形状再据此分配窗口。具体来说势能变化平缓的区域窗口可以稀疏一些比如每0.1 nm一个窗口。势垒顶部或曲线陡峭的区域窗口要加密可能每0.02到0.05 nm就要一个窗口。相邻窗口的伞势要有充分重叠一般要求相邻窗口的直方图重叠区域不少于总采样范围的20%到30%。窗口数量方面一维反应坐标通常需要15到40个窗口具体取决于反应坐标的总跨度和曲线复杂度。跨度大、曲线复杂就多设窗口。我做过的一个配体解离体系反应坐标跨度约1.2 nm最后用了32个窗口势垒区域每0.03 nm一个两端每0.08 nm一个效果比较理想。2.3 力常数怎么定太大太小都是坑伞形势的力常数pull-coord-k直接决定了每个窗口的采样范围。力常数太大体系被死死摁在一个点上直方图很窄相邻窗口重叠不够力常数太小体系在窗口里乱跑甚至翻过势垒跑到别的窗口区域采样效率极低。一个实用的经验公式是力常数大致取 2 到 5 倍的热能除以你期望的窗口宽度平方。在Gromacs单位制下300 K时kT约2.5 kJ/mol。如果你希望每个窗口的采样宽度约0.05 nm那力常数大约在 2×2.5/0.05² 2000 kJ/mol/nm² 到 5×2.5/0.05² 5000 kJ/mol/nm² 之间。我通常先用2000到3000 kJ/mol/nm²试跑看直方图宽度和重叠情况再调整。如果直方图太窄宽度小于0.02 nm就减小力常数如果直方图太宽甚至出现双峰说明力常数太小体系在窗口内跑到了不该去的地方。2.4 采样时间不是跑得越久越好而是要跑够每个窗口要跑多久取决于体系的复杂度和反应坐标的扩散特性。一个粗略的判断标准是每个窗口的采样时间至少要让体系在窗口内来回穿越几次。实际操作中我一般先每个窗口跑1到2 ns做测试看直方图是否平滑、相邻窗口是否重叠。如果直方图有明显毛刺或者重叠不足要么加时间要么调窗口。确认没问题后再延长到5到10 ns做正式采样。对于复杂的生物体系每个窗口10到20 ns也不罕见。这里有个容易被忽略的点采样时间不是均匀分配就合理。势垒附近的窗口往往需要更长的采样时间因为那里的体系动力学更慢需要更久才能充分探索局部空间。我通常会给势垒附近的窗口多分配30%到50%的时间。3. 核心细节解析与实操要点3.1 pull code 的几种模式到底该用哪个Gromacs的pull code支持多种模式伞形采样常用的是pull-coord-type umbrella配合pull-coord-geometry指定几何类型。但很多人不清楚pull-coord-geometry的几个选项到底有什么区别几何类型含义适用场景distance两组原子质心间的距离配体解离、离子穿透direction沿某方向的投影距离需要指定方向的拉伸cylinder圆柱半径限制在圆柱内的径向距离angle角度角度类反应坐标dihedral二面角构象翻转大部分伞形采样用distance就够了。但如果你做的是膜通道穿透可能需要用direction配合pull-coord-vec指定膜法线方向这样反应坐标才是沿通道方向的位移而不是简单的质心距离。还有一个关键参数是pull-coord-start和pull-coord-init。前者是参考位置后者是初始位置。在伞形采样中每个窗口的pull-coord-init应该设成该窗口的目标位置pull-coord-start也设成同一个值这样伞势的中心就在目标位置上。3.2 别让旋转和平移毁掉你的采样这是我最想强调的一个坑伞形采样中如果你不对体系做位置约束整个分子可能在模拟过程中整体漂移或旋转导致反应坐标的定义失效。具体来说如果你用两个原子组的质心距离做反应坐标但其中一个组所在的分子整体在盒子里乱跑那这个距离的变化就不完全是你要研究的构象变化还混入了整体平移。解决办法是在mdp里加上comm-mode linear去掉质心平移并且对参考组加位置约束pull-coord-group配合pull-group的设置。对于旋转如果反应坐标涉及方向还需要考虑用pull-coord-vec显式指定方向向量而不是依赖体系自身的取向。我见过太多案例PMF跑出来形状诡异最后发现是体系在模拟中整体翻转了。3.3 初始构型怎么生成别用同一个构型复制粘贴每个窗口的初始构型不能简单地把同一个构型复制过来改一下pull-coord-init就完事。如果所有窗口都从同一个构型出发势垒附近的窗口可能因为初始构型离目标位置太远在伞势的拉扯下产生很大的初始应力导致体系局部变形甚至崩溃。正确的做法是用 steered MDSMD或缓慢拉动的方式生成一条从初态到末态的连续轨迹然后从这条轨迹上按窗口位置截取构型作为各窗口的初始结构。Gromacs里可以用pull-coord-rate设置一个较小的拉动速度比如0.001 nm/ps跑一条SMD轨迹然后用gmx trjconv按时间截取。如果不想跑SMD至少也要用能量最小化加短时间平衡的方式让每个窗口的体系在伞势下先弛豫一段时间再开始正式采样。我通常每个窗口先跑500 ps到1 ns的平衡确认反应坐标值稳定在目标位置附近后再开始计时。3.4 WHAM的输入文件别小看直方图的数量WHAM需要每个窗口的直方图数据作为输入。Gromacs跑完伞形采样后用gmx wham可以直接读取pullx或pulll文件。但这里有几个细节确保每个窗口的采样都达到了足够的统计量。如果某个窗口的直方图只有几百个点WHAM会给出很大的统计误差。检查直方图的重叠。可以用gmx wham -hist生成直方图文件画出来看相邻窗口是否有重叠。没有重叠的窗口对WHAM来说就是断开的PMF会出现不连续的跳变。注意单位。Gromacs的pull输出默认单位是nmWHAM输出的PMF单位是kJ/mol。如果你用的力常数单位是kJ/mol/nm²那PMF的单位就是对的不需要额外转换。WHAM的收敛判断不能只看PMF曲线是否平滑还要看bootstrap误差。gmx wham -bs-method b-hist可以做bootstrap分析给出PMF的误差棒。如果误差棒在势垒区域超过2到3 kJ/mol说明采样还不够需要加时间或加窗口。4. 实操过程与核心环节实现4.1 体系准备与平衡地基没打好后面全白搭伞形采样的体系准备和普通MD没有本质区别但有几个额外注意点盒子要足够大。反应坐标跨度加上两端各留1 nm以上的缓冲避免周期性边界条件下分子和自己镜像相互作用。我一般让盒子边长至少是反应坐标跨度的2.5倍。溶剂化要充分。用gmx solvate加水后检查盒子边缘是否有真空层。如果有用gmx editconf调整盒子尺寸或补加水。离子中和。加Na或Cl-中和体系电荷浓度控制在0.15 mol/L左右模拟生理条件。能量最小化。用最陡下降法跑5000步以上确保最大力小于1000 kJ/mol/nm。NVT和NPT平衡。先NVT平衡100 ps稳住温度再NPT平衡100 ps稳住压力。这里可以用Parrinello-Rahman或Berendsen控压但正式采样建议用Parrinello-Rahman或C-rescale控压更准确。4.2 生成伞形采样窗口的初始构型假设你已经跑完平衡现在要生成各窗口的初始构型。我的标准流程是第一步跑一条短SMD轨迹。在mdp里设置pull yes pull-ngroups 2 pull-ncoords 1 pull-coord1-type umbrella pull-coord1-geometry distance pull-coord1-groups 1 2 pull-coord1-dim N N N pull-coord1-start yes pull-coord1-init 0 pull-coord1-rate 0.001 pull-coord1-k 1000这里pull-coord1-rate 0.001表示每ps拉动0.001 nm跑1 ns就能拉1 nm。力常数用1000 kJ/mol/nm²比正式采样小一些让体系有足够时间弛豫。第二步从SMD轨迹截取构型。用gmx trjconv按时间截取每个窗口对应一个时间点。比如你要32个窗口反应坐标从0到1.2 nm那每个窗口间隔约0.0387 nm对应SMD轨迹上约38.7 ps的时间点。第三步检查截取构型。用gmx check确认每个构型的原子数和盒子尺寸一致用VMD或PyMOL目视检查反应坐标值是否接近目标位置。4.3 各窗口的mdp配置与运行每个窗口的mdp文件除了pull相关参数不同其他都一样。关键参数pull yes pull-ngroups 2 pull-ncoords 1 pull-coord1-type umbrella pull-coord1-geometry distance pull-coord1-groups 1 2 pull-coord1-dim N N N pull-coord1-start yes pull-coord1-init 窗口目标位置 pull-coord1-rate 0 pull-coord1-k 2500 pull-coord1-vec 0 0 0注意pull-coord1-rate 0因为伞形采样是固定窗口位置采样不需要继续拉动。pull-coord1-k用2500 kJ/mol/nm²这是我常用的起始值。运行命令gmx grompp -f umbrella.mdp -c window_XX.gro -p topol.top -o umbrella_XX.tpr gmx mdrun -deffnm umbrella_XX -pf pullf_XX.xvg -px pullx_XX.xvg-pf输出力文件-px输出位置文件WHAM两个都需要。4.4 WHAM分析从数据到PMF跑完所有窗口后用gmx wham做分析gmx wham -it tpr_files.dat -if pullf_files.dat -o pmf.xvg -hist histo.xvg -bs-method b-hist -bsres bs_result.xvg其中tpr_files.dat列出所有窗口的tpr文件路径pullf_files.dat列出所有pullf文件路径。-hist输出直方图-bs-method b-hist做bootstrap误差分析。跑完后用xmgrace或Python画PMF曲线和误差棒。检查几个点PMF曲线是否平滑有没有不连续跳变。误差棒在势垒区域是否小于3 kJ/mol。直方图重叠是否充分。如果PMF在某个区域出现不连续通常是那个区域的窗口重叠不够需要补窗口或加采样时间。5. 常见问题与排查技巧实录5.1 直方图不重叠怎么办这是最常见的问题。表现是WHAM跑出来的PMF在某些位置出现台阶或跳变。排查步骤用gmx wham -hist生成直方图画出来看哪两个相邻窗口没有重叠。检查这两个窗口的力常数是否太大导致直方图太窄。如果是减小力常数重跑。检查窗口间距是否太大。如果是在中间补窗口。检查采样时间是否足够。如果直方图毛刺很多加时间。我遇到过一个案例两个相邻窗口的直方图完全不重叠中间隔了0.1 nm的空白。最后发现是力常数用了5000 kJ/mol/nm²直方图宽度只有0.015 nm而窗口间距是0.05 nm。把力常数降到2000后问题解决。5.2 PMF曲线出现假势垒假势垒的表现是PMF在某个位置出现一个不该有的峰但物理上那里不应该有势垒。常见原因反应坐标选得不好。比如用距离做坐标但路径上分子会绕路导致距离先减小后增大再减小PMF就会出现假峰。窗口采样不充分。势垒附近的窗口如果采样不够WHAM会给出偏高的自由能。体系在窗口内发生了非目标构象变化。比如配体在窗口内翻转了取向导致反应坐标值相同但物理状态不同。解决办法先检查反应坐标是否单调再检查势垒附近窗口的采样时间和直方图质量。如果都正常可能需要换反应坐标。5.3 WHAM报错“No data in histogram”这个报错通常是因为pullf或pullx文件为空或者格式不对。排查确认mdrun运行时确实加了-pf和-px参数。检查pullf文件是否有内容用wc -l看行数。确认tpr文件和pullf文件对应的是同一个窗口。检查mdp里pull-coord1-start是否设为yes如果设为nopull不会在第一步就生效。5.4 采样过程中体系崩溃表现是mdrun报错“LINCS warnings”或“Segmentation fault”。原因通常是初始构型应力太大或者伞势拉扯太猛。解决办法用SMD生成的构型而不是直接改init值。先跑短时间的平衡力常数从500开始逐步加到目标值。检查盒子是否足够大避免分子和镜像相互作用。用gmx energy检查温度和压力是否稳定。5.5 常见问题速查表问题可能原因解决办法直方图不重叠力常数太大、窗口间距太大减小力常数、补窗口PMF有台阶窗口重叠不足补窗口、加采样时间PMF有假势垒反应坐标不单调、采样不足换坐标、加时间WHAM报错无数据pullf为空、tpr不匹配检查运行参数、文件对应体系崩溃初始构型应力大、盒子太小用SMD构型、加大盒子PMF误差棒大采样时间不够加时间、加窗口5.6 几个我踩过的坑和独家技巧坑一忘了设pull-coord1-start yes。这个参数如果不设pull在第一步不会生效体系会先自由跑一段导致初始构型和目标位置偏离。我因为这个参数浪费过整整两天。坑二用同一个构型跑所有窗口。早期图省事把所有窗口的gro文件都复制同一个只改init值。结果势垒附近的窗口在伞势拉扯下局部变形PMF完全不可信。后来老老实实跑SMD生成构型问题再没出现过。技巧一用gmx wham -nBootstrap 200做更多次bootstrap。默认是100次增加到200次误差估计更稳定代价是计算时间稍长。技巧二势垒附近窗口多给50%时间。这个前面提过但值得再强调。势垒附近的慢动力学是PMF误差的主要来源多给时间比全局加时间更划算。技巧三用pull-coord1-vec显式指定方向。如果你的反应坐标涉及方向不要依赖体系自身取向显式指定向量可以避免整体旋转带来的噪声。技巧四定期检查pullx文件里的反应坐标值。如果某个窗口的反应坐标值在采样过程中持续漂移说明伞势没拉住需要检查力常数或初始构型。6. 从数据到结论怎么判断你的PMF可信跑完WHAM拿到PMF曲线怎么判断结果可不可信我通常从四个维度检查第一看误差棒。势垒区域的误差棒如果小于2 kJ/mol基本可信2到4 kJ/mol需要谨慎解读大于4 kJ/mol说明采样严重不足结论不能要。第二看直方图重叠。所有相邻窗口的直方图重叠区域应该不少于20%。如果某个位置重叠只有5%那里的PMF就是不可信的。第三看PMF的对称性。如果你做的是两个方向的可逆过程正反两个方向的PMF应该基本一致。如果差异很大说明采样有系统性偏差。第四看物理合理性。PMF的势垒高度应该和实验值或其他计算结果大致吻合。如果差了一个数量级那肯定是哪里出了问题。我个人的习惯是正式跑之前先用短时间低精度的参数跑一遍全流程确认流程通、结果合理再上长时间高精度。这样即使有问题也能快速定位不至于跑了一周才发现参数错了。伞形采样这个技术入门不难做对不容易。核心就三件事反应坐标选对、窗口和力常数调好、采样跑够。剩下的就是耐心和细心。我见过太多人跑完WHAM看到一条曲线就以为完事了结果曲线上的假势垒和巨大误差棒全被忽略。多花半小时检查直方图和误差棒比多跑三天模拟更有价值。