Неожиданная сходимость расчетных значений, вы можете помочь мне понять, почему?

Я пишу программу молекулярной динамики, чтобы создать решетку и заполнить ее атомами / молекулами. Затем они получают случайные скорости, и система инициализируется. Затем в течение времени молекулы взаимодействуют друг с другом и оказывают силы друг на друга.

Я пытался сделать свою программу максимально читабельной.

Моя проблема в том, что когда код выполняется, значения всего, что меня интересует, в основном векторов _tempreture, _pressure, _totalEnergy, _interactionEnergy, _kineticEnergy, все сходятся (после 5 запусков) к определенным значениям. Эти значения затем не изменяются в течение остальных циклов, а иногда будут случайным образом всплывать, а затем уменьшаться.

Я не могу понять неожиданное поведение моего кода, и я смотрю на это часами. Я надеюсь, что один из вас умных парней сможет мне помочь.

//in MD.h

#ifndef MD_H
#define MD_H

#include <iostream>
#include <random>
#include <time.h>
#include <vector>
#include <string>
#include <fstream>

class MD{
public:
MD();
~MD();
void initLatice();
void simulate(int _runs);
void wrtieToFile(std::string fileName);
double randomValue(double upper, double lower);
private:
double _nMolecules;
double _estKineticEnergy;
double _dt;
double _density;
bool _tempscale;
double _energyCorrection;
double _pressureCorrection;
double _nFCC = 4;
double _truncationSeperation;
double _sidelength;
double _pi = 3.142;
std::vector<double> x_data, y_data, z_data;
std::vector<double> _vX, _vY, _vZ;
std::vector<double> _xPos, _yPos, _zPos;
std::vector<double> _totalEnergy, _kineticEnergy, _interactionEnergy, _pressure, _tempreture;
std::vector<double> _radialDist, _radialDistCoord;
};

#endif

в случае, если кому-то интересно, что следующий файл не содержит определения для writeToFile, я решил не вставлять его сюда, чтобы сэкономить место.

//in MD.cpp
#include "MD.h"
MD::MD() : _radialDist(201, 0.0), _radialDistCoord(201, 0.0) {
_dt = 0.005;
_density = 1.0;
_tempscale = true;
_nMolecules = pow(_nFCC, 3) * 4;
x_data = { 0.25, 0.75, 0.75, 0.25 };
y_data = { 0.25, 0.75, 0.25, 0.75 };
z_data = { 0.25, 0.25, 0.75, 0.75 };
_sidelength = pow((_nFCC / _density), 1. / 3.);
_truncationSeperation = 2.5;
if (_truncationSeperation > _sidelength / 2.0)
_truncationSeperation = _sidelength / 2.0;
_energyCorrection = (8.0 * _pi * _density)*(1.0 / (9.0 * pow(_truncationSeperation, 9)) - 1 / (3.0 * pow(_truncationSeperation, 3)));
_pressureCorrection = (16.0 * _pi * pow(_density, 2))*(2.0 / (9.0 * pow(_truncationSeperation, 9)) - 1 / (3.0 * pow(_truncationSeperation, 3)));
}

MD::~MD(){}

double MD::randomValue(double upper, double lower)
{
double returnValue;
do{
returnValue = upper + (rand() / (RAND_MAX / (lower - upper)));
} while (returnValue > upper || returnValue < lower);
return returnValue;
}

void MD::initLatice(){
int count = 0;
double x_init = 0, y_init = 0, z_init = 0;
double _tempreture = 1.0;
_estKineticEnergy = 0.5 * (3.0 * _nMolecules - 4.0) * _tempreture;
srand(time(NULL));

for (int i = 0; i < 4; i++){
for (int a = 1; a <= _nFCC; a++){
for (int b = 1; b <= _nFCC; b++){
for (int c = 1; c <= _nFCC; c++){
_xPos.push_back(_sidelength*(a - 1 + x_data[i]) / _nFCC);
_vX.push_back(randomValue(0.5, -0.5));
x_init += _vX[count] / _nMolecules;

_yPos.push_back(_sidelength*(b - 1 + y_data[i]) / _nFCC);
_vY.push_back(randomValue(0.5, -0.5));
y_init += _vY[count] / _nMolecules;

_zPos.push_back(_sidelength*(c - 1 + z_data[i]) / _nFCC);
_vZ.push_back(randomValue(0.5, -0.5));
z_init += _vZ[count] / _nMolecules;

count++;
}
}
}
}

double velocityScale = 0, kineticEnergy = 0;

for (int i = 0; i < _nMolecules; i++){
_vX[i] -= x_init;
_vY[i] -= y_init;
_vZ[i] -= z_init;
kineticEnergy += 0.5 * (pow(_vX[i], 2) + pow(_vY[i], 2) + pow(_vZ[i], 2));
}

velocityScale = sqrt(_estKineticEnergy / kineticEnergy);

for (int i = 0; i < _nMolecules; i++){
_vX[i] = _vX[i] * velocityScale;
_vY[i] = _vY[i] * velocityScale;
_vZ[i] = _vZ[i] * velocityScale;
}
}

