Passa al tema normale
Discussioni su argomenti di Informatica

Regole del forum

Consulta il nostro regolamento e la guida per scrivere le formule
Rispondi al messaggio

[C++] Dinamica molecolare Lennard-Jones

12/09/2018, 21:22

Sto cercando di scrivere un programma di dinamica molecolare seguendo il Frenkel-Smit

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

Re: [C++] Dinamica molecolare Lennard-Jones

13/09/2018, 22:50

Un primo errore che vedo è nella funzione forza. Nel ciclo azzeri il valore della forza corrente, ma poi modifichi il valore della forza di particelle che devono ancora essere visitate. Perché il codice funzioni hai bisogno di separare l'inizializzazione dai due cicli (facendo quindi un ciclo a parte) o considerare il contributo delle particelle precedenti invece che seguenti. Non avendo idea di quali equazioni stai cercando di risolvere non ti saprei dire se il codice è giusto. Ma inizierei certamente dal correggere le forze..

Re: [C++] Dinamica molecolare Lennard-Jones

14/09/2018, 08:27

Ti ringrazio! con questa modifica il codice va meglio, per lo meno le forze non esplodono. In realtà lo pseudocodice sul libro fa esattamente come dici tu, un ciclo a parte, ma io ho sempre ignorato questa cosa perchè pensavo (evidentemente stupidamente e senza riflettere per bene sulla questione) che non cambiasse nulla.

Ho modificato anche le condizioni periodiche al contorno che erano sbagliate, da così $dx=dx-L*\text{int}(dx/L)$ a così $dx=dx-L*\text{round}(dx/L)$

Le equazioni che uso sono queste(nel mio codice $\epsilon=\sigma=1$):

Potenziale di Lennard-Jones $4\epsilon\left[(\frac{\sigma}{r})^12-(\frac{\sigma}{r})^6\right]$

la componente $x$ della forza $f_x(r)=48\frac{x}{r^2}\left[(\frac{\sigma}{r})^12-\frac{1}{2}(\frac{\sigma}{r})^6\right]$

il viriale, e le correzioni ( en_corr(), P_corr() ) dovrebbero essere corretti perchè ho fatto un'altro programma con l'algoritmo di Metropolis e funziona bene.

Le velocità e le posizioni le aggiorno con l'algoritmo di Verlet $r(t+dt) ~~ 2r(t)-r(t-dt)+\frac{f(t)}{m}dt^2$

$v(t)~~\frac{r(t+dt)-r(t-dt)}{2dt}$

nel mio codice $p[i]->x = r(t) \quad,\quad p[i]->xm = r(t-dt)$

riporto anche gli pseudocidici del libro (non so perchè non ci avevo pensato prima)

Testo nascosto, fai click qui per vederlo
Immagine



Immagine




Immagine

Re: [C++] Dinamica molecolare Lennard-Jones

14/09/2018, 08:56

Mi sono dimenticato di dire che il potenziale che uso, Lennard Jones cutted shifted, è fatto così

$u_{l-j}^{cs}(r)={(u_{l-j}(r)-u_{l-j}(r_c),if r<r_c),(0,if r>r_c):}$

dove $u_{l-j}(r)$ è il potenziale di Lennard-Jones che ho scritto sopra

Re: [C++] Dinamica molecolare Lennard-Jones

14/09/2018, 10:44

Guardando il codice ho notato alcune cose:
  1. Usi molte variabili globali, alcune non costanti, che rendono più difficile seguire cosa succede nel codice.
  2. Non sembra che tu stia facendo alcuno sforzo per semplificare o ridurre le operazioni che fai. Per esempio \(\displaystyle 4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right] = 4\epsilon \left( \frac{\sigma}{r} \right)^{6} \left[\left(\frac{\sigma}{r}\right)^{6}-1\right] \)
    e molti valori li ricalcoli molte volte tra le varie funzioni.
  3. Non usi alcune funzioni che potresti usare. Per esempio dal C++11 puoi usare cbrt per la radice terza. E potevi definire npart come std::ceil( std::cbrt( static_cast<double>( N ) ) ) seppur non sia una buona idea comparare interi come variabili di tipo double. Considerando che N è una costante, avresti potuto calcularti npart direttamente (è uguale a 5). Tra l'altro sono personalmente contrario all'uso di pow per potenze piccole e conosciute.
  4. Non comprendo il senso di avere p come un array fisso di puntatori, a seconda delle ragioni per cui hai fatto qualcosa del genere vedo più sensato usare un std::vector o semplicemente un array di dimensione fissa.

