Clase # 10
Simulación de Monte Carlo
Prof. Ramón Garduño Juárez
Modelado Molecular
Diseño de Fármacos
Contorno
•
•
•
•
Introducción
Ensamble Canónico
Simple Monte Carlo
Método de Metrópolis
– Cadena de Markov
– Implementación
– Generación de números aleatorios
• MC para moléculas
– Moléculas Rígidas
– Moléculas Flexibles
• Modelos usados en Polímeros
– Modelo Lattice
– Modelo Continuo
Introducción
• Ensamble Canónico (NVT) función de partición Z i.e. Q
Q NVT 
1
N! h
1

3N
N
H (p ,r ) 
N
N

i 1
Q NVT 
1
N! h
1
3N
N
dp dr
pi
N
 H (pN ,rN )
exp  

k
T
B


2
 V (r )
N
2m
 dp
N
exp[ 
p
2
2 mk B T
N
]  dr
N
exp[ 
V (r )
k BT
]
Introducción
La separación es posible sólo si la función de potencial V(rN) no es dependiente
de las velocidades (es seguro así suponerlo para casi todas las funciones de
potencial comúnmente usadas en biología). Así la integral sobre el momento
puede ser llevada a cabo analíticamente:
 dp
N
p
exp[ 
2
2 mk B T
]   2 mk B T 
3N /2
La función de partición puede ser descrita como:
Q NVT
1  2 mk B T 


  dr
2
N!
h

N
N
exp( 
V (r )
)
k BT
La integral sobre las posiciones es a menudo referida como la integral
configuracional, ZNVT
N
Z NVT 
 dr
N
exp( 
V (r )
k BT
)
Calculando las propiedades por integración
Cálculo del potencial promedio:
V (r ) 
N
 dr
N
V (r )  (r )
N
N
exp[  V ( r ) / k B T ]
N
 (r ) 
N
Z
El denominador Z, la integral configuracional no puede ser evaluada analíticamente, pero los
métodos numéricos pueden dar la solución..
Para una función de dos variables (f(x,y)), es necesario elevar al cuadrado el número de
evaluaciones requeridas de la función. Para una integral 3N-dimensional el número total de
evaluaciones requeridas de la función deberá ser m3N, donde m es el número de puntos
necesarios para determinar la integral en cada dimensión. Este número es enorme aún
para un número muy pequeño de partículas (para 50 partículas y tres puntos por
dimensión, un total de 3150 ~ 1071 evaluaciones serán requeridas).
Método Azaroso como una posible alternativa
Una área bajo la curva puede ser estimada como la relación de los puntos bajo la curva al número total de
puntos generados.
El cálculo de la función de partición para un sistema de N átomos usando esta simple integración de Monte
Carlo que involucra los siguientes pasos:
1.Obtener una configuración al generar azarosamente 3N coordenadas Cartesianas que se asignan a las
partículas.
2.Calcular la energía potencial V(r) de la configuración y forma de esta energía calcular el factor de
Boltzmann exp(-V(r)/kBT).
3.Añadir el factor de Boltzmann a la suma acumulada de los factores de Boltzmann y la contribución de
energía potencial a la suma acumulada y regresar al paso 1.
4.Después de un número de iteraciones Ntrial, el valor promedio de la energía potencial deberá calcularse
usando:
V (r ) 
N

N tria l
i 1

V i ( r ) exp[ V i ( r ) / k B T ]
N
N tria l
i 1
N
exp[ V i ( r ) / k B T ]
N
Método Azaroso como una posible alternativa
• Éste no es un enfoque factible para calcular propiedades termodinámica debido al gran
número de configuraciones que tienen factores de Boltzmann extremadamente pequeños
(efectivamente cero) debido a la alta energía (principalmente causada por el traslape entre
las partículas).
• Esto refleja la naturaleza del espacio fase, en su mayoría corresponde a configuraciones no
físicas con energías muy altas. Sólo una proporción muy pequeña del espacio fase
corresponde a configuraciones de baja energía donde no hay partículas traslapadas y
donde el factor de Boltzmann tiene un valor apreciable (el exponente en el factor de
Boltzmann es el valor negativo de la energía).
• Una forma de resolver este problema es el generar configuraciones que hagan la mayor
contribución a la integral:
V (r ) 
N

