Simulação - Lennard-Jones
!--------------------------------------------------------------------!
MODULE globais IMPLICIT NONE INTEGER :: Lx,Ly,N,Nf,nhis REAL(8) :: pi,eps,sigma,sigma2,rho END MODULE
!--------------------------------------------------------------------
PROGRAM lennard_jones USE globais IMPLICIT NONE INTEGER :: i,t,k,j,m,ngr INTEGER :: q,semente,seed(8) INTEGER,ALLOCATABLE :: hist_g(:) REAL(4) :: t_inicial,t_final,tempo REAL(8) :: d,r,dt,tmax,p,delg,u REAL(8),ALLOCATABLE :: x(:),y(:),Vx(:),Vy(:),Fx(:),Fy(:),dist(:),g(:) REAL(8) :: e_cin,e_pot,Temp,pressao,sum_Temp,sum_pressao
!
CALL CPU_time(t_inicial)
!-------------------------------------------------------------------- ! Semente para números aleatórios!
semente = 1234 DO q=0,7 seed(q+1) = semente + q*50 END DO CALL random_seed(put = seed)
!--------------------------------------------------------------------
pi = 4.0*atan(1.0) Lx = 25 Ly = 25 N = 144 !Número de partículas! Nf = 2*N-2 !Nuúmero de graus de liberdade! nhis = 200 !Número de bins! d = 2.0 !Distância entre as partículas! r = d*0.5 !Raio das partículas! eps = 1.0 sigma = 1.0 sigma2 = sigma*sigma rho = 1.0*N/(Lx*Ly) !Densidade numérica! Temp = 1.0 !Temperatura desejada, usada para! !reescalar as velocidades!
!-------------------------------------------------------------------- ! Aloca os vetores
ALLOCATE(x(N),y(N),Vx(N),Vy(N),Fx(N),Fy(N)) ALLOCATE(g(nhis),dist(nhis),hist_g(nhis))
!-------------------------------------------------------------------- ! Zera variáveis
Fx(:) = 0.0 Fy(:) = 0.0 x(:) = 0.0 y(:) = 0.0 Vx(:) = 0.0 Vy(:) = 0.0 g(:) = 0.0 hist_g(:) = 0 dist(:) = 0.0 sum_Temp = 0.0 sum_pressao = 0.0
!-------------------------------------------------------------------- ! Inicializa posições, velocidades, forças e fdr
CALL inicializacao(x,y,Vx,Vy,r,Temp,e_cin) CALL posicoes(x,y,r) CALL forcas(x,y,Fx,Fy,e_pot,p) CALL fdr(0,x,y,ngr,delg,dist,g)
!-------------------------------------------------------------------- ! Loop para integrar as equações de movimento
dt = 0.01 tmax = 100.0/dt OPEN(10,FILE='dados_particulas.dat',ACTION='write',STATUS='unknown') OPEN(100,FILE='dados_termo.dat',ACTION='write',STATUS='unknown') OPEN(200,FILE='dados_fdr.dat',ACTION='write',STATUS='unknown') DO t=0,int(tmax) Temp = e_cin/Nf pressao = rho*Temp + 0.5/(Lx*Ly)*p
! Retira o valor inicial da T e p do cálculo da média, pois é muito diferente
IF (t == 0) THEN sum_Temp = -Temp sum_pressao = -pressao GO TO 20 ELSE GO TO 20 END IF
20 sum_Temp = sum_Temp + Temp
sum_pressao = sum_pressao + pressao
! Escreve os dados nos arquivos
DO i=1,N WRITE(10,200) t,i,x(i),y(i),Vx(i),Vy(i),Fx(i),Fy(i) END DO WRITE(10,*) WRITE(10,*)
200 FORMAT(i5,2x,i3,2x,f10.5,2x,f10.5,2x,f10.5,2x,f10.5,2x,f10.5,2x,f10.5)
WRITE(100,300) t,e_cin,e_pot,Temp,pressao
300 FORMAT(i5,2x,f10.5,2x,f10.5,2x,f10.5,2x,f10.5) ! Algoritmo Velocity-Verlet
DO k=1,10 CALL verlet(x,y,Vx,Vy,Fx,Fy,dt,e_cin,e_pot,p) END DO CALL fdr(1,x,y,ngr,delg,dist,g) END DO CALL fdr(2,x,y,ngr,delg,dist,g)
! Imprime a fdr
DO m=1,nhis WRITE(200,400) dist(m),g(m) END DO
400 FORMAT(f10.5,2x,f10.5)
CLOSE(10) CLOSE(100) CLOSE(200)
! Calcula a temperatura média e a pressão média
sum_Temp = sum_Temp/tmax sum_pressao = sum_pressao/tmax WRITE(*,600) 'Densidade:',rho WRITE(*,700) 'Temperatura média:',sum_Temp WRITE(*,800) 'Pressão média:',sum_pressao
600 FORMAT(a11,1x,f7.4) 700 FORMAT(a19,1x,f7.4) 800 FORMAT(a16,1x,f7.4) !
CALL CPU_time(t_final) tempo = t_final - t_inicial WRITE(*,900) 'Tempo de execução:',tempo,'s'
900 FORMAT(a20,1x,f6.2,1x,a1) !
END PROGRAM
!--------------------------------------------------------------------- ! Inicialização
SUBROUTINE inicializacao(rx,ry,vx,vy,r,Temp,e_cin) USE globais IMPLICIT NONE INTEGER :: num,i REAL(8) :: d,r,x,y,Temp,e_cin,rand1,rand2 REAL(8) :: sumv2,sum_vx,sum_vy,fs REAL(8) :: rx(N),ry(N),vx(N),vy(N)
! Inicializa as posições
num = 1 x = r y = r DO WHILE (y .lt. Ly-r+0.05) DO WHILE (x .lt. Lx-r+0.05) IF (num .le. N) THEN rx(num) = x ry(num) = y num = num + 1 END IF x = x + 2*r END DO x = r y = y + 2*r END DO
! Inicializa as velocidades
sumv2 = 0.0 sum_vx = 0.0 ; sum_vy = 0.0 DO i=1,N CALL random_number(rand1) CALL random_number(rand2) vx(i) = rand1 - 0.5 vy(i) = rand2 - 0.5 sumv2 = sumv2 + vx(i)*vx(i) + vy(i)*vy(i) !energia cinética! sum_vx = sum_vx + vx(i) !veloc x do CM! sum_vy = sum_vy + vy(i) !veloc y do CM! END DO e_cin = 0.5*sumv2 sumv2 = sumv2/N fs = sqrt(2.0*Temp/sumv2) !Reescala as velocidades! sum_vx = sum_vx/N !veloc x do CM! sum_vy = sum_vy/N !veloc y do CM! DO i=1,N vx(i) = (vx(i) - sum_vx)*fs vy(i) = (vy(i) - sum_vy)*fs END DO END SUBROUTINE
!---------------------------------------------------------------------- ! Implime a configuração inicial do sistema
SUBROUTINE posicoes(x,y,r) USE globais IMPLICIT NONE INTEGER :: k REAL(8) :: r REAL(8) :: x(N),y(N) OPEN(9,FILE='posicoes.dat',ACTION='write',STATUS='unknown') DO k=1,N WRITE(9,*) x(k),y(k),r END DO CLOSE(9) END SUBROUTINE
!---------------------------------------------------------------------- ! Calcula as forças em um instante de tempo
SUBROUTINE forcas(x,y,fx,fy,e_pot,p) USE globais IMPLICIT NONE INTEGER :: i,j REAL(8) :: a,r2,dx,dy,e_pot,p REAL(8) :: x(N),y(N),fx(N),fy(N) fx(:) = 0.0 fy(:) = 0.0 e_pot = 0.0 p = 0.0 DO i=1,N-1 DO j=i+1,N
! Convenção da imagem mínima
dx = x(i)-x(j) dy = y(i)-y(j) dx = dx - Lx*nint(dx/Lx) dy = dy - Ly*nint(dy/Ly) r2 = dx*dx + dy*dy a = 48*eps/sigma2*(sigma2/r2)**4*((sigma2/r2)**3 - 0.5) fx(i) = fx(i) + a*dx fy(i) = fy(i) + a*dy fx(j) = fx(j) - a*dx fy(j) = fy(j) - a*dy e_pot = e_pot + 4.0*eps*(sigma2/r2)**3*((sigma2/r2)**3-1.0) !Energia potencial! p = p + a*(dx*dx + dy*dy) !Pressão! END DO END DO END SUBROUTINE
!--------------------------------------------------------------------- ! Implementa a condição de contorno
SUBROUTINE condicao_contorno(x,y) USE globais IMPLICIT NONE INTEGER :: i REAL(8) :: x(N),y(N) DO i=1,N IF (x(i) .gt. Lx) THEN x(i) = x(i) - Lx ELSE IF (x(i) .lt. 0) THEN x(i) = x(i) + Lx END IF IF (y(i) .gt. Ly) THEN y(i) = y(i) - Ly ELSE IF (y(i) .lt. 0) THEN y(i) = y(i) + Ly END IF END DO END SUBROUTINE
!---------------------------------------------------------------------- ! ALgoritmo Velocity-Verlet
SUBROUTINE verlet(rx,ry,vx,vy,fx,fy,deltat,e_cin,e_pot,p) USE globais IMPLICIT NONE INTEGER :: i REAL(8) :: deltat,e_cin,e_pot,p REAL(8) :: rx(N),ry(N),vx(N),vy(N),fx(N),fy(N),g(nhis) e_cin = 0.0 DO i=1,N vx(i) = vx(i) + 0.5*fx(i)*deltat vy(i) = vy(i) + 0.5*fy(i)*deltat rx(i) = rx(i) + vx(i)*deltat ry(i) = ry(i) + vy(i)*deltat END DO CALL condicao_contorno(rx,ry) CALL forcas(rx,ry,fx,fy,e_pot,p) DO i=1,N vx(i) = vx(i) + 0.5*fx(i)*deltat vy(i) = vy(i) + 0.5*fy(i)*deltat e_cin = e_cin + vx(i)*vx(i) + vy(i)*vy(i) END DO e_cin = 0.5*e_cin !Energia cinética! END SUBROUTINE
!--------------------------------------------------------------------- ! Calcula a função distribuição radial
SUBROUTINE fdr(switch,x,y,ngr,delg,dist,g) USE globais IMPLICIT NONE INTEGER :: i,j,k,ig,ngr,switch ! INTEGER :: hist_g(nhis) REAL(8) :: dx,dy,r,r2,delg,vb,nid REAL(8) :: x(N),y(N),dist(nhis),g(nhis) IF (switch == 0) THEN !Inicializa! ngr = 0 delg = Lx/(2.0*nhis) !Tamanho do bin! dist(:) = 0.0 ! hist_g(:) = 0 g(:) = 0.0 ELSE IF (switch == 1) THEN !amostra! ngr = ngr + 1 DO i=1,N-1 DO j=i+1,N ! Convenção da imagem mínima dx = x(i)-x(j) dy = y(i)-y(j) dx = dx - Lx*nint(dx/Lx) dy = dy - Ly*nint(dy/Ly) r2 = dx*dx + dy*dy r = sqrt(r2) ig = int(r/delg) g(ig) = g(ig) + 2 END DO END DO ELSE IF (switch == 2) THEN !calcula g(r)! DO k=1,nhis dist(k) = delg*(k + 0.5) vb = (2*k+1)*delg*delg nid = pi*vb*rho g(k) = g(k)/(ngr*N*nid) !Normaliza g(r)! END DO END IF END SUBROUTINE