program demonm * *------------------------------------------------------------------------------ * * Este programa ilustra o algoritimo "demon" para o ensemble microcanonico * Sistema considerado: gas ideal * * Autor: Sergio R. Souza * e-mail: srsouza@if.ufrj.br * *------------------------------------------------------------------------------ implicit none * logical ter_sis * * * common/infsim/ --> guarda informacao sobre parametros da simulacao * integer*4 n_mc,n_term real*8 e_part,dv_max common/infsim/e_part,dv_max,n_mc,n_term * *------------------------------------------------------------------------------ * * ---------------- * Entrada de dados * ---------------- * call ent_dados * * ------------------------ * inicializacao do sistema * ------------------------ * call inicia * * * --------------------- * termalizacao do demon * --------------------- * open(file='term.dat',unit=10,status='new') * ter_sis=.true. call mcsamp(ter_sis) * close(10) * * * ------------------ * calculo da energia * ------------------ * if(n_mc.gt.0)then ter_sis=.false. call mcsamp(ter_sis) call saida end if * end *------------------------------------------------------------------------------ * * Esta subrotina le dados de entrada do programa * *------------------------------------------------------------------------------ * ==================== subroutine ent_dados * ==================== implicit none * include 'demon.inp' * *------------------------------------------------------------------------------ * variaveis de saida *------------------------------------------------------------------------------ * * * common/infsim/ --> guarda informacao sobre parametros da simulacao * integer*4 n_mc,n_term real*8 e_part,dv_max common/infsim/e_part,dv_max,n_mc,n_term * *------------------------------------------------------------------------------ * open(file='demon.dat',unit=10,status='old') * read(10,*)e_part ! energia por particula read(10,*)dv_max ! maior variacao da velocidade read(10,*)n_mc ! numero de sorteios Monte Carlo read(10,*)n_term ! numero de passos para termalizacao * close(10) * end *------------------------------------------------------------------------------ * * Esta subrotina inicializa o sistema * *------------------------------------------------------------------------------ * ================= subroutine inicia * ================= implicit none * include 'demon.inp' * *------------------------------------------------------------------------------ * variaveis de entrada *------------------------------------------------------------------------------ * * * common/infsim/ --> guarda informacao sobre parametros da simulacao * integer*4 n_mc,n_term real*8 e_part,dv_max common/infsim/e_part,dv_max,n_mc,n_term * *------------------------------------------------------------------------------ * variaveis de saida *------------------------------------------------------------------------------ * * * common/psiste/ -> guarda informacoes sobre as propriedades do sistema * real*8 vel(n_part),e_demon,e_sist,n_aceit common/psiste/vel,e_demon,e_sist,n_aceit * * * common/semente/ -> guarda semente para gerador de numeros aleatorios * integer seed common/semente/seed * * * common/caixas/ -> guarda informacoes sobre as caixas * real*8 de,n_e(0:n_box) common/caixas/de,n_e * *------------------------------------------------------------------------------ * variaveis locais *------------------------------------------------------------------------------ * integer*4 i real*8 e_tot * *------------------------------------------------------------------------------ * * * ----------------------------- * inicializacao das velocidades * ----------------------------- * do i=1,n_part vel(i)=sqrt(2.0*e_part) end do * * -------------------------------------------- * semente para o gerador de numeros aleatorios * -------------------------------------------- * seed=-1 * * --------------------------------------------- * caixas para o histograma com energia do demon * --------------------------------------------- * e_tot=n_part*e_part ! energia total de=e_tot/n_box ! tamanho das caixas * do i=0,n_box n_e(i)=0.d0 ! numero de eventos com energia entre end do ! de*i e de*(i+1) * end *------------------------------------------------------------------------------ * * * Esta funcao calcula grandezas fisicas relevantes utilizando o metodo * de Monte Carlo. * * .true. - apenas a termalizacao do demon e feita * ter_sis -> * .false. - calculo Monte Carlo e realmente feito * *------------------------------------------------------------------------------ * ========================== subroutine mcsamp(ter_sis) * ========================== implicit none * include 'demon.inp' * *------------------------------------------------------------------------------ * variaveis de entrada/saida *------------------------------------------------------------------------------ * logical ter_sis * * * common/infsim/ --> guarda informacao sobre parametros da simulacao * integer*4 n_mc,n_term real*8 e_part,dv_max common/infsim/e_part,dv_max,n_mc,n_term * * * common/psiste/ -> guarda informacoes sobre as propriedades do sistema * real*8 vel(n_part),e_demon,e_sist,n_aceit common/psiste/vel,e_demon,e_sist,n_aceit * * * common/semente/ -> guarda semente para gerador de numeros aleatorios * integer seed common/semente/seed * * * common/caixas/ -> guarda informacoes sobre as caixas * real*8 de,n_e(0:n_box) common/caixas/de,n_e * real ran3 ! funcao externa * *------------------------------------------------------------------------------ * variaveis locais *------------------------------------------------------------------------------ * save n_cont integer*4 i,j,n,i_part,n_cont,i_box real*8 e_dsoma,e_ssoma,v_test,dv,delta_e,n_ev * *------------------------------------------------------------------------------ * * * ------------- * inicializacao * ------------- * n_aceit=0.d0 * e_dsoma=0.d0 e_ssoma=0.d0 * if(ter_sis)then n=n_term ! apenas termalizar o sistema e_sist=n_part*e_part e_demon=0.d0 n_cont=0 else n=n_mc ! fazer simulacao end if * * --------- * simulacao * --------- * n_ev=0.d0 * do i=1,n ! loop em eventos do j=1,n_part ! loop em particulas n_ev=n_ev+1.d0 i_part=int(n_part*ran3(seed)+1) ! sorteio da particula * dv=dv_max*(2.d0*ran3(seed)-1.d0) ! sorteio do incremento em v * v_test=vel(i_part)+dv ! incrementando v delta_e=0.5*(v_test*v_test-vel(i_part)*vel(i_part)) * if(delta_e.le.e_demon)then vel(i_part)=vel(i_part)+dv e_demon=e_demon-delta_e e_sist=e_sist+delta_e n_aceit=n_aceit+1.d0 end if * e_dsoma=e_dsoma+e_demon e_ssoma=e_ssoma+e_sist * if(ter_sis)then n_cont=n_cont+1 if(n_cont.eq.500)then write(10,*)sngl(n_ev),sngl(e_dsoma/n_ev) n_cont=0 end if else i_box=int(e_demon/de) n_e(i_box)=n_e(i_box)+1.d0 end if end do end do * if(n_aceit.eq.0)then write(*,*) write(*,*)' e necessario aumentar o numero de passos' write(*,*) stop end if * e_sist=e_ssoma/n_ev ! energia media e_demon=e_dsoma/n_ev ! energia media * if(ter_sis)then n_aceit=0.d0 else e_sist=e_sist/n_part ! energia por particula end if * end *------------------------------------------------------------------------------ * * Esta subrotina faz a saida de dados * *------------------------------------------------------------------------------ * ================ subroutine saida * ================ implicit none * include 'demon.inp' * *------------------------------------------------------------------------------ * variaveis de entrada *------------------------------------------------------------------------------ * * * common/psiste/ -> guarda informacoes sobre as propriedades do sistema * real*8 vel(n_part),e_demon,e_sist,n_aceit common/psiste/vel,e_demon,e_sist,n_aceit * * * common/caixas/ -> guarda informacoes sobre as caixas * real*8 de,n_e(0:n_box) common/caixas/de,n_e * *------------------------------------------------------------------------------ * variaveis locais *------------------------------------------------------------------------------ * integer*4 i_box real*8 e_d * *------------------------------------------------------------------------------ * write(*,*)'< e_sist > /n_part = ',e_sist write(*,*)'< e_demon > = ',e_demon * open(file='hist.dat',unit=10,status='new') * do i_box=0,n_box-1 if(n_e(i_box).gt.0.5d0)then ! escrever apenas se > 0 e_d=(i_box+0.5d0)*de ! energia no meio da caixa write(10,*)sngl(e_d),sngl(n_e(i_box)) end if end do * close(10) * end *------------------------------------------------------------------------------ * * Esta funcao gera um numero aleatorio entre zero e um com distribuicao * uniforme. * * E necessario que idum seja menor ou igual a zero na primeira vez * que a funcao seja chamada para que variaveis internas sejam inicializadas * *------------------------------------------------------------------------------ FUNCTION ran3(idum) INTEGER idum INTEGER MBIG,MSEED,MZ C REAL MBIG,MSEED,MZ REAL ran3,FAC PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1./MBIG) C PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG) INTEGER i,iff,ii,inext,inextp,k INTEGER mj,mk,ma(55) C REAL mj,mk,ma(55) SAVE iff,inext,inextp,ma DATA iff /0/ if(idum.lt.0.or.iff.eq.0)then iff=1 mj=MSEED-iabs(idum) mj=mod(mj,MBIG) ma(55)=mj mk=1 do 11 i=1,54 ii=mod(21*i,55) ma(ii)=mk mk=mj-mk if(mk.lt.MZ)mk=mk+MBIG mj=ma(ii) 11 continue do 13 k=1,4 do 12 i=1,55 ma(i)=ma(i)-ma(1+mod(i+30,55)) if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG 12 continue 13 continue inext=0 inextp=31 idum=1 endif inext=inext+1 if(inext.eq.56)inext=1 inextp=inextp+1 if(inextp.eq.56)inextp=1 mj=ma(inext)-ma(inextp) if(mj.lt.MZ)mj=mj+MBIG ma(inext)=mj ran3=mj*FAC return END C (C) Copr. 1986-92 Numerical Recipes Software '(*+('.