Mudanças entre as edições de "Termostato de Nosé-Hoover"

De Física Computacional
Ir para: navegação, pesquisa
(Termostato de Nosé-Hoover)
Linha 9: Linha 9:
  
 
== Método ==
 
== Método ==
 +
 +
=== Termostato de Nose ===
  
 
Para entender o termostado de Nóse-Hoover, primeiramente será mostrado o termostato de Nosé<ref name=nose> NOSÉ, Shuichi, '''A molecular dynamics method for simulations in the canonical ensemble''', Molecular Physics, 1984, Vol. 52, No. 2, 255-268 </ref>.  
 
Para entender o termostado de Nóse-Hoover, primeiramente será mostrado o termostato de Nosé<ref name=nose> NOSÉ, Shuichi, '''A molecular dynamics method for simulations in the canonical ensemble''', Molecular Physics, 1984, Vol. 52, No. 2, 255-268 </ref>.  
Linha 63: Linha 65:
 
#include <stdio.h>
 
#include <stdio.h>
 
#include <stdlib.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;   
 
return 0;   

Edição das 19h17min de 25 de maio de 2021

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[1]. 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 E TAMBÉM INFOS SOBRE LJ

Método

Termostato de Nose

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

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

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

e

onde é 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:

Como não é explicitamente dependente do tempo:

Como se conserva, esse sistema é numericamente estável [3]

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

onde é o número de graus de liberdade do sistema;

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, 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 , 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 , o tempo também será escalado com , o que não acontece em sistemas reais e extendidos. [3]

Para contornar esses problemas, Hoover utilizou uma parametrização diferente, sem o termo [4]. O parâmetro 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

Resultados

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>
  
return 0;  
}

Referências

  1. https://www2.ph.ed.ac.uk/~dmarendu/MVP/MVP03.pdf
  2. NOSÉ, Shuichi, A molecular dynamics method for simulations in the canonical ensemble, Molecular Physics, 1984, Vol. 52, No. 2, 255-268
  3. 3,0 3,1 http://www.courses.physics.helsinki.fi/fys/moldyn/lectures/L5.pdf
  4. William G. Hoover, Canonical Dynamics: Equilibrium phase-space distributions, Physical Review A, 1985, Vol. 31, No. 3.