问分子动力学模拟中的势能计算

EN

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/*GENERAL PARAMETERS*/
int dim=3; //number of dimensions
int np=5; //number of particles
double L=36.413; //box length (A)
double invL=1/36.413; //inverse of box length

/*ARGON CHARACTERISTICS*/
double sig=3.4; // Angstroms (A)
double e=0.001; // eV


double distSQ(double array[]){
    /*calculates the squared distance given the array x=[dx,dy,dz]*/

    int i;
    double r2=0;
    for(i=0;i<dim;i++) r2+=array[i]*array[i];
    return r2;
}//distSQ

void MIC(double dr[],double L, int dim){
    /* MINIMUM IMAGE CONVENTION: dr[] is the array dr = [dx,dy,dz] describing relative
       positions of two particles, L is the box length, dim the number
       of dimensions */

    int i;
    for(i=0;i<dim;i++) dr[i]-=round(dr[i]*invL)*L;
}//MIC

void POT(double x1[],double* potential){
    /*given the positions of each particle in the form x=[x0,y0,z0;x1,y1,z1;...;xn-1,yn-1,zn-1],
      the number of dimensions dim and particles np, it calculates the potential energy of the configuration*/

    //variables for potential calculation
    int i,j,k;
    double *x2;
    double r2inv; // 1/r^2
    double foo,bar;
    double dr[dim];
    *potential=0; // set potential energy to zero

    //main part of POT
    for(i=0;i<np-1;i++){
        x2=x1+dim;
        for(j=i+1;j<np;j++){
                //calculate relative distances between particles i & j
                //apply periodic BCs and then calculate squared distance
                //and the potential energy between them.
                for(k=0;k<dim;k++) dr[k] = x2[k]-x1[k];
                MIC(dr,L,dim); //periodic boundary conditions
                r2inv=1/distSQ(dr);
                //calculate potential energy
                foo = sig*sig*r2inv;
                bar = foo*foo*foo;
                *potential+=bar*(bar-1);

        }//for j
    x1+=dim;
    }//for i
    *potential*=4*e; //scale and give energy units
}//POT

void POT_ONE(int particle,double pos[],double* potential){

    *potential=0;
    int i,k;
    double dr[dim];
    double r2inv,foo,bar;
    double par_pos[dim];
    int index=particle*dim;
    par_pos[0]=pos[index];
    par_pos[1]=pos[index+1];
    par_pos[2]=pos[index+2];

    for(i=0;i<np;i++){
        if(i!=particle){
            for(k=0;k<dim;k++) dr[k]=pos[k]-par_pos[k];
            MIC(dr,L,dim);
            r2inv=1/distSQ(dr);
            foo=sig*sig*r2inv;
            bar=foo*foo*foo;
            *potential+=bar*(bar-1);
        }
        pos+=dim;
    }
    *potential*=4*e; //scale and give energy units
}//POT_ONE

int main(){

    int D=np*dim;
    double* pos=malloc(D*sizeof(double));
    double potential=0; //calculated with POT
    double U=0; ////calculated with POT_ONE
    double tempU=0;
    pos[0]=0;pos[1]=0;pos[2]=0;
    pos[3]=4;pos[4]=0;pos[5]=0;
    pos[6]=0;pos[7]=4;pos[8]=0;
    pos[9]=0;pos[10]=0;pos[11]=4;
    pos[12]=1;pos[13]=1;pos[14]=1;

    POT(pos,&potential);
    printf("POT: %f\n", potential);

    int i,j;
    for(i=0;i<np;i++){
        POT_ONE(i,pos,&tempU);
        U+=tempU;
    }
    U=U/2;
    printf("POT_ONE: %f\n\n", U);

    return 0;
}