1. 为什么伞形采样值得你花时间折腾做分子动力学模拟的人迟早会碰到一个绕不开的需求想知道某个分子穿过某个通道、或者两个分子互相靠近时自由能到底怎么变。普通平衡态模拟里这种事件要么发生得太快根本抓不住要么发生得太慢跑几百纳秒都看不到一次。伞形采样就是专门用来解决这类问题的增强采样方法它通过在反应坐标上施加一系列简谐势也就是“伞”把体系强行约束在坐标轴的不同位置上从而把整条自由能曲线一段一段地“拼”出来。Gromacs 作为分子动力学领域使用最广的开源软件之一对伞形采样的支持相当成熟配套的gmx wham工具可以直接把多窗口的采样结果重加权成 PMFPotential of Mean Force平均力势。但问题在于从体系搭建到最终出图中间有太多容易踩坑的地方pull code 怎么配、窗口怎么分布、采样多久算够、WHAM 报错怎么排查……每一步都可能让你卡上几天。这篇内容适合已经跑过基础分子动力学、想上手伞形采样但被各种报错劝退的人也适合跑过几次但结果总是不太对、想系统梳理一遍流程的人。我会把整个流程拆成设计思路、核心配置、实操步骤、问题排查几个部分尽量把每个参数背后的逻辑讲清楚而不是只丢一堆命令让你照抄。2. 伞形采样的整体设计与思路拆解2.1 伞形采样到底在算什么先把这个方法的物理图像说清楚。假设你想算一个配体从蛋白质口袋往外拉出来的自由能变化反应坐标就是配体与口袋之间的距离。在自然状态下配体大部分时间待在口袋里你几乎采不到它跑到外面的构象。伞形采样的做法是在距离轴上取一系列位置比如 0.3、0.4、0.5……一直到 2.0 纳米每个位置放一个独立的模拟窗口每个窗口里加一根“弹簧”把配体往这个位置拽。弹簧的劲度系数决定了约束的松紧太软了配体会跑偏太硬了采样效率低。每个窗口跑完之后你会得到该窗口内反应坐标的分布以及弹簧施加的力。WHAM 的核心思想就是利用这些有偏采样数据反推出无偏的自由能曲线。它本质上是一个自洽迭代的直方图重加权方法把所有窗口的数据放在一起找到一组权重因子使得各个窗口的分布能够拼成一条连续的曲线。这里有个关键点很多人一开始不理解为什么不能只跑一个窗口把弹簧慢慢移动那样做确实存在叫“牵引模拟”或者“steered MD”但它得到的是非平衡过程的功需要额外做 Jarzynski 平均才能转成自由能收敛慢且方差大。伞形采样是平衡态方法每个窗口单独平衡统计上更稳健这也是它成为主流的原因。2.2 反应坐标的选择逻辑反应坐标选得好不好直接决定整个项目能不能出结果。好的反应坐标应该满足几个条件能区分反应物和产物状态、在物理上合理、采样过程中不会出现其他慢自由度干扰。以配体拉出为例常用的反应坐标有配体与结合位点残基的质心距离、配体沿某条通道的投影距离、或者配体与蛋白质整体质心的距离。质心距离最直观但如果通道是弯曲的直线距离就不能真实反映路径这时候可能需要用路径投影或者多个集体变量组合。我个人的经验是第一次做某个体系时先用牵引模拟快速拉一遍看看配体到底沿什么路径出去、中间有没有卡住的地方。这一步花不了多少机时但能帮你判断反应坐标选得对不对。如果牵引过程中配体撞墙或者绕路那说明你选的坐标有问题直接上伞形采样只会浪费更多时间。2.3 窗口分布与弹簧常数的取舍窗口怎么排、弹簧多硬这两个参数是耦合的。基本原则是相邻窗口的采样分布要有重叠一般建议重叠区域至少占每个窗口分布的 20% 到 30%。如果窗口间距太大分布不重叠WHAM 就无法把相邻窗口连起来PMF 会出现断裂或者跳变。弹簧常数 k 的选择要结合热涨落来考虑。在温度 T 下反应坐标的涨落方差大约是 kT/k。如果你希望涨落标准差在 0.05 纳米左右室温下 kT 约 2.5 kJ/mol那么 k 大约在 1000 kJ/mol/nm² 量级。实际用的时候我一般先用 1000 到 2000 之间的值试跑看直方图分布宽度再决定要不要调整。窗口间距可以按 0.05 到 0.1 纳米来设具体取决于你关心的坐标范围。如果总跨度是 2 纳米用 0.1 纳米间距就是 20 个窗口用 0.05 就是 40 个。窗口越多每个窗口的采样负担越轻但总机时也越多。我的建议是先用粗间距跑一遍看整体形状再在关键区域加密。3. 核心配置与实操要点3.1 拓扑与坐标文件的准备伞形采样对拓扑文件的要求和普通模拟基本一致但有一点要特别注意如果你用的是联合原子力场或者粗粒化力场反应坐标的定义要确保在拓扑里有对应的原子组。比如你要拉一个配体配体的原子索引必须明确不能混在蛋白质组里。坐标文件方面起始结构最好来自平衡过的体系。我见过有人直接拿能量最小化后的结构就开始拉结果配体还没出口袋蛋白质就散了。正确的做法是先跑一段 NPT 平衡让体系密度和盒子尺寸稳定下来再作为伞形采样的起点。另外盒子尺寸要足够大。如果你要拉的距离是 2 纳米盒子在对应方向上至少要有 4 纳米以上的长度否则配体拉到一半就撞到周期性镜像了。这一点在设置盒子时就要考虑进去不要等跑完了才发现。3.2 pull code 的关键参数逐条拆解Gromacs 的 pull code 写在 mdp 文件里核心参数有这么几个pull yes pull-ngroups 2 pull-ncoords 1 pull-group1-name Protein pull-group2-name Ligand pull-coord1-type umbrella pull-coord1-geometry distance pull-coord1-groups 1 2 pull-coord1-dim N N Y pull-coord1-init 0.5 pull-coord1-rate 0 pull-coord1-k 1000pull-coord1-geometry决定了反应坐标怎么算。distance是两组质心距离direction是沿某方向的投影cylinder是柱坐标距离。选哪个取决于你的体系配体拉出一般用distance或direction。pull-coord1-init是窗口的初始位置每个窗口这个值不同。pull-coord1-rate 0表示不做牵引只做静态约束这是伞形采样的标准做法。如果你设了非零速率那就变成牵引模拟了。pull-coord1-k就是弹簧常数。单位是 kJ/mol/nm²。这个值不能太小否则约束不住也不能太大否则采样效率低。前面说的 1000 到 2000 是常用范围。pull-coord1-dim控制哪些方向参与计算。比如你只想沿 z 轴拉就设N N Y。这个参数很容易被忽略但如果设错了反应坐标会包含你不想要的方向分量结果就偏了。3.3 各窗口的初始化与平衡策略每个窗口的起始结构可以从同一个平衡结构生成但初始位置不同。做法是用gmx grompp配合不同的 mdp 文件每个 mdp 里改pull-coord1-init的值。或者更高效的方式是用gmx genrestr之类的工具批量生成但最稳妥的还是脚本化处理。每个窗口开始采样前建议先做一段短平衡让体系适应该窗口的约束。平衡时间不用太长100 到 500 皮秒通常够用但具体要看体系。判断标准是反应坐标的均值稳定在设定值附近且涨落不再有漂移趋势。采样阶段每个窗口一般跑 5 到 20 纳秒。总机时等于窗口数乘以单窗口时间。如果窗口多可以考虑用 GPU 并行跑Gromacs 对多窗口并行支持很好每个窗口一个 GPU 或者多个窗口共享一个 GPU 都可以。注意每个窗口的采样必须保存反应坐标的轨迹也就是 pull 相关的输出。在 mdp 里要确保pull-nstxout和pull-nstfout设置合理一般每 500 到 1000 步输出一次就够了太密会拖慢速度太疏会影响 WHAM 的统计精度。3.4 WHAM 输入文件的生成WHAM 需要每个窗口的 pull 输出文件通常是pullx.xvg或pullf.xvg。gmx wham的命令行大致是这样gmx wham -it tpr_files.dat -if pullf_files.dat -o pmf.xvg -hist histo.xvg其中tpr_files.dat是每行一个 tpr 文件路径pullf_files.dat是每行一个 pullf 文件路径。这两个文件的顺序必须一一对应不能错位。WHAM 还有一些可选参数比如-b和-e指定起止时间-temp指定温度-bins指定直方图分箱数。分箱数一般设 100 到 200太少会平滑掉细节太多会有噪声。如果要做 bootstrap 误差分析可以加-nBootstrap参数WHAM 会自动做重采样并给出误差估计。这个功能很实用建议每次都开。4. 完整实操流程与关键环节4.1 从平衡结构到窗口结构生成假设你已经有一个平衡好的eq.gro和eq.tpr接下来要生成 20 个窗口的起始结构。最直接的方式是写一个 bash 脚本循环修改 mdp 里的pull-coord1-init然后依次跑grompp。for i in $(seq 0 19); do init$(echo 0.3 $i * 0.1 | bc) sed s/pull-coord1-init .*/pull-coord1-init $init/ template.mdp window_$i.mdp gmx grompp -f window_$i.mdp -c eq.gro -p topol.top -o window_$i.tpr -maxwarn 1 done这里-maxwarn 1是为了跳过一些无害的警告比如盒子尺寸和截断半径的关系。但不要滥用这个参数如果警告涉及物理设置错误还是要认真处理。生成 tpr 之后每个窗口先跑一段平衡。平衡的 mdp 和采样 mdp 可以共用只是把nsteps设短一些。平衡完之后用-t参数从平衡的 cpt 文件继续跑采样。4.2 采样阶段的并行与监控采样阶段最耗机时合理利用并行很重要。如果机器上有多个 GPU可以用-multidir或者简单的 shell 并行。比如for i in $(seq 0 19); do (gmx mdrun -deffnm window_$i -s window_$i.tpr -gpu_id $((i % 4)) ) done wait这段脚本把 20 个窗口分配到 4 块 GPU 上并行跑。注意gpu_id的分配要均匀不要让某块 GPU 过载。跑的过程中要监控反应坐标是否稳定。可以用gmx energy或者直接看 pullx 文件。如果某个窗口的坐标均值偏离设定值超过 0.1 纳米说明弹簧太软或者平衡不够需要调整。4.3 WHAM 执行与 PMF 出图采样完成后收集所有窗口的 pullf 文件生成列表ls window_*/pullf.xvg pullf_files.dat ls window_*/window_*.tpr tpr_files.dat然后跑 WHAMgmx wham -it tpr_files.dat -if pullf_files.dat -o pmf.xvg -hist histo.xvg -nBootstrap 100 -bins 200跑完之后pmf.xvg就是自由能曲线。可以用 Grace 或者 Python 的 matplotlib 画图。检查曲线是否平滑、两端是否收敛、误差棒是否合理。如果曲线中间有断裂或者跳变大概率是窗口重叠不够。这时候要么增加窗口要么增大弹簧常数让分布更宽。如果两端翘起可能是采样不够需要延长端部窗口的时间。4.4 结果验证与收敛性检查PMF 跑出来不代表就完事了还要做收敛性检查。常用的方法有几种一是把总采样时间分成前后两半分别做 WHAM看两条曲线是否重合二是看 bootstrap 误差是否在可接受范围三是检查每个窗口的直方图是否有足够的重叠。我一般会画一张直方图重叠图把所有窗口的分布画在一起直观地看相邻窗口是否有交集。如果某两个窗口之间明显断开那一段的 PMF 就不可信。另外反应坐标的选取是否合理也要回头验证。如果 PMF 曲线出现不物理的振荡或者平台可能是坐标里混入了其他自由度。这时候可能需要换坐标或者加限制。5. 常见问题与排查技巧实录5.1 WHAM 报错与直方图问题WHAM 最常见的报错是“no overlap between windows”或者“histogram is empty”。前者说明相邻窗口分布不重叠后者说明某个窗口根本没采到数据。排查思路先看 pullx 文件里反应坐标的范围确认每个窗口是否真的在设定位置附近采样。如果某个窗口的坐标跑到了很远的地方说明弹簧太软或者初始结构没平衡好。如果坐标范围正常但直方图空可能是分箱数设得不对或者输出频率太低导致数据点太少。解决办法增大弹簧常数、增加窗口密度、延长采样时间、调整分箱数。这几个手段可以组合使用。5.2 反应坐标漂移与约束失效有时候跑着跑着反应坐标慢慢偏离设定值最后完全失控。这种情况通常是弹簧太软加上体系里有其他力在干扰。比如配体在口袋里受到蛋白质的吸引力如果弹簧不够硬配体就会被拉回口袋。解决办法是增大 k 值或者在窗口初始化时把配体放在更靠近设定位置的地方。另外检查pull-coord1-dim是否设对如果方向设错了约束力会作用在错误的方向上。还有一种可能是周期性边界条件导致的。如果配体拉到了盒子边缘和镜像发生了相互作用坐标就会跳变。这时候需要加大盒子或者用pull-coord1-geometry direction配合合适的参考向量。5.3 PMF 曲线不收敛的典型表现PMF 不收敛的表现有很多种曲线两端翘起、中间有深谷、误差棒很大、前后半段不一致。每种表现对应的原因不同。两端翘起通常是端部窗口采样不足因为端部窗口的构象空间更大需要更长时间平衡。中间有深谷可能是反应坐标经过了一个真实的能垒但也可能是采样不够导致的假象。误差棒大说明统计量不够需要延长采样或者增加窗口。我自己的经验是先跑短时间看整体形状确定形状合理后再延长采样。如果形状本身就不对延长时间也没用得先改设置。5.4 常见问题速查表问题现象可能原因排查方法解决手段WHAM 报无重叠窗口间距太大或弹簧太软画直方图重叠图加密窗口或增大 k反应坐标漂移弹簧太软或方向设错检查 pullx 轨迹增大 k 或修正 dimPMF 两端翘起端部采样不足分半收敛检查延长端部窗口时间直方图有空窗输出频率太低检查 pullx 数据点提高输出频率PMF 噪声大分箱数太多或采样不足调整 bins 看变化减少 bins 或延长采样约束力异常大初始结构不合理看 pullf 数值重新平衡起始结构5.5 几个容易被忽略的实操细节第一个细节是温度耦合和压力耦合的选择。伞形采样过程中体系被约束在特定构象如果压力耦合太强盒子尺寸会波动可能影响反应坐标的计算。我一般建议用较弱的压力耦合或者只在平衡阶段用采样阶段可以关掉压力耦合改用 NVT。第二个细节是输出频率和文件大小。pullx 和 pullf 文件如果输出太密会占用大量磁盘空间而且 WHAM 读起来也慢。一般每 500 到 1000 步输出一次足够。第三个细节是随机数种子。如果多个窗口用同一个种子可能会引入系统性偏差。Gromacs 默认会自动生成种子但如果你手动设了gen-seed要确保每个窗口不同。第四个细节是 WHAM 的-b和-e参数。如果你在采样初期有不平衡的部分可以用-b跳过。但不要跳太多否则统计量不够。一般跳过前 1 到 2 纳秒比较稳妥。6. 提升效率的进阶技巧6.1 用哈密顿量交换加速采样如果体系比较复杂单个窗口的采样可能很慢。这时候可以考虑用哈密顿量交换Hamiltonian replica exchange结合伞形采样也就是所谓的“伞形采样 副本交换”。Gromacs 本身不直接支持这个组合但可以通过外挂脚本实现或者用 PLUMED 这样的增强采样库。PLUMED 对伞形采样的支持更灵活可以定义更复杂的集体变量也支持多 walker 并行。如果你经常做伞形采样花点时间学 PLUMED 是值得的。6.2 自适应窗口分布传统的伞形采样窗口是均匀分布的但实际自由能曲线往往不均匀有些区域变化快有些区域平坦。自适应方法可以在采样过程中动态调整窗口位置把更多窗口放在自由能变化快的区域。Gromacs 本身没有内置自适应伞形采样但可以通过脚本实现简单的自适应先跑一轮粗采样根据 PMF 的斜率重新分布窗口再跑第二轮。这样可以用更少的窗口达到同样的精度。6.3 结果的可视化与报告PMF 出图之后怎么呈现也很重要。我一般会画三张图第一张是 PMF 曲线加误差棒第二张是直方图重叠图第三张是收敛性检查图。这三张图放在一起审稿人或者合作者一眼就能看出结果可不可靠。画图工具方面Python 的 matplotlib 最灵活可以完全自定义。如果不想写代码Grace 也可以但样式调整比较麻烦。xvg 文件可以直接用xmgrace打开也可以转成 csv 再用其他工具画。提示PMF 的零点可以任意选一般把反应物状态设为零。但报告的时候要说明零点在哪里否则别人看不懂。6.4 机时估算与资源分配最后说一下机时估算。假设你有 20 个窗口每个窗口跑 10 纳秒体系有 5 万个原子用一块中端 GPU 大概每纳秒 2 到 3 小时那么总机时大约是 20 乘以 10 乘以 2.5等于 500 GPU 小时。如果只有一块 GPU要跑将近一个月。所以实际项目中要么减少窗口数要么缩短单窗口时间要么用多 GPU 并行。我的建议是先用 10 个窗口、每个 2 纳秒跑一轮快速测试看看 PMF 形状和收敛趋势。如果形状合理再决定要不要加窗口和延长时间。这样可以用最小的代价判断方案是否可行避免一上来就投入大量机时结果发现方向错了。这个流程我前后跑过十几个体系从简单的离子穿膜到复杂的配体-蛋白质解离踩过的坑基本都在这了。最深的体会是伞形采样的成功与否八成取决于前期设计两成取决于后期调参。反应坐标选对了、窗口分布合理了、弹簧常数合适了后面基本就是等结果。反过来如果前期设计有问题后面怎么调都是事倍功半。所以动手之前多花点时间想清楚物理图像比急着跑命令重要得多。
Gromacs伞形采样实战：从反应坐标设计到WHAM自由能计算全流程
1. 为什么伞形采样值得你花时间折腾做分子动力学模拟的人迟早会碰到一个绕不开的需求想知道某个分子穿过某个通道、或者两个分子互相靠近时自由能到底怎么变。普通平衡态模拟里这种事件要么发生得太快根本抓不住要么发生得太慢跑几百纳秒都看不到一次。伞形采样就是专门用来解决这类问题的增强采样方法它通过在反应坐标上施加一系列简谐势也就是“伞”把体系强行约束在坐标轴的不同位置上从而把整条自由能曲线一段一段地“拼”出来。Gromacs 作为分子动力学领域使用最广的开源软件之一对伞形采样的支持相当成熟配套的gmx wham工具可以直接把多窗口的采样结果重加权成 PMFPotential of Mean Force平均力势。但问题在于从体系搭建到最终出图中间有太多容易踩坑的地方pull code 怎么配、窗口怎么分布、采样多久算够、WHAM 报错怎么排查……每一步都可能让你卡上几天。这篇内容适合已经跑过基础分子动力学、想上手伞形采样但被各种报错劝退的人也适合跑过几次但结果总是不太对、想系统梳理一遍流程的人。我会把整个流程拆成设计思路、核心配置、实操步骤、问题排查几个部分尽量把每个参数背后的逻辑讲清楚而不是只丢一堆命令让你照抄。2. 伞形采样的整体设计与思路拆解2.1 伞形采样到底在算什么先把这个方法的物理图像说清楚。假设你想算一个配体从蛋白质口袋往外拉出来的自由能变化反应坐标就是配体与口袋之间的距离。在自然状态下配体大部分时间待在口袋里你几乎采不到它跑到外面的构象。伞形采样的做法是在距离轴上取一系列位置比如 0.3、0.4、0.5……一直到 2.0 纳米每个位置放一个独立的模拟窗口每个窗口里加一根“弹簧”把配体往这个位置拽。弹簧的劲度系数决定了约束的松紧太软了配体会跑偏太硬了采样效率低。每个窗口跑完之后你会得到该窗口内反应坐标的分布以及弹簧施加的力。WHAM 的核心思想就是利用这些有偏采样数据反推出无偏的自由能曲线。它本质上是一个自洽迭代的直方图重加权方法把所有窗口的数据放在一起找到一组权重因子使得各个窗口的分布能够拼成一条连续的曲线。这里有个关键点很多人一开始不理解为什么不能只跑一个窗口把弹簧慢慢移动那样做确实存在叫“牵引模拟”或者“steered MD”但它得到的是非平衡过程的功需要额外做 Jarzynski 平均才能转成自由能收敛慢且方差大。伞形采样是平衡态方法每个窗口单独平衡统计上更稳健这也是它成为主流的原因。2.2 反应坐标的选择逻辑反应坐标选得好不好直接决定整个项目能不能出结果。好的反应坐标应该满足几个条件能区分反应物和产物状态、在物理上合理、采样过程中不会出现其他慢自由度干扰。以配体拉出为例常用的反应坐标有配体与结合位点残基的质心距离、配体沿某条通道的投影距离、或者配体与蛋白质整体质心的距离。质心距离最直观但如果通道是弯曲的直线距离就不能真实反映路径这时候可能需要用路径投影或者多个集体变量组合。我个人的经验是第一次做某个体系时先用牵引模拟快速拉一遍看看配体到底沿什么路径出去、中间有没有卡住的地方。这一步花不了多少机时但能帮你判断反应坐标选得对不对。如果牵引过程中配体撞墙或者绕路那说明你选的坐标有问题直接上伞形采样只会浪费更多时间。2.3 窗口分布与弹簧常数的取舍窗口怎么排、弹簧多硬这两个参数是耦合的。基本原则是相邻窗口的采样分布要有重叠一般建议重叠区域至少占每个窗口分布的 20% 到 30%。如果窗口间距太大分布不重叠WHAM 就无法把相邻窗口连起来PMF 会出现断裂或者跳变。弹簧常数 k 的选择要结合热涨落来考虑。在温度 T 下反应坐标的涨落方差大约是 kT/k。如果你希望涨落标准差在 0.05 纳米左右室温下 kT 约 2.5 kJ/mol那么 k 大约在 1000 kJ/mol/nm² 量级。实际用的时候我一般先用 1000 到 2000 之间的值试跑看直方图分布宽度再决定要不要调整。窗口间距可以按 0.05 到 0.1 纳米来设具体取决于你关心的坐标范围。如果总跨度是 2 纳米用 0.1 纳米间距就是 20 个窗口用 0.05 就是 40 个。窗口越多每个窗口的采样负担越轻但总机时也越多。我的建议是先用粗间距跑一遍看整体形状再在关键区域加密。3. 核心配置与实操要点3.1 拓扑与坐标文件的准备伞形采样对拓扑文件的要求和普通模拟基本一致但有一点要特别注意如果你用的是联合原子力场或者粗粒化力场反应坐标的定义要确保在拓扑里有对应的原子组。比如你要拉一个配体配体的原子索引必须明确不能混在蛋白质组里。坐标文件方面起始结构最好来自平衡过的体系。我见过有人直接拿能量最小化后的结构就开始拉结果配体还没出口袋蛋白质就散了。正确的做法是先跑一段 NPT 平衡让体系密度和盒子尺寸稳定下来再作为伞形采样的起点。另外盒子尺寸要足够大。如果你要拉的距离是 2 纳米盒子在对应方向上至少要有 4 纳米以上的长度否则配体拉到一半就撞到周期性镜像了。这一点在设置盒子时就要考虑进去不要等跑完了才发现。3.2 pull code 的关键参数逐条拆解Gromacs 的 pull code 写在 mdp 文件里核心参数有这么几个pull yes pull-ngroups 2 pull-ncoords 1 pull-group1-name Protein pull-group2-name Ligand pull-coord1-type umbrella pull-coord1-geometry distance pull-coord1-groups 1 2 pull-coord1-dim N N Y pull-coord1-init 0.5 pull-coord1-rate 0 pull-coord1-k 1000pull-coord1-geometry决定了反应坐标怎么算。distance是两组质心距离direction是沿某方向的投影cylinder是柱坐标距离。选哪个取决于你的体系配体拉出一般用distance或direction。pull-coord1-init是窗口的初始位置每个窗口这个值不同。pull-coord1-rate 0表示不做牵引只做静态约束这是伞形采样的标准做法。如果你设了非零速率那就变成牵引模拟了。pull-coord1-k就是弹簧常数。单位是 kJ/mol/nm²。这个值不能太小否则约束不住也不能太大否则采样效率低。前面说的 1000 到 2000 是常用范围。pull-coord1-dim控制哪些方向参与计算。比如你只想沿 z 轴拉就设N N Y。这个参数很容易被忽略但如果设错了反应坐标会包含你不想要的方向分量结果就偏了。3.3 各窗口的初始化与平衡策略每个窗口的起始结构可以从同一个平衡结构生成但初始位置不同。做法是用gmx grompp配合不同的 mdp 文件每个 mdp 里改pull-coord1-init的值。或者更高效的方式是用gmx genrestr之类的工具批量生成但最稳妥的还是脚本化处理。每个窗口开始采样前建议先做一段短平衡让体系适应该窗口的约束。平衡时间不用太长100 到 500 皮秒通常够用但具体要看体系。判断标准是反应坐标的均值稳定在设定值附近且涨落不再有漂移趋势。采样阶段每个窗口一般跑 5 到 20 纳秒。总机时等于窗口数乘以单窗口时间。如果窗口多可以考虑用 GPU 并行跑Gromacs 对多窗口并行支持很好每个窗口一个 GPU 或者多个窗口共享一个 GPU 都可以。注意每个窗口的采样必须保存反应坐标的轨迹也就是 pull 相关的输出。在 mdp 里要确保pull-nstxout和pull-nstfout设置合理一般每 500 到 1000 步输出一次就够了太密会拖慢速度太疏会影响 WHAM 的统计精度。3.4 WHAM 输入文件的生成WHAM 需要每个窗口的 pull 输出文件通常是pullx.xvg或pullf.xvg。gmx wham的命令行大致是这样gmx wham -it tpr_files.dat -if pullf_files.dat -o pmf.xvg -hist histo.xvg其中tpr_files.dat是每行一个 tpr 文件路径pullf_files.dat是每行一个 pullf 文件路径。这两个文件的顺序必须一一对应不能错位。WHAM 还有一些可选参数比如-b和-e指定起止时间-temp指定温度-bins指定直方图分箱数。分箱数一般设 100 到 200太少会平滑掉细节太多会有噪声。如果要做 bootstrap 误差分析可以加-nBootstrap参数WHAM 会自动做重采样并给出误差估计。这个功能很实用建议每次都开。4. 完整实操流程与关键环节4.1 从平衡结构到窗口结构生成假设你已经有一个平衡好的eq.gro和eq.tpr接下来要生成 20 个窗口的起始结构。最直接的方式是写一个 bash 脚本循环修改 mdp 里的pull-coord1-init然后依次跑grompp。for i in $(seq 0 19); do init$(echo 0.3 $i * 0.1 | bc) sed s/pull-coord1-init .*/pull-coord1-init $init/ template.mdp window_$i.mdp gmx grompp -f window_$i.mdp -c eq.gro -p topol.top -o window_$i.tpr -maxwarn 1 done这里-maxwarn 1是为了跳过一些无害的警告比如盒子尺寸和截断半径的关系。但不要滥用这个参数如果警告涉及物理设置错误还是要认真处理。生成 tpr 之后每个窗口先跑一段平衡。平衡的 mdp 和采样 mdp 可以共用只是把nsteps设短一些。平衡完之后用-t参数从平衡的 cpt 文件继续跑采样。4.2 采样阶段的并行与监控采样阶段最耗机时合理利用并行很重要。如果机器上有多个 GPU可以用-multidir或者简单的 shell 并行。比如for i in $(seq 0 19); do (gmx mdrun -deffnm window_$i -s window_$i.tpr -gpu_id $((i % 4)) ) done wait这段脚本把 20 个窗口分配到 4 块 GPU 上并行跑。注意gpu_id的分配要均匀不要让某块 GPU 过载。跑的过程中要监控反应坐标是否稳定。可以用gmx energy或者直接看 pullx 文件。如果某个窗口的坐标均值偏离设定值超过 0.1 纳米说明弹簧太软或者平衡不够需要调整。4.3 WHAM 执行与 PMF 出图采样完成后收集所有窗口的 pullf 文件生成列表ls window_*/pullf.xvg pullf_files.dat ls window_*/window_*.tpr tpr_files.dat然后跑 WHAMgmx wham -it tpr_files.dat -if pullf_files.dat -o pmf.xvg -hist histo.xvg -nBootstrap 100 -bins 200跑完之后pmf.xvg就是自由能曲线。可以用 Grace 或者 Python 的 matplotlib 画图。检查曲线是否平滑、两端是否收敛、误差棒是否合理。如果曲线中间有断裂或者跳变大概率是窗口重叠不够。这时候要么增加窗口要么增大弹簧常数让分布更宽。如果两端翘起可能是采样不够需要延长端部窗口的时间。4.4 结果验证与收敛性检查PMF 跑出来不代表就完事了还要做收敛性检查。常用的方法有几种一是把总采样时间分成前后两半分别做 WHAM看两条曲线是否重合二是看 bootstrap 误差是否在可接受范围三是检查每个窗口的直方图是否有足够的重叠。我一般会画一张直方图重叠图把所有窗口的分布画在一起直观地看相邻窗口是否有交集。如果某两个窗口之间明显断开那一段的 PMF 就不可信。另外反应坐标的选取是否合理也要回头验证。如果 PMF 曲线出现不物理的振荡或者平台可能是坐标里混入了其他自由度。这时候可能需要换坐标或者加限制。5. 常见问题与排查技巧实录5.1 WHAM 报错与直方图问题WHAM 最常见的报错是“no overlap between windows”或者“histogram is empty”。前者说明相邻窗口分布不重叠后者说明某个窗口根本没采到数据。排查思路先看 pullx 文件里反应坐标的范围确认每个窗口是否真的在设定位置附近采样。如果某个窗口的坐标跑到了很远的地方说明弹簧太软或者初始结构没平衡好。如果坐标范围正常但直方图空可能是分箱数设得不对或者输出频率太低导致数据点太少。解决办法增大弹簧常数、增加窗口密度、延长采样时间、调整分箱数。这几个手段可以组合使用。5.2 反应坐标漂移与约束失效有时候跑着跑着反应坐标慢慢偏离设定值最后完全失控。这种情况通常是弹簧太软加上体系里有其他力在干扰。比如配体在口袋里受到蛋白质的吸引力如果弹簧不够硬配体就会被拉回口袋。解决办法是增大 k 值或者在窗口初始化时把配体放在更靠近设定位置的地方。另外检查pull-coord1-dim是否设对如果方向设错了约束力会作用在错误的方向上。还有一种可能是周期性边界条件导致的。如果配体拉到了盒子边缘和镜像发生了相互作用坐标就会跳变。这时候需要加大盒子或者用pull-coord1-geometry direction配合合适的参考向量。5.3 PMF 曲线不收敛的典型表现PMF 不收敛的表现有很多种曲线两端翘起、中间有深谷、误差棒很大、前后半段不一致。每种表现对应的原因不同。两端翘起通常是端部窗口采样不足因为端部窗口的构象空间更大需要更长时间平衡。中间有深谷可能是反应坐标经过了一个真实的能垒但也可能是采样不够导致的假象。误差棒大说明统计量不够需要延长采样或者增加窗口。我自己的经验是先跑短时间看整体形状确定形状合理后再延长采样。如果形状本身就不对延长时间也没用得先改设置。5.4 常见问题速查表问题现象可能原因排查方法解决手段WHAM 报无重叠窗口间距太大或弹簧太软画直方图重叠图加密窗口或增大 k反应坐标漂移弹簧太软或方向设错检查 pullx 轨迹增大 k 或修正 dimPMF 两端翘起端部采样不足分半收敛检查延长端部窗口时间直方图有空窗输出频率太低检查 pullx 数据点提高输出频率PMF 噪声大分箱数太多或采样不足调整 bins 看变化减少 bins 或延长采样约束力异常大初始结构不合理看 pullf 数值重新平衡起始结构5.5 几个容易被忽略的实操细节第一个细节是温度耦合和压力耦合的选择。伞形采样过程中体系被约束在特定构象如果压力耦合太强盒子尺寸会波动可能影响反应坐标的计算。我一般建议用较弱的压力耦合或者只在平衡阶段用采样阶段可以关掉压力耦合改用 NVT。第二个细节是输出频率和文件大小。pullx 和 pullf 文件如果输出太密会占用大量磁盘空间而且 WHAM 读起来也慢。一般每 500 到 1000 步输出一次足够。第三个细节是随机数种子。如果多个窗口用同一个种子可能会引入系统性偏差。Gromacs 默认会自动生成种子但如果你手动设了gen-seed要确保每个窗口不同。第四个细节是 WHAM 的-b和-e参数。如果你在采样初期有不平衡的部分可以用-b跳过。但不要跳太多否则统计量不够。一般跳过前 1 到 2 纳秒比较稳妥。6. 提升效率的进阶技巧6.1 用哈密顿量交换加速采样如果体系比较复杂单个窗口的采样可能很慢。这时候可以考虑用哈密顿量交换Hamiltonian replica exchange结合伞形采样也就是所谓的“伞形采样 副本交换”。Gromacs 本身不直接支持这个组合但可以通过外挂脚本实现或者用 PLUMED 这样的增强采样库。PLUMED 对伞形采样的支持更灵活可以定义更复杂的集体变量也支持多 walker 并行。如果你经常做伞形采样花点时间学 PLUMED 是值得的。6.2 自适应窗口分布传统的伞形采样窗口是均匀分布的但实际自由能曲线往往不均匀有些区域变化快有些区域平坦。自适应方法可以在采样过程中动态调整窗口位置把更多窗口放在自由能变化快的区域。Gromacs 本身没有内置自适应伞形采样但可以通过脚本实现简单的自适应先跑一轮粗采样根据 PMF 的斜率重新分布窗口再跑第二轮。这样可以用更少的窗口达到同样的精度。6.3 结果的可视化与报告PMF 出图之后怎么呈现也很重要。我一般会画三张图第一张是 PMF 曲线加误差棒第二张是直方图重叠图第三张是收敛性检查图。这三张图放在一起审稿人或者合作者一眼就能看出结果可不可靠。画图工具方面Python 的 matplotlib 最灵活可以完全自定义。如果不想写代码Grace 也可以但样式调整比较麻烦。xvg 文件可以直接用xmgrace打开也可以转成 csv 再用其他工具画。提示PMF 的零点可以任意选一般把反应物状态设为零。但报告的时候要说明零点在哪里否则别人看不懂。6.4 机时估算与资源分配最后说一下机时估算。假设你有 20 个窗口每个窗口跑 10 纳秒体系有 5 万个原子用一块中端 GPU 大概每纳秒 2 到 3 小时那么总机时大约是 20 乘以 10 乘以 2.5等于 500 GPU 小时。如果只有一块 GPU要跑将近一个月。所以实际项目中要么减少窗口数要么缩短单窗口时间要么用多 GPU 并行。我的建议是先用 10 个窗口、每个 2 纳秒跑一轮快速测试看看 PMF 形状和收敛趋势。如果形状合理再决定要不要加窗口和延长时间。这样可以用最小的代价判断方案是否可行避免一上来就投入大量机时结果发现方向错了。这个流程我前后跑过十几个体系从简单的离子穿膜到复杂的配体-蛋白质解离踩过的坑基本都在这了。最深的体会是伞形采样的成功与否八成取决于前期设计两成取决于后期调参。反应坐标选对了、窗口分布合理了、弹簧常数合适了后面基本就是等结果。反过来如果前期设计有问题后面怎么调都是事倍功半。所以动手之前多花点时间想清楚物理图像比急着跑命令重要得多。
