- #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σ.
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?
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?