基于Lennard-Jones系统的Leapfrog积分实现正确性验证问询

AI 小北

2026-7-11

验证Leapfrog积分实现与结果正确性的核心要点

一、标准Leapfrog积分的核心规则

先明确Leapfrog的半步更新逻辑，这是判断实现是否合规的核心：

初始状态：已知所有粒子的位置x(t)、速度v(t)，并基于x(t)计算加速度a(t) = F(x(t))/m
第一步：更新半步速度
v(t + dt/2) = v(t) + 0.5 * dt * a(t)
第二步：用半步速度更新全步位置
x(t + dt) = x(t) + dt * v(t + dt/2)
（若用周期性边界，需在此步后对x(t+dt)做边界折叠，确保粒子在模拟盒内）
第三步：基于新位置x(t+dt)重新计算所有粒子的加速度a(t+dt)
第四步：更新全步速度
v(t + dt) = v(t + dt/2) + 0.5 * dt * a(t+dt)
循环以上步骤完成模拟。

合规性判断：
- 检查代码是否严格遵循半步更新顺序：先算半速→更新位置→算新加速度→更新全速。如果是直接用v(t)更新x(t+dt)再调整v(t)，则不符合Leapfrog规则。
- 单原子情况下，若无外势/边界作用，LJ力为0，加速度恒为0，此时Leapfrog退化为匀速运动——若代码输出的位置随时间线性增长、速度恒定，则逻辑正确。
结果验证：给单原子加一个LJ势阱（比如固定一个虚拟粒子作为作用源），观察粒子是否做符合LJ势的周期性运动，总能量（动能+势能）是否近似守恒。

合规性判断：
- 核心检查点：计算t+dt时刻的力时，所有粒子的位置必须是t+dt时刻的更新后值，不能在更新位置的过程中混用t和t+dt时刻的位置（比如先更新粒子i的位置，再用这个新位置和粒子j的旧位置算力）。
- 若代码中先完成所有粒子的位置更新，再统一遍历所有粒子对计算力，则符合要求；反之则存在逻辑错误。
结果验证：模拟小体系（比如8个粒子），观察总能量的波动——Leapfrog是辛积分，NVE系综下总能量应无明显漂移（数值误差范围内的小幅波动是正常的）。

合规性判断：
- 改进点通常是利用牛顿第三定律（F_ij = -F_ji）减少重复计算，遍历粒子对时只计算i<j的情况，然后同时更新F[i]和F[j]。只要这一步逻辑正确（比如给F[i]加力矢量，给F[j]减同一个力矢量），就不会破坏Leapfrog的积分逻辑。
- 必须验证LJ力的公式是否正确：
  正确的LJ力矢量实现逻辑示例：
```
vec r = x[j] - x[i];
// 处理周期性边界的最小镜像（若使用）
for (int d=0; d<3; d++) {
    if (r[d] > box[d]/2) r[d] -= box[d];
    else if (r[d] < -box[d]/2) r[d] += box[d];
}
double r_sq = dot(r, r);
double sigma_sq = sigma*sigma;
double r6 = pow(r_sq, 3);
double r12 = r6*r6;
double force_scalar = 48 * epsilon * (pow(sigma_sq,6)/r12 - 0.5*pow(sigma_sq,3)/r6) / r_sq;
vec force = r * force_scalar;
F[i] += force;
F[j] -= force;
```
  若代码中的力计算符合上述逻辑（包括最小镜像处理），则力的部分正确。
结果验证：对比改进前后的总能量和粒子运动轨迹，若结果一致（或误差在数值精度内），说明优化有效且未引入错误；同时检查计算效率是否提升（比如粒子数N=100时，力计算的循环次数从N²变为N(N-1)/2）。

能量守恒测试：NVE系综下，总能量（动能+势能）的相对误差应控制在1e-5量级以内（取决于步长dt，dt越小误差越小）。若出现能量随时间线性漂移，大概率是积分顺序错误或力计算错误。
径向分布函数（RDF）对比：模拟LJ液态（比如N=100，密度ρ=0.8442，温度T=0.722），计算RDF并与标准结果对比——RDF的第一个峰位置应在r≈1.12σ，峰值约为2.8。
步长收敛测试：分别用dt=0.005、dt=0.001模拟，观察总能量波动和RDF的变化——步长越小，结果应越稳定，RDF越平滑。

内容的提问来源于stack exchange，提问作者user366312