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