我正在使用Lennard-Jones 势和Velocity Verlet算法进行 MD 模拟。代码编译运行，但每个时间步的能量不守恒；它在运行时不断减小。势能的大小远大于动能的大小，但我一直无法检测到错误在哪里。

感谢您提供的任何帮助！

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