1. No category

MINIMIZACION DE ENERGIA
HIPERFICIE DE ENERGÍA POTENCIAL: describe como
varía la energía con respecto a las coordenadas.
Ej: variación energética al rotar diedros en n-pentano. La
energía en este caso es un a función de solamente dos variables
y puede visualizarse de la siguiente forma.
En Modelado Molecular hay especial interés en los puntos
mínimos de la superficie de energía potencial.
1
1
En general existe más de un mínimo de energía en una superficie
de energía potencial de una molécula:
Mínimo global: mínimo con la menor energía posible.
Mínimo locales: los restantes mínimos.
Algoritmo de minimización: identifica geometrías del sistema
que corresponden a mínimos en la superficie de energía
potencial.
Otros aspectos:
• variación en posiciones atómicas durante una reacción.
• adopción de diferentes mínimos.
• diferencias de energías relativas.
• otros puntos estacionarios: estructuras de transición.
CONCEPTO DE MÍNIMO DE UNA FUNCIÓN
Dada una función f que depende de un conjunto de variables
independientes: x1, x2, …..xi,….., xn
Minimizar la función equivale a encontrar el conjunto de valores
de las variables independientes (x*) tal que:
f(x*) = min(f(x))
2
2
En el mínimo se cumple:
∂f
∂2 f
= 0;
>0
∂xi
∂xi 2
En Modelado Molecular pretendemos minimizar energías.
La diferencia radica en la expresiones utilizadas para la
energía.
Mecánica Cuántica: Hamiltoniano
Mecánica Molecular: Campo de fuerza
Variables independientes = coordenadas de los átomos:
- cartesianas
-internas
La minimización exacta (analítica) no es posible en general
Æ Métodos numéricos
ALGORITMOS DE MINIMIZACIÓN
No derivables: no utilizan derivadas de la energía
Derivables: calculan derivadas primeras y/o segundas
3
3
En general, los algoritmos de minimización localizan el mínimo
más cercano (“downhill”) al punto de partida.
DERIVADAS DE LA FUNCIÓN DE ENERGÍA
POTENCIAL (MM)
Bond stretching
∂V
= k i (li − l i , 0 )
∂li
armónico
∂V
= 2Dea{1−exp[−a(l −l 0)]} × exp[−a(l − l 0)]
∂li
Morse
Angle Bending
∂θ i
∂V
= ki (θ i − θ i , 0 )
∂rj
∂θi
Torsión
∂ω i
∂V
Vn
= −
nsen ( n ω i )
∂ω i
2
∂rj
4
4
No enlazantes
12
6
⎡
⎤
⎛
⎞
⎛
⎞
σ
ij
q iq j
∂V 24ε ⎢ ⎜ σ ij ⎟
⎜
⎟
⎥
−
=
−2
+
⎜
⎟
⎜
⎟
rij ⎢ ⎝ rij ⎠
rij ⎠ ⎥ 4πε o rij 2
∂rij
⎝
⎣
⎦
DERIVADAS NUMÉRICAS
∂E
∂xi
Se hace la diferencia de E entre dos conformaciones
donde se varió la coordenada xi y luego el cociente
entre ambas.
MÉTODOS NO DERIVATIVOS
Simplex
Figura geométrica con M + 1 vértices interconectados, con M la
dimensionalidad de la función de energía.
Movimiento “ameboide”:
• Reflexión.
• Expansión.
• Contracción.
5
5
Procedimiento
1) Generación de vértices
•Un vértice corresponde a la configuración inicial del sistema.
•El resto se genera incrementando en una constante a una
coordenada por vez.
2) Cálculo de los valores de la energía para cada “vértice”.
3) Aplicación de movimiento hasta llegar a un valor de energía
aceptable.
Método muy costoso computacionalmente para sistemas de
muchas coordenadas (MM).
Util cuando se está muy lejos del mínimo
6
6
Secuencial univariado
Cicla a través de cada coordenada.
Procedimiento
1) Para cada coordenada se generan dos nuevas estructuras:
1) xi +δxi
2) xi + 2δxi
2) Se calculan las energías para las nuevas estructuras.
3) Ajuste de parábola a los tres puntos y determinación de mínimo
de esta función.
4) Cambio de la coordenada a la posición del mínimo.
5) Repite el procedimiento con las siguientes coordenadas hasta
que el cambio en todas las coordenadas es suficientemente
pequeño.
7
7
Más rápido que simplex (útil en mecánica cuántica).
Problemas de convergencia cuando hay acoplamiento fuerte entre
coordenadas.
MÉTODOS DERIVATIVOS
Son los más utilizados en minimización energética.
Derivadas proveen con información útil:
1) Derivada primera: gradiente.
•
•
dirección indica la ubicación del mínimo.
magnitud indica la pendiente de la curva
2) Derivada segunda:
• Indica la curvatura de la función (importante para
determinar la naturaleza de los puntos estacionarios).
∂V
= −Fi
∂xi
La fuerza sobre cada átomo es igual y
contraria al gradiente evaluado en el mismo
átomo
La energía del sistema puede ser disminuida moviendo los
átomos en respuesta a la fuerza que actúa sobre ellos.
8
8
Los métodos derivativos consideran la función de energía
potencial en la forma de polinomio de Taylor (evaluado en xk)
V ( x) = V ( xk ) + ( x − xk )V ´(xk ) + ( x − xk ) 2 V ´´(xk ) / 2 + .....
Se suponen funciones cuadráticas y con comportamiento
armónico.
Para el caso multidimensional:
V´(xk) -Æ vector 3N dimensional :
V´´ (xk) Æ matriz 3N X 3N:
(Hessiana)
∂V
∂xi
∂ 2V
∂xi∂xj
Un método derivativo se clasifica por su orden:
• Orden cero (no derivativos).
• Orden uno (utiliza el gradiente).
• Orden dos (utiliza el gradiente y la hessiana).
9
9
MÉTODOS DE ORDEN 1 O PRIMER ORDEN
-Descenso más pronunciado (steepest descents)
-Gradientes conjugados (conjugate gradients)
Steepest descent
Mueve el sistema en la dirección paralela a la fuerza neta:
sk = -gk/|gk|
gk = gradiente
sk = vector multidimensilnal de
norma 1 que representa la
dirección del movimiento.
Luego de elegida la dirección, hace falta asignar la magnitud del
movimiento:
xk+1 = xk +λksk
λk es el tamaño del movimiento
10
10
La dirección del gradiente está determinada en gran medida por
las fuerzas interatómicas de mayor magnitud:
Steepest descent es eficiente en eliminar los conflictos estéricos
más importantes.
Funciona bien lejos del mínimo.
Gradientes conjugados
Direcciones son conjugadas y no ortogonales como en steepest
descents.
Si la función es cuadrática, el mínimo se alcanza en N pasos
N = número de variables
Como en steepest descent, se calcula el gradiente en cada paso
La dirección en cambio se calcula como:
vk = -gk + γkvk-1
gk = gradiente en la iteración k
vk-1 = dirección del paso k-1
γk es un escalar dado por:
γk = |gk.gk|/|gk-1.gk-1|
Fletcher-Reeves
γk = |(gk-gk-1).gk|/|gk-1.gk-1
Polak-Ribiere
Gradientes conjugados localiza más fácilmente el mínimo que
steepest descent.
11
11
MÉTODOS DE ORDEN 2 O SEGUNDO ORDEN
Cálculan las derivadas primera y segunda para obtener el mínimo.
Newton-Raphson
Una dimensión:
x∗ = xk −V´(xk) / V´´(xk)
Generalización para una función multidimensional:
x∗ = xk − V´(xk) V´´− 1(xk)
Xk: coordenadas de partida
X*= coordenadas en el mínimo
Ventajas
•Determina inequívocamente el carácter de un punto estacionario.
•Para funciones cuadráticas, el mínimo se halla exactamente en
un paso desde cualquier punto de partida
Desventajas
•Necesidad de calcular (e invertir!) la matriz Hessiana.
•Lejos del mínimo, la aproximación harmónica puede no
cumplirse Æ el método falla. Suele usarse otro método hasta estar
cerca del mínimo y luego Newton-Raphson.
12
12
ALTERNATIVAS AL MÉTODO NEWTON-RAPHSON
Se han descrito una serie de métodos para agilizar la evaluación
computacional de la matriz Hessiana, principal inconveniente del
método Newton-Raphson.
• Cálculo de la matriz cada n pasos.
• Método N-R de diaginal por bloques.
• Métodos cuasi-newtonianos (cuasi segundo orden).
Diagonal por bloques
Supone que sólo un átomo se mueva a la vez (un átomo por
iteración). Como consecuencia, todos los términos de derivadas
segundas que involucran coordenadas de dos átomos distintos
valdrán cero.
∂ 2V = 0
∂xi∂xj
Si xi y xj no pertenecen al mismo
átomo
En cada iteración sólo tenemos que invertir una matriz de 3x3!!!
13
13
Cuasi-Newton
Esta categoría engloba a una serie de métodos los cuales van
construyendo gradualmente el inverso de la matriz Hessiana en
iteraciones sucesivas.
Se genera una secuencia de matrices Hk tal que:
lim Hk = V´´ -1
kÆ ∞
En cada iteración, las nuevas posiciones se obtienen según:
xk+1 = xk –Hkgk
gk = gradiente
Hk = aproximación al inverso de la
Hessiana.
y se deriva un nueva valor de H:
(xk + 1 − xk ) ⊗ (xk − 1 − xk )
−
Hk + 1 = Hk +
(xk + 1 − xk ).(gk + 1 − gk )
[Hk.(gk + 1 − gk )] ⊗ [Hk.(gk + 1 − gk )]
(gk + 1 − gk ).Hk.(gk + 1 − gk )
Davidon-Fletcher-Powell (DFP)
La matriz H se inicializa como I (la matriz unidad) aunque puede
mejorarse la performance con un “initial guess” de H calculado
con MM o semiempricos.
14
14
ELECCIÓN DEL MÉTODO DE MINIMIZACIÓN
Influyen diversos factores:
1) Capacidad de almacenamiento de datos.
2) Velocidad requerida en el cálculo.
3) Disponibilidad de derivadas analíticas.
4) Tamaño del sistema a modelar.
MECÁNICA MOLECULAR
• Steepest Descents y Gradientes Conjugados para sistemas de
porte mediano a grande.
• Newton-Raphson para moléculas más pequeñas y cercanas al
mínimo.
Punto crítico: evaluación de derivadas primeras y segundas (hay
muchos átomos)
QUÍMICA CUÁNTICA
• Newton-Raphson para niveles bajos de teoría.
• Cuasi-Newton para niveles más altos.
• Steepest Descents y Gradientes Conjugados para semiempíricos
Punto crítico: cálculo de la energía.
15
15
DISTINCIÓN ENTRE MÍNIMO, MÁXIMO Y PUNTOS
DE ENSILLADURA
∂f
=0
∂x i
en los tres casos
Mínimo: todos los valores propios de la Hessiana son positivos.
Máximo: todos los valores de la Hessiana son negativos.
Punto de ensilladura de orden n: existen n valores propios
negativos en la Hessiana.
(Un ESTADO DE TRANSICIÓN es un punto de ensilladura de
orden 1: máximo en la coordenada que conecta dos mínimos pero
mínimo en las demás).
Criterios de convergencia en minimización de energía
Imposibilidad de llegar al mínimo exacto en la práctica -Æ
Aproximación al mínimo
Monitoreo en:
•Diferencia de energía.
•Cambio en las coordenadas.
•Cambio en el gradiente (rms= sqrt(sum(gi2/3N))).
•Valor máximo del gradiente.
16
16
APLICACIONES DEL MÉTODO DE MINIMIZACIÓN DE
ENERGÍA
1) Preparación de sistemas para Monte Carlo.
2) Preparación de sistemas para Dinámica Molecular.
3) Análisis de los modos normales de vibración.
4) Determinación de estructuras de transición y caminos de
reacción.
DETERMINACIÓN DE ESTRUCTURAS DE TRANSICIÓN
Y CAMINOS DE REACCIÓN
“CAMINO DE REACCIÓN”: camino entre dos mínimos
cualesquiera. El uso de la palabra “reacción” no necesariamente
implica ruptura o formación de enlaces.
No es lo mismo una estructura de transición que un estado de
transición.
ESTRUCTURA DE TRANSICIÓN: punto de mayor energía
potencial a lo largo del camino.
ESTADO DE TRANSICIÓN: geometría en el pico de la curva
de energía libre.
17
17
MÉTODOS PARA LOCALIZAR PUNTOS DE
ENSILLADURA
Grid search
• Se utiliza una grilla o malla para escanear la superficie de
energía potencial de modo de localizar la posición aproximada de
la estructura de transición.
• Las coordenadas son variadas sistemáticamente para generar un
set de estructuras a las cuales se les mide la energía.
• Se ajustan los puntos a una expresión analítica y se predice el
punto de ensilladura por métodos de cálculo estandar.
• Restringido a sistemas con pocos grados de libertad.
• No provee directamente con la estructura de transición.
Mapeo adiabático
Método de tipo “deriva de coordenadas”
Estimación de barreras energéticas y estructuras de transición
durante transiciones conformacionales (“torsion angle driving”)
Ej: rotación de anillo de Phe o Tyr en proteína
18
18
Ciclo de mapeo adiabático
1) Se varía la coordenada relevante en pequeños incrementos.
2) Se permite la relajación (minimización energética) del entorno.
Reaction coordinate distance method
• Se parte de las estructuras del mínimo (reactivo y producto). A y
B son las especies en cuestión.
• Se define R (distancia de coordenada de reacción).
• Se minimiza el valor de R modificando las estructuras. Se supone
que cuando R es lo duficientemente pequeña, ambas estructuras
deben ser una aproximación razonable a la estructura de transición.
CAMINO DE REACCIÓN
Intrinsic reaction coordinate (IRC)
Es el camino que seguiría una partícula moviéndose a lo largo del
camino de descenso más pronunciado usando un paso
infinitesimal desde la estructura de transición hacia cada mínimo.
Las direcciones iniciales hacia los mínimos se obtienen a partir
del vector propio asociado a la estructura de transición (vector
propio asociado al único valor propio negativo de la Hessiana).
19
19
Geometría vieja
Función de onda
Ciclo
SCF
Ciclo de
optimización
de geometría
Energía + gradiente
Geometría nueva
20
20