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