Time step independence of Molecular Dynamics code

181 views Asked by At

I am writing a basic MD code in C++ using LJ potential for an NVE system. The starting configuration is FCC and the starting velocities are randomly generated.

I am facing a strange problem in that the evolution of the system seems to be independent of the time step I implement, it is my understanding that the energy losses are smaller for small time steps and larger for larger time steps. However I am getting the same result at the end of the simulation in terms of energy whether I run (0.0001step)*(10000steps) or 0.001*1000 and so on.

The entire code is to big for me to post here, so I am posting what I think is relevant and leaving out binning etc., kindly let me know if any additional information is required. I have been through countless codes available online and though they look similar to mine I just am not able to figure out what the difference is and where I am going wrong.

The main cpp contains the following loop

for (int i=0; i<t;i++)
{
    potential_calc(neighlist,fromfile, run_parameters,i);//calculating the force fields
    velverlet(neighlist,fromfile, run_parameters, bin, dt);//calculating the velocities
}

The declarations of the 2 cpp files for potential calculation & verlet integration are

void potential_calc(neighborlist_type *neighlist, config_type *fromfile, potential *run_parameters, int t)

void velverlet(neighborlist_type *neighlist, config_type *fromfile, potential *run_parameters, bin_type *bin, double dt)

The code for calculating the force - potential_calc.cpp is below

for (long i=0; i<fromfile->N; i++)
{
long atom_p=i;


for (long j=0; j<neighlist[i].countsn; j++)
    {

    long atom_s=neighlist[i].numb[j];

    for (int k=0; k<Dim; k++)
        {
        dist[k]= fromfile->r[atom_p][k] - (fromfile->r[atom_s][k] + neighlist[atom_p].xyz[j][k]*fromfile->L[k]);
//the .xyz indicates the image being considered real or mirror(if mirror then in which direction)
        }
    disp2 = pow(dist[0],2)+pow(dist[1],2)+pow(dist[2],2);
    if (disp2<rb2)
        {
        int c1=fromfile->c[atom_p];
        int c2=fromfile->c[atom_s];
        double long force_temp;
        disp=pow(disp2,0.5);

        sig_r6=pow(run_parameters->sigma[c1-1][c2-1]/disp,6);//(sigma/r)^6
        sig_r8=pow(run_parameters->sigma[c1-1][c2-1]/disp,8);//(sigma/r)^8

        run_parameters->pe[atom_p] += (4*run_parameters->epsilon[c1-1][c2-1]*((sig_r6*sig_r6)-sig_r6)) - potential_correction[c1-1][c2-1];

        force_temp=(-1*((48*run_parameters->epsilon[c1-1][c2-1])/pow(run_parameters->sigma[c1-1][c2-1],2)*((sig_r6*sig_r8)-((sig_r8)*0.5))));
        for (int k=0; k<Dim;k++)
            {
            run_parameters->force[atom_p][k]+=force_temp*(-1*dist[k]);
            }
        }
    }
//calculating kinetic energy
run_parameters->ke[atom_p] = 0.5*(pow(fromfile->v[atom_p][0],2)+pow(fromfile->v[atom_p][1],2)+pow(fromfile->v[atom_p][2],2));

}

Once the force calculation is done it goes to the updation of velocity and position in the velverlet.cpp

for (long i=0; i<fromfile->N; i++)
{
for (int j=0; j<Dim; j++)
    {
    fromfile->v[i][j] += (dt*run_parameters->force[i][j]);
    }
}

for (long i=0; i<fromfile->N; i++)
{
for (int j=0; j<Dim; j++)
    {
    fromfile->r[i][j] += dt*fromfile->v[i][j];
    }
}

There may be slight differences in how velocity verlet is implemented by different people but I can't figure out how I am getting time step independent results.

Please help. Any input is appreciated Sorry if any formatting/tagging is wrong, this is the first time I am posting here

0

There are 0 answers