Termostato de Nosé-Hoover

Grupo: Gabriel Azevedo, Rafael Abel e Thierre F. Conceição.

Termostato de Nosé-Hoover

O termostato de Nosé-Hoover é um algoritmo utilizado para simulação de dinâmica molecular. Esse ensemble é relevante quando o sistema em estudo está em contato com um banho térmico, para manter a temperatura constante. A maneira que o algoritmo de Nosé-Hoover mantém a temperatura constante é a partir da adição de uma variável dinâmica fictícia (um "agente" externo), que atua sobre as velocidades das partículas no sistema, as acelerando ou desacelerando até que estas atinjam a temperatura desejada.

ADICIONAR O RESTO DAS INFORMAÇÕES

Método

Para entender o termostado de Nóse-Hoover, primeiramente será mostrado o termostato de Nosé.

Este termostato atribui coordenadas generalizados adicionais $s$ e o seu momento conjugado $p_{s}$ ao banho térmico. O fator $s$ é definido como um fator de escala das velocidades, onde:

${\dot {\mathbf {v}}}=s{\dot {\mathbf {r}}}=s{\mathbf {p}}/m$

E também são definidas as energia potenciais e cinética associadas a $s$ como:

${\mathcal {U}}_{s}=(N_{f}+1)k_{B}Tln(s)$ e ${\mathcal {K}}_{s}={\frac {1}{2}}Q{\dot {s}}^{2}={\frac {p_{s}^{2}}{2Q}}$

onde $Q$ é entendido como a "inércia térmica", ele determina a escala do tempo da flutuação de temperatura.

O Lagrangiano do sistema extendida (consistente das partículas e do banho térmico) então é postulado como:

${\mathcal {L}}={\mathcal {K}}+{\mathcal {K}}_{s}-{\mathcal {U}}-{\mathcal {U}}_{s}=\sum _{i}{\frac {{\mathbf {p}}_{i}^{2}}{2m_{i}s^{2}}}+{\frac {p_{s}^{2}}{2Q}}-{\mathcal {U}}({\mathbf {r}})-(N_{f}+1)k_{B}Tln(s)$

Como não é explicitamente dependente do tempo:

${\mathcal {H}}_{N}={\mathcal {K}}+{\mathcal {K}}_{s}+{\mathcal {U}}+{\mathcal {U}}_{s}=\sum _{i}{\frac {{\mathbf {p}}_{i}^{2}}{2m_{i}s^{2}}}+{\frac {p_{s}^{2}}{2Q}}+{\mathcal {U}}({\mathbf {r}})+(N_{f}+1)k_{B}Tln(s)$

Como ${\mathcal {H}}_{N}$ se conserva, esse sistema é numericamente estável

Assim, as equações de movimento podem ser deduzidas:

${\mathbf {\dot {r}}}_{i}={\frac {\partial {\mathcal {H}}_{N}}{\partial {\mathbf {p}}_{i}}}={\mathbf {p}}_{i}/(m_{i}s^{2})$

${\mathbf {\dot {p}}}=-{\frac {\partial {\mathcal {H}}_{N}}{\partial {\mathbf {r}}_{i}}}={\mathbf {f}}_{i}$ onde $f$ é o número de graus de liberdade do sistema;

${\dot {s}}={\frac {\partial {\mathcal {H}}_{N}}{\partial p_{s}}}=p_{s}/Q$

${\dot {p}}_{s}=-{\frac {\partial {\mathcal {H}}_{N}}{\partial s}}=\sum _{i}{\frac {{\mathbf {p}}_{i}^{2}}{m_{i}s^{3}}}-(N_{f}+1)k_{B}T/s$

Assim, o termostato de Nose pode ser tratado como um sistema de partículas junto a um banho térmico como um ensemble NVE. Entretanto, neste sistema, $Q$ precisa ser determinado por tentativa e erro. Caso o valor escolhido seja muito pequeno, o sistema possuirá muitas oscilações, logo é necessário aumentar o valor de $Q$, porém caso o valor escolhido seja muito alto, o tempo para atingir equilíbrio térmico será demasiadamente longo. Outro problema do termostato de Nose é o fato de que, por as velocidades serem escaladas com o $s$, o tempo também será escalado com $s$, o que não acontece em sistemas reais e extendidos.

Para contornar esses problemas, Hoover utilizou uma parametrização diferente, sem o termo $s.$ . O parâmetro $s$ pode ser removido das equações reescrevendo , Assim, as equações de movimento do termostato de Nosé-Hoover são:

TERMINAR EQUAÇÕES

ADICIONAR EQUAÇÕES DAS CADEIAS DE NOSE HOOVER

Programas Utilizados

