Simulação - Lennard-Jones

De Física Computacional
Ir para navegação Ir para pesquisar

!--------------------------------------------------------------------!

      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