2
$\begingroup$

I'm working on a MD simulation using the Lennard-Jones potential and the Velocity Verlet algorithm. The code compiles and runs, but the energy for each time step is not conserved; it decreases constantly when run. The magnitude of the potential energy is much larger than that of the kinetic energy, but I haven't been able to detect where the error is.

Thank you for any help you can offer!

#include <stdio.H>    
#include <math.H>
#include <stdlib.H>
struct particle
{
    double x[3],v[3],a[3];
};

double sigma, epsilon, m, dt;
#define N 1000
struct particle particles[N];

double rfinder(int id_a, int id_b)
{
    int i;
    double rtot,dist,dist2;
    dist2 = 0;
    for(i=0;i<3;i++)
    {
        dist = particles[id_a].x[i] - particles[id_b].x[i];
        dist2 += dist * dist;
    }
    rtot = sqrt(dist2);
    return rtot;
}

double calcacceleration()
{
    int b,c,i;
    double dist,r,r7,r13,sigma6,sigma12, r12, sor6, sor12, pe;
    pe = 0;
    for(b=0;b<N;b++)
    {
        for(i=0;i<3;i++)
        {
            particles[b].a[i] = 0;
        }
    }

    for(b=0;b<N;b++)
    {
        for(c=0;c<N;c++)
        {
            if(b!=c)
            {
                double dist2;
                dist2 = 0;
                for(i=0;i<3;i++)
                {
                    dist = particles[b].x[i] - particles[c].x[i];
                    dist2 += dist * dist;
                }
                if (dist2<8*8)
                {
                    r = sqrt(dist2);
                    //printf("%d %d %f\n",b,c,r);
                    double rinv = 1.0/r;
                    double r2 = rinv*rinv;
                    double r4 = r2*r2;
                    double r6 = r4 * r2;
                    r7 = r4*r2*rinv;
                    r13 = r7*r4*r2;
                    r12 = r6 * r6;
                    sigma6 = sigma*sigma*sigma*sigma*sigma*sigma;
                    sigma12 = sigma6*sigma6;
                    sor6 = sigma6*r6;
                    sor12 = sigma12*r12;
                    for(i=0;i<3;i++)
                    {
                        dist = particles[b].x[i] - particles[c].x[i];
                        particles[b].a[i] += -(0.5 * (4 * epsilon *((-12 * sigma12 * r13) +(6 * sigma6 * r7))) * ((dist/r))/m);
                        particles[c].a[i] +=  (0.5 * (4 * epsilon *((-12 * sigma12 * r13) +(6 * sigma6 * r7))) * ((dist/r))/m);
                        pe += 0.5*(4*epsilon*(sor12 - sor6));
                    }
                }
            }
        }
    }
    return pe;
}

void position_mover()
{
    double x, x0, v;
    int p, i;
    for(p=0;p<N;p++)
    {
        for(i=0;i<3;i++)
        {
            particles[p].x[i] = particles[p].x[i] + (particles[p].v[i] * dt);
        }
    }
    return;
}

double velocity_mover()
{
    double ke;
    int p, i;
    ke = 0;
    for(p=0;p<N;p++)
    {
        for(i=0;i<3;i++)
        {
            particles[p].v[i] = particles[p].v[i] + (0.5 * particles[p].a[i] * dt);
            ke += 0.5 * m * particles[p].v[i] * particles[p].v[i];
        }
    }
    return ke;
}

void output_to_file()
{
int p;
FILE * fp;
fp = fopen("data.xyz", "a");
fprintf(fp,"%d \n\n",N);
for(p=0;p<N;p++)
{
    fprintf(fp,"H %f %f %f\n",particles[p].x[0],particles[p].x[1],particles[p].x[2]);
}
fclose(fp);
}

int main()
{
int p, i, j, k;
double ft,t,pe, ke,te;
sigma = 1;
epsilon = 1;
FILE * fp;
fp = fopen("data.xyz", "w");
fclose(fp);
printf("What timestep would you like?\n");
scanf("%lf", &dt);
printf("In seconds, what is the final time?\n");
scanf("%lf", &ft);
printf("What mass would you like?\n");
scanf("%lf", &m);
p = 0;
//sets initial positions of atoms to be evenly spaced in a 10 * 10 * 10 cube
for(i=0;i<10;i++)
{
    for(j=0;j<10;j++)
    {
        for(k=0;k<10;k++)
        {
            particles[p].x[0] = i*2;
            particles[p].x[1] = j*2;
            particles[p].x[2] = k*2;
            p++;
        }
    }
}

/*sets initial velocities and accelerations of atoms to 0 to get rid of any
garbage in memory */
for(i=0;i<3;i++)
{
    for(p=0;p<N;p++)
    {
        particles[p].v[i] = 0;
        particles[p].a[i] = 0;
    }
}


for(t=0;t<ft;t+=dt)
{
output_to_file();
calcacceleration();
velocity_mover();
position_mover();
pe = calcacceleration();
ke = velocity_mover();
te = ke + pe;
printf("Time: %f\n PE: %f\n KE: %f\n TE: %f\n",t,pe, ke, te);
}
printf("Done!\n");
return 0;
}
$\endgroup$
  • $\begingroup$ Are you using a time-step that is small enough? $\endgroup$ – nicoguaro Jun 30 '16 at 21:58
  • $\begingroup$ @nicoguaro I've been using .005. $\endgroup$ – L L Jun 30 '16 at 22:04
  • $\begingroup$ What are the units? You are not showing the system, so I can't tell that 0.005 is small or big. How do you determine that value? What happens when you increase it or decrease it? $\endgroup$ – nicoguaro Jun 30 '16 at 22:13
  • $\begingroup$ @nicoguaro I'm using reduced units. A timestep of 1 results in the PE going to 0 after 18 seconds while the kinetic energy blows up. A tilmestep of of .000005 results in the same constant decrease. $\endgroup$ – L L Jun 30 '16 at 22:16
  • $\begingroup$ There are infinitely many reduce units systems. I was trying to get the information for people to help you. But if you don't want to provide it, it's your choice. You are the one asking the question... $\endgroup$ – nicoguaro Jun 30 '16 at 22:18

Your Answer

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

Browse other questions tagged or ask your own question.