/*Simulação de DM de um fluido de Lennard-Jones com termostato Nose-Hoover Compile usando "gcc -o NVT_NH NVT_NH.c -lm -lgsl" */
/*********************************************/  
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
/*********CONSTANTES GLOBAIS******************/
#define rho 0.75 //Densidade 
#define V N*1./rho //Volume
#define N 256 //Numeros de particulas
#define T0 0.0
#define M 2 //Numero de correntes Nosé-Hoover
#define TEMPO_MAX 1000 //Tempo de simulação
/*********************MSD*********************/
#define BLOCKS     100    // num part para calc

/*********VARIAVEIS GLOBAIS******************/
double rx[N],ry[N],rz[N];
double vx[N],vy[N],vz[N];
double xi1[N],xi2[N],vxi1[N],vxi2[N];  
double fx[N],fy[N],fz[N];
double Q1[M],Q2[M];// Massa de resistencia do Termostato
double uk,e; // Energia potencial e cinetica
double T;
double l=3*N;// Variaveis reais - Frenkel --> L=3N 
double L = pow(V,1./3.);//Comprimento da aresta da caixa
double delt=0.001, delt2, delt4, delt8;
double To = 1.0;//Temperatura desejada
int tempo_p=1200; // Tempo do pulo de temperatura
double Temp_p=1.0;// Temperatura do pulo
double pe, et, et0,r2;
/********************MSD*********************/
double msd[TEMPO_MAX/BLOCKS], msdx[TEMPO_MAX/BLOCKS],
  msdy[TEMPO_MAX/BLOCKS],msdz[TEMPO_MAX/BLOCKS];
double vacf[TEMPO_MAX/BLOCKS],vacfx[TEMPO_MAX/BLOCKS],
  vacfy[TEMPO_MAX/BLOCKS],vacfz[TEMPO_MAX/BLOCKS];
int xpbc[N],ypbc[N],zpbc[N],MSD_BLOCK_SIZE;
double rx0[N],ry0[N],rz0[N],vx0[N],vy0[N],vz0[N];
/********************************************/



gsl_rng * r;// Gerador global
/****INCIAR AS POSIÇÕES E VELOCIDADES********/
void init ( void) 
/*******************************************/
{
int i;
double cmvx=0.0,cmvy=0.0,cmvz=0.0;
double fac,mu=1.0;
double dx,dy,dz;
int M1,n1;
unsigned long int Seed = 23410981;

gsl_rng_set(r,Seed);
uk=0; 
n1=ceil(pow(N,1./3.));
M1=pow(n1,3);

dx=pow(V/M1,(1./3.));
dy=pow(V/M1,(1./3.));
dz=pow(V/M1,(1./3.));
  
for (i=0;i<N;i++)
 {
rx[i] = (i%n1)*dx;
ry[i] = ((i/n1)%n1)*dy;
rz[i] = ((i/(n1*n1))%n1)*dz;
// Velocidade aleatoria de distribuição normal 

vx[i]=gsl_ran_gaussian(r,1.0);
vy[i]=gsl_ran_gaussian(r,1.0);
vz[i]=gsl_ran_gaussian(r,1.0);

cmvx+=vx[i];
cmvy+=vy[i];
cmvz+=vz[i];

vx[i]-=cmvx/N;
vy[i]-=cmvy/N;
vz[i]-=cmvz/N;
   
uk+=(vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i])/2;// Energia cinética média
 }
/********************************************/
if (T0>0.0)
 {

T = (2/3)*(uk/N);// Temperatura do sistema
fac=sqrt(T0/T);
    
for (i=0;i<N;i++)
  {
  
vx[i]*=fac;
vy[i]*=fac;
vz[i]*=fac;
      
uk+= (vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i])/2;
  }
  
 }    

