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.
//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
in 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.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);
}
}
and finally...
//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;
}
Many thanks!