Se ho tempo ricontrollerò usando gli pseudocodici.

Re: [C++] Dinamica molecolare Lennard-Jones

14/09/2018, 12:12

Si lo so, il codice è un po' un disastro dal punto di vista dei punti (scusa la ripetizione) 1 e 2 del tuo messaggio. Questo codice l'ho scritto in fretta, ed è un copia-incolla da un'altro codice di dinamica molecolare scritto da me, nella quale uso le linked-list per ridurre il costo computazionale (ho 'riscritto' questo codice più semplice per capire qual'era l'errore di fondo).
Per quanto riguarda gli altri punti, certe cose non le faccio perchè proprio non le conosco (non ero a conoscenza di cbrt, per esempio), non so se si nota ma non è che sappia programmare così tanto :D...diciamo che questo è il terzo programma (un po' più serio di 'hello world') che scrivo nella mia vita, e si vede!

Re: [C++] Dinamica molecolare Lennard-Jones

14/09/2018, 21:20

Guardando il codice noto che ci sono diverse modifiche che hai apportato agli algoritmi (come gli if intramezzati ai cicli nella inizializzazione per esempio). Siccome non sei particolarmente esperto di programmazione o con questi algoritmi ti consiglio di cercare di seguire gli algoritmi alla lettera prima di cercare di fare ottimizzazione. Come per l'inizializzazione a zero per le forze, ci sono spesso delle ragioni per cui gli algoritmi sono scritti in un modo invece che un altro. Eventuali ottimizzazioni o modifiche vanno in genere fatte dopo che si è arrivati ad un algoritmo funzionante.. Non ti so dire se il problema è legato a queste modifiche, ma avere qualcosa di confrontabile direttamente con il codice di riferimento è molto utile in quanto aiuta a trovare un eventuale discrepanza.

Re: [C++] Dinamica molecolare Lennard-Jones

14/09/2018, 21:28

L'inizializzazione è corretta (almeno quella!), e lo dico perchè nell'altro programma col metodo Metropolis funziona, e plottando il risultato dell'inizializzazione viene esattamente quello che mi aspetto.
Ci sono anche dei codici già scritti da Frenkel online, ma io non conoscco per nulla fortran...
Un problema che ho notato che appare è che l'energia potenziale fluttua pochissimo: a qualsiasi temperatura parto, il sistema (per quanto riguarda l'energia potenziale) sembra sempre già equilibrato, ed è molto strano.
Domani provo a rimetterci mano!
Ma oltre a quella dell'inizializzazione che modifiche grosse hai notato? Perchè ho cercato di rispettare l'algoritmo per quanto possibile

Re: [C++] Dinamica molecolare Lennard-Jones

15/09/2018, 17:38

Con le modifiche che hai già detto di aver fatto alla forza gli algoritmi mi sembrano corretti. Non ho controllato se le formule fossero effettivamente corrette e il calcolo dell'energia. Un metodo utile per fare il debug in problemi di questo tipo è esportare i valori dopo ogni iterazione (o ad intervalli regolari) per vedere come i tuoi valori cambiano nel tempo.

Re: [C++] Dinamica molecolare Lennard-Jones

15/09/2018, 22:16

Ad ora non sono riuscito a risolvere i problemi che ha il programma, escono valori sbagliati e non riesco a trovare gli errori.
Mi spiace che il professore non può ricevermi prima dell'esame e quindi non posso completare il programma (che non è obbligatorio per fare l'esame, ma mi sarebbe piaciuto fare qualcosa fatta bene).
Nonostante non sia facile discutere di queste cose su un forum, un po' siete riusciti ad aiutarmi, quindi grazie :D
Rispondi al messaggio


Skuola.net News è una testata giornalistica iscritta al Registro degli Operatori della Comunicazione.
Registrazione: n° 20792 del 23/12/2010.
©2000— Skuola Network s.r.l. Tutti i diritti riservati. — P.I. 10404470014.