DINAMICA MOLECULAR Métodos de simulación -Permitir el estudio de propiedades de sistemas complejos

1. DINAMICA MOLECULAR • Métodos de simulación • -Permitir el estudio de propiedades de sistemas complejos • -Generación de conjunto de configuraciones distintas para un mismo sistema • -Predicción del comportamiento temporal de un sistema (MD) • Dos métodos: • -Dinámica molecular - promedio temporal • -Monte Carlo - promedio de ensamblado • Propiedad A: dependiente de las posiciones y momentos de todas las partículas: • A(pN(t), rN(t)) • pN(t) = momentos de las partículas en el tiempo t • rN(t) = posiciones de las partículas en el tiempo t • Momentos y posiciones definen el ESPACIO DE FASE

2. Valor exp. = promedio durante el tiempo del exp. (promedio temporal) Aave = lim t ∞ • Se puede calcular el valor promedio de una propiedad simulando el comportamiento de un sistema, i.e. determinando los valores instantáneos de la propiedad • Alternativa: promedio de ensamble  mecánica estadística <A> = valor esperado de la propiedad r(pN, rN) = densidad de probabilidad Hipótesis ergódica <A> = Aave

3. r corresponde a la distribución de Boltzmann: E(pN, rN) = energía Q = función de partición Q para el ensamble canónico = N, V y T ctes. Propiedades termodinámicas Mecánicas -Energía -Capacidad calorífica

4. A partir de una sola simulación: CV = {<E2> - <E>2}/kBT2 = <(E -<E>)2>/kBT2 -Presión Teorema del virial xi = coordenada i p´xi = derivada de p a lo largo de coordenada i (fuerza) Para un gas real:

5. Temperatura NC = número de “constricciones” del sistema • Térmicas • -Entropía • -Energía libre • -Potencial químico • No son fácilmente derivables  técnicas especiales • Propiedades mecánicas  dependen de la derivada primera de la función de partición • Propiedades térmicas - dependen de la función de partición misma

6. Método de Dinámica Molecular (DM o MD) Concepto básico: generación de configuraciones sucesivas de un sistema a partir de la integración de las leyes de movimiento de Newton: Fxi = fuerza sobre una partícula en la dirección xi • Tres situaciones • Conjunto de moléculas que sufren colisiones en ausencia de fuerzas externas entre colisiones • Cada partícula experimenta una fuerza constante entre las colisiones (ej: carga en un campo eléctrico) • La fuerza sobre cada partícula depende de la posición relativa de las demás partículas Primeras simulaciones fueron del tipo 1 (modelos de esferas duras) Primera simulación “real”: Argon (1964)

7. Simulaciones con potenciales continuos y variables • Para resolver las ecuaciones diferenciales (Newton) se aplican métodos de “diferencia finita” • Idea básica: la integración de las ecuaciones se realiza en etapas de tiempo muy pequeñas (Dt) • Para cada Dt: • 1) Evaluar la energía potencial de sistema (FF) a partir de las posiciones actuales de las partículas • 2) Derivar la fuerza para cada partícula • 3) Obtener la aceleración para cada partícula • 4) Calcular nuevas posiciones a partir de aceleración, velocidad y posiciones actuales • Se asume fuerzas ctes. durante cada Dt

8. Algortimos de integración Verlet Ventajas: -Simple y expeditivo -Combinable con otros algoritmos Desventajas: -Poca precisión -Cálculo de velocidades 2) Leap-frog Utiliza la aceleración y velocidad para el cálculo de la posición:

9. Ventajas: -Simple -Calculas las velocidades explícitamente Desventajas: -Posiciones y velocidades se obtienen desfasadas en el tiempo 3) Velocity Verlet

10. Ventajas: -Calcula velocidades explícitamente y en fase con posiciones 4) Beeman Asume una variación lineal de la aceleración en el intervalo Dt Ventajas: -Mayor precisión -Mayor tamaño de Dt -Cálculo explícito de velocidades Desventajas: -Mayor costo computacional

11. 5) Gear (predictor-corrector) • 3 pasos: • 1) Predicción de r(t + Dt), v(t + Dt), a(t + Dt), etc. a partir de expansiones de Taylor • 2) Cálculo de la fuerza en las nuevas posiciones y obtención de a(t + Dt) • 3) Comparación de a (t + Dt) calculadas y predichas • La diferencia entre los valores predichos y calculados se usan para corregir las posiciones, velocidades, etc., en el paso de corrección: aceleración teórica

