Lennard-Jones势分子动力学模拟势能计算异常排查求助
2026-8-12
问题：WCA势分子动力学模拟势能始终为0排查
正在实现基于Lennard-Jones势（WCA截断）的分子动力学模拟，现有99个配置文件（n=0到99）存储各时间步粒子的位置与速度。已完成配置文件读取、基于粒子间距的势能/受力计算（利用牛顿第三定律避免重复计算）、周期性边界条件（最小镜像约定），但测试时输出的势能始终为0（仅当r_ij>截断半径r_cut时势能返回0），需要排查原因。
相关C++代码
#include <iostream> #include <math.h> #include <fstream> #include <stdlib.h> #include <vector> #include <string> #include <utility> #include <stdexcept> #include <sstream> using namespace std; // 读取配置文件中的粒子状态 void read_input(int num_file, vector<double>& n, vector<double>& x, vector<double>& y, vector<double>& vx, vector<double>& vy, double& lx, double& ly) { ifstream file("configurations/config_" + to_string(num_file) + ".dat"); double num, temp_x, temp_y, temp_vx, temp_vy; if (file.is_open()) { string line; while (getline(file, line)) { // 处理头部的模拟盒尺寸（假设长度为5，如"14 14"） if (line.length() == 5) { stringstream dim_str(line); dim_str >> lx >> ly; continue; } // 处理粒子数据行 stringstream temp_str(line); temp_str >> num >> temp_x >> temp_y >> temp_vx >> temp_vy; n.push_back(num); x.push_back(temp_x); y.push_back(temp_y); vx.push_back(temp_vx); vy.push_back(temp_vy); } file.close(); } } // 计算WCA势能 double calc_pot(double r) { double sigma = 1.0; double epsilon = 1.0; if (r <= pow(2, (1.0 / 6)) * sigma) { double res = 4.0 * epsilon * (pow(sigma / r, 12) - pow(sigma / r, 6)) + epsilon; return res; } else { return 0; } } // 计算受力大小（基于WCA势的导数） double calc_force(double r) { double sigma = 1.0; double epsilon = 1.0; if (r <= pow(2, (1.0 / 6)) * sigma) { double res = (48.0 * epsilon / pow(r, 13)) - (24 * epsilon / pow(r, 7)); return res; } else { return 0; } } // 计算二维距离 double dist(double rx, double ry) { return sqrt(rx * rx + ry * ry); } int main() { int N = 144; double mass = 1.0; double sigma = 1.0; double epsilon = 1.0; double tau = 1.0; double dt = 0.02; vector<double> n, x, y, vx, vy; double lx, ly; vector<double> epot(N), f_x(N), f_y(N); for (int i = 0; i < N; i++) { epot[i], f_x[i], f_y[i] = 0; } for (int k = 0; k <= 99; k++) { string fname = "output/epot_" + to_string(k) + ".txt"; ofstream output(fname); read_input(k, n, x, y, vx, vy, lx, ly); for (int i = 0; i < N - 1; i++) { double rix = x[i]; double riy = y[i]; for (int j = i + 1; j < N; j++) { double rjx = x[j]; double rjy = y[j]; // 周期性边界修正粒子位置 if (rix > lx) { rix -= lx; } if (riy > ly) { riy -= ly; } if (rjx > lx) { rjx -= lx; } if (rjy > ly) { rjy -= ly; } if (rix < 0) { rix += lx; } if (riy < 0) { riy += ly; } if (rjx < 0) { rjx += lx; } if (rjy < 0) { rjy += ly; } // 计算相对位置分量 double dist_x = rix - rjx; double dist_y = riy - rjy; // 最小镜像约定修正相对距离 if (abs(dist_x) > lx / 2) { dist_x = (lx - abs(dist_x)) * (-dist_x) / abs(dist_x); } if (abs(dist_y) > ly / 2) { dist_y = (ly - abs(dist_y)) * (-dist_y) / abs(dist_y); } // 计算受力分量（此处存在错误） f_x[i] += calc_force(dist_x) * (1 / dist(dist_x, dist_y)); f_y[i] += calc_force(dist_y) * (1 / dist(dist_x, dist_y)); f_y[j] += -calc_force(dist_x) * (1 / dist(dist_x, dist_y)); f_y[j] += -calc_force(dist_y) * (1 / dist(dist_x, dist_y)); // 累加势能 epot[i] += calc_pot(dist(dist_x, dist_y)); } output << fixed << std::setprecision(4) << epot[i] / (N) << endl; } } }
配置文件示例
14 14 0 0 0 1.0292605474705 0.394157727758591 1 0 1.16666666666667 1.05721528014223 1.9850461002085 2 0 2.33333333333333 1.18385526103892 0.143930912297367 3 0 3.5 -0.938850340823852 1.71993225409788 4 0 4.66666666666667 1.99468650405917 0.952210892864475 5 0 5.83333333333333 -0.985361963654284 3.05201529674118 6 0 7 2.84071317501321 0.0689241023507716 7 0 8.16666666666667 3.56152464385237 2.88858201933488 8 0 9.33333333333333 0.147896423269195 1.40592679110988
配置文件头部为模拟盒尺寸，后续每行存储粒子编号、x坐标、y坐标、x方向速度、y方向速度。
WCA势能公式：当r≤2(1/6)σ时，U(r)=4ε[(σ/r)12-(σ/r)^6]+ε，否则为0；受力由F=-dU/dr计算。
问题排查与修复
1. 粒子数据未重置导致索引错位
外层循环处理每个配置文件时，n、x、y等存储粒子位置的向量会不断追加数据，导致后续循环处理的是多个配置文件的混合数据，粒子索引完全错位，实际计算的粒子对间距远大于截断半径，所以势能始终为0。
修复：在每个时间步循环开始时，清空这些向量：
for (int k = 0; k <= 99; k++) { n.clear(); x.clear(); y.clear(); vx.clear(); vy.clear(); // 后续读取文件、计算逻辑 }
2. 势能与受力向量未重置
每个时间步的势能和受力计算应该基于当前配置文件的粒子状态，但代码中epot、f_x、f_y仅初始化一次，后续循环会累加之前的结果，且初始赋值语句无效：
// 原错误语句：仅f_y[i]被赋值0，epot和f_x为随机值 epot[i], f_x[i], f_y[i] = 0;
修复：在每个时间步循环内，用fill重置向量：
for (int k = 0; k <= 99; k++) { fill(epot.begin(), epot.end(), 0.0); fill(f_x.begin(), f_x.end(), 0.0); fill(f_y.begin(), f_y.end(), 0.0); // 后续读取文件、计算逻辑 }
3. calc_force参数误用
calc_force设计为接收粒子间的欧氏距离r，但代码中错误传入了x/y方向的相对距离分量，导致受力计算逻辑混乱，同时也会间接影响势能计算的上下文（粒子对判断错误）。
修复：先计算粒子间欧氏距离，再传入calc_force，并正确拆分受力分量：
double dist_x = rix - rjx; double dist_y = riy - rjy; // 最小镜像修正... double r = dist(dist_x, dist_y); double force_mag = calc_force(r); // 拆分x/y方向受力分量 double fx = force_mag * (dist_x / r); double fy = force_mag * (dist_y / r); f_x[i] += fx; f_y[i] += fy; // 牛顿第三定律：j受到的力与i大小相等方向相反 f_x[j] += -fx; f_y[j] += -fy;
4. 受力分量赋值错误
原代码中把j的x/y方向受力都加到了f_y[j]上，属于低级语法错误，导致受力计算完全错误。
修复：参考上面的代码，分别给f_x[j]和f_y[j]赋值反向受力。
5. 头部判断逻辑不可靠
依赖line.length() == 5判断配置文件头部，当模拟盒尺寸不是两位数时（如"10 10"长度为4），会导致头部被当作粒子数据读取，进而引发所有粒子位置错误。
修复：改为判断是否是第一行，或者判断行内是否仅包含两个数值：
bool is_first_line = true; while (getline(file, line)) { if (is_first_line) { stringstream dim_str(line); dim_str >> lx >> ly; is_first_line = false; continue; } // 处理粒子数据... }
内容的提问来源于stack exchange，提问作者ugur
相关产品推荐
相关产品推荐
