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