N tria l
i 1

V i ( r ) exp[ V i ( r ) / k B T ]
N
N tria l
i 1
N
exp[ V i ( r ) / k B T ]
N
• La extensión de Metrópolis al método azaroso genera configuraciones preferenciadas con la
probabilidad exp(-V(rN)/kBT) que es diferente a la del Monte Carlo simple con igual probabilidad
y luego les asigna un peso exp(-V(rN)/kBT).
Metrópolis
• “EL Método de Monte Carlo” en la comunidad
de modelado biomolecular.
• Tendencia hacia aquellas configuraciones que
hacen grandes contribuciones a la integral.
• Éste genera estados con una probabilidad
exp(-V/kT) y las cuenta de manera igual.
• Metrópolis genera una cadena de Markov
– El resultado de las pruebas depende solamente de
la prueba previa, pero no de todas las pruebas
anteriores.
– Cada prueba pertenece a un grupo finito de
posibles resultados.
• Matriz de Transición
• Aceptación del estado nuevo
Implementación de Metrópolis
• Procedimiento en la simulación de un fluido
atómico
– Energía
– Movimiento
– Coeficiente de aceptación (50%)
• Selección/ajuste de Rmax
• Selección de partícula : random vs secuencial
• Mover varias partículas a la vez
El algoritmo de Metrópolis genera una cadena de estados Markov
1. El resultado de cada prueba depende solamente de la prueba previa y no de todas las demás pruebas anteriores.
2. Cada prueba pertenece a un grupo finito de posibles resultados.
Condición 1 provee una clara distinción entre MC y dinámica molecular.
Formalización:
– Una matriz estocástica es una matriz n x n de valores no negativos, donde cada hilera suma a 1.
– Hagamos Q = {1, . . . , n} un grupo finito de estados, y considérese la distribución de probabilidad inicial
  = (p1, …, pn), considerado como un vector hilera, y la matriz estocástica n x n P= (pmn).