12. Ventajas: • -Gran precisión • -Uso de Dt mayor • Desventajas: • -Mayor espacio de almacenamiento • -Evaluación doble de las fuerzas en cada Dt • SHAKE • Algoritmo de constricción • Permite el uso de Dt significativamente mayores  simulaciones más largas • Constraints del tipo: • rij2 – dij2 = 0 constricción holonómica • rij = distancia entre átomos i y j • dij= valor fijo de distancia entre i y j

13. SHAKE utiliza un procedimiento iterativo que en cada Dt ajusta todas las distancias de enlace a valores pre-establecidos • Elección del algoritmo de DM • Factores: • 1) Esfuerzo computacional en cada paso • 2) Tamaño de Dt • 3) Conservación de la energía (ens. microcanónico) • 4) Requerimientos de memoria y espacio en disco • 5) Precisión • En general, el algoritmo óptimo debe demandar poco costo computacional (cálculo, memoria y espacio), debe ser preciso (resultados similares a los experimentales), posibilitar trayectorias largas en tiempos razonables, ser moldeable (sencillo de modificar) y adaptarse adecuadamente a diferentes condiciones de simulación

14. Dinámica Molecular a temperatura y presión constantes -Ensamble microcanónico: N, V y E ctes. -Ensamble canónico: N, V y T ctes. -Ensamble isotermo-isobárico: N, T y P ctes. • Temperatura: • <K>NVT = 3/2NkBT • Escalamiento de velocidades • Se multiplican las velocidades por un cierto factor (l) T0 = temperatura de referencia T(t) = temperatura instantánea DT = (l2 -1)T(t)

15. 2) Acople térmico • Se acopla el sistema a un baño térmico (T) con el cual intercambia energía • Tba = temperatura del baño • = constante de acoplamiento DT = Dt/t(Tba – T(t)) Dt = tiempo de integración Ventaja -Permite la fluctuación de temperatura alrededor del valor de referencia

16. 3) Colisiones estocásticas • Se reasignan velocidades a partir de una distribución de Maxwell-Boltzmann a la temperatura de referencia • Baño térmico emitiendo partículas térmicas! • La frecuencia de “colisiones” para una partícula es: k = conductividad térmica h = densidad de partículas N = número de partículas o n =nc/N2/3nc = frecuencia de colisión intermolecular Desventaja -Introduce discontinuidades en la trayectoria

17. Presión • Escalamiento del volumen • Acople a baño de presión k = compresibilidad isotérmica r´i = l1/3ri Escalamiento: -isotrópico -anisotrópico

18. Algoritmo de Dinámica Molecular • En cada Dt de integración, el algoritmo tipo (basado en leap-frog) de DM puede resumirse en los siguientes pasos: • 0)Las posiciones (ri(tn)) y velocidades (vi(tn – Dt/2) para todos los átomos y el largo (Rbox(tn)) y volumen (Vbox(tn) ) de la caja son conocidos en t = tn • Cálculo de las fuerzas (no constreñidas): • 2) Cálculo de la energía cinética, el virial y la presión:

19. Algoritmo de Dinámica Molecular • 3) Cálculo de la temperatura: • 4) Cálculo de las nuevas velocidades no constreñidas: • 5) En caso de ajuste de temperatura, multiplicación de velocidades por factor: • 6) Cálculo de las nuevas posiciones no constreñidas:

20. Algoritmo de Dinámica Molecular • 7) En caso de aplicación de constricciones (SHAKE), almacenamiento de posiciones no constreñidas: • 8) En caso de aplicación de constricciones (SHAKE), cálculo de nuevas posiciones constreñidas: • 9) En caso de aplicación de constricciones (SHAKE) cálculo de nuevas velocidades constreñidas:

21. Algoritmo de Dinámica Molecular • 10) En caso de aplicación de constricciones (SHAKE), cálculo de fuerzas no constreñidas: • 11) En caso de ajuste de presión, multiplicación de posiciones por factor: • y cálculo del nuevo largo y volumen de la caja:

22. Fases de una simulación (MD o MC) • 1) Inicio • Elección del modelo energético (FF) • Elección de la configuración de partida y entorno • Asignación de velocidades • 2) Equilibramiento • Evolución del sistema desde la configuración inicial hasta el equilibrio (incluye pre-calentamiento) • 3) Producción • Recolección de datos y cálculo de propiedades simples • 4) Análisis • Cálculos más complejos de propiedades • Examen estructural • Detección de problemas • Estimación de errores

23. Condiciones de borde o frontera • Aspecto relevante en una simulación • Tipos de borde • 1) Vacío • El sistema se encuentra totalmente aislado • Aproximación pobre: • -Distorsión de la forma molecular en la interfase • -Cambios en la forma global por efecto de interacciones intramoleculares • 2) Seudovacío • Sin solvente explícito pero se incluyen efectos del solvente • Ej: campo de fuerza NIS, dinámica estocástica • Mejor aproximación: algunos efectos del solvente son tenidos en cuenta: apantallamiento de cargas e interacciones de van der Waals.

