问题描述
我正在尝试模拟50个粒子系统(氩气m = 40amu)的Lennard Jones电位,我的代码可以编译并运行良好,但是在计算动能时存在问题。我是C ++编码的新手,并且希望对问题出在哪里有一些了解。在运行时打印值对查找错误没有任何帮助。我正在尝试使用速度Verlet算法进行积分。
#include <iostream>
#include <math.h>
#include <fstream>
int main()
{
//std::ofstream fout;
//fout.open("HW2.txt");
std::ofstream out("HW2.txt");
static const double m = 40; // 40 amu Convert later
static const double e = 0.103; // electron-volts
static const double s = 3.405; // angstroms
double t = 0.0;// in ps
static const double dt = 0.05;
static const double stop = 50;// 50 ps
double atoms_x [50];
double atoms_y [50];
double atoms_vx [50];
double atoms_vy [50];
double atoms_fx [50];
double atoms_fy [50];
int i = 0; // iteration var (atoms initial)
while (i < 50) {
atoms_vx[i] = 0.0;
atoms_vy[i] = 0.0;
atoms_x[i] = i % 5;
atoms_y[i] = i / 5;
i++;
}
int q = 0;
while (q < 50) {
int r = 0;
while (r < 50) {
if (q == r) {
r++;
continue; // break for the i=j case
}
if (r < q) {
r++;
continue; // break for double count
}
double dx = atoms_x[r] - atoms_x[q];
double dy = atoms_y[r] - atoms_y[q];
double sqr_dx = pow(dx,2.0);
double sqr_dy = pow(dy,2.0);
double inner = sqr_dx + sqr_dy;
double rij = pow(inner,0.5);
double powsix = s / rij;
double fxadd = (((24 * e * pow(s,6) * dx) / (pow(rij,8))) * (1 - 2 * pow(powsix,6)));
double fyadd = (((24 * e * pow(s,6) * dy) / (pow(rij,6)));
atoms_fx[q] -= fxadd;
atoms_fy[q] -= fyadd;
atoms_fx[r] += fxadd;
atoms_fy[r] += fyadd;
r++;
}
q++;
}
while (t < stop) {
double atoms_x1[50];
double atoms_y1[50];
double atoms_vx1[50];
double atoms_vy1[50];
double atoms_fx1[50];
double atoms_fy1[50];
double utot = 0;
double ktot = 0;
int a = 0; // iteration var (atoms main)
while (a < 50) {
atoms_x1[a] = atoms_x[a] + dt * atoms_vx[a] - (dt * dt * atoms_fx[a]) / (2 * m);
atoms_y1[a] = atoms_y[a] + dt * atoms_vy[a] - (dt * dt * atoms_fy[a]) / (2 * m);
a++;
}
int q = 0;
while (q < 50) {
int r = 0;
while (r < 50) {
if (q == r) {
r++;
continue; // break for the i=j case
}
if (r < q) {
r++;
continue; // break for double count
}
double dx = atoms_x1[q] - atoms_x1[r];
double dy = atoms_y1[q] - atoms_y1[r];
double sqr_dx = pow(dx,2.0);
double sqr_dy = pow(dy,2.0);
double inner = sqr_dx + sqr_dy;
double rij = pow(inner,0.5);
double powsix = s / rij;
double fxadd = (((24 * e * pow(s,6)));
double fyadd = (((24 * e * pow(s,6)));
atoms_fx1[q] -= fxadd;
atoms_fy1[q] -= fyadd;
atoms_fx1[r] += fxadd;
atoms_fy1[r] += fyadd;
double uadd = 4 * e * (pow(powsix,12) - pow(powsix,6));
utot += uadd;
r++;
}
q++;
}
int va = 0;
while (va < 50) {
atoms_vx1[va] = atoms_vx[va] + (dt * (atoms_fx[va] + atoms_fx1[va]) / (2 * m));
atoms_vy1[va] = atoms_vy[va] + (dt * (atoms_fy[va] + atoms_fy1[va]) / (2 * m));
double kadd = 0.5 * m * (pow(atoms_vy1[va],2) + pow(atoms_vx1[va],2));
ktot += kadd;
va++;
}
int n = 0;
while (n < 50) {
atoms_x[n] = atoms_x1[n];
atoms_y[n] = atoms_y1[n];
atoms_vx[n] = atoms_vx1[n];
atoms_vy[n] = atoms_vy1[n];
atoms_fx[n] = atoms_fx1[n];
atoms_fy[n] = atoms_fy1[n];
n++;
}
t += dt;
std::cout << utot << "\t" << ktot << "\t" << t << "\n" ;
}
for (int z = 50 - 1; z >= 0; z--)
std::cout << atoms_x[z] << "\t" << atoms_y[z] << "\t" << z << "\n" ;
}
下面列出的一些输出是总势能-总动能-时间:
3.9917e+07 inf 0.05
7.80208e+07 inf 0.1
5.6369e+06 inf 0.15
2.78139e+09 inf 0.2
-7.5547e-07 inf 0.25
-8.20938e-10 inf 0.3
-4.56581e-08 inf 0.35
-5.83223e-12 inf 0.4
-8.9669e-13 inf 0.45
-1.78096e-13 inf 0.5
-4.28183e-14 inf 0.55
-1.19697e-14 inf 0.6
-3.78067e-15 inf 0.65
-1.32112e-15 inf 0.7
-5.02862e-16 inf 0.75
-2.06359e-16 inf 0.8
-9.13939e-17 inf 0.85
-4.63778e-17 inf 0.9
-4.14284e-17 inf 0.95
-2.06414e-16 inf 1
.
.
.
-2.32717e-35 inf 49
-2.30239e-35 inf 49.05
-2.2779e-35 inf 49.1
-2.25368e-35 inf 49.15
-2.22975e-35 inf 49.2
-2.20609e-35 inf 49.25
-2.1827e-35 inf 49.3
-2.15957e-35 inf 49.35
-2.13671e-35 inf 49.4
-2.11412e-35 inf 49.45
-2.09177e-35 inf 49.5
-2.06969e-35 inf 49.55
-2.04785e-35 inf 49.6
-2.02627e-35 inf 49.65
-2.00493e-35 inf 49.7
-1.98383e-35 inf 49.75
-1.96297e-35 inf 49.8
-1.94235e-35 inf 49.85
-1.92196e-35 inf 49.9
-1.9018e-35 inf 49.95
-1.88187e-35 inf 50
-1.86216e-35 inf 50.05
解决方法
变量atoms_fx
和atoms_fy
未初始化(因为它们是局部变量,所以它们的初始值是Undefined)。因此,对它们的首次访问(atoms_fx[q] -= fxadd;
导致未定义行为。
您可以使用循环或使用
进行初始化double atoms_fx [50] = {0};
如果您使用的不是旧的编译器。