Monte-Carlo Lennard-Jones Potential

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
“”“”"

1 Like

The code you posted shows no loops as far as I can tell. Is it possible to show more of your program? It is difficult to see what is taking the computational effort and you also do not show the calculation of the Lennard-Jones potential as such.

The obvious first step to optimize the LJ potential evaluation for Metropolis MC moves, is to note that as long as you stick to single-particle updates, there is no need to re-evaluate the total energy in full with the double loop. One way of achieving this is to note that the evaluation of the change in energy Delta-E required by the detailed balance algorithm requires only a single loop to compute the change.

To improve readability it would help if you enclosed your code between

$``` fortran

$```

For example

module algoritmo_metropolis
use def_precision
use variables_comunes ! Usa las variables comunes necesarias para la subrutina
use random_module
use subrutinas