循环展开优化Lennard-Jones粒子加速度计算时输出异常求助
2026-7-8
循环展开优化Lennard-Jones分子动力学函数的输出异常排查
我尝试用步长为2的循环展开技术优化计算Lennard-Jones粒子实际气体性质的C++分子动力学加速度函数，但修改后的代码改变了原函数的输出矩阵a，请帮忙排查问题。
原函数
void computeAccelerations() { int i, j, k; double f, rSqd, rSqd4, rSqd7; double result1, result2, result3; double rij[3]; // position of i relative to j for (i = 0; i < N; i++) { a[i][0] = 0; a[i][1] = 0; a[i][2] = 0; } for (i = 0; i < N - 1; i++) { for (j = i + 1; j < N; j++) { result1 = r[i][0] - r[j][0]; result2 = r[i][1] - r[j][1]; result3 = r[i][2] - r[j][2]; rSqd = result1 * result1 + result2 * result2 + result3 * result3; rSqd4 = rSqd * rSqd * rSqd * rSqd; rSqd7 = rSqd4 * rSqd * rSqd * rSqd; f = 24 * (2 / rSqd7 - 1 / rSqd4); a[i][0] += result1 * f; a[j][0] -= result1 * f; a[i][1] += result2 * f; a[j][1] -= result2 * f; a[i][2] += result3 * f; a[j][2] -= result3 * f; } } }
我的实现
void computeAccelerations() { int i, j, k; double f, rSqd, rSqd4, rSqd7; double result1, result2, result3, result4, result5, result6; double rij[6]; // position of i relative to j for (i = 0; i < N; i++) { a[i][0] = 0; a[i][1] = 0; a[i][2] = 0; } for (i = 0; i < N - 1; i += 2) { for (j = i + 1; j < N; j += 2) { result1 = r[i][0] - r[j][0]; result2 = r[i][1] - r[j][1]; result3 = r[i][2] - r[j][2]; rSqd = result1 * result1 + result2 * result2 + result3 * result3; rSqd4 = rSqd * rSqd * rSqd * rSqd; rSqd7 = rSqd4 * rSqd * rSqd * rSqd; f = 24 * (2 / rSqd7 - 1 / rSqd4); a[i][0] += result1 * f; a[j][0] -= result1 * f; a[i][1] += result2 * f; a[j][1] -= result2 * f; a[i][2] += result3 * f; a[j][2] -= result3 * f; result4 = r[i + 1][0] - r[j + 1][0]; result5 = r[i + 1][1] - r[j + 1][1]; result6 = r[i + 1][2] - r[j + 1][2]; rSqd = result4 * result4 + result5 * result5 + result6 * result6; rSqd4 = rSqd * rSqd * rSqd * rSqd; rSqd7 = rSqd4 * rSqd * rSqd * rSqd; f = 24 * (2 / rSqd7 - 1 / rSqd4); a[i + 1][0] += result4 * f; a[j + 1][0] -= result4 * f; a[i + 1][1] += result5 * f; a[j + 1][1] -= result5 * f; a[i + 1][2] += result6 * f; a[j + 1][2] -= result6 * f; } } }
问题根源
- 大量粒子对未被处理：原函数遍历了所有满足
i < j的粒子对，而你的代码只处理了(i,j)和(i+1,j+1)这两类配对，完全遗漏了(i,j+1)、(i+1,j)这些关键粒子对，直接导致加速度计算不完整。 - 边界处理错误：当粒子总数
N为奇数时，最后一个粒子（索引N-1）完全没有参与任何计算；同时内层循环j += 2会跳过中间的粒子，进一步减少了计算的粒子对数量。 - 循环展开思路错误：循环展开的核心是在单次循环迭代内处理多个连续的循环变量，而不是修改粒子的配对逻辑。你错误地将i和j都按步长2跳跃，破坏了原有的所有粒子对遍历逻辑。
修正方案示例
正确的循环展开应该保留原有的粒子对遍历逻辑，仅在单个循环内多处理几次迭代。以下是两种可行的修正方式：
方式1：对i循环展开步长2
void computeAccelerations() { int i, j; double f, rSqd, rSqd4, rSqd7; double result1, result2, result3; // 初始化加速度矩阵 for (i = 0; i < N; i++) { a[i][0] = 0; a[i][1] = 0; a[i][2] = 0; } // 步长2展开i循环 for (i = 0; i < N - 1; i += 2) { // 处理第i个粒子与所有j>i的粒子对 for (j = i + 1; j < N; j++) { result1 = r[i][0] - r[j][0]; result2 = r[i][1] - r[j][1]; result3 = r[i][2] - r[j][2]; rSqd = result1 * result1 + result2 * result2 + result3 * result3; rSqd4 = rSqd * rSqd * rSqd * rSqd; rSqd7 = rSqd4 * rSqd * rSqd * rSqd; f = 24 * (2 / rSqd7 - 1 / rSqd4); a[i][0] += result1 * f; a[j][0] -= result1 * f; a[i][1] += result2 * f; a[j][1] -= result2 * f; a[i][2] += result3 * f; a[j][2] -= result3 * f; } // 处理第i+1个粒子（需确保i+1不越界） if (i + 1 < N - 1) { for (j = i + 2; j < N; j++) { result1 = r[i+1][0] - r[j][0]; result2 = r[i+1][1] - r[j][1]; result3 = r[i+1][2] - r[j][2]; rSqd = result1 * result1 + result2 * result2 + result3 * result3; rSqd4 = rSqd * rSqd * rSqd * rSqd; rSqd7 = rSqd4 * rSqd * rSqd * rSqd; f = 24 * (2 / rSqd7 - 1 / rSqd4); a[i+1][0] += result1 * f; a[j][0] -= result1 * f; a[i+1][1] += result2 * f; a[j][1] -= result2 * f; a[i+1][2] += result3 * f; a[j][2] -= result3 * f; } } } }
方式2：对j循环展开步长2（更高效）
void computeAccelerations() { int i, j; double f, rSqd, rSqd4, rSqd7; double dx1, dy1, dz1, dx2, dy2, dz2; // 初始化加速度矩阵 for (i = 0; i < N; i++) { a[i][0] = 0; a[i][1] = 0; a[i][2] = 0; } for (i = 0; i < N - 1; i++) { // 步长2展开j循环 for (j = i + 1; j < N - 1; j += 2) { // 处理j dx1 = r[i][0] - r[j][0]; dy1 = r[i][1] - r[j][1]; dz1 = r[i][2] - r[j][2]; rSqd = dx1*dx1 + dy1*dy1 + dz1*dz1; rSqd4 = rSqd*rSqd*rSqd*rSqd; rSqd7 = rSqd4*rSqd*rSqd*rSqd; f = 24 * (2 / rSqd7 - 1 / rSqd4); a[i][0] += dx1 * f; a[j][0] -= dx1 * f; a[i][1] += dy1 * f; a[j][1] -= dy1 * f; a[i][2] += dz1 * f; a[j][2] -= dz1 * f; // 处理j+1 dx2 = r[i][0] - r[j+1][0]; dy2 = r[i][1] - r[j+1][1]; dz2 = r[i][2] - r[j+1][2]; rSqd = dx2*dx2 + dy2*dy2 + dz2*dz2; rSqd4 = rSqd*rSqd*rSqd*rSqd; rSqd7 = rSqd4*rSqd*rSqd*rSqd; f = 24 * (2 / rSqd7 - 1 / rSqd4); a[i][0] += dx2 * f; a[j+1][0] -= dx2 * f; a[i][1] += dy2 * f; a[j+1][1] -= dy2 * f; a[i][2] += dz2 * f; a[j+1][2] -= dz2 * f; } // 处理N为奇数时的最后一个j if (j < N) { dx1 = r[i][0] - r[j][0]; dy1 = r[i][1] - r[j][1]; dz1 = r[i][2] - r[j][2]; rSqd = dx1*dx1 + dy1*dy1 + dz1*dz1; rSqd4 = rSqd*rSqd*rSqd*rSqd; rSqd7 = rSqd4*rSqd*rSqd*rSqd; f = 24 * (2 / rSqd7 - 1 / rSqd4); a[i][0] += dx1 * f; a[j][0] -= dx1 * f; a[i][1] += dy1 * f; a[j][1] -= dy1 * f; a[i][2] += dz1 * f; a[j][2] -= dz1 * f; } } }
内容的提问来源于stack exchange，提问作者PeRaDi
相关产品推荐
相关产品推荐
