I am writing a molecular dynamics program to create an lattice and populate it with atoms/molecules. These then are given random velocities and the system is initialized. Then throughout time the molecules interact with each other and exert forces on each other. I've tried to make my program as readable as possible.My issue is when the code is run, the values of everything I'm interested in, mainly the vectors _tempreture, _pressure, _totalEnergy, _interactionEnergy, _kineticEnergy all seem to converge (after 5 runs) to certain values. These values then do not change throughout the rest of the runs and sometimes will randomly spike, then drop back down. I can't figure out the unexpected behavior of my code, and I've been looking at it hours upon hours. I hoping one of you bright guys will be able to help me. The code should produce values which fluctuate, they should be similar values but definitely not the same exact values, which is what I get.The code that does most the work is initLatice, this creates a lattice, with random positions, and velocities of the atoms inside. Then simulate is used to calculate interactions, change positions and velocities of all particles, then calculate the energies of that timestep.//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;};#endifin case anyone is wondering the next file doesn't contain the definition for writeToFile, I decided not to paste it on here to save space.//in MD.cpp#include MD.hMD::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);    }}and finally...//in main.cpp#include MD.hint main(){    MD myMD;    myMD.initLatice();    myMD.simulate(400);    std::string fileName = myFile.txt;    myMD.wrtieToFile(fileName);    return 0;}Many thanks!
Unexpected behavior, values tend to converge instead of fluctuate. (MD)
c++;computational physics;molecular dynamics;oop