时间戳

首页 > 解决方案 > Lennard Jones 势能返回 inf 的动能

c++ - Lennard Jones 势能返回 inf 的动能

问题描述

我正在尝试模拟 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, 8))) * (1 - 2 * pow(powsix, 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) * dx) / (pow(rij, 8))) * (1 - 2 * pow(powsix, 6)));

            double fyadd = (((24 * e * pow(s, 6) * dy) / (pow(rij, 8))) * (1 - 2 * pow(powsix, 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

标签： c++simulation

解决方案

变量atoms_fx和atoms_fy是未初始化的（因为它们是局部变量，它们的初始值是未定义的）。因此，第一次访问它们 ( atoms_fx[q] -= fxadd;) 会导致未定义的行为。

您可以使用循环初始化它们，或者使用

double atoms_fx [50] = {0};

如果您不使用旧的编译器。

推荐阅读