企业官网建设流程全解析

1. 为什么伞形采样值得你花时间啃下来

做分子动力学模拟的人，迟早会碰到一个绕不开的需求：某个分子要从一个位置移动到另一个位置，中间隔着一道能量势垒，常规模拟在有限时间内根本翻不过去。比如配体从蛋白口袋中解离、离子穿过膜通道、构象从一种稳态翻到另一种稳态，这些过程的自由能变化才是真正有价值的信息，但直接跑模拟你只能看到分子在原地抖动。

伞形采样就是解决这类问题的经典手段。它的核心思路并不复杂：沿着你关心的反应坐标，人为施加一系列简谐势（也就是“伞”），把体系限制在反应坐标的不同窗口里，让每个窗口都能充分采样，最后再用WHAM（Weighted Histogram Analysis Method）或类似方法把各窗口的数据拼起来，还原出整条自由能曲线（PMF，Potential 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配合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 MD（SMD）或缓慢拉动的方式，生成一条从初态到末态的连续轨迹，然后从这条轨迹上按窗口位置截取构型作为各窗口的初始结构。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输出默认单位是nm，WHAM输出的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. 盒子要足够大。反应坐标跨度加上两端各留1 nm以上的缓冲，避免周期性边界条件下分子和自己镜像相互作用。我一般让盒子边长至少是反应坐标跨度的2.5倍。
2. 溶剂化要充分。用gmx solvate加水后，检查盒子边缘是否有真空层。如果有，用gmx editconf调整盒子尺寸或补加水。
3. 离子中和。加Na+或Cl-中和体系电荷，浓度控制在0.15 mol/L左右模拟生理条件。
4. 能量最小化。用最陡下降法跑5000步以上，确保最大力小于1000 kJ/mol/nm。
5. 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在某些位置出现台阶或跳变。排查步骤：

1. 用gmx wham -hist生成直方图，画出来看哪两个相邻窗口没有重叠。
2. 检查这两个窗口的力常数是否太大，导致直方图太窄。如果是，减小力常数重跑。
3. 检查窗口间距是否太大。如果是，在中间补窗口。
4. 检查采样时间是否足够。如果直方图毛刺很多，加时间。

我遇到过一个案例，两个相邻窗口的直方图完全不重叠，中间隔了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文件为空，或者格式不对。排查：

1. 确认mdrun运行时确实加了-pf和-px参数。
2. 检查pullf文件是否有内容，用wc -l看行数。
3. 确认tpr文件和pullf文件对应的是同一个窗口。
4. 检查mdp里pull-coord1-start是否设为yes，如果设为no，pull不会在第一步就生效。

5.4 采样过程中体系崩溃

表现是mdrun报错“LINCS warnings”或“Segmentation fault”。原因通常是初始构型应力太大，或者伞势拉扯太猛。解决办法：

• 用SMD生成的构型而不是直接改init值。
• 先跑短时间的平衡，力常数从500开始逐步加到目标值。
• 检查盒子是否足够大，避免分子和镜像相互作用。
• 用gmx energy检查温度和压力是否稳定。

5.5 常见问题速查表

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看到一条曲线就以为完事了，结果曲线上的假势垒和巨大误差棒全被忽略。多花半小时检查直方图和误差棒，比多跑三天模拟更有价值。