24. 3) Condiciones periódicas (PBC) • Copias del modelo repetidas periódicamente

25. Ventaja: -Permite simular sistemas pequeños atenuando efectos de vacío: los átomos siempre están rodeados de otros átomos Desventajas: -Pérdida de fluctuaciones con longitud de onda mayor al largo de la celda -Pérdida de interacciones no enlazantes 4) Solvente -Condiciones periódicas -Pared extendida (“extended wall”)

26. Tratamiento de las interacciones no enlazantes • Evaluación costosa - pares de interacciones prop. a N2 • Potenciales no enlazantes dependen inversamente de la distancia: • Uso de radios de cut off: • 1) Twin range method • Dos radios de cut off • -rc1 Lennard Jones + Coulomb • -rc2 Coulomb • rc2 > rc1 (ej: rc1 = 8Å y rc2 = 15Å) • En PBC  convención de mínima imagen: la interacción se calcula con respecto al átomo mismo o a su imagen más cercana (una sólo vez) • rc2 < l/2 (caja cúbica) • Combinación de radios cut off con lista de vecinos no enlazantes (non bonded neighbour list)

27. Los vecinos se determinan cada n Dts y se almacenan en una lista: • En cada Dt, se calculan las distancias con respecto a la lista de vecinos • rcnei >= rc2 • Necesidad de actualizar la lista de vecinos periódicamente • 2) Twin range + neighbour pair list • r < rc1 todas las interacciones en cada Dt • rc1 < r < rc2 - sólo interacciones electrostáticas cada vez que se actualiza la lista de vecinos

28. 3) Radios de cut off entre grupos • Cálculo de interacciones no enlazantes entre grupos de átomos en vez de átomos aislados. • Concepto de grupos de carga (charge groups): • Carga total = 0  dipolo-dipolo término dominante • Radios de cut off entre grupos: • -Medido entre los centros de masa • -Considerando un átomo marcador • Ventajas • -Menor alcance de interacciones electrostáticas (1/r3) • -Funciona mejor con radios de cut off Desventajas del uso de radios de cut off 1) Aproximación 2) Introducción de discontinuidades en el potencial y la fuerza

29. Mejoras en la discontinuidad del potencial 1) Potenciales corridos: “shifted potentials” V´(r) = V(r) – Vc r <= rc V´(r) = 0 r > rc Vc = valor del potencial en rc 2) Funciones de cambio: “switching functions” S = 1 rij < rc1 S = (rc2 –rij)/(rc2 –rc1) rc1<=rij<=rc2 S = 0 rc2 < rij V´(r) = V(r)S(r) Tratamiento de las fuerzas de largo alcance -Sumatoria de Ewald -Campo de reacción -Celdas multipolo

30. Campo de reacción • Término adicional al potencial electrostático (calculado dentro de rc) • Vi = Ei.mi • mi: dipolo de molécula i • mj: dipolo de molécula j dentro de rc • es: constante dieléctrica del medio fuera de rc • Se modelan los átomos fuera del radio de cut off como un dieléctrico continuo (buena aproximación para fluidos) • Ventajas • -Simple conceptualmente • -Fácil de implementar • Desventajas: • -Necesidad de conocer es

31. Análisis de la simulación 1) Energético: • Evolución de la energía total (cinética + potencial) • Evolución de la energía potencial • Evaluación de componentes de energía potencial: -Energía covalente: enlaces y ángulos -Energía no covalente: diedros e interacciones no enlazantes -Energías de interacción no enlazante entre grupos (soluto-solvente, ligando- proteína, etc.) • Promedios de energía total y potencial y cálculo de energías libres 2) Estructural: RMSD (desviaciones cuadráticas medias) • Valor de la desviación promedio entre todos los átomos o coordenadas de conformaciones

32. Factores B (temperatura) • Valor de la desviación promedio de cada átomo por separado en las distintas conformaciones m: masa del átomo i ri(tk): posición de átomo i en el tiempo tk <r>i: posición promedio de átomo i

33. P.M.I. (momentos principales de inercia) • Medida de la distribución de masa del sistema entre los ejes cartesianos • Noción de volumen y forma del sistema Centro de masas para fragmento k: ri: posición de átomo i ; mi: masa del átomo i nk: número de átomos de fragmento k Radio de giro para fragmento k: Matriz de momentos de inercia: 3 valores propios de Mk --- momentos de inercia del fragmento K 1: matriz identidad