Molecular dynamics Lennard Jones

[cuppls]

Following Frekel-Smit book: https://www.amazon.it/dp/0122673514/
I am trying to write a program based on the algorithms 3,4,5 and 6 in the book.
I use a cutted-shifted Lennard Jones potential, with cutoff radius at 2σ.

#include <iostream>
#include <cstdlib>
#include <cmath>
#include<fstream>
#include<string>
#include<iomanip>
#include<ctime>
#include<chrono>

using namespace std;const double pi=3.141592653589793;
const double kb=pow(1.380649,-23);      //Cost. Boltz. J/K
//const double kb=8.61673324*pow(10,-5);    //Cost. Boltz. eV/K
const int N=108; // Number of particles
double sigma=1,epsilon=1; // Lennard jones parameter
double rc; // Cutoff distance
double L;  // L box lenght
double rhos,rho,Ts,T; //rho*,rho,T*,T
double dt; //time step reduced units dt*=dt/sqrt(m*sigma*sigma/epsilon)
double massa=1; // mass of particles

struct particella{ //class particle
double x,y,z;
double xm,ym,zm;
double vx,vy,vz;
double fx,fy,fz;
};

double lj (double r1) {
     double s12=pow(sigma/r1,12);
     double s6=pow(sigma/r1,6);
     return 4*epsilon*(s12-s6);
}

double viriale(double d){ //virial
  double a=(pow(sigma/d,12)-0.5*pow(sigma/d,6));
return a;
}

double en_corr(){  //energy tail correction
double a=(8./3)*pi*epsilon*rho*pow(sigma,3);
double b=(1./3)*pow(sigma/rc,9)-pow(sigma/rc,3);
return a*b;
}

void forza(particella *p[N],double &en,double &vir){ //force calculation
en=0,vir=0;
double ecut=4*(pow(rc,-12)-pow(rc,-6));
for(int i=0;i<N;++i){
  p[i]->fx=0;
  p[i]->fy=0;
  p[i]->fz=0;
for(int j=i+1;j<N;++j){
  double dx=p[i]->x-p[j]->x;
  double dy=p[i]->y-p[j]->y;
  double dz=p[i]->z-p[j]->z;
  dx=dx-L*round(dx/L);  //periodic buondary conditions
  dy=dy-L*round(dy/L);
  dz=dz-L*round(dz/L);
  double r2=dx*dx+dy*dy+dz*dz;
  if(r2<rc*rc){
      double r2i=1/r2;
      double r6i=r2i*r2i*r2i;
      double ff=48*r2i*r6i*(r6i-0.5);
      p[i]->fx+=ff*dx;
      p[i]->fy+=ff*dy;
      p[i]->fz+=ff*dz;
      p[j]->fx-=ff*dx;
      p[j]->fy-=ff*dy;
      p[j]->fz-=ff*dz;
      en+=4*r6i*(r6i-1)-ecut;
      vir+=r6i*(r6i-0.5);
   }
  }
}
}

void verlet (particella *p[N],double &vx_cm,double &vy_cm,double &vz_cm,double &v2){  //verlet algorithm
vx_cm=0,vy_cm=0,vz_cm=0,v2=0;
for(int i=0;i<N;++i){
  double px,py,pz;

  px=2*(p[i]->x)-p[i]->xm+(p[i]->fx)*dt*dt;
  p[i]->vx=(px-p[i]->xm)/(2*dt);

  py=2*(p[i]->y)-p[i]->ym+(p[i]->fy)*dt*dt;
  p[i]->vy=(py-p[i]->ym)/(2*dt);

  pz=2*(p[i]->z)-p[i]->zm+(p[i]->fz)*dt*dt;
  p[i]->vz=(pz-p[i]->zm)/(2*dt);
  vx_cm+=p[i]->vx;
  vy_cm+=p[i]->vy;
  vz_cm+=p[i]->vz;
  v2+=vx_cm*vx_cm+vy_cm*vy_cm+vz_cm*vz_cm;
 
  p[i]->xm=p[i]->x;
  p[i]->ym=p[i]->y;
  p[i]->zm=p[i]->z;
  p[i]->x=px;
  p[i]->y=py;
  p[i]->z=pz;
}
}

double P_corr (){
double a=(16./3)*pi*pow(rho,2)*epsilon*pow(sigma,3);
double b=(2./3)*pow(sigma/rc,9)-pow(sigma/rc,3);
return a*b;
}

