Simulación de Dinámica Molecular Microcanónica
Abra el archivo md-simple.jl, que contiene el algoritmo de simulación. La simulación empieza con velocidades aleatorias, ajustadas para la media termodinámica deseada de 0.6 unidades/átomo ($kT=0.6$ en un sistema bidimensional). A esta energía cinética media le llamaremos "temperatura". El algoritmo de integración es Velocity-Verlet, que consiste en propagar la posiciones con
\[\vec{x}(t+\Delta t) = \vec{x}(t) + \vec{v}(t)\Delta t + \frac{1}{2}\vec{a}(t)\Delta t^2\]
siendo $\vec{a}(t)=\vec{F}(t)/m$, donde $\vec{F}(t)$ es la fuerza en el tiempo corriente. La fuerza en seguida es calculada en un tiempo posterior de tiempo con
\[\vec{F}(t+\Delta t) = -\nabla V\left[\vec{x}(t)\right]\]
y entonces las velocidades en el instante siguiente son calculadas con
\[\vec{v}(t+\Delta t) = \vec{v}(t) + \frac{1}{2}\left[ \frac{\vec{F}(t)}{m}+\frac{\vec{F}(t+\Delta t)}{m}\right]\]
completando el ciclo. En este ejemplo las masas son consideradas unitarias, para simplificar. La simulación es ejecutada por nsteps
pasos, con paso de integración $\Delta t$, este siendo un parámetro de entrada, dt
, definido en Options
.
La simulación no tiene control de temperatura o de presión. Es una propagación de la trayectoria según las leyes de Newton, que deberían conservar la energía. A esto se le llama una simulación "microcanónica", o "NVE" (que conserva, en principio, el número de partículas, el volumen y la energía total).
3.1. Paso de integración
Para realizar una MD simple, con un paso de integración de dt=1.0
, ejecute el comando:
julia> out = md(sys,Options(dt=0.1));
En principio, está previsto realizar 2000 pasos de integración de las ecuaciones de movimiento. Pruebe pasos de integración entre 1.0
y 0.01
. Note que pasa con la energía. Note que pasa con la energía cinética media, la cual fue inicializada en 0.6 unidades/átomo. Discuta la elección del paso de integración, y los valores de energía cinética obtenidos. Las simulaciones que siguen van a usar un paso de integración dt = 0.05
.
Es posible controlar la frecuencia de impresión y el número de puntos salvos en el archivo de trayectoria con las opciones iprint
y iprintxyz
:
julia> out = md(sys,Options(dt=0.1,iprint=1,iprintxyz=5))
El número total de pasos se controla con el parámetro nsteps
.
La variable out
que sale de la simulación es una matriz con las energías y la temperatura en cada paso de la simulación. Es probable que la simulación "explote" con pasos de tiempo grandes. Para visualizar este proceso, podemos hacer:
julia> using Plots
julia> plot(
out,ylim=[-100,100],
label=["Potential" "Kinetic" "Total" "Temperature"],
xlabel="step"
)
Y obtendremos un gráfico similar a:
Para pasos de tiempo menores la simulación debe conseguir llegar hasta el fin. Podemos ver el resultado nuevamente, y debe ser algo similar a:
Observe y trate de entender las amplitudes de las oscilaciones de las energías cinética y potencial, y las amplitudes de las oscilaciones de la energía total. A que se deben cada una de las oscilaciones? Observe como estas oscilaciones dependen del paso de integración.
3.2. Visualización de la trayectoria
Por fin, abra la trayectoria usando VMD. No es necesario salir de la sección de Julia
. Al apretar ;
(punto y coma) aparecerá un prompt shell>
, desde el cual es posible ejecutar a VMD
, si este está instalado correctamente y disponible en el path
:
shell> vmd traj.xyz
Dentro de VMD, elija la representación de VDW
, en
Graphics -> Representations -> Drawing Method -> VDW
y dé el comando
vmd> pbc set { 100. 100. 100. } -all
para indicar la periodicidad del sistema. Para representar explícitamente el sistema periódico, elija +X;+Y;-X;-Y
en
Graphics -> Representations -> Periodic
Para salir de VMD
use el comando exit
, y para volver al prompt de Julia
desde shell>
, use backspace
.
3.3. Código completo resumido
using FundamentosDMC, Plots
sys = System(n=100)
minimize!(sys)
out = md(sys,Options(dt=0.05))
plot(out,ylim=[-100,100],
label=["Potential" "Kinetic" "Total" "Temperature"],
xlabel="step"
)
Es generado el archivo traj.xyz
que puede ser visualizado en VMD
.