Dinámica molecular

Transcript Dinámica molecular

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(

p

N (t),

r

N (t))

p

N (t) = momentos de las partículas en el tiempo t

r

N (t) = posiciones de las partículas en el tiempo t Momentos y posiciones definen el

ESPACIO DE FASE

 Valor exp. = promedio durante el tiempo del exp. (promedio temporal) A ave = lim t  ∞ 1 t

t   0

(

))

 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 

  

N d

N A

(

)  (

) = valor esperado de la propiedad  (

p

N ,

r

N ) = densidad de probabilidad

Hipótesis ergódica

= A ave

 corresponde a la distribución de Boltzmann:  (

)  exp( 

(

) /

k B T

) /

E(p N , r N ) = energía Q = función de partición

Q NVT

 1 1

 

N d

exp[ 

(

) /

k B T

)] Q para el ensamble canónico = N, V y T ctes.

Propiedades termodinámicas Mecánicas

-Energía U

 1 /

M M

 

E i -Capacidad calorífica C V

 ( 



)

A partir de una sola simulación: C V = { - 2 }/k B T 2 = <(E -) 2 >/k B T 2

-Presión

Teorema del virial

 

x i p

x i

  3

Nk B T

x i = coordenada i p´ xi = derivada de p a lo largo de coordenada i (fuerza) Para un gas real:

  3



i N

 1

j N

  

 1

r ij dV

(

r ij

)

dr ij

  3

Nk B T P

 1 /

 

Nk B T

 1 3

k B T N N

   

r ij f ij

 

Temperatura K



i N

  1 |

m i

| 2 

k B T

2 ( 3



N c

) N C = 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  partición dependen de la derivada primera de la función de  Propiedades térmicas  dependen de la función de partición misma

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:

x i dt

2 

F xi m i

F xi = fuerza sobre una partícula en la dirección x i

Tres situaciones

1) Conjunto de moléculas que sufren colisiones en ausencia de fuerzas externas entre colisiones 2) Cada partícula experimenta una fuerza constante entre las colisiones (ej: carga en un campo eléctrico) 3) 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)

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 ( D t) Para cada D t: 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 D t

Algortimos de integración

Verlet

(

 D

)  D

) 

(

)  D

(

)  1 / 2 D

(

) 

(

)  D

(

)  1 / 2 D

(

)

(

 D

)  2

(

) 

(

 D

)  D

(

)

(

)  [

(

 D

) 

(

 D

)] / 2 Δ

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:

(

 D

) 

(

)  D

(

 1 / 2 D

)

(

 1 / 2 D

) 

(

 1 / 2 D

)  D

(

)

Ventajas: -Simple -Calculas las velocidades explícitamente Desventajas: -Posiciones y velocidades se obtienen desfasadas en el tiempo

3) Velocity Verlet

(

 D

) 

(

)  D

(

)  1 / 2 D

(

)

(

 D

) 

(

)  1 / 2 D

[

(

) 

(

 D

)]

Ventajas: -Calcula velocidades explícitamente y en fase con posiciones

4) Beeman

Asume una variación lineal de la aceleración en el intervalo D t

(

 D

) 

(

) 

[

(

) 

(

 D

)] / D

(

 D

) 

(

)  D

(

)  2 / 3 D

(

)  1 / 6 D

(

 D

)

(

 1 /  6 D

) D

( 

 (

D )

)  1 / 3 D

(

 D

)  5 / 6 D

(

) Ventajas: -Mayor precisión -Mayor tamaño de D t -Cálculo explícito de velocidades Desventajas: -Mayor costo computacional

5) Gear (predictor-corrector)

 3 pasos: 1) Predicción de r(t + D t), v(t + D t), a(t + D t), etc. a partir de expansiones de Taylor 2) Cálculo de la fuerza en las nuevas posiciones y obtención de a(t + D t) 3) Comparación de a (t + D t) 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: D

(

 D

) 

(

 D

) 

(

 D

) aceleración teórica

(

 D

) 

(

 D

) 

0 D

(

 D

)

(

 D

) 

(

 D

) 

1 D

(

 D

)

(

 D

) / 2 

(

 D

) / 2 

2 D

(

 D

)

(

 D

) / 6 

(

 D

) / 6 

3 D

(

 D

)

Ventajas: -Gran precisión -Uso de D t mayor Desventajas: -Mayor espacio de almacenamiento -Evaluación doble de las fuerzas en cada D t

SHAKE

 Algoritmo de constricción  Permite el uso de D t significativamente mayores  simulaciones más largas  Constraints del tipo: r ij 2 – d ij 2 = 0 constricción holonómica r ij = distancia entre átomos i y j d ij = valor fijo de distancia entre i y j

 SHAKE utiliza un procedimiento iterativo que en cada D t 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 D t 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

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:

NVT

= 3/2Nk

T

1) Escalamiento de velocidades

Se multiplican las velocidades por un cierto factor (l ) l  (

T o

(

)) T 0 = temperatura de referencia T(t) = temperatura instantánea D

 1 / 2

i N

  1 2 / 3    

m i

( l

v i

) 2

Nk B

     1 / 2

i N

  1 2 / 3    

m i

(

v i

) 2

Nk B

    D T = ( l 2 -1)T(t)

2) Acople térmico

 Se acopla el sistema a un baño térmico (T) con el cual intercambia energía

