GROMACS SMD拉伸动力学模拟指南：弹簧常数、拉速与力谱分析

蛋白-配体复合物的结合自由能可以用FEP算，但结合/解离路径(它们是怎么分开的?要克服多大的力垒?)需要另一种工具。拉伸分子动力学(Steered Molecular Dynamics，SMD)就是为此设计的——它给系统施加一个人为的方向约束和拉力，观察它如何响应。

SMD在GROMACS中通过pull代码实现。它的应用场景远比想象中广泛：蛋白-配体解离路径、蛋白去折叠力学响应、高分子链拉伸弹性模量、以及配合Jarzynski等式反推平衡自由能(PMF)。

恒速拉伸 vs 恒力拉伸：两种模式的选择

GROMACS的pull代码支持两种基本模式：

恒速模式（pull=umbrella）：假想的弹簧一端固定在参考组，另一端连在拉动组，弹簧的平衡位置以恒定速度移动。体系感受到一个不断增大的回复力，直到克服能垒完成跃迁——跃迁时刻对应力-位移曲线的峰值。这是最常用的SMD模式，因为它和单分子力谱实验(AFM、光镊)的逻辑最接近。

恒力模式（pull=constant-force）：对拉动组施加恒定大小的外力。体系在恒定偏压下演化——如果外力足够大，能垒被压低到kT量级以下，体系自发跃迁。恒力模式适合研究临界解离力的阈值效应，但取样效率低于恒速，目前在单分子生物物理中用得较少。

实际项目中，恒速SMD占据绝对主流。

弹簧常数与拉速：SMD的两个核心参数

SMD结果对弹簧常数(pull-coord1-k)和拉速(pull-coord1-rate)极其敏感，选错任何一个，你的力谱看的是弹簧响应而非分子事件的信号。

弹簧常数（k）： 太软(<100 kJ/mol·nm²)，弹簧的形变吸收了大部分位移信号，蛋白-配体界面的真实解离过程被掩埋在弹簧噪声里。太硬(>5000 kJ/mol·nm²)，体系感受到的是急剧增加的应力，可能在非物理路径上”撕开”界面。对于蛋白-配体解离，1000-2000 kJ/mol·nm²是经验范围，上下限对应需要更高灵敏度还是更高稳定性。

拉速（v）： 影响的是过程的可逆程度和力谱的噪声水平。低拉速(<0.001 nm/ps = 1 m/s)接近准静态，但要在微秒量级模拟时间内观察到解离事件，拉速不能低于这个值的10倍以上。高拉速(>0.01 nm/ps = 10 m/s)会在力谱中引入显著的粘性阻尼贡献——”力峰值”中有一部分是溶剂摩擦而非真实能垒。

典型的折中拉速是0.005-0.01 nm/ps(5-10 m/s)。注意这是MD的拉速，比单分子实验(μm/s量级)快6-9个数量级——这是时间尺度限制的结果，也是SMD结果不能直接与实验力值一一对应的原因。

pull代码的具体配置

以蛋白-配体解离为例，沿Z轴向上拉配体：

; pull code

pull = yes

pull-ngroups = 2

pull-ncoords = 1

; 组定义

pull-group1-name = protein_ref ; 蛋白骨架(固定参考)

pull-group2-name = ligand ; 配体质心(被拉)

; 拉动坐标

pull-coord1-type = umbrella ; 恒速

pull-coord1-geometry = direction ; 沿向量方向

pull-coord1-vec = 0.0 0.0 1.0 ; Z轴方向

pull-coord1-groups = 1 2 ; 组1和组2的距离

pull-coord1-dim = N N Y ; 只在Z方向有约束

pull-coord1-k = 1000 ; kJ/(mol·nm²)

pull-coord1-rate = 0.01 ; nm/ps

pull-coord1-start = yes ; 从初始距离开始

pull-coord1-dim = N N Y是关键设置——只在Z方向施加力约束。如果写成了Y Y Y(三个方向)，相当于全向约束，蛋白-配体在XY方向的自由相对运动也被限制，SMD变成拖拽而非解离。

初始距离(pull-coord1-start=yes对应的距离)应该略大于复合物中蛋白-配体质心的实际距离。常用蛋白-配体距离在2.0-3.0 nm左右，从这个初始值拉开到5.0-6.0 nm完成解离——全程约200-400 nm的弹簧平衡位置移动，对应0.01 nm/ps拉速下20-40 ns模拟。

力-位移曲线的分析

SMD最常见的输出是力-时间曲线，从拉速换算到位移即得力-位移曲线。典型的蛋白-配体解离SMD力谱呈四个区段：

弹性区（小幅增长）： 配体和蛋白骨架之间的非键相互作用在弹性形变范围内响应拉力，力近似线性增长。这是弹簧力和界面的范德华/静电作用达到准平衡的阶段。

力垒区（快速增长→峰值）： 配体开始从结合口袋中被拉出，疏水作用和氢键逐个断裂，力达到峰值。峰值力的大小取决于拉速——拉速越高，水分子来不及重新排布到暴露的界面，疏水粘附贡献更大，峰值力被高估。

力下降区（指数衰减）： 配体离开结合口袋进入溶剂，非键相互作用迅速衰减。力减小的速率反映口袋几何形状的”陡峭度”——深而窄的口袋延缓下降，浅而宽的口袋下降快。

完全解离区（基线附近波动）： 配体进入体相溶液，与蛋白只在偶然碰撞中产生微弱的随机力信号，力在零附近波动。

多次独立SMD轨迹(至少5-10条，从不同初始速度分布开始)的统计平均可以得到峰值力的均值和方差。单条SMD轨迹的力峰值有20-30%的随机性——不要基于一次计算就下结论。

PMF反推：从非平衡力谱到平衡自由能

Jarzynski等式(JE)是连接非平衡SMD力谱和平衡PMF的桥梁。它的核心思想是：对无数次非平衡做功求指数平均，等于平衡态的自由能差。在SMD语境中：多条恒定拉速的力-位移曲线做指数加权平均，可反推出沿反应坐标的PMF。

实际操作中的困难：JE的统计收敛需要大量(>100条)SMD轨迹，因为指数平均在小样本下被稀有的大波动事件主导。第二累积展开(2nd cumulant expansion，即功的均值和方差近似自由能)对接近平衡的慢拉速更稳健，但需要低速下额外的轨迹数。

更实用的PMF计算方法还是Umbrella Sampling(US)——沿着反应坐标设多个窗口，每个窗口独立约束采样，用WHAM整合成PMF。US在SMD的力谱标定基础上有序布置窗口位置(窗口间距可以参考力垒位置)，效率远高于从头网格搜索。

SMD虽然看似简单，但弹簧常数和拉速的选择决定了结果的可信度。科研学术网(https://www.keyanxueshu.com）-配体解离、蛋白去折叠和高分子拉伸等场景。