ello! I am currently working on an assignment where I have to simulate an ideal gas using fortran 95 and the Lennard-Jones potential.
After programming a lattice of 500 particles (using a fcc) and shifting them a bit (to simulate a gas) I want to reach an equilibrium point using verlet to forward the time. My issue has encountered when for some reason, teh equilibration code doesn’t calculate properly the energies. I was expecting a small dent in the total energy that I was going to correct, but after fixing my energy at -640, I have that by applying this algorithm I get an evolution of the system at energy -300. I have searched and tried everything but can’t seem to get what the error is. I belive both the crea red and potljsubroutine calculate the enrgy correctly, so I believe now that the problem can be when I grab the different positions and velocities ? I am going to paste my current code and would really appreciate some help on this topic. Thank you!!!
PD: Sorry for mistakes in writting, english isn’t my first language.
Equilibration program:
program equilibrio
use def_precision
use variables_comunes
use verlet_mod
implicit none
integer :: np, kcuenta, kk, kpaso, ktotal
real(kind=doblep), allocatable :: rx(:), ry(:), rz(:), vx(:), vy(:), vz(:), ax(:), ay(:), az(:)
real(kind=doblep):: ecin, epot, etot, factor, dfiv, d2fiv, tiempo, xnp
! Abrir el archivo 'variables.dat' para lectura
open(unit=10, file='variables.dat')
! Leer datos
read(10, '(i8)') np
read(10, '(f12.6)') pl
read(10, '(f12.6)') rc
read(10, '(f12.6, 2x, f12.6)') vol, dens
read(10, '(f12.6, 2x, f12.6, 2x, f12.6)') etot, ecin, epot
close(10)
! Asignar tamaño dinámico a los arreglos basado en np
allocate(rx(np), ry(np), rz(np), vx(np), vy(np), vz(np), ax(np), ay(np), az(np))
! Comprobar si los datos leĂdos son correctos
print *, "EnergĂa total: ", etot, "EnergĂa cinĂ©tica: ", ecin, "EnergĂa potencial: ", epot
! Leer datos de posiciones, velocidades y aceleraciones
open(20, file='resultados.dat', form='unformatted')
read(20) rx, ry, rz, vx, vy, vz, ax, ay, az
close(20)
! InicializaciĂłn de variables de tiempo y factores
dt = 1.0d-04
dt12 = dt / 2.0d0
dt2 = dt * dt / 2.0d0
rc2 = rc * rc
xnp = dble(np)
factor = pi * xnp * xnp / (vol * rc * rc * rc)
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)
! Bucle principal de integraciĂłn
kpaso = 1
ktotal = 50
kcuenta = 0
do kk = 1, ktotal
call verlet(np, rx, ry, rz, vx, vy, vz, ax, ay, az, epot, ecin, dfiv, d2fiv)
etot = ecin + epot
! Grabar datos cada kpaso pasos
if (mod(kk, kpaso) == 0) then
kcuenta = kcuenta + 1
tiempo = dble((kcuenta - 1) * kpaso) * dt
open(unit=10, file='simulacion.dat', access='append')
write(10, 9001) tiempo, etot, ecin, epot
close(10)
end if
end do
write(*,*) 'grabados', kcuenta, 'pasos'
! Escribir la configuraciĂłn final
open(20, file='resultados.dat', form='unformatted')
write(20) rx, ry, rz, vx, vy, vz, ax, ay, az
close(20)
9001 format(1pe13.6, 2x, e13.6, 2x, e13.6, 2x, e13.6)
end program equilibrio
Verlet program:
module verlet_mod
use def_precision
use variables_comunes
use subrutinas
contains
subroutine verlet(np, rx, ry, rz, vx, vy, vz, ax, ay, az, epot, ecin, dfiv, d2fiv)
implicit none
!Algoritmo velocidad de Verlet
!Preparado para 500 particulas con el parametro npmax
!The Verlet algorithm is the numerical method we used to calculate updates to particle
!position based on the total force of each particle. The algorithm is generally used to
!integrate Newton’s equations of motion
!DeclaraciĂłn de variables:
!Variables de entrada:
! np : numero de particulas en el sistema
! rx, ry, rz : posiciones de las particulas
! vx, vy, vz : velocidades de las particulas
! ax, ay, az : aceleraciones de las particulas
! dfiv :: dervada primera del potencial respecto el volumen
! d2fiv :: derivada segunda del potencial respecto del volumen
integer(kind=entero), intent(in) :: np
real(kind=doblep), intent(inout) :: epot
real(kind=doblep), dimension(npmax), intent(inout):: rx, ry, rz, vx, vy, vz, ax, ay, az
real(kind=doblep), intent(out) :: dfiv, d2fiv, ecin
rx=rx+vx*dt+ax*dt2
ry=ry+vy*dt+ay*dt2
rz=rz+vz*dt+az*dt2
vx=vx+ax*dt12
vy=vy+ay*dt12
vz=vz+az*dt12
!print*, 'epot dins verlet', epot
call potlj(np, rx, ry, rz, ax, ay, az, epot, dfiv, d2fiv)
vx=vx+ax*dt12
vy=vy+ay*dt12
vz=vz+az*dt12
ecin=0.5d00*(sum(vx*vx)+sum(vy*vy)+sum(vz*vz))
return
end subroutine verlet
end module verlet_mod
The module where the subroutine of teh Lennard-Jones potential is kept:
module subrutinas
use def_precision
use variables_comunes ! Usa las variables comunes necesarias para la subrutina
implicit none
contains
!SUBRUTINA DEL POTENCIAL DE LENNARD-JONES
subroutine potlj(np, rx, ry, rz, ax, ay, az, 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), dimension(:), intent(out) :: ax, ay, az
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 ! Primera derivada inicializada a 0
d2fiv = 0.d00 ! Segunda derivada inicializada a 0
ax = 0.d00 ! Aceleraciones inicializadas a 0
ay = 0.d00
az = 0.d00
!print*, rc2
! Bucle doble para recorrer todos los pares de partĂculas i-j
do i = 1, np-1
rrx = rx(i) ! Coordenadas de la partĂcula i
rry = ry(i)
rrz = rz(i)
do j = i + 1, np
! Cálculo de las distancias relativas entre las partĂculas i y j
rijx = rrx - rx(j)
rijy = rry - ry(j)
rijz = rrz - rz(j)
! AplicaciĂłn de condiciones periĂłdicas de contorno para evitar efectos de borde:
! Se ajustan las distancias para que las partĂculas interactĂşen con sus imágenes más cercanas.
rijx = rijx - pl * dnint(rijx * pli)
rijy = rijy - pl * dnint(rijy * pli)
rijz = rijz - pl * dnint(rijz * pli)
! Cálculo de la distancia al cuadrado entre las partĂculas i y j
dis2 = (rijx * rijx) + (rijy * rijy) + (rijz * rijz)
! Solo calculamos la interacciĂłn si la distancia es menor o igual al radio de corte (rc)
if (dis2 <= rc2) then
a2 = 1.d00 / dis2 ! Inverso de la distancia al cuadrado
a6 = a2 * a2 * a2 ! a6 = (1/r^2)^3 = 1/r^6
a12 = a6 * a6 ! a12 = (1/r^6)^2 = 1/r^12
! Actualizamos la energĂa potencial usando el potencial de Lennard-Jones (epot):
epot = epot + (a12 - a6)
! Auxiliar para el cálculo de las derivadas y las fuerzas:
aux1 = (2.d00 * a12) + a6
! Suma de las derivadas del potencial respecto al volumen
dfiv = dfiv - aux1
d2fiv = d2fiv + 26.d00 * a12 - 7.d00 * a6
! Magnitud de la fuerza (fmod) entre las partĂculas i y j
fmod = aux1 * a2
! Actualizamos las aceleraciones de las partĂculas i y j
ax(i) = ax(i) + fmod * rijx
ay(i) = ay(i) + fmod * rijy
az(i) = az(i) + fmod * rijz
ax(j) = ax(j) - fmod * rijx
ay(j) = ay(j) - fmod * rijy
az(j) = az(j) - fmod * rijz
end if
end do
end do
! Una vez calculada la interacciĂłn entre todos los pares de partĂculas, aplicamos factores comunes:
epot = 4.d00 * epot + corr_ener ! Factor de correcciĂłn de energĂa potencial
dfiv = 24.d00 * dfiv + corr_suma_rvp ! CorrecciĂłn de la primera derivada
d2fiv = 24.d00 * d2fiv + corr_sum_r2vpp ! CorrecciĂłn de la segunda derivada
! Escalamos las aceleraciones calculadas:
ax = 24.d00 * ax
ay = 24.d00 * ay
az = 24.d00 * az
! Finalmente, calculamos las derivadas del potencial respecto al volumen:
d2fiv = (d2fiv - 2.d00 * dfiv) / (9.d00 * vol * vol)
dfiv = dfiv / (3.d00 * vol)
return
end subroutine potlj
end module subrutinas
Creation of the lattice (and where i store the different values of r, v, a)
PROGRAM crea_red
! Este programa genera un sistema de particulas en una red cubica, aplica un random shift a sus posiciones,
! asigna unas velocidades aleatorias y computa su energia. Se escriben los resultados en un fichero txt.
use def_precision ! Modulo con las definiciones de precisiĂłn
use variables_comunes ! Modulo con variables comunas
use subrutinas ! Modulo que contiene subrutinas: en este caso contiene la subrutina del potencial de Lennard-Jones
use random_module ! Modulo para generar numeros aleatorios
implicit none
! DEFINICION DE LAS VARIABLES USADAS:
integer(kind=entero):: np, i, j, k, iflag, idum
real(kind=doblep), dimension(npmax):: rx, ry, rz, vx, vy, vz, ax, ay, az
!Arrays para almacenar posiciones (rx, ry, rz), velocidades (vx, vy, vz) y aceleraciones (ax, ay, az)
real(kind=doblep):: xnp, pa, pma, factor, px, py, pz, ecin, etotq, ecinq, epot, dfiv, d2fiv
! xnp: nĂşmero total de partĂculas (como double)
! pa: parámetro de red (espaciado entre partĂculas en la red cĂşbica)
! pma: mitad del parámetro de red
! factor: factor de escala usado en la correcciĂłn de energĂa
! ecin: energĂa cinĂ©tica
! etotq: energĂa total establecida (en mi caso es de -640)
! ecinq: energĂa cinĂ©tica deseada tras ajustar la energĂa total
! dfiv :: dervada primera del potencial respecto el volumen
! d2fiv :: derivada segunda del potencial respecto del volumen
!Valores necesarios para el calculo
pl=10.d0
rc=5.d0
idum = -1 !Para que inicialice el mi_random
xnp = dble(npmax)
pli=1.d00/pl
vol = pl*pl*pl
dens = xnp/vol
rc2 = rc*rc
pa = pl/dble(numk)
pma = pa/2.d00
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)
print*, corr_ener, corr_sum_r2vpp, corr_suma_rvp
!! POSICIONES:
!Colocamos las particulas y las voy contando
np=0 !Inicializamos
do i = 0, numk-1
do j = 0, numk-1
do k = 0, numk-1
if (np < npmax) then
np = np + 1
rx(np) = 0.d00 + dble(i) * pa
ry(np) = 0.d00 + dble(j) * pa
rz(np) = 0.d00 + dble(k) * pa
end if
if (np < npmax) then
np = np + 1
rx(np) = pma + dble(i) * pa
ry(np) = 0.d00 + dble(j) * pa
rz(np) = pma + dble(k) * pa
end if
if (np < npmax) then
np = np + 1
rx(np) = 0.d00 + dble(i) * pa
ry(np) = pma + dble(j) * pa
rz(np) = pma + dble(k) * pa
end if
if (np < npmax) then
np = np + 1
rx(np) = pma + dble(i) * pa
ry(np) = pma + dble(j) * pa
rz(np) = 0.d00 + dble(k) * pa
end if
end do
end do
end do
write(*,*) "Se han colocado en la red", np, "particulas"
iflag = 64 ! Establezco aquĂ un numero para iniciar el generador de numeros aleatorios
! TambiĂ©n podrĂa introducirlo por consola pero de esta manera es más rápido para ver los resultados
! Call the mi_random function
px = generator_random(idum)
! Aplicar perturbaciones aleatorias a las posiciones de las partĂculas
do i=1, np
rx(i) = rx(i)+(2.d00*generator_random(iflag)-1.d00)*0.1d0*pma
ry(i) = ry(i)+(2.d00*generator_random(iflag)-1.d00)*0.1d0*pma
rz(i) = rz(i)+(2.d00*generator_random(iflag)-1.d00)*0.1d0*pma
enddo
! Calcular la energĂa potencial del sistema usando un potencial Lennard-Jones
call potlj (np,rx,ry,rz,ax,ay,az,epot,dfiv,d2fiv)
!! VELOCIDADES:
! Asignar velocidades aleatorias entre -1 y 1 para cada partĂcula
do i = 1, np
vx(i) = 2.d00*generator_random(iflag)-1.d00
vy(i) = 2.d00*generator_random(iflag)-1.d00
vz(i) = 2.d00*generator_random(iflag)-1.d00
enddo
!Calculamos los momentos
px = sum(vx)
py = sum(vy)
pz = sum(vz)
!TendrĂan que ser 0 (Velocidad CM = 0) -> Corregimos para que sea nulo
px = px/xnp !Velocidad promedio
py = py/xnp
pz = pz/xnp
vx = vx-px !Restamos esta velocidad promedio para equilibrar el sistema
vy = vy-py
vz = vz-pz
! Recalculamos con los ajustes
px = sum(vx)
py = sum(vy)
pz = sum(vz)
ecin = 0.5d00 * (sum(vx*vx)+sum(vy*vy)+sum(vz*vz))
!Imprimimos valores por consola:
write(*,*)
write(*,*) 'Momentos:', px, py, pz
write(*,*) 'EnergĂa cinĂ©tica:', ecin, "EnergĂa potencial:", epot, "EnergĂa total:", ecin+epot
! Ajustamos la energĂa total al valor definido por el usuario
!write(*,*) "ÂżQuĂ© energĂa total desea fijar?"
!read(*,*) etotq
etotq = -640.d0
ecinq = etotq - epot
if (ecinq.le.0.d00) then
write(*,*) "Error: energĂa cinĂ©tica negativa" ! Nos aseguramos que la energĂa cinĂ©tica no puede ser negativa
stop
end if
!Ajustamos las velocidades en consecuencia de esta energĂa
factor = dsqrt(ecinq / ecin)
vx = factor * vx
vy = factor * vy
vz = factor * vz
!Recalculamos para verificar
px = sum(vx)
py = sum(vy)
pz = sum(vz)
write(*,*)
write(*,*) 'Despues de modificar la energĂa:'
write(*,*) 'Momentos:', px, py, pz
write(*,*) 'EnergĂa cinĂ©tica:', ecinq, "EnergĂa potencial:", epot, "EnergĂa total:", ecinq+epot
! FICHEROS:
open (10, file='variables.dat')
! Abrimos el archivo 'variables.dat' para escribir las variables principales en formato de texto
write(10, '(i8)') np ! NĂşmero de partĂculas
write(10, '(f12.6)') pl ! Lado de la caja
write(10, '(f12.6)') rc ! Radio de corte del potencial
write(10, '(f12.6, 2x, f12.6)') vol, dens ! Volumen y Densidad
write(10, '(f12.6, 2x, f12.6, 2x, f12.6)') ecinq + epot, ecinq, epot ! EnergĂas
close(10)
! Abrimos otro archivo llamado 'resultados.dat' para guardar los datos en formato binario
open(20, file='resultados.dat', form='unformatted')
write(20) rx,ry,rz,vx,vy,vz,ax,ay,az
close(20)
stop
END PROGRAM crea_red
Common variables:
module variables_comunes
use def_precision
!preparado para 500 partĂculas
!podemos modificar npmax con un numero k adecuado :: npmax = 4*k**3
! k = 1, 2, 3, 4, 5, 6...
! npmax = 4, 32, 108, 256, 500, 864 ...
implicit none
integer(kind=entero), parameter :: npmax=500
integer(kind=entero), parameter :: numk=5
real(kind=doblep), parameter :: pi=3.141592653589d00
real(kind=doblep):: pl, pli, vol, dens, rc, rc2, dt, dt12, dt2
real(kind=doblep):: corr_ener=0.d00, corr_suma_rvp=0.d00, corr_sum_r2vpp=0.d00
end module variables_comunes