void MD::simulate(int _runs){
double force = 0;
double seperationSquared = 0;
double interactionEnergy = 0;
double kineticEnergy = 0;
double pressure = 0;

for (int run = 0; run <= _runs; run++){
interactionEnergy = 0;
pressure = 0;
std::vector<double> force_x(256, 0.0);
std::vector<double> force_y(256, 0.0);
std::vector<double> force_z(256, 0.0);

for (int i = 0; i < _nMolecules - 1; i++){
double xI = _xPos[i], yI = _yPos[i], zI = _zPos[i];
for (int j = i + 1; j < _nMolecules; j++){
double x = xI - _xPos[j];
double y = yI - _yPos[j];
double z = zI - _zPos[j];

if (x > _sidelength / 2.0)
x -= _sidelength;
if (y > _sidelength / 2.0)
y -= _sidelength;
if (z > _sidelength / 2.0)
z -= _sidelength;
if (x < -_sidelength / 2.0)
x += _sidelength;
if (y < -_sidelength / 2.0)
y += _sidelength;
if (z < -_sidelength / 2.0)
z += _sidelength;

seperationSquared = pow(x, 2) + pow(y, 2) + pow(z, 2);

if (seperationSquared <= (pow(_truncationSeperation, 2))){
interactionEnergy += 4. * ((1. / pow(seperationSquared, 6)) - (1. / pow(seperationSquared, 3)));
force = 24. * ((2. / pow(seperationSquared, 7)) - (1. / pow(seperationSquared, 4)));
force_x[i] += x * force;
force_y[i] += y * force;
force_z[i] += z * force;
force_x[j] -= x * force;
force_y[j] -= y * force;
force_z[j] -= z * force;
pressure += force * seperationSquared;

int histCounter = ceil(sqrt(seperationSquared) * (200.0/_truncationSeperation));
_radialDist[histCounter] += 1.0;
}
}
}

//thermostating
double velocityScale = 0;

if (_tempscale == true){
kineticEnergy = 0;
for (int i = 0; i < _nMolecules; i++)
kineticEnergy += 0.5 * (pow(_vX[i], 2) + pow(_vY[i], 2) + pow(_vZ[i], 2));
velocityScale = sqrt(_estKineticEnergy / kineticEnergy);
}
else { velocityScale = 1.0; }

kineticEnergy = 0;

//applying verlets leapfrog algorithm
for (int i = 0; i < _nMolecules; i++){
_vX[i] = _vX[i] * velocityScale + force_x[i] * _dt;
_xPos[i] += _vX[i] * _dt;

_vY[i] = _vY[i] * velocityScale + force_y[i] * _dt;
_yPos[i] += _vY[i] * _dt;

_vZ[i] = _vZ[i] * velocityScale + force_z[i] * _dt;
_zPos[i] += _vZ[i] * _dt;

kineticEnergy += 0.5 * (pow(_vX[i], 2) + pow(_vY[i], 2) + pow(_vZ[i], 2));

//check if the particle is still within the box
if (_xPos[i] > _sidelength)
_xPos[i] -= _sidelength;
if (_yPos[i] > _sidelength)
_yPos[i] -= _sidelength;
if (_zPos[i] > _sidelength)
_zPos[i] -= _sidelength;
if (_xPos[i] < 0.0)
_xPos[i] += _sidelength;
if (_yPos[i] < 0.0)
_yPos[i] += _sidelength;
if (_zPos[i] < 0.0)
_zPos[i] += _sidelength;
}

_kineticEnergy.push_back(kineticEnergy / _nMolecules);
_interactionEnergy.push_back((interactionEnergy / _nMolecules) + _energyCorrection);
_totalEnergy.push_back(_kineticEnergy[run] + _interactionEnergy[run]);

if (_tempscale == true)
_tempreture.push_back(_nMolecules * 2. * _kineticEnergy[run] / (3.0 * _nMolecules - 4.0));
else
_tempreture.push_back(_nMolecules * 2. * _kineticEnergy[run] / (3.0 * _nMolecules - 3.0));

_pressure.push_back((2.0 * _kineticEnergy[run] * _nMolecules + pressure) / (3.0*(pow(_sidelength, 3) + _pressureCorrection)));
}

for (int i = 1; i <= 200; i++){
double r = i * _truncationSeperation / 200.0;
_radialDistCoord[i] = _radialDistCoord[i - 1] + 2.0 * _radialDist[i] * 10 / _runs / _nMolecules;
_radialDist[i] = _radialDist[i] / (2.0*_pi*r*r*(_truncationSeperation/200.0)*_density*_runs*_nMolecules);
}}

и наконец…

//in main.cpp
#include "MD.h"
int main(){
MD myMD;
myMD.initLatice();
myMD.simulate(400);
std::string fileName = "myFile.txt";
myMD.wrtieToFile(fileName);
return 0;
}

Большое спасибо!