时间戳

首页 > 解决方案 > 在分子动力学模拟中计算势能

c - 在分子动力学模拟中计算势能

问题描述

背景

想象一下，我们在一个长度为 L 的盒子里有 N 个粒子，它们相互作用（通过 Lennard Jones 势）。

我想计算系统的总势能。我实现了函数 POT，它计算所有粒子的所有贡献并给出正确的结果（这是经过测试的，可以假设为真）。

我还写了一个函数 POT_ONE，它只计算一个粒子相对于所有其他粒子的势能。这意味着如果我想计算总势能，我将不得不调用这个函数 N 次（确保粒子不与自身相互作用）然后除以 2，因为我重复计算了相互作用。

目标

我的目标是使第二个函数产生与第一个函数相同的结果。

问题

发生了一些非常奇怪的事情：如果我放 4 个粒子，这两个函数给出相同的结果。如果我放第五个，那就有偏差。然后对于 6,7,8 个粒子，它再次给出正确的结果，然后对于 N=9，我得到不同的结果。在 N=1000 的情况下，我从 POT_ONE 得到的结果类似于 113383820348202024。

我的 N=5 的结果是：-0.003911 使用 POT 和 12.864234 使用 POT_ONE

如果有人尝试运行代码并想要检查 N=4 的情况，他/她应该更改定义为全局变量的粒子数 (np)，然后注释该行pos[12]=1;pos[13]=1;pos[14]=1;。

代码

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

标签： c

解决方案

您的错误在POT，您忘记x2在内部循环结束时进行更新。

for (i = 0; i < np - 1; i++) {
    double *x2 = x1 + dim;

    for (j = i + 1; j < np; j++) {
        // ... calculate stuff ..

        x2 += dim;
    }

    x1 += dim;
}

一个更容易且可以说更具可读性的变体是完全放弃指针运算并使用无聊的旧索引：

    for (k = 0; k < dim; k++) {
        dr[k] = x[j * dim + k] - x[i * dim + k];
    }

进一步的观察：

  • 请使您的变量在使用它们的范围内成为本地变量。函数顶部的大量未初始化变量使得跟踪变量变得非常困难，即使在像你这样的短函数中也是如此。
  • 请考虑从函数返回单个值，而不是传入指针。在我看来，这使得像距离平方这样的函数更具可读性。
  • 你的代码结构很难看出来，因为一切都非常紧密地运行在一起，甚至是注释。

推荐阅读