简介这份资源是经典教材《The Art of Molecular Dynamics Simulation》D.C. Rapaport 著Cambridge University Press 出版随书附带的 C 语言程序代码面向物理、化学、材料与生物等方向的学习者和科研人员帮助读者把分子动力学理论落到可编译运行的实践层面。压缩包共 127 个文件约 264KB其中 66 个 .c 源文件构成核心算法与示例程序52 个 .in 输入文件用于配置模拟参数另有 4 个 .h 头文件、2 个 .sh 脚本及 readme、copying、errata 等说明文档结构清晰、便于按章节检索。已有 50 人学习。代码覆盖分子模型构建、能量最小化、系统时间演化等关键环节并涉及并行计算与向量化等加速思路读者可据此理解科学计算程序的框架设计与优化技巧也可将其作为教学辅助或科研工作的起点进行算法改进与特定问题的模拟验证。1. 分子动力学模拟的 C 语言实现从理论公式到可运行代码很多人第一次接触分子动力学MD模拟是从一本经典教材开始的但真正动手把书里的公式变成能跑的 C 代码时往往会卡在几个地方Verlet 列表怎么建、周期性边界怎么处理、能量为什么一直漂移。这个标题指向的正是这条路径——用 C 语言从零实现一套分子动力学模拟的核心框架。它适合两类人一是学过统计物理但没写过 MD 代码的研究生二是需要把 MD 嵌入自己仿真流程的工程师。C 语言在这里不是随便选的MD 每步要算几十万次粒子间作用力C 的指针操作和内存布局控制能直接把性能拉满而 Python 或 MATLAB 在这个场景下往往差一个数量级。下面从最小可运行系统开始把每一步的参数和坑讲清楚。2. 最小可运行 MD 系统Lennard-Jones 流体与 Velocity-Verlet 积分2.1 为什么选 Lennard-Jones 势作为起点分子动力学模拟的核心是求解牛顿方程。对于最简单的流体系统Lennard-JonesLJ势是最常用的模型势[ U(r) 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right] ]其中 (\epsilon) 是势阱深度(\sigma) 是零势距离。选它作为起点有三个理由第一公式简单求导后力表达式干净第二氩气等简单流体的性质已经被研究得很透方便验证结果第三截断处理有成熟方案不会一上来就陷入长程力的复杂性。常见做法是把 (\epsilon) 和 (\sigma) 作为约化单位即令 (\epsilon1)、(\sigma1)质量 (m1)。这样所有物理量都变成无量纲数代码里不需要反复乘常数。实际模拟氩气时再按 (\epsilon0.0104) eV、(\sigma3.40) Å 换算回真实单位。2.2 初始化位置、速度与周期性边界初始化分三步放置粒子、赋初速度、处理边界。放置粒子最忌讳的是让两个粒子靠得太近因为 LJ 势在 (r \to 0) 时趋于无穷直接导致力爆炸。我一般用面心立方FCC晶格排布再给每个粒子加一点随机扰动这样既保证初始无重叠又不会让系统一直保持晶格结构。// 初始化FCC 晶格 随机扰动 #include math.h #include stdlib.h #define N 256 // 粒子数取 4 的倍数便于 FCC 排布 #define RHO 0.8 // 约化密度 #define BOX pow(N / RHO, 1.0/3.0) // 立方盒子边长 typedef struct { double x, y, z; // 位置 double vx, vy, vz; // 速度 double fx, fy, fz; // 受力 } Particle; void init_positions(Particle *p) { int n (int)cbrt(N / 4.0); // 每边晶胞数 double a BOX / n; // 晶胞边长 int idx 0; // FCC 基元4 个原子 per 晶胞 double basis[4][3] { {0, 0, 0}, {0.5, 0.5, 0}, {0.5, 0, 0.5}, {0, 0.5, 0.5} }; for (int i 0; i n; i) for (int j 0; j n; j) for (int k 0; k n; k) for (int b 0; b 4; b) { if (idx N) break; p[idx].x (i basis[b][0]) * a 0.01 * (rand()/(double)RAND_MAX - 0.5); p[idx].y (j basis[b][1]) * a 0.01 * (rand()/(double)RAND_MAX - 0.5); p[idx].z (k basis[b][2]) * a 0.01 * (rand()/(double)RAND_MAX - 0.5); idx; } }这段代码里RHO0.8是约化密度对应液态氩的典型值。扰动幅度取 0.01 倍晶胞边长太大可能造成近邻重叠太小则系统长时间保持晶格特征。BOX用pow计算保证密度精确。初速度一般从麦克斯韦-玻尔兹曼分布采样然后减去质心速度防止整个系统漂移。温度通过速度平方和定义[ T \frac{1}{3N} \sum_i m v_i^2 ]在约化单位下 (m1)所以直接累加速度平方即可。周期性边界是 MD 的命门。粒子从盒子右边出去要从左边回来。位置更新后立刻做环绕// 周期性边界处理 void apply_pbc(Particle *p) { for (int i 0; i N; i) { p[i].x - BOX * floor(p[i].x / BOX); p[i].y - BOX * floor(p[i].y / BOX); p[i].z - BOX * floor(p[i].z / BOX); } }用floor而不是if判断是因为粒子一步可能跨过多个盒子长度虽然正常模拟不会但温度失控时会。这个写法比取模更稳负数也能正确处理。2.3 Velocity-Verlet 积分一步更新位置、速度、受力积分算法选 Velocity-Verlet因为它时间可逆、能量守恒好而且每步只算一次力。公式如下[ r(t\Delta t) r(t) v(t)\Delta t \frac{1}{2}a(t)\Delta t^2 ] [ v(t\Delta t) v(t) \frac{1}{2}[a(t) a(t\Delta t)]\Delta t ]实现时先更新位置再算新受力最后更新速度// Velocity-Verlet 积分一步 void integrate(Particle *p, double dt) { // 1. 更新位置 for (int i 0; i N; i) { p[i].x p[i].vx * dt 0.5 * p[i].fx * dt * dt; p[i].y p[i].vy * dt 0.5 * p[i].fy * dt * dt; p[i].z p[i].vz * dt 0.5 * p[i].fz * dt * dt; } apply_pbc(p); // 2. 保存旧力 double old_fx[N], old_fy[N], old_fz[N]; for (int i 0; i N; i) { old_fx[i] p[i].fx; old_fy[i] p[i].fy; old_fz[i] p[i].fz; } // 3. 计算新力 compute_forces(p); // 4. 更新速度 for (int i 0; i N; i) { p[i].vx 0.5 * (old_fx[i] p[i].fx) * dt; p[i].vy 0.5 * (old_fy[i] p[i].fy) * dt; p[i].vz 0.5 * (old_fz[i] p[i].fz) * dt; } }dt的选取很关键。LJ 系统典型振动周期在 (10^{-2}) 约化时间量级dt一般取 0.001 到 0.005。取 0.005 能跑得快但能量漂移明显取 0.001 稳但同样步数覆盖的物理时间短。我一般先用 0.001 验证能量守恒确认无误后再放宽到 0.002 做生产运行。2.4 力计算与截断别让长尾力拖垮性能LJ 力公式[ \mathbf{F}{ij} 24\epsilon \left[ 2\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right] \frac{\mathbf{r}{ij}}{r^2} ]直接对每对粒子算力是 (O(N^2))256 个粒子还能忍上千个就必须用 Verlet 列表或格点法。但即使不优化截断也必须做当 (r r_c) 时力已经很小忽略它。常见截断半径 (r_c 2.5\sigma)。// 计算所有粒子受力含截断 #define RCUT 2.5 #define RCUT2 (RCUT * RCUT) void compute_forces(Particle *p) { // 清零 for (int i 0; i N; i) p[i].fx p[i].fy p[i].fz 0.0; // 双重循环 for (int i 0; i N - 1; i) { for (int j i 1; j N; j) { double dx p[i].x - p[j].x; double dy p[i].y - p[j].y; double dz p[i].z - p[j].z; // 最小镜像约定 dx - BOX * round(dx / BOX); dy - BOX * round(dy / BOX); dz - BOX * round(dz / BOX); double r2 dx*dx dy*dy dz*dz; if (r2 RCUT2 || r2 1e-12) continue; double inv_r2 1.0 / r2; double inv_r6 inv_r2 * inv_r2 * inv_r2; double inv_r12 inv_r6 * inv_r6; double f 24.0 * (2.0 * inv_r12 - inv_r6) * inv_r2; p[i].fx f * dx; p[i].fy f * dy; p[i].fz f * dz; p[j].fx - f * dx; p[j].fy - f * dy; p[j].fz - f * dz; } } }这里用round实现最小镜像约定比if判断更简洁。r2 1e-12是防止两个粒子完全重合导致除零。截断半径 2.5 是 LJ 系统的经验值此时势能约为 (-0.016\epsilon)力约为 (0.039\epsilon/\sigma)对大多数热力学量影响在 1% 以内。如果做相变或界面研究需要加长程修正或换更复杂的势。3. 让模拟跑得稳温度控制、能量诊断与性能优化3.1 温度耦合Berendsen 热浴的参数怎么设刚初始化的系统温度往往偏离目标值需要热浴调节。Berendsen 热浴是最简单的每步把速度乘以一个因子 (\lambda)让温度向目标 (T_0) 弛豫。[ \lambda \sqrt{1 \frac{\Delta t}{\tau_T}\left(\frac{T_0}{T} - 1\right)} ]// Berendsen 热浴 void berendsen_thermostat(Particle *p, double dt, double T0, double tau_T) { double sum_v2 0.0; for (int i 0; i N; i) sum_v2 p[i].vx*p[i].vx p[i].vy*p[i].vy p[i].vz*p[i].vz; double T sum_v2 / (3.0 * N); // 约化单位下 m1, kB1 double lambda sqrt(1.0 dt / tau_T * (T0 / T - 1.0)); for (int i 0; i N; i) { p[i].vx * lambda; p[i].vy * lambda; p[i].vz * lambda; } }tau_T是耦合时间常数一般取 0.1 到 1.0 约化时间。取太小如 0.01会让速度分布被强行拉成窄峰动力学性质失真取太大如 10则温度弛豫慢浪费机时。我一般用 0.5平衡阶段用 0.1 加速收敛生产阶段换回 0.5。注意 Berendsen 热浴不产生正确的正则系综它只是让温度“看起来对”。如果要做涨落相关的计算得换 Nosé-Hoover 或 Andersen 热浴。但作为入门和快速平衡Berendsen 足够。3.2 能量诊断总能量漂移多少算正常每步输出总能量动能势能看它是否守恒。Velocity-Verlet 在dt0.001时能量漂移应该在 (10^{-4}) 每粒子每步量级。如果漂移达到 (10^{-2})说明dt太大或力计算出错。// 计算总能量 double total_energy(Particle *p) { double ke 0.0, pe 0.0; for (int i 0; i N; i) ke 0.5 * (p[i].vx*p[i].vx p[i].vy*p[i].vy p[i].vz*p[i].vz); for (int i 0; i N - 1; i) { for (int j i 1; j N; j) { double dx p[i].x - p[j].x; double dy p[i].y - p[j].y; double dz p[i].z - p[j].z; dx - BOX * round(dx / BOX); dy - BOX * round(dy / BOX); dz - BOX * round(dz / BOX); double r2 dx*dx dy*dy dz*dz; if (r2 RCUT2 || r2 1e-12) continue; double inv_r6 1.0 / (r2*r2*r2); pe 4.0 * (inv_r6 * inv_r6 - inv_r6); } } return ke pe; }势能求和时注意截断否则会重复计算。输出频率一般每 100 步一次太频繁拖慢速度太稀疏看不出漂移趋势。3.3 Verlet 列表与格点法什么时候必须上(O(N^2)) 在 (N256) 时每步约 3 万次对算还能接受。但 (N1000) 时是 50 万次(N10000) 时是 5000 万次必须优化。Verlet 列表的思路是对每个粒子记录截断半径加缓冲区内所有邻居每 10 到 20 步重建一次列表中间步只遍历列表。// Verlet 列表结构简化版 typedef struct { int count; int *neighbors; // 动态数组 } VerletList; // 重建列表缓冲半径 RCUT SKIN #define SKIN 0.3 void build_verlet_list(Particle *p, VerletList *vl) { double r_skin2 (RCUT SKIN) * (RCUT SKIN); for (int i 0; i N; i) { vl[i].count 0; for (int j 0; j N; j) { if (i j) continue; double dx p[i].x - p[j].x; double dy p[i].y - p[j].y; double dz p[i].z - p[j].z; dx - BOX * round(dx / BOX); dy - BOX * round(dy / BOX); dz - BOX * round(dz / BOX); if (dx*dx dy*dy dz*dz r_skin2) vl[i].neighbors[vl[i].count] j; } } }SKIN取 0.3 是经验值太小则列表频繁重建太大则列表变长、遍历变慢。重建频率取决于粒子最大位移一般每 10 步重建一次。格点法Cell List更适合超大系统把盒子切成边长大于截断半径的格子只检查相邻格子复杂度降到 (O(N))。但格点法实现复杂容易在边界处理上翻车建议先用 Verlet 列表把 (N) 跑到 5000 左右再考虑格点。4. 避坑与排查MD 代码最容易翻车的五个地方4.1 能量突然爆炸温度冲到 (10^6)现象跑了几百步后总能量指数上升粒子飞出盒子。原因两个粒子距离过近LJ 势的 (r^{-12}) 项产生巨大斥力一步就把速度打到天上。解决初始化时检查最小粒子间距确保大于 (0.8\sigma)积分步长dt降到 0.0005 先跑 1000 步加一个力上限保护当 (r 0.5\sigma) 时直接跳过或设力为固定值。4.2 温度一直往下掉系统“冻结”现象Berendsen 热浴开着但温度缓慢下降。原因截断半径处势能不为零能量不守恒Velocity-Verlet 会慢慢损失能量。解决加截断位移修正或者换用带长程修正的势。更简单的办法是每步都调热浴但这样动力学性质不准。我一般接受小幅漂移只要在生产运行的时间尺度内温度变化不超过 5%。4.3 最小镜像约定写错力方向反了现象两个粒子明明靠近却互相吸引而不是排斥。原因round(dx/BOX)用错或者dx没有先做减法。解决记住最小镜像的公式是 (dx dx - L \cdot \text{round}(dx/L))其中 (L) 是盒子边长。测试方法放两个粒子在盒子两边距离应该是 (L - r) 而不是 (r)。4.4 随机数没播种每次运行结果一样现象想跑不同初始条件但输出完全一致。原因rand()默认种子是 1。解决在main开头加srand(time(NULL))或者从命令行读入种子。但注意做对照实验时反而要固定种子否则无法复现。4.5 输出文件太大硬盘写满现象跑了 10 万步轨迹文件几十 GB。原因每步都写所有粒子坐标。解决每 100 或 1000 步写一次只写需要分析的粒子或只写能量。如果要做可视化可以写 XYZ 格式但控制帧数。我一般每 500 步写一帧10 万步也就 200 帧足够看结构演化。5. 从能跑到好用径向分布函数与扩散系数的计算技巧5.1 径向分布函数 (g(r))分箱宽度与归一化径向分布函数是 MD 最常输出的结构量。它描述以某个粒子为中心距离 (r) 处找到另一个粒子的概率密度相对于均匀分布的比例。计算时先做直方图再归一化[ g(r) \frac{V}{N^2} \frac{1}{4\pi r^2 \Delta r} \sum_{i} \sum_{j \neq i} \delta(r - r_{ij}) ]// 计算径向分布函数 #define NBINS 200 #define RMAX 5.0 void compute_rdf(Particle *p, double *g_r, int *count) { double dr RMAX / NBINS; for (int b 0; b NBINS; b) { g_r[b] 0.0; count[b] 0; } for (int i 0; i N - 1; i) { for (int j i 1; j N; j) { double dx p[i].x - p[j].x; double dy p[i].y - p[j].y; double dz p[i].z - p[j].z; dx - BOX * round(dx / BOX); dy - BOX * round(dy / BOX); dz - BOX * round(dz / BOX); double r sqrt(dx*dx dy*dy dz*dz); if (r RMAX) { int bin (int)(r / dr); count[bin] 2; // i-j 和 j-i 都算 } } } for (int b 0; b NBINS; b) { double r (b 0.5) * dr; double shell_vol 4.0 * M_PI * r * r * dr; double rho N / (BOX * BOX * BOX); g_r[b] count[b] / (N * rho * shell_vol); } }分箱宽度dr取 0.025即 RMAX/NBINS太宽则峰形模糊太窄则统计噪声大。归一化时注意除以 (N \times \rho \times 4\pi r^2 dr)其中 (\rho) 是数密度。跑 1000 步后取平均(g(r)) 的第一个峰应该在 (r \approx 1.1\sigma)高度约 2.8这是 LJ 液体的特征。5.2 扩散系数从均方位移斜率提取扩散系数 (D) 通过均方位移MSD的长时间斜率得到[ D \lim_{t \to \infty} \frac{1}{6t} \langle |\mathbf{r}(t) - \mathbf{r}(0)|^2 \rangle ]实现时记录每个粒子相对于初始位置的位移平方和对时间做平均。注意要用未折叠坐标否则粒子穿过边界时位移会跳变。// 均方位移计算需保存初始位置和未折叠坐标 void compute_msd(Particle *p, Particle *p0, double *msd, int step) { double sum 0.0; for (int i 0; i N; i) { double dx p[i].x - p0[i].x; double dy p[i].y - p0[i].y; double dz p[i].z - p0[i].z; sum dx*dx dy*dy dz*dz; } msd[step] sum / N; }实际代码里要维护未折叠坐标每次 PBC 环绕时记录粒子移动了多少个盒子长度累加到未折叠坐标上。MSD 曲线先经历弹道区(t^2)再进入扩散区线性。取线性段做最小二乘拟合斜率除以 6 就是 (D)。LJ 液体在 (T0.8)、(\rho0.8) 时(D) 约 0.05 到 0.1 约化单位。5.3 一个容易忽略的细节相关时间与统计误差MSD 和 (g(r)) 都需要对多个时间起点做平均否则统计误差很大。我一般把生产运行分成 10 个独立段每段重新初始化速度从同一分布采样然后对 10 段结果取平均。这样得到的误差棒才是可信的。如果只跑一条长轨迹相邻帧高度相关误差会被严重低估。另外输出频率影响相关时间。每步都输出相邻帧的 MSD 差很小拟合时权重不均。我一般每 10 步输出一次保证相邻帧位移足够大又不至于丢失短时信息。写 MD 代码这些年最大的教训是别急着优化性能先把能量守恒跑对。我见过太多人一上来就上格点法、并行化结果能量漂移得没法看回头查 bug 花的时间比省下来的机时多得多。先用 256 个粒子、dt0.001、Berendsen 热浴把 (g(r)) 和 MSD 跑出合理结果再逐步加粒子数、换算法。这个顺序不能反。希望帮到你。本文还有配套的精品资源点击获取
简介这份资源是经典教材《The Art of Molecular Dynamics Simulation》D.C. Rapaport 著Cambridge University Press 出版随书附带的 C 语言程序代码面向物理、化学、材料与生物等方向的学习者和科研人员帮助读者把分子动力学理论落到可编译运行的实践层面。压缩包共 127 个文件约 264KB其中 66 个 .c 源文件构成核心算法与示例程序52 个 .in 输入文件用于配置模拟参数另有 4 个 .h 头文件、2 个 .sh 脚本及 readme、copying、errata 等说明文档结构清晰、便于按章节检索。已有 50 人学习。代码覆盖分子模型构建、能量最小化、系统时间演化等关键环节并涉及并行计算与向量化等加速思路读者可据此理解科学计算程序的框架设计与优化技巧也可将其作为教学辅助或科研工作的起点进行算法改进与特定问题的模拟验证。1. 分子动力学模拟的 C 语言实现从理论公式到可运行代码很多人第一次接触分子动力学MD模拟是从一本经典教材开始的但真正动手把书里的公式变成能跑的 C 代码时往往会卡在几个地方Verlet 列表怎么建、周期性边界怎么处理、能量为什么一直漂移。这个标题指向的正是这条路径——用 C 语言从零实现一套分子动力学模拟的核心框架。它适合两类人一是学过统计物理但没写过 MD 代码的研究生二是需要把 MD 嵌入自己仿真流程的工程师。C 语言在这里不是随便选的MD 每步要算几十万次粒子间作用力C 的指针操作和内存布局控制能直接把性能拉满而 Python 或 MATLAB 在这个场景下往往差一个数量级。下面从最小可运行系统开始把每一步的参数和坑讲清楚。2. 最小可运行 MD 系统Lennard-Jones 流体与 Velocity-Verlet 积分2.1 为什么选 Lennard-Jones 势作为起点分子动力学模拟的核心是求解牛顿方程。对于最简单的流体系统Lennard-JonesLJ势是最常用的模型势[ U(r) 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right] ]其中 (\epsilon) 是势阱深度(\sigma) 是零势距离。选它作为起点有三个理由第一公式简单求导后力表达式干净第二氩气等简单流体的性质已经被研究得很透方便验证结果第三截断处理有成熟方案不会一上来就陷入长程力的复杂性。常见做法是把 (\epsilon) 和 (\sigma) 作为约化单位即令 (\epsilon1)、(\sigma1)质量 (m1)。这样所有物理量都变成无量纲数代码里不需要反复乘常数。实际模拟氩气时再按 (\epsilon0.0104) eV、(\sigma3.40) Å 换算回真实单位。2.2 初始化位置、速度与周期性边界初始化分三步放置粒子、赋初速度、处理边界。放置粒子最忌讳的是让两个粒子靠得太近因为 LJ 势在 (r \to 0) 时趋于无穷直接导致力爆炸。我一般用面心立方FCC晶格排布再给每个粒子加一点随机扰动这样既保证初始无重叠又不会让系统一直保持晶格结构。// 初始化FCC 晶格 随机扰动 #include math.h #include stdlib.h #define N 256 // 粒子数取 4 的倍数便于 FCC 排布 #define RHO 0.8 // 约化密度 #define BOX pow(N / RHO, 1.0/3.0) // 立方盒子边长 typedef struct { double x, y, z; // 位置 double vx, vy, vz; // 速度 double fx, fy, fz; // 受力 } Particle; void init_positions(Particle *p) { int n (int)cbrt(N / 4.0); // 每边晶胞数 double a BOX / n; // 晶胞边长 int idx 0; // FCC 基元4 个原子 per 晶胞 double basis[4][3] { {0, 0, 0}, {0.5, 0.5, 0}, {0.5, 0, 0.5}, {0, 0.5, 0.5} }; for (int i 0; i n; i) for (int j 0; j n; j) for (int k 0; k n; k) for (int b 0; b 4; b) { if (idx N) break; p[idx].x (i basis[b][0]) * a 0.01 * (rand()/(double)RAND_MAX - 0.5); p[idx].y (j basis[b][1]) * a 0.01 * (rand()/(double)RAND_MAX - 0.5); p[idx].z (k basis[b][2]) * a 0.01 * (rand()/(double)RAND_MAX - 0.5); idx; } }这段代码里RHO0.8是约化密度对应液态氩的典型值。扰动幅度取 0.01 倍晶胞边长太大可能造成近邻重叠太小则系统长时间保持晶格特征。BOX用pow计算保证密度精确。初速度一般从麦克斯韦-玻尔兹曼分布采样然后减去质心速度防止整个系统漂移。温度通过速度平方和定义[ T \frac{1}{3N} \sum_i m v_i^2 ]在约化单位下 (m1)所以直接累加速度平方即可。周期性边界是 MD 的命门。粒子从盒子右边出去要从左边回来。位置更新后立刻做环绕// 周期性边界处理 void apply_pbc(Particle *p) { for (int i 0; i N; i) { p[i].x - BOX * floor(p[i].x / BOX); p[i].y - BOX * floor(p[i].y / BOX); p[i].z - BOX * floor(p[i].z / BOX); } }用floor而不是if判断是因为粒子一步可能跨过多个盒子长度虽然正常模拟不会但温度失控时会。这个写法比取模更稳负数也能正确处理。2.3 Velocity-Verlet 积分一步更新位置、速度、受力积分算法选 Velocity-Verlet因为它时间可逆、能量守恒好而且每步只算一次力。公式如下[ r(t\Delta t) r(t) v(t)\Delta t \frac{1}{2}a(t)\Delta t^2 ] [ v(t\Delta t) v(t) \frac{1}{2}[a(t) a(t\Delta t)]\Delta t ]实现时先更新位置再算新受力最后更新速度// Velocity-Verlet 积分一步 void integrate(Particle *p, double dt) { // 1. 更新位置 for (int i 0; i N; i) { p[i].x p[i].vx * dt 0.5 * p[i].fx * dt * dt; p[i].y p[i].vy * dt 0.5 * p[i].fy * dt * dt; p[i].z p[i].vz * dt 0.5 * p[i].fz * dt * dt; } apply_pbc(p); // 2. 保存旧力 double old_fx[N], old_fy[N], old_fz[N]; for (int i 0; i N; i) { old_fx[i] p[i].fx; old_fy[i] p[i].fy; old_fz[i] p[i].fz; } // 3. 计算新力 compute_forces(p); // 4. 更新速度 for (int i 0; i N; i) { p[i].vx 0.5 * (old_fx[i] p[i].fx) * dt; p[i].vy 0.5 * (old_fy[i] p[i].fy) * dt; p[i].vz 0.5 * (old_fz[i] p[i].fz) * dt; } }dt的选取很关键。LJ 系统典型振动周期在 (10^{-2}) 约化时间量级dt一般取 0.001 到 0.005。取 0.005 能跑得快但能量漂移明显取 0.001 稳但同样步数覆盖的物理时间短。我一般先用 0.001 验证能量守恒确认无误后再放宽到 0.002 做生产运行。2.4 力计算与截断别让长尾力拖垮性能LJ 力公式[ \mathbf{F}{ij} 24\epsilon \left[ 2\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right] \frac{\mathbf{r}{ij}}{r^2} ]直接对每对粒子算力是 (O(N^2))256 个粒子还能忍上千个就必须用 Verlet 列表或格点法。但即使不优化截断也必须做当 (r r_c) 时力已经很小忽略它。常见截断半径 (r_c 2.5\sigma)。// 计算所有粒子受力含截断 #define RCUT 2.5 #define RCUT2 (RCUT * RCUT) void compute_forces(Particle *p) { // 清零 for (int i 0; i N; i) p[i].fx p[i].fy p[i].fz 0.0; // 双重循环 for (int i 0; i N - 1; i) { for (int j i 1; j N; j) { double dx p[i].x - p[j].x; double dy p[i].y - p[j].y; double dz p[i].z - p[j].z; // 最小镜像约定 dx - BOX * round(dx / BOX); dy - BOX * round(dy / BOX); dz - BOX * round(dz / BOX); double r2 dx*dx dy*dy dz*dz; if (r2 RCUT2 || r2 1e-12) continue; double inv_r2 1.0 / r2; double inv_r6 inv_r2 * inv_r2 * inv_r2; double inv_r12 inv_r6 * inv_r6; double f 24.0 * (2.0 * inv_r12 - inv_r6) * inv_r2; p[i].fx f * dx; p[i].fy f * dy; p[i].fz f * dz; p[j].fx - f * dx; p[j].fy - f * dy; p[j].fz - f * dz; } } }这里用round实现最小镜像约定比if判断更简洁。r2 1e-12是防止两个粒子完全重合导致除零。截断半径 2.5 是 LJ 系统的经验值此时势能约为 (-0.016\epsilon)力约为 (0.039\epsilon/\sigma)对大多数热力学量影响在 1% 以内。如果做相变或界面研究需要加长程修正或换更复杂的势。3. 让模拟跑得稳温度控制、能量诊断与性能优化3.1 温度耦合Berendsen 热浴的参数怎么设刚初始化的系统温度往往偏离目标值需要热浴调节。Berendsen 热浴是最简单的每步把速度乘以一个因子 (\lambda)让温度向目标 (T_0) 弛豫。[ \lambda \sqrt{1 \frac{\Delta t}{\tau_T}\left(\frac{T_0}{T} - 1\right)} ]// Berendsen 热浴 void berendsen_thermostat(Particle *p, double dt, double T0, double tau_T) { double sum_v2 0.0; for (int i 0; i N; i) sum_v2 p[i].vx*p[i].vx p[i].vy*p[i].vy p[i].vz*p[i].vz; double T sum_v2 / (3.0 * N); // 约化单位下 m1, kB1 double lambda sqrt(1.0 dt / tau_T * (T0 / T - 1.0)); for (int i 0; i N; i) { p[i].vx * lambda; p[i].vy * lambda; p[i].vz * lambda; } }tau_T是耦合时间常数一般取 0.1 到 1.0 约化时间。取太小如 0.01会让速度分布被强行拉成窄峰动力学性质失真取太大如 10则温度弛豫慢浪费机时。我一般用 0.5平衡阶段用 0.1 加速收敛生产阶段换回 0.5。注意 Berendsen 热浴不产生正确的正则系综它只是让温度“看起来对”。如果要做涨落相关的计算得换 Nosé-Hoover 或 Andersen 热浴。但作为入门和快速平衡Berendsen 足够。3.2 能量诊断总能量漂移多少算正常每步输出总能量动能势能看它是否守恒。Velocity-Verlet 在dt0.001时能量漂移应该在 (10^{-4}) 每粒子每步量级。如果漂移达到 (10^{-2})说明dt太大或力计算出错。// 计算总能量 double total_energy(Particle *p) { double ke 0.0, pe 0.0; for (int i 0; i N; i) ke 0.5 * (p[i].vx*p[i].vx p[i].vy*p[i].vy p[i].vz*p[i].vz); for (int i 0; i N - 1; i) { for (int j i 1; j N; j) { double dx p[i].x - p[j].x; double dy p[i].y - p[j].y; double dz p[i].z - p[j].z; dx - BOX * round(dx / BOX); dy - BOX * round(dy / BOX); dz - BOX * round(dz / BOX); double r2 dx*dx dy*dy dz*dz; if (r2 RCUT2 || r2 1e-12) continue; double inv_r6 1.0 / (r2*r2*r2); pe 4.0 * (inv_r6 * inv_r6 - inv_r6); } } return ke pe; }势能求和时注意截断否则会重复计算。输出频率一般每 100 步一次太频繁拖慢速度太稀疏看不出漂移趋势。3.3 Verlet 列表与格点法什么时候必须上(O(N^2)) 在 (N256) 时每步约 3 万次对算还能接受。但 (N1000) 时是 50 万次(N10000) 时是 5000 万次必须优化。Verlet 列表的思路是对每个粒子记录截断半径加缓冲区内所有邻居每 10 到 20 步重建一次列表中间步只遍历列表。// Verlet 列表结构简化版 typedef struct { int count; int *neighbors; // 动态数组 } VerletList; // 重建列表缓冲半径 RCUT SKIN #define SKIN 0.3 void build_verlet_list(Particle *p, VerletList *vl) { double r_skin2 (RCUT SKIN) * (RCUT SKIN); for (int i 0; i N; i) { vl[i].count 0; for (int j 0; j N; j) { if (i j) continue; double dx p[i].x - p[j].x; double dy p[i].y - p[j].y; double dz p[i].z - p[j].z; dx - BOX * round(dx / BOX); dy - BOX * round(dy / BOX); dz - BOX * round(dz / BOX); if (dx*dx dy*dy dz*dz r_skin2) vl[i].neighbors[vl[i].count] j; } } }SKIN取 0.3 是经验值太小则列表频繁重建太大则列表变长、遍历变慢。重建频率取决于粒子最大位移一般每 10 步重建一次。格点法Cell List更适合超大系统把盒子切成边长大于截断半径的格子只检查相邻格子复杂度降到 (O(N))。但格点法实现复杂容易在边界处理上翻车建议先用 Verlet 列表把 (N) 跑到 5000 左右再考虑格点。4. 避坑与排查MD 代码最容易翻车的五个地方4.1 能量突然爆炸温度冲到 (10^6)现象跑了几百步后总能量指数上升粒子飞出盒子。原因两个粒子距离过近LJ 势的 (r^{-12}) 项产生巨大斥力一步就把速度打到天上。解决初始化时检查最小粒子间距确保大于 (0.8\sigma)积分步长dt降到 0.0005 先跑 1000 步加一个力上限保护当 (r 0.5\sigma) 时直接跳过或设力为固定值。4.2 温度一直往下掉系统“冻结”现象Berendsen 热浴开着但温度缓慢下降。原因截断半径处势能不为零能量不守恒Velocity-Verlet 会慢慢损失能量。解决加截断位移修正或者换用带长程修正的势。更简单的办法是每步都调热浴但这样动力学性质不准。我一般接受小幅漂移只要在生产运行的时间尺度内温度变化不超过 5%。4.3 最小镜像约定写错力方向反了现象两个粒子明明靠近却互相吸引而不是排斥。原因round(dx/BOX)用错或者dx没有先做减法。解决记住最小镜像的公式是 (dx dx - L \cdot \text{round}(dx/L))其中 (L) 是盒子边长。测试方法放两个粒子在盒子两边距离应该是 (L - r) 而不是 (r)。4.4 随机数没播种每次运行结果一样现象想跑不同初始条件但输出完全一致。原因rand()默认种子是 1。解决在main开头加srand(time(NULL))或者从命令行读入种子。但注意做对照实验时反而要固定种子否则无法复现。4.5 输出文件太大硬盘写满现象跑了 10 万步轨迹文件几十 GB。原因每步都写所有粒子坐标。解决每 100 或 1000 步写一次只写需要分析的粒子或只写能量。如果要做可视化可以写 XYZ 格式但控制帧数。我一般每 500 步写一帧10 万步也就 200 帧足够看结构演化。5. 从能跑到好用径向分布函数与扩散系数的计算技巧5.1 径向分布函数 (g(r))分箱宽度与归一化径向分布函数是 MD 最常输出的结构量。它描述以某个粒子为中心距离 (r) 处找到另一个粒子的概率密度相对于均匀分布的比例。计算时先做直方图再归一化[ g(r) \frac{V}{N^2} \frac{1}{4\pi r^2 \Delta r} \sum_{i} \sum_{j \neq i} \delta(r - r_{ij}) ]// 计算径向分布函数 #define NBINS 200 #define RMAX 5.0 void compute_rdf(Particle *p, double *g_r, int *count) { double dr RMAX / NBINS; for (int b 0; b NBINS; b) { g_r[b] 0.0; count[b] 0; } for (int i 0; i N - 1; i) { for (int j i 1; j N; j) { double dx p[i].x - p[j].x; double dy p[i].y - p[j].y; double dz p[i].z - p[j].z; dx - BOX * round(dx / BOX); dy - BOX * round(dy / BOX); dz - BOX * round(dz / BOX); double r sqrt(dx*dx dy*dy dz*dz); if (r RMAX) { int bin (int)(r / dr); count[bin] 2; // i-j 和 j-i 都算 } } } for (int b 0; b NBINS; b) { double r (b 0.5) * dr; double shell_vol 4.0 * M_PI * r * r * dr; double rho N / (BOX * BOX * BOX); g_r[b] count[b] / (N * rho * shell_vol); } }分箱宽度dr取 0.025即 RMAX/NBINS太宽则峰形模糊太窄则统计噪声大。归一化时注意除以 (N \times \rho \times 4\pi r^2 dr)其中 (\rho) 是数密度。跑 1000 步后取平均(g(r)) 的第一个峰应该在 (r \approx 1.1\sigma)高度约 2.8这是 LJ 液体的特征。5.2 扩散系数从均方位移斜率提取扩散系数 (D) 通过均方位移MSD的长时间斜率得到[ D \lim_{t \to \infty} \frac{1}{6t} \langle |\mathbf{r}(t) - \mathbf{r}(0)|^2 \rangle ]实现时记录每个粒子相对于初始位置的位移平方和对时间做平均。注意要用未折叠坐标否则粒子穿过边界时位移会跳变。// 均方位移计算需保存初始位置和未折叠坐标 void compute_msd(Particle *p, Particle *p0, double *msd, int step) { double sum 0.0; for (int i 0; i N; i) { double dx p[i].x - p0[i].x; double dy p[i].y - p0[i].y; double dz p[i].z - p0[i].z; sum dx*dx dy*dy dz*dz; } msd[step] sum / N; }实际代码里要维护未折叠坐标每次 PBC 环绕时记录粒子移动了多少个盒子长度累加到未折叠坐标上。MSD 曲线先经历弹道区(t^2)再进入扩散区线性。取线性段做最小二乘拟合斜率除以 6 就是 (D)。LJ 液体在 (T0.8)、(\rho0.8) 时(D) 约 0.05 到 0.1 约化单位。5.3 一个容易忽略的细节相关时间与统计误差MSD 和 (g(r)) 都需要对多个时间起点做平均否则统计误差很大。我一般把生产运行分成 10 个独立段每段重新初始化速度从同一分布采样然后对 10 段结果取平均。这样得到的误差棒才是可信的。如果只跑一条长轨迹相邻帧高度相关误差会被严重低估。另外输出频率影响相关时间。每步都输出相邻帧的 MSD 差很小拟合时权重不均。我一般每 10 步输出一次保证相邻帧位移足够大又不至于丢失短时信息。写 MD 代码这些年最大的教训是别急着优化性能先把能量守恒跑对。我见过太多人一上来就上格点法、并行化结果能量漂移得没法看回头查 bug 花的时间比省下来的机时多得多。先用 256 个粒子、dt0.001、Berendsen 热浴把 (g(r)) 和 MSD 跑出合理结果再逐步加粒子数、换算法。这个顺序不能反。希望帮到你。本文还有配套的精品资源点击获取
