📚 课程目录第3章能量极小化

能量极小化：最速下降法与共轭梯度法

能量极小化，说白了就是帮你的分子结构找个「坑」蹲着。这个坑不是随便哪个坑，是势能面上最舒服的那个局部最低点。我刚开始做分子模拟时，总觉得优化就是跑几步的事，结果结构乱七八糟，能量满天飞。后来才明白，选不对极小化方法，后面的一切采样都是白搭。

这章我们就聊聊最常用的两种优化算法：最速下降法和共轭梯度法。它们各有脾气，用对了事半功倍，用错了…嗯，你会花一整天在服务器前发呆。

最速下降法：简单但容易「打转」

最速下降法的思路非常直观：沿着梯度方向往下走，能量下降最快的方向就是它选的路。每走一步，重新算一次梯度，再沿着新方向继续走。

你可以想象一个人蒙着眼站在山坡上，每一步都朝感觉最陡的方向迈。好处是简单、稳定、占用内存少。坏处呢？越靠近谷底，步子就越容易「打转」，左右晃动就是下不去。

我的经验：最速下降法最适合处理初始结构离极小点很远的情况。比如刚搭建的蛋白-配体复合物，原子之间可能重叠严重，这时候用别的算法容易直接炸掉。

最速下降法的核心更新公式很简单：

x_{k+1} = x_k - α_k * ∇f(x_k)

其中 α_k 是步长，∇f(x_k) 是当前梯度。步长的选择很关键。太大，结构会跳过极小点；太小，收敛慢得让你怀疑人生。我在项目中习惯先令步长在 0.01-0.1 nm 之间，观察能量下降曲线再调整。

一个常见的配置片段长这样：

; 最速下降法参数
nsteps          = 5000
nstcgmin        = 100   ; 共轭梯度前先跑100步最速下降
emstep          = 0.01  ; 初始步长 (nm)
emtol           = 10.0  ; 能量收敛判据 (kJ/mol)
nstxout         = 100   ; 每100步输出坐标

避坑：我曾经接过一个膜蛋白项目，初始结构离平衡态太远，直接上共轭梯度法，结果跑了没几步力场就炸了。后来我学乖了，不管三七二十一，先用 500 步最速下降法「降温」，再切到共轭梯度法。这个习惯保留到现在。

共轭梯度法：走直线，不回头

共轭梯度法比最速下降法聪明一点。它会记住上一次的搜索方向，在保证当前方向与之前所有方向「共轭」的前提下选择下一步。说白了，它不会像最速下降法那样在每个拐弯处重新选择方向，而是试图走一条更直接的路径。

共轭梯度的更新公式：

d_{k+1} = -∇f(x_{k+1}) + β_k * d_k
x_{k+1} = x_k + α_k * d_k

β_k 的计算方式有好几种，最常见的包括 Fletcher-Reeves 和 Polak-Ribiere。我个人的习惯是：用 Polak-Ribiere 公式，它在处理非二次型势能面时表现得更好。

特性	最速下降法	共轭梯度法
收敛速度	慢，线性收敛	快，超线性收敛
内存占用	低	中等
初始结构适应性	好，稳定	差，容易发散
局部精细收敛	吃力	高效
适用阶段	预优化、消除原子重叠	精调、逼近极小点

那这两种方法怎么搭配使？我个人建议的顺序是：先用最速下降法跑几百步「开荒」，等能量下降趋势平缓了，再换共轭梯度法收割。你想想看，这不就跟开荒种地一样吗？先刨地，再细作。

  调试要点：

  观察每一步的能量是否单调下降。如果不降反升，说明步长太大或结构发散。
最速下降法阶段，能量下降曲线大致呈指数衰减。如果衰减太慢，可以稍微提高步长。
共轭梯度法阶段，能量应该快速收敛。如果出现震荡，检查是否在切换算法时梯度没清零。

收敛判据怎么选？别被表面数字骗了

收敛判据是决定你何时停止优化的「尺子」。选得太松，结构没完全稳定；选得太严，计算资源浪费不说，可能永远收敛不到。

最常用的两个判据：

能量判据：设定一个能量变化阈值，比如连续几步能量变化小于 tol。优点是直观，缺点是你可能只盯着能量看。
力判据：计算所有原子受力的最大值或均方根，当力小于某个值（如 10 kJ/mol/nm）时认为收敛。这是更物理的判据，因为力直接反映了结构离平衡点的距离。

我之前在一个项目中只用能量判据，优化后看着能量挺稳的，一跑动力学发现结构立刻崩了。后来检查才发现，虽然能量没怎么变，但某些原子的受力还很大，只是能量对局部应变不敏感。从那以后，我坚持能量判据和力判据一起用，而且是后者优先。

我曾经踩过的坑：有一回优化一个环肽，能量收敛到小数点后两位，力判据也过了，结果跑 NVT 没到 10 ps 就二面角翻转了。为什么？因为环状结构的应力分布非常复杂，全局能量收敛不代表每个键角都合理。后来我加了一个分步优化：先放开所有键角约束粗优化，再收紧约束精细调整。之后再跑动力学就稳多了。

常见的收敛参数示例：

; GROMACS .mdp 中的收敛判据设置
emtol           = 100.0   ; 最大力 (kJ/mol/nm)
emstep          = 0.01    ; 初始步长 (nm)
nsteps          = 10000   ; 最大步数
nstcgmin        = 500     ; 最速下降最小步数

我一般会把 emtol 设成 100、50、10 分三档，先粗跑再细调。一来节约时间，二来可以分阶段检查结构是否合理。如果你要跑自由能计算或增强采样，收敛判据得严一些，力判据建议下探到 5 kJ/mol/nm 以下。

算法选择的整体思路

说了这么多，我画一张简图帮大家理清整个极小化流程。从初始结构出发，根据结构质量选择算法，逐步逼近收敛。其中每一步的判据都值得你停下来想一想：这个力正常吗？能量下降的速度合理吗？

好了，能量极小化的核心内容就这些。最速下降法和共轭梯度法，一个像大砍刀，一个像手术刀。知道什么时候该用哪一把，比会背公式重要得多。调试的时候多留个心眼，看看能量曲线，查查受力分布，别只盯着收敛判据的数值过不过。