Метод Монте-Карло для систем многих частиц - C++ - Ответ 9006997

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <iomanip>
 
using namespace std;
 
const int ntot = 20000;
FILE *out = stdout;
 
 
class MD
{
    private:
        const int N = 5; // количество частиц
        double dt;
    double rn[N][3], L[3], L2[3];
        double r[N][3];
        double f[N][3]; // координаты, сила
        double utot; 
        double E_kin;
    public:
        double v[N][3];
 
        MD()
        {
            N = 5;
            dt = 0.0005;
            utot = 0;
            E_kin = 0;
 
            for(int i = 0; i < 3; i++)
            {
                L[i] = 3;
                L2[i] = L[i] / 2;
            }
 
            out = fopen("result.csv","ab+");
        }
 
 
    void nearest_image(); // приведение к основной ячейке
        void ClearForces();// f = 0;
        void ReadCoordinates(); // чтение координат из файла
        void CalcForces(); //Расчет сил
        void EqMotion();
        void KineticEnergy();
        void SaveState(int n);
};
 
void MD::ClearForces()
{
 
    for (int i = 0; i < N; i ++)
        for (int j = 0; j < 3; j ++)
            f[i][j] =0;
}
 
void MD::ReadCoordinates()
{
    FILE *a = stdout;
    a = fopen("coordinates.csv","r");
 
    for (int i = 0; i < N; i ++)
        for (int j = 0; j < 3; j ++)
            fscanf(a,"%lf;", &r[i][j]);
 
    fclose(a);
}
 
void MD::CalcForces() //Расчет сил
{
 
 
    double rij[3]; // расстояние между частицами
    double Force = 0; // сила
    double r1, u1, u6, u12;
    utot = 0;
 
    for(int i = 1; i < 5; i ++)
        for(int j = 0; j < i; j ++)
        {
            double r2 = 0;
            for(int k=0; k<3; k++)
            {
                rij[k] = rn[i][k] - rn[j][k];
                if(rij[k] > L2[k]) rij[k] -= L[k];
                else if(rij[k] < -L2[k]) rij[k] += L[k];
                r2 += rij[k] * rij[k];
            }
 
            r1 = sqrt(r2);
            u1 = r1*r1*r1;
            u6 = u1*u1;
            u12 = u6*u6;
            utot += 4.*(1./u12 - 1./u6); // расчет потенциала
            Force = -24.*(u6-2.)/(u12 *r1*r1);
 
            for(int k = 0; k < 3; k ++)
            {
                f[i][k] += Force * rij[k];
                f[j][k] -= Force * rij[k];
            }
        }
}
 
void MD::EqMotion() //Решение ур.движения
{
    for(int i = 0; i < 5; i ++)
        for(int k = 0; k < 3; k ++)
        {
            v[i][k] += f[i][k] * dt;
            r[i][k] += v[i][k] * dt;
        }
}
 
void MD::KineticEnergy()
{
    double velocity = 0;
 
    for(int i = 0; i < 5; i ++)
        for(int j = 0; j < 3; j ++)
            velocity += v[i][j]*v[i][j];
 
        E_kin = velocity/2.;
}
 
void MD::nearest_image()
{
    for(int i=0; i<5; i++)
        for(int k=0; k<3; k++)
        {
            if(r[i][k] > 0)
                rn[i][k] = fmod(r[i][k] + L2[k], L[k]) - L2[k];
            else
                rn[i][k] = fmod(r[i][k] - L2[k], L[k]) + L2[k];
        }
 
}
 
void MD::SaveState(int n) // сохранение результата
{
    fprintf(out, "%.3lf; %.10lf; %.10lf; %.10lf;\n", n*dt, utot, E_kin, utot+E_kin);
}
 
int main()
{
 
    MD A;
    for (int i = 0; i < 5; i ++)
        for (int j = 0; j < 3; j ++)
            A.v[i][j] = 0;
 
    A.ReadCoordinates();
    
 
    for(int n = 0; n < ntot; n ++)
    {
    A.nearest_image();
        A.ClearForces();
        A.CalcForces();
        A.KineticEnergy();
        A.EqMotion();
 
        if(n%10 == 0) A.SaveState(n);
    }
    fclose(out);
    return 0;
}