Hello, I am currently making a simulation of a gas in the NVT ensamble. To do so I am using de Monte-Carlo method. Because I have to do a lot of metropolis steps, I was looking for a way to optimize the Lennard-Jones potential so it still computes the total energy of the system but maybe only taking into account that I move each iteration only one particle. I’ll add the current code without any optimization. I have the positions to initialize in equilibrium as well as the potential energy. Does anyone have any ideas on how I can optimize this code?
“”“”"
module algoritmo_metropolis
use def_precision
use variables_comunes ! Usa las variables comunes necesarias para la subrutina
use random_module
use subrutinas
contains
subroutine metropolis(np, rx, ry, rz, epot, dfiv, d2fiv, iflag, temp)
implicit none
integer, intent(in) :: np
real(kind=doblep), intent(inout) :: rx(np), ry(np), rz(np)
real(kind=doblep) :: old_rx, old_ry, old_rz
real(kind=doblep), intent(out) :: epot, dfiv, d2fiv
integer, intent(inout) :: iflag
real(kind=doblep), intent(in) :: temp
real(kind=doblep) :: old_epot, new_epot, deltaE, factor, pli, xnp
integer :: i, part
real(kind=doblep) :: random_val
!Valores necesarios para el calculo
xnp = dble(npmax)
pli=1.d00/pl
vol = pl*pl*pl
rc2 = rc*rc
factor = pi*xnp*xnp/(vol*rc*rc*rc)
!Correcciones de la energía relacionadas con el potencial Lennard-Jones
corr_ener=8.d00*factor*(1.d00/(3.d00*rc**6)-1.d00)/3.d00
corr_suma_rvp = 16.d00*factor*(-2.d00/(3.d00*rc**6)+1.d00)
corr_sum_r2vpp = 16.d00 * factor*(26.d00/(3.d00*rc**6)-7.d00)
! Seleccionar una partícula aleatoria
part= int( generator_random(iflag)* 500.d0)+1 !Generamos un numero entero aleatorio entre 1 y 500
! Guardar las posiciones originales de la partícula
old_rx = rx(part)
old_ry = ry(part)
old_rz = rz(part)
old_epot = epot
!Movemos la particula aleatoriamente
rx(part) = rx(part)+(2.d00*generator_random(iflag)-1.d00)*0.08d0
ry(part) = ry(part)+(2.d00*generator_random(iflag)-1.d00)*0.08d0
rz(part) = rz(part)+(2.d00*generator_random(iflag)-1.d00)*0.08d0
!Calculamos la nueva energía
call potlj(np, rx, ry, rz,epot, dfiv, d2fiv) ! Calcular la energía actual
new_epot = epot
! Calcular la diferencia de energía
deltaE = new_epot - old_epot
! Aceptar o rechazar el cambio según la regla de Metropolis
if (deltaE < 0.d0) then
! Aceptamos el cambio si la energía disminuye siempre -> el sistema va hacia energías más bajas
epot = new_epot
else
! Aceptamos el cambio con una probabilidad dependiendo de la temperatura
if (exp(-deltaE/temp) >= generator_random(iflag)) then
epot = new_epot
else
! Revertimos el cambio si se rechaza
rx(part) = old_rx
ry(part) = old_ry
rz(part) = old_rz
epot = old_epot
end if
end if
end subroutine metropolis
end module algoritmo_metropolis
“”“”
“”“”"
module subrutinas
use def_precision
use variables_comunes ! Usa las variables comunes necesarias para la subrutina
contains
!SUBRUTINA DEL POTENCIAL DE LENNARD-JONES
subroutine potlj(np, rx, ry, rz,epot, dfiv, d2fiv)
implicit none
! Esta subrutina calcula la energía potencial total del sistema (epot), las fuerzas entre partículas
! (almacenadas en las aceleraciones ax, ay, az), y las derivadas primera (dfiv) y segunda (d2fiv) del
! potencial de Lennard-Jones respecto al volumen.
integer (kind=entero), intent(in) :: np
real (kind=doblep), dimension(:), intent(in) :: rx, ry, rz
real (kind=doblep), intent(out) :: dfiv, d2fiv
real(kind=doblep), intent(out) :: epot
! Variables locales:
integer (kind=entero) :: i, j ! Índices para recorrer pares de partículas
real (kind=doblep) :: rrx, rry, rrz, rijx, rijy, rijz ! Variables para almacenar posiciones relativas
real (kind=doblep) :: dis2, fmod, a2, a6, a12, aux1 ! Variables para cálculos de distancias, fuerzas, etc.
! Inicialización de acumuladores:
epot = 0.d00
dfiv = 0.d00
d2fiv = 0.d00
rc2=rc*rc
pli=1.d0/pl
! Bucle doble para recorrer todos los pares de partículas i-j
do i = 1, np-1
rrx = rx(i)
rry = ry(i)
rrz = rz(i)
do j = i + 1, np
rijx = rrx - rx(j)
rijy = rry - ry(j)
rijz = rrz - rz(j)
rijx = rijx - pl * dnint(rijx * pli)
rijy = rijy - pl * dnint(rijy * pli)
rijz = rijz - pl * dnint(rijz * pli)
dis2 = (rijx * rijx) + (rijy * rijy) + (rijz * rijz)
if (dis2 <= rc2) then
a2 = 1.d00 / dis2
a6 = a2 * a2 * a2
a12 = a6 * a6
epot = epot + (a12 - a6)
aux1 = (2.d00 * a12) - a6
dfiv = dfiv - aux1
d2fiv = d2fiv + 26.d00 * a12 - 7.d00 * a6
end if
end do
end do
epot = 4.d00 * epot + corr_ener
dfiv = 24.d00 * dfiv + corr_suma_rvp
d2fiv = 24.d00 * d2fiv + corr_sum_r2vpp
d2fiv = (d2fiv - 2.d00 * dfiv) / (9.d00 * vol * vol)
dfiv = dfiv / (3.d00 * vol)
return
end subroutine potlj
end module subrutinas
“”“”"