分子动力学的计算：系综选择、时间步长与恒温器对比

分子动力学的计算中，项目组遇到过一个诡异现象：用Nosé-Hoover恒温器跑蛋白质在水中的NVT模拟，温度曲线在300K附近振荡幅度±25K——看起来在”恒温”，但蛋白的RMSD波动异常大，二级结构在模拟过程中部分解折叠。换成Langevin恒温器后，温度振荡缩到±3K，蛋白结构稳定。同一个体系、同样的力场参数，恒温器不同导致模拟结论完全不同。这之后项目组建立了一套严格的参数选择流程。

系综选择：NVT/NPT/NVE的物理边界

系综选择由模拟目标决定，不是个人偏好：

NVE（微正则系综）：粒子数N、体积V、能量E恒定。这是最”纯”的MD——没有外部热浴干预，系统完全遵循牛顿方程。适用于：

• 验证力场能量守恒（NVE中总能量漂移应<0.01%）
• 冲击波模拟（非平衡过程，不能加恒温器）
• 自由能计算中的伞形采样（恒温器会干扰势能面采样）

NVT（正则系综）：粒子数N、体积V、温度T恒定。适用于：

• 固定晶胞的结构性质计算（RDF、MSD）
• 表面吸附模拟（体积不变）
• 温度梯度的初始平衡阶段

NPT（等温等压系综）：粒子数N、压力P、温度T恒定。适用于：

• 液体密度预测（NVT需要预先知道密度）
• 相变模拟（体积随相变变化）
• 热膨胀系数计算

项目组的标准流程是：NVT预平衡（10万步）→ NPT平衡（20万步）→ NVE或NVT采样（100万步以上）。先用NVT让结构稳定，再用NPT调整到正确的密度，最后在目标系综下做生产模拟。

时间步长：不是越小越好

时间步长的选择取决于体系中最高频振动模式。含氢体系（水、有机分子、蛋白质）的O-H键振动频率约3400-3700 cm⁻¹，对应周期约9-10 fs；C-H键振动频率约2900-3100 cm⁻¹，对应周期约11-12 fs。时间步长必须小于最高频周期的1/10——否则积分误差会导致能量漂移。

项目组测试了不同步长下水分子体系（SPC/E水模型，4096个水分子）的能量守恒：

0.5 fs是含氢体系的黄金标准——能量漂移可忽略。1.0 fs在短时模拟中可用，但超过1 ns后能量漂移开始影响统计。2.0 fs在NVE中已不可靠，但如果使用SHAKE/LINCS约束算法固定C-H和O-H键，步长可以安全地提升到2.0 fs——因为约束后最高频振动变成了C-C键扭转（周期~40 fs）。

金属体系不含氢，最高频是声子振动（周期~100 fs），时间步长可以用1-2 fs。项目组在LAMMPS中跑铜熔体用2 fs步长，NVE中1 ns能量漂移仅0.008%。

恒温器对比：三种主流方法的优劣

Nosé-Hoover (NH)：

NH通过引入一个额外的”热浴”变量来控制温度。它的优点是确定性（非随机），可以正确生成正则系综的分布。但NH有一个致命弱点——当系统远离平衡态时（如初始结构很差、力场参数不匹配），NH的反馈机制可能产生温度的周期性振荡（”hot-cold”循环），即所谓的”flying ice cube”问题。

项目组在蛋白质体系中的经验：NH适合平衡良好的体系做长时间采样。恒温器耦合常数（Q参数，LAMMPS中的Tdamp）需要谨慎选择——太短（<10 fs）会导致温度振荡，太长（>500 fs）温度控制不灵敏。推荐值：Tdamp=100 fs。

Langevin：

Langevin恒温器在牛顿方程上加了随机力（模拟热浴的碰撞）和阻尼力（模拟耗散）：

m·a = F - mγv + R(t)

其中γ是阻尼系数，R(t)是随机力（满足涨落-耗散定理）。

Langevin的优势是鲁棒性——即使初始结构很差，也不会产生温度振荡。缺点是随机力破坏了动量守恒，在某些场景下会影响动力学性质（如扩散系数的计算）。项目组在蛋白质模拟中偏好Langevin，因为生物分子力场（AMBER、CHARMM）的初始结构通常远离平衡态，Langevin能平稳过渡。

阻尼系数γ的选择：γ太大会过阻尼（运动被”冻结”），γ太小会欠阻尼（温度控制弱）。项目组用γ=1 ps⁻¹（对应LAMMPS中damping=100），这是生物分子模拟的标准值。

CSVR (Bussi-Donadio-Parrinello)：

CSVR（Canonical Sampling through Velocity Rescaling）是一种基于速度重标定的随机恒温器。它结合了NH的确定性框架和Langevin的随机性，但通过精巧的设计确保了正确的正则分布。

CSVR的优势：温度控制精确（振荡最小），动量守恒比Langevin好，且不产生NH的周期振荡。CSVR在动力学性质（如扩散系数）的计算上偏差通常小于Langevin。

三种恒温器的总结：

反思：参数选择是一个整体

系综、步长、恒温器不是独立选择的——它们相互影响。NPT模拟中，恒压器和恒温器的耦合常数需要匹配（恒压器通常比恒温器慢3-5倍），否则会出现压力-温度耦合振荡。使用SHAKE约束时，步长可以加大，但LINCS的阶数（order参数）也要相应提高（从4阶升到8阶），否则约束精度不够。

分子动力学的计算不是”设好参数跑就完事”，每一步选择都需要物理依据。更多分子动力学的实战经验，可以参考分子动力学栏目，或返回科研学术网首页。