(

)

 1 / t (

T ba



(

)) T ba = temperatura del baño t = constante de acoplamiento D T = D t/ t (T ba – T(t)) D t = tiempo de integración l 2  1  D

T ba

(

)  1 )

Ventaja

-Permite la fluctuación de temperatura alrededor del valor de referencia

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:   2

 3

k B

 1 / 3

2 / 3 o  =  c /N 2/3  = conductividad térmica  = densidad de partículas N = número de partículas  c = frecuencia de colisión intermolecular Desventaja -Introduce discontinuidades en la trayectoria

Presión

1) Escalamiento del volumen 2) Acople a baño de presión

dP

(

t

)  1 / t

(

P



P

(

t

))

dt

l  1   D

(



P ba

)    1

( 



)

r´ i = l 1/3 r i

Escalamiento

: -isotrópico -anisotrópico  = compresibilidad isotérmica

Algoritmo de Dinámica Molecular

 En cada D t de integración, el algoritmo tipo (basado en leap-frog) de DM puede resumirse en los siguientes pasos: 0) Las posiciones (r i (t n )) y velocidades (v i (t n – D t/2) para todos los átomos y el largo (R box (t n )) y volumen (V box (t n ) ) de la caja son conocidos en t = t n 1) Cálculo de las fuerzas (no constreñidas):

F

(

t

)   

V

(

r



r

(

t

)) 2) Cálculo de la energía cinética, el virial y la presión:

E kin

(

t n

)  (

t n

)   

i N

  1 1 1 / / 2 (

m i

2 (

 

j v i

2 (

 D

t r ij

(

t n

F ij

(

t n

) ) / 2 ))

(

t n

)  2 (

E kin

(

t n

 D

/ 2 3

V box

(

t n

) )   (

t n

))

Algoritmo de Dinámica Molecular

3) Cálculo de la temperatura:

(

t n

 D

/ 2 )  2

E kin

(

t n

 D

/ 2 )

N df k B

c

k

t

T

(

t

t

Algoritmo de Dinámica Molecular

7) En caso de aplicación de constricciones (SHAKE), almacenamiento de posiciones no constreñidas:

(

t n

 D

) 

r i

(

t n

) 

m

     1  

(

0 

(

t n

))     1 / 3 y cálculo del nuevo largo y volumen de la caja:

R

box

(

)  l

V

box

(

t

)

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

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.

3) Condiciones periódicas (PBC)

 Copias del modelo repetidas periódicamente

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”)

Tratamiento de las interacciones no enlazantes

 Evaluación costosa  pares de interacciones prop. a N 2  Potenciales no enlazantes dependen inversamente de la distancia:  Uso de radios de

cut off: 1) Twin range method

 Dos radios de cut off -rc 1 Lennard Jones + Coulomb -rc 2 Coulomb rc 2 > rc 1 (ej: rc 1 = 8Å y rc 2 = 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) rc 2 < l/ 2 (caja cúbica)

Combinación de radios cut off con lista de vecinos no enlazantes (non bonded neighbour list)

 Los vecinos se determinan cada n D ts y se almacenan en una lista:  En cada D t, se calculan las distancias con respecto a la lista de vecinos rc nei >= rc 2 Necesidad de actualizar la lista de vecinos periódicamente 2)

Twin range + neighbour pair list

r < rc 1  todas las interacciones en cada D t rc 1 < r < rc 2  sólo interacciones electrostáticas cada vez que se actualiza la lista de vecinos

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/r 3 ) -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

Mejoras en la discontinuidad del potencial

1) Potenciales corridos: “shifted potentials” V´(r) = V(r) – V c V´(r) = 0 V c = valor del potencial en r c r <= r r > r c c 2) Funciones de cambio: “switching functions” S = 1 S = (r c2 –r ij )/(r c2 S = 0 –r c1 ) r ij < r c1 r c1 <=r ij <=r c2 r c2 < r ij V´(r) = V(r)S(r)

Tratamiento de las fuerzas de largo alcance

-Sumatoria de Ewald -Campo de reacción -Celdas multipolo

Campo de reacción

 Término adicional al potencial electrostático (calculado dentro de r c ) Vi =

E

m i

 2 (  

s s

  1 ) 1 1 (

r c

2 )

 

;

r ij

j c

m i : dipolo de molécula i m j : dipolo de molécula j dentro de r c  s : constante dieléctrica del medio fuera de r c  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  s

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

coord



N at

  ( 1

X i



X i

ref

) 2  (

Y i



Y i

ref

) 2  (

Z i



Z i

ref

)  2 3

N at R



N at

  ( 1

X i



X i

ref

) 2  (

Y i



Y i

ref

) 2  (

Z i



Z i

ref

)  2

N at Factores B (temperatura)

 Valor de la desviación promedio de cada átomo por separado en las distintas conformaciones

B

 8  3

m

k n

  1 (

r i

(

t k

)  



 2 m: masa del átomo i r i (t k ): posición de átomo i en el tiempo tk i : posición promedio de átomo i

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:

R gk



i n k

  1

m i r i i n k

  1

m i

r i : posición de átomo i; mi: masa del átomo i n k : número de átomos de fragmento k Radio de giro para fragmento k:

I k



i nk

  1

m i

(

r i



R gk

) 2 Matriz de momentos de inercia:

M k

 1

I K



n k i

  1

m i

(

r i



R gk

) 3 valores propios de Mk --   (

r i



R gk

) 1: matriz identidad

momentos de inercia

del fragmento K