第1章：分子动力学基本原理

大家好，我是你们这门课的老朋友。今天咱们来聊聊分子动力学模拟的根基——那些最基础、最核心的原理。说实话，我刚开始接触这个领域时，也被一堆公式和概念搞得头大。但后来我发现，只要把几个关键点吃透了，整个框架就清晰了。

1.1 牛顿运动方程：模拟的发动机

分子动力学模拟，说白了就是让原子动起来。怎么动？靠牛顿第二定律。

公式很简单：F = ma

但这里有个坑——我们面对的是成千上万个原子，每个原子都受到周围原子的力。你得同时解出所有原子的运动轨迹。我当年第一次写模拟程序时，就天真地以为一个循环就能搞定，结果跑出来的结果完全不对。

嗯，这里要注意：每个原子的加速度由它受到的合力决定。而合力，又取决于所有其他原子的位置。所以这是个耦合问题，不能单独处理。

m_i * a_i = F_i = -∇U(r_i)

1.2 势能函数：原子间的"社交规则"

原子之间怎么相互作用？靠势能函数。你可以把它想象成原子间的"社交规则"——距离近了排斥，远了吸引。

最常见的势能函数是Lennard-Jones势：

U(r) = 4ε[(σ/r)^12 - (σ/r)^6]

这个公式里，ε 代表势阱深度（说白了就是结合能），σ 是原子直径。r^12 项描述排斥，r^6 项描述吸引。

我做过一个项目，模拟液态氩的相变。当时用的就是LJ势。结果发现，如果σ设得不对，模拟出来的沸点能差出几十度。所以参数标定真的很重要。

1.3 力场概念：势能函数的"组合拳"

单个势能函数往往不够用。真实材料中，原子间作用很复杂——有键伸缩、键角弯曲、二面角扭转，还有非键相互作用（范德华力、静电作用）。把这些打包在一起，就是力场。

一个典型的力场长这样：

U_total = U_bond + U_angle + U_dihedral + U_nonbond

其中：

• U_bond：键伸缩能，类似弹簧模型
• U_angle：键角弯曲能
• U_dihedral：二面角扭转能
• U_nonbond：非键相互作用（LJ + 库仑）

常用的力场有AMBER、CHARMM、OPLS等。选哪个？取决于你模拟什么。我习惯用OPLS做有机小分子，用CHARMM做蛋白质。但这不是绝对的，得看具体问题。

1.4 时间积分算法：Verlet算法

有了力，有了势能，怎么让原子动起来？我们需要时间积分算法。说白了，就是解微分方程。

最经典的是Verlet算法。它基于泰勒展开：

r(t+Δt) = 2r(t) - r(t-Δt) + a(t) * Δt²

这个公式妙在哪里？它不需要速度就能更新位置。但代价是——速度得另外算。

我常用的变体是速度Verlet算法：

1. r(t+Δt) = r(t) + v(t) * Δt + 0.5 * a(t) * Δt²
2. 计算 a(t+Δt) 基于新位置
3. v(t+Δt) = v(t) + 0.5 * [a(t) + a(t+Δt)] * Δt

为什么推荐这个版本？因为它同时更新位置和速度，而且数值稳定性好。我做过对比测试，同样的模拟，Verlet算法跑100万步，能量漂移只有0.1%。而用欧拉法，10万步就崩了。

1.5 知识体系总览

下面这张图，是我自己总结的分子动力学核心逻辑。你看一遍，应该能对整个框架有个直观认识。

1.6 小结

这一章我们讲了四个核心概念：

• 牛顿方程：F=ma，驱动原子运动
• 势能函数：描述原子间相互作用
• 力场：多种势能函数的组合
• Verlet算法：数值积分，更新原子位置

说实话，这些内容看起来多，但真正上手跑一个模拟后，你会发现它们其实是一体的。我建议你先把Verlet算法的代码写一遍，哪怕只是模拟两个原子。跑通了，后面的大体系自然就理解了。