return;  
}
/*********CALCULAR FORÇAS******************/
double force(void)
/*******************************************/ 
{
int i,j;
double dx, dy, dz,r2i, r6i;
double ff,rc,rc2;
double ecut;

e = 0.0;
rc = L/2;
rc2 = rc*rc;
ecut = 4*(pow(rc2,-6)-pow(rc2,-3));

for (i=0;i<N;i++)
  {
  
fx[i]=0.0;
fy[i]=0.0;
fz[i]=0.0;
  }
   

for (i=0;i<(N-1);i++)
{
for (j=i+1;j<N;j++)
 {
 
dx  = (rx[i]-rx[j]);
dy  = (ry[i]-ry[j]);
dz  = (rz[i]-rz[j]);
	
dx=dx-rint(dx/L)*L;
dy=dy-rint(dy/L)*L;
dz=dz-rint(dz/L)*L;
	
r2 = dx*dx + dy*dy + dz*dz;
	
if (r2<rc2) 
{
	
r2i = 1.0/r2;	
r6i  = r2i*r2i*r2i;
ff = 48.0*r6i*r2i*(r6i-0.5);
	  
fx[i] += dx*ff;
fx[j] -= dx*ff;
	  
fy[i] += dy*ff;
fy[j] -= dy*ff;
	  
fz[i] += dz*ff;
fz[j] -= dz*ff;
	   
e+= 4.0*r6i*(r6i - 1.0) - ecut;//Energia potencial
	   

}
   }
      }
      
 return e;
}
/*******************************************/
void aplica_PBC_msd(void) 
/*******************************************/
{

  int i;
  //
  for(i=0;i<N;i++) {
    // PBC - x
    if(rx[i]>L) {
      rx[i] = rx[i] - L;
      xpbc[i]++;
    }
    else if(rx[i]<0) {
      rx[i] = rx[i] + L;
      xpbc[i]--;
    }
    // y cord
    // PBC - Y
    if(ry[i]>L) {
      ry[i] = ry[i] - L;
      ypbc[i]++;
    }
    else if(ry[i]<0) {
      ry[i] = ry[i] + L;
      ypbc[i]--;
    }
    // z cord
    // PBC - Z
    if(rz[i]>L) {
      rz[i] = rz[i] - L;
      zpbc[i]++;
    }
    else if(rz[i]<0) {
      rz[i] = rz[i] + L;
      zpbc[i]--;
     }
    
  }
  //
  return;
}
/*******************************************/
void measure_msd(long int steps)
/*******************************************/
{

  int i,local_step;
  double mmsdx,mmsdy,mmsdz;
  double mvacfx,mvacfy,mvacfz;
  double xreal,yreal,zreal;
  //
  // calculo o tempo dentro de um bloco
  local_step = steps%MSD_BLOCK_SIZE;
  //
  // testo se tenho que atualizar as pos e vel para
  // novos r0 e v0
  if(local_step==0)
   {
    // redefino os pontosiniciais para msd
    for(i=0;i<N;i++) {
      rx0[i] = rx[i] + xpbc[i]*L;
      ry0[i] = ry[i] + ypbc[i]*L;
      rz0[i] = rz[i] + zpbc[i]*L;
      vx0[i] = vx[i];
      vy0[i] = vy[i];
      vz0[i] = vz[i];
    }
  //  printf("aqui: atualizei conf MSD -- step: %ld local: %d  msdblock: %d\n",steps,local_step,MSD_BLOCK_SIZE);
  }
  // 
  //
  mmsdx = 0.0; mvacfx = 0.0;
  mmsdy = 0.0; mvacfy = 0.0;
  mmsdz = 0.0; mvacfz = 0.0;
  //
  for(i=0;i<N;i++)
   {
    //
    xreal = rx[i] + xpbc[i]*L;
    yreal = ry[i] + ypbc[i]*L;
    zreal = rz[i] + zpbc[i]*L;
    
    //
    //mmsdx = mmsdx + (rx[i]-rx0[i])*(rx[i]-rx0[i]);
    //mmsdy = mmsdy + (ry[i]-ry0[i])*(ry[i]-ry0[i]);
    //mmsdz = mmsdz + (rz[i]-rz0[i])*(rz[i]-rz0[i]);
    
    mmsdx = mmsdx + (xreal-rx0[i])*(xreal-rx0[i]);
    mmsdy = mmsdy + (yreal-ry0[i])*(yreal-ry0[i]);
    mmsdz = mmsdz + (zreal-rz0[i])*(zreal-rz0[i]);
    
    mvacfx = mvacfx + vx[i]*vx0[i];
    mvacfy = mvacfy + vy[i]*vy0[i];
    mvacfz = mvacfz + vz[i]*vz0[i];
  }
  mmsdx = mmsdx/(1.0*N);
  mmsdy = mmsdy/(1.0*N);
  mmsdz = mmsdz/(1.0*N);
  
  msdx[local_step] = msdx[local_step] + mmsdx;
  msdy[local_step] = msdy[local_step] + mmsdy;
  msdz[local_step] = msdz[local_step] + mmsdz;  
    
  msd[local_step] = msd[local_step] + mmsdx + mmsdy + mmsdz;
  //
  mvacfx = mvacfx/(1.0*N);
  mvacfy = mvacfy/(1.0*N);
  mvacfz = mvacfz/(1.0*N);
  
  vacfx[local_step] = vacfx[local_step] + mvacfx;
  vacfy[local_step] = vacfy[local_step] + mvacfy;  
  vacfz[local_step] = vacfz[local_step] + mvacfz;  
  vacf[local_step] = vacf[local_step] + mvacfx + mvacfy +mvacfz;
  
  //printf("aqui: %ld %d %f %f %f %f\n",steps,local_step,msd,msdx,msdy);
  //
  return;
}
/*******************************************/
void inicializa_msd(void) 
/*******************************************/
{

int i,block_size;

  // guardo os vetores das posicoes iniciais
for(i=0;i<N;i++)
 {
    rx0[i] = rx[i];
    ry0[i] = ry[i];
    rz0[i] = rz[i];
    vx0[i] = vx[i];
    vy0[i] = vy[i];
    vz0[i] = vz[i];
    xpbc[i] = 0;
    ypbc[i] = 0;  
    zpbc[i] = 0;  
  }

  block_size = TEMPO_MAX/BLOCKS;
  MSD_BLOCK_SIZE = block_size;
  //printf("\n\n m: %d  b: %d   bs: %d\n\n",MEASURES,BLOCKS,block_size);

for(i=0;i<block_size;i++)
 {
    msd[i] = 0.0;
    msdx[i] = 0.0;
    msdy[i] = 0.0;
    msdz[i] = 0.0;
    vacf[i] = 0.0;
    vacfx[i] = 0.0;
    vacfy[i] = 0.0;
    vacfz[i] = 0.0;
  }
  
  return;
}
/*******************************************/
void imprime_msd(void)
/*******************************************/
 {

FILE *fpmsd;
double f,diff;
int i;


fpmsd = fopen("fmsd.dat","w");


fprintf(fpmsd,"#part msd vcaf diff\n");
f = 1.0/(1.0*BLOCKS);
diff = 0.0;

for(i=0;i<MSD_BLOCK_SIZE;i++)
{
    
    fprintf(fpmsd,"%f\t%f\t%f\t%f\n",i*delt,f*msd[i],f*vacf[i],diff);
    

diff = diff + f*vacf[i]*delt;
  
  
  }
fprintf(fpmsd,"#diff = %f\n",diff);


   fclose(fpmsd);
 
  return;
}
/*******************************************/
void chain ( double T) 
/*******************************************/
{

double G1, G2, s;
int i;

Q1[0] = 0.1;
Q2[1] = 0.1;

delt2 = delt/2;
delt4 = delt/4;
delt8 = delt/8;

G2= (Q1[0]*vxi1[0]*vxi1[0]-T);
vxi2[1]= vxi2[1] + G2*delt4;
vxi1[0]=vxi1[0]*exp(-vxi2[1]*delt8);

G1=(2*uk-l*T)/Q1[0];
vxi1[0]=vxi1[0] + G1*delt4;
vxi1[0]=vxi1[0]*exp(-vxi2[1]*delt8); 
 
xi1[0]=xi2[0] + vxi1[0]*delt2;
xi2[1]=xi2[1] + vxi2[1]*delt2;

s=exp(-vxi1[0]*delt2);

for (i=0;i<N;i++)
 {
vx[i]=s*vx[i]; 
vy[i]=s*vy[i]; 
vz[i]=s*vz[i]; 

 }
uk = uk*s*s;

vxi1[0]=vxi1[0]*exp(-vxi2[1]*delt8);
G1=(2*uk-l*T)/Q1[0];
vxi1[0]=vxi1[0] + G1*delt4;

vxi1[0]=vxi1[0]*exp(-vxi2[1]*delt8);
G2=(Q1[0]*vxi1[0]*vxi1[0]-T)/Q2[1];
vxi2[1]=vxi2[1] + G2*delt4;

return;
}
/*******************************************/
int main ( void ) 
/*******************************************/
{
int i,t;
const gsl_rng_type * A;

gsl_rng_env_setup();
A = gsl_rng_default;
r = gsl_rng_alloc (A);

init();
/*******************************************/
for (t=0;t<TEMPO_MAX;t++)
   {
/*******************************************/
if (t==tempo_p)
    
To = Temp_p;

chain(To);
/******************************************/
uk = 0.0;
    
    
    
for (i=0;i<N;i++)
   {
rx[i]+=vx[i]*delt2;
ry[i]+=vy[i]*delt2;
rz[i]+=vz[i]*delt2;
  
   }
/************Calculo das Forças**************/
pe = force();
      
/********************************************/
for (i=0;i<N;i++)
  {
  
vx[i]+=delt*fx[i];
vy[i]+=delt*fy[i];
vz[i]+=delt*fz[i];
      
rx[i]+=vx[i]*delt2;
ry[i]+=vy[i]*delt2;
rz[i]+=vz[i]*delt2;
      
uk+=(vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i])/2;

  }

chain(To);

T = (2/3.)*(uk/N);
et = pe + uk;

/*Tempo,Potencial,Cinetica,Total,Temperatura*/    
printf("%d\t%lf\t%lf\t%lf\t%lf\n", t,pe,uk,et,T);
/********************************************/  
  }
  
return 0;  
}