Ho scritto un programma praticamente identico a quello che indica Frenkel (in particolare se avete modo di procurarvi il libro, ho 'copiato' gli algoritmi 3,4,5,6) sul suo libro.
Il problema è che per qualche motivo che non sono riuscito ad individuare le particelle si avvicinano troppo, di conseguenza il potenziale di Lennard-Jones esplode per piccoli $r$ e i risultati mi vengono sempre dell'ordine di $10^9$ o giù di lì.
Purtroppo il professore non può ricevermi prima dell'esame,e sono tre giorni che cerco un maledetto errore che non riesco a trovare, quindi spero che qualcuno di voi possa aiutarmi!
- Codice:
#include <iostream>
#include <cstdlib>
#include <cmath>
#include<fstream>
#include<string>
#include<iomanip>
#include<ctime>
#include<chrono>
using namespace std;
//Parametri
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; // Numero particelle
double sigma=1,epsilon=1; // Parametri lennard-jones
double rc; // Distanza di cutoff
double L; // L box simulazione
double rhos,rho,Ts,T; //rho*,rho,T*,T
double dt; //time step unità ridotte dt*=dt/sqrt(m*sigma*sigma/epsilon)
double massa=1; // massa particelle
struct particella{
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){
double a=(pow(sigma/d,12)-0.5*pow(sigma/d,6));
return a;
}
double en_corr(){
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){
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*int(dx/L);
dy=dy-L*int(dy/L);
dz=dz-L*int(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+=lj(sqrt(r2))-ecut;
vir+=viriale(sqrt(r2));
}
}
}
}
void verlet (particella *p[N],double &vx_cm,double &vy_cm,double &vz_cm,double &v2){
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;
}
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};//velocità del centro di massa
double v_quadro=0;
for(int i=0; i<npart;++i){//piazzo le particelle in un scc, le assegno velocità a caso
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);//scalo le v di questo fattore (dall'equipartizione dell'energia), così assegno la T voluta
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<1000;++j){
double en,vir,vx_cm,vy_cm,vz_cm,v2;
forza(p,en,vir);
verlet(p,vx_cm,vy_cm,vz_cm,v2);
}
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';
}//fine main
Uso un potenziale cutted-shifted. Pongo inizialmente le particelle in un reticolo cubico, poi faccio 1000 step per equilibrare, e con altri 500 step vorrei calcolare l'energia media per particella e la pressione