0
$\begingroup$

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;
};

#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!

$\endgroup$
  • 1
    $\begingroup$ Have you tried with higher temperature? I recommend you either use LAMMPS or have a look at this code github.com/akohlmey/ljmd-c it's a simple code written by a professional for educational purposes (Axel Kohlmeyer is renown im MD) $\endgroup$ – Azad May 6 '15 at 13:18
  • $\begingroup$ Have you tried to do the simulation of a spring-mass system? This problem is 1D, but the integration is the same. $\endgroup$ – nicoguaro May 13 '15 at 13:48
2
$\begingroup$

Your algorithm is a first-order Euler and not the Verlet algorithm.

What is the actual box size when you run the code? If you made a mistake and your particles end up being well separated (i.e. further apart than your cutoff) then the force is zero and the kinetic energy will remain constant.

This is just one suggestion to get started.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.