void stampapos (particella *p[N],string name_file){
       ofstream os;                    
      os.open(name_file);
      os<<"# X Y Z\n";
     for(int i=0; i<N;++i){
    os<<setprecision(6)<< p[i]->x << ' ' <<
    p[i]->y << ' '<<p[i]->z << '\n';
    }
    os.close();
}
int main(){
srand48(time(NULL));

cout << "Densità (in unità ridotte): ";
cin >> rhos;
cout << "Temperatura (in unità ridotte): ";
cin >> Ts;
rc=2*sigma;
L=pow(N/rhos,1./3)*sigma;
rho=rhos/pow(sigma,3);
T=Ts*epsilon/kb;
dt=0.001;

cout << "# particelle (N): " << N << '\n';
cout << "Lato box (L): " << L << '\n';
cout << "Raggio cutoff (rc): " << rc << '\n';
cout << "Densità (reale): " << rho << '\n';
cout << "Time step : "<<dt<<'\n';
cout << "Temperatura (reale): "<< T << " K " << '\n';particella  *p[N];double npart=(pow(N,1./3)-int(pow(N,1./3)))==0?int(pow(N,1./3)):int(pow(N,1./3))+1;
const double dist=L/npart; //distanza tra particelle nel scc
const double e=dist/2;
static int count=0;

double v_cm[3]={0,0,0};//center of mass velocity
double v_quadro=0;
for(int i=0; i<npart;++i){//initialization, put particles in cubic lattice with random velocities
      if(i*dist+e<L && count < N){
  for(int j=0;j<npart;++j && count < N){
      if(j*dist+e<L){
    for(int k=0;k<npart;++k){
      if(k*dist+e<L && count < N){
         double a=i*dist+e;
         double b=j*dist+e;
         double c=k*dist+e;
         p[count]=new particella;
         p[count]->x=a;
         p[count]->y=b;
         p[count]->z=c;
         p[count]->vx=drand48()-0.5;
         p[count]->vy=drand48()-0.5;
         p[count]->vz=drand48()-0.5;
         double v0=p[count]->vx;  
         double v1=p[count]->vy;
         double v2=p[count]->vz;
         v_cm[0]+=v0;  
         v_cm[1]+=v1;
         v_cm[2]+=v2;
         v_quadro+=v0*v0+v1*v1+v2*v2;
         count+=1;                 
          }
          else break;
         }
        }
        else break;
       }        
      }
       else break;
    }v_cm[0]/=N;  
v_cm[1]/=N;
v_cm[2]/=N;
v_quadro/=N;

double fs=sqrt(3*Ts*epsilon/v_quadro);//scale factor for temperature

for(int i=0;i<N;++i){
p[i]->vx=(p[i]->vx-v_cm[0])*fs;
p[i]->vy=(p[i]->vy-v_cm[1])*fs;
p[i]->vz=(p[i]->vz-v_cm[2])*fs;
p[i]->xm=p[i]->x-(p[i]->vx)*dt;
p[i]->ym=p[i]->y-(p[i]->vy)*dt;
p[i]->zm=p[i]->z-(p[i]->vz)*dt;
}

double beta=1/(T*kb); //beta*=1/T*
double U=0,V=0,V2=0;

ofstream os,os1,os2;
os.open("energia.dat");
os1.open("k.dat");
os2.open("u.dat");
os<<"# X Y\n";
os1<<"# X Y\n";
os2<<"# X Y\n";for(int j=0;j<100000;++j){
double en,vir,vx_cm,vy_cm,vz_cm,v2;
forza(p,en,vir);
verlet(p,vx_cm,vy_cm,vz_cm,v2);
}

stampapos(p,"pos.dat");

for(int j=0;j<500;++j){
double en,vir,vx_cm,vy_cm,vz_cm,v2;
forza(p,en,vir);
verlet(p,vx_cm,vy_cm,vz_cm,v2);
U+=en;
V+=vir;
V2+=v2;
double etot=(U+0.5*V2)/N;
os<<setprecision(6)<<' '<<j<<' '<<etot<<'\n';
os1<<setprecision(6)<<' '<<j<<' '<<0.5*v2/N<<'\n';
os2<<setprecision(6)<<' '<<j<<' '<<en/N<<'\n';
}
os.close();

cout << "\nu*="<<U/(N*500*epsilon)+en_corr()<<'\n';
cout << "\nP*="<<(48*epsilon*V*pow(sigma/L,3))/(3*500)+rho/beta+P_corr()<<'\n';}

This is the code. My problem is that particles get too close,the Lennard Jones force blow up, and the energy and pressure that result are always nearly 10⁹.
Can someone help me?

 

[Mark44]

My problem is that particles get too close,the Lennard Jones force blow up, and the energy and pressure that result are always nearly 10⁹.
Can someone help me?

You haven't given us much to go on. Do you know how to use a debugger?

 

[cuppls]

Thanks for reply. I used sometimes gdb to find segmentation faults, but not more than that. I have followed exactly the pseudocodes in the Frenkel's book, and I have much difficulties to find where is the error..
What informations should I give to allow you to help me?

 

[Mark44]

What informations should I give to allow you to help me?

Which specific line/s of code is/are producing results that aren't correct.

 

[DrClaude]

My problem is that particles get too close,the Lennard Jones force blow up, and the energy and pressure that result are always nearly 10⁹.

Check your time step. It must small enough (compare with velocity) that the particles can't move artificially close to each other.