Molecular dynamics Lennard Jones

In summary: Z\n";for(int i=0; i<1000000;++i){ double en,vir; forza(p,en,vir); verlet(p,v_cm[0],v_cm[1],v_cm[2],v_quadro); U+=en; V+=vir; V2+=en*en; if(i%100==0){ os<<i*dt*sqrt(massa*sigma*sigma/epsilon)<<'\t'<<en<<endl; os1<<i*dt*sqrt(massa*sigma*sigma/epsilon)<<'\t'<<v_quadro*N*0.5<<endl; os2<<i*dt*sqrt(massa*sigma*sigma
  • #1
cuppls
9
0
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σ.

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?
 
Technology news on Phys.org
  • #2
cuppls said:
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
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
cuppls said:
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
cuppls said:
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 to Molecular dynamics Lennard Jones

1. What is Molecular Dynamics Lennard Jones?

Molecular Dynamics Lennard Jones is a computational method used to simulate the interactions and movements of particles in a molecular system. It is based on the Lennard-Jones potential, which describes the attractive and repulsive forces between particles, and allows for the prediction of the system's behavior over time.

2. How does Molecular Dynamics Lennard Jones work?

Molecular Dynamics Lennard Jones works by using Newton's laws of motion to calculate the forces acting on each particle in a system. These forces are then used to update the positions and velocities of the particles, allowing for the simulation of the system's behavior over time. This process is repeated for each time step, resulting in a trajectory of the particles' movements.

3. What type of systems can be simulated using Molecular Dynamics Lennard Jones?

Molecular Dynamics Lennard Jones can be used to simulate a wide range of systems, including liquids, solids, and gases. It can also be applied to biological systems, such as proteins and DNA, as well as more complex systems, such as nanoparticles and polymers.

4. What are the limitations of Molecular Dynamics Lennard Jones?

One limitation of Molecular Dynamics Lennard Jones is that it assumes that the particles in the system are spherical and have a fixed size. This means that it may not accurately model systems with irregularly shaped particles or those with significant size differences. Additionally, the method does not account for quantum effects, so it may not be suitable for systems where these effects are important.

5. What are the applications of Molecular Dynamics Lennard Jones?

Molecular Dynamics Lennard Jones has a wide range of applications in various fields, including chemistry, physics, materials science, and biology. It can be used to study the behavior of molecules and materials, predict their properties, and design new compounds and materials with specific properties. It is also commonly used in drug discovery and development, as well as in understanding biological processes at the molecular level.

Similar threads

  • Programming and Computer Science
Replies
1
Views
693
  • Programming and Computer Science
Replies
15
Views
2K
  • Programming and Computer Science
Replies
23
Views
1K
  • Programming and Computer Science
Replies
2
Views
2K
  • Programming and Computer Science
Replies
17
Views
1K
  • Programming and Computer Science
Replies
5
Views
1K
  • Programming and Computer Science
Replies
2
Views
2K
  • Programming and Computer Science
Replies
5
Views
1K
  • Programming and Computer Science
Replies
14
Views
5K
  • Programming and Computer Science
Replies
2
Views
3K
Back
Top