Molecular dynamics Lennard Jones

  • C/++/#
  • Thread starter cuppls
  • Start date
  • #1
9
0

Main Question or Discussion Point

Following Frekel-Smit book: https://www.amazon.it/Understanding-Molecular-Simulation-Algorithms-Applications/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σ.

Code:
#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 109.
Can someone help me?
 

Answers and Replies

  • #2
33,075
4,776
My problem is that particles get too close,the Lennard Jones force blow up, and the energy and pressure that result are always nearly 109.
Can someone help me?
You haven't given us much to go on. Do you know how to use a debugger?
 
  • #3
9
0
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?
 
  • #4
33,075
4,776
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.
 
  • #5
DrClaude
Mentor
7,049
3,199
My problem is that particles get too close,the Lennard Jones force blow up, and the energy and pressure that result are always nearly 109.
Check your time step. It must small enough (compare with velocity) that the particles can't move artificially close to each other.
 

Related Threads for: Molecular dynamics Lennard Jones

Replies
0
Views
316
Replies
0
Views
1K
Replies
3
Views
2K
Replies
1
Views
770
Replies
1
Views
918
  • Last Post
Replies
1
Views
550
  • Last Post
Replies
3
Views
3K
Top