– Una cadena de Markov de primer orden, tiempo-homogéneo M=(Q, , P) es un proceso estocástico, cuyo
estado qt a tiempo t es una variable azarosa determinada por:
Pr [ q 0  i ]   i ,
Pr [ q t 1 | q t  i ]  p i , j
Que define
Pr [ q 0  i ]   i ,
Pr [ q t 1 | q t  i ]  p i , j
Es obvio que p i ( t )  Pr [ q t  i ]
La entrada
y
p i , j  Pr [ q t  j | q 0  i ]
(t )
(t )
(i,j)th de la tth potencia de Pt de P es igual a p i , j
La homogeneidad del tiempo da p i , j  Pr [ q t 0  t  j | q t 0  i ] para todo t0.
(t )
La probabilidad del tercer estado es
p 3 ( t )   2 ( t  1) Pr [ q t 1  2 | q t  3 ]   1 ( t  2 ) Pr [ q t  2  1 | q t  2  2 ] Pr [ q t 1  2 | q t  3 ]
Implementación de Metrópolis
Si el generador de números azarosos producen números (x en el intervalo de 0
a 1, se mueve en ambas, las direcciones positivas y negativas son cambiadas de
acuerdo a las siguientes reglas:
x new  x old  ( 2 x  1) rmax
y new  y old  ( 2 x  1) rmax
y new  y old  ( 2 x  1) rmax
• Un valor único azaroso se genera para cada una de las tres coordenadas.
• La energía de la nueva configuración se calcula y es comparada con el valor
de la anterior. Si la configuración es una de más baja energía es aceptada.
Si la energía es más grande, el factor de Boltzmann de la configuración
exp(-V(rN)/kBT) es comparado con el número azaroso entre 0 y 1. Si el
número azaroso es más grande que el factor de Boltzmann la configuración
es aceptada; en cualquier otro caso la configuración es rechazada. Esta
regla de aceptación se puede expresar como:
rand(0,1) ≤ exp(-V(rN)/kBT)
Generación de números aleatorios
• Generador de números aleatorios
– No es verdaderamente azaroso
– Algunos números deberán ser generados con la
misma semilla
• Propiedades estadísticas
– Distribución azarosa en el espacio k-dimensional
• Congruente Lineal
–
e[1] = semilla
–
e[i] = MOD{(e[i-1]*a+b), m}
Monte Carlo de moléculas
• Más difícil que en sistemas atómicos
contra uno solamente necesita
considerar los grados de libertad de
traslación.
• El conjunto de movimientos debe incluir
el cambio de los grados internos de
libertad.
• Rotación
– Ángulos de Euler
– Muestreo un informe
Moléculas Rígidas
Reglas de rotación
 x'  1
  
 y '   0
 z'  0
  
0
cos 
 sin 


sin   
cos   
0
x

y
z 
• Muestreo un informe de los ángulos de Euler
– Ángulos de Euler f q y
– Muestrear cosq más que q
Moléculas Rígidas
Cuaterniones: q1, q2, q3, q4, un vector 4-D.
Revisar el coeficiente de aceptación deseado para los componentes de
traslación y rotación
q 0  cos
1
q 1  sin
1
q 2  sin
1
q 3  cos
1
q cos
2
1
f
2
q cos
2
1
y

(f  y )
2
q sin
2
1
(f  y )
2
q sin
2
1
(f  y )
2
La matriz de rotación de los ángulos de Euler se puede expresar
usando estas relaciones:
 q 02  q 12  q 22  q 32

A   2 ( q1 q 2  q 0 q 3 )

 2 ( q1 q 3  q 0 q 2 )
2 ( q1 q 2  q 0 q 3 )
q 0  q1  q 2  q 3
2
2
2
2 ( q 2 q 3  q 0 q1 )
2
2 ( q1 q 3  q 0 q 2 ) 

2 ( q 2 q 3  q 0 q1 ) 
2
2
2
2 
q 0  q1  q 2  q 3 
Moléculas Flexibles
• Difícil
• A menos que el sistema sea
pequeño, o con algunos
grados de libertad internos
congelados, o cuando se usan
modelos especiales
• La manera más simple es el
cambio aleatorio de las
coordenadas Cartesianas
• Grado de cobertura es lento
• Pequeñas rotaciones a
menudo causante
movimientos muy grandes
Modelos usados en Polímeros
• Un polímero es una macromolecular está construida
al ligar químicamente otra secuencia semejante o
otros fragmentos moleculares.
• Diferente escala de tiempo
– Desde fs a segundos u horas (estiramiento de la unión al
plegado de proteínas/desnaturalización a asociación)
• Diferente escala de longitud
– Desde ~1 to ~1000 Ángstrom
• Modelos de energía empírica
Modelo Lattice de polímeros
•
•
•
•
•
•
Los modelos energía son simples
Cálculo rápido de la energía
Diferentes longitudes de unión son posibles
Camino azaroso
Camino a auto-evitante
Conjunto de movimientos
–
–
–
–
Crankshaft
Kink jump
End rotation
Slithering snake
Modelo Continuo de polímero
• Cadena de cuentas conectadas con varillas rígidas o
por resortes armónicos
• Cadena que rota libremente
• El coeficiente de aceptación es muy pequeño
• Modelo del estado rotacional isomérico (RIS)
– Cada unión se supone que adopta alguna de un número
pequeño de estados rotacionales discretos, el cual
generalmente corresponde al mínimo de energía.
• Un generador de matrices se usa para establecer
cierta propiedades dependientes de la conformación.
• Una matriz de peso estadístico: trata con la influencia
que una unión tiene sobre sus vecinos.
Descargar

Monte Carlo Metropolis simulation