Minimización de energía

En el campo de la química computacional , la minimización de energía (también llamada optimización de energía , minimización de geometría u optimización de geometría ) es el proceso de encontrar una disposición en el espacio de una colección de átomos donde, según algún modelo computacional de enlace químico, la fuerza interatómica neta sobre cada átomo sea aceptablemente cercana a cero y la posición en la superficie de energía potencial (PES) sea un punto estacionario (descrito más adelante). La colección de átomos podría ser una sola molécula , un ion , una fase condensada , un estado de transición o incluso una colección de cualquiera de estos. El modelo computacional de enlace químico podría, por ejemplo, ser la mecánica cuántica.

A modo de ejemplo, al optimizar la geometría de una molécula de agua , se intenta obtener las longitudes de enlace hidrógeno-oxígeno y el ángulo de enlace hidrógeno-oxígeno-hidrógeno que minimizan las fuerzas que de otro modo unirían o separarían los átomos.

La motivación para realizar una optimización geométrica es el significado físico de la estructura obtenida: las estructuras optimizadas a menudo corresponden a una sustancia tal como se encuentra en la naturaleza y la geometría de dicha estructura se puede utilizar en una variedad de investigaciones experimentales y teóricas en los campos de la estructura química , termodinámica , cinética química , espectroscopia y otros.

Por lo general, pero no siempre, el proceso busca encontrar la geometría de una disposición particular de los átomos que representa un mínimo de energía local o global. En lugar de buscar el mínimo de energía global, podría ser deseable optimizar hasta un estado de transición , es decir, un punto de silla en la superficie de energía potencial. ^[1] Además, ciertas coordenadas (como la longitud de un enlace químico) podrían fijarse durante la optimización.

Geometría molecular e interpretación matemática

La geometría de un conjunto de átomos se puede describir mediante un vector de posiciones de los átomos. Este podría ser el conjunto de coordenadas cartesianas de los átomos o, si se consideran moléculas, podrían ser las llamadas coordenadas internas formadas a partir de un conjunto de longitudes de enlace, ángulos de enlace y ángulos diedros.

Dado un conjunto de átomos y un vector, $r$ , que describe las posiciones de los átomos, se puede introducir el concepto de energía como función de las posiciones, $E$ $($ $r$ $)$ . La optimización geométrica es entonces un problema de optimización matemática , en el que se desea encontrar el valor de $r$ para el cual $E$ $($ $r$ $)$ está en un mínimo local , es decir, la derivada de la energía con respecto a la posición de los átomos, $\partial$ $E$ $/\partial$ $r$ , es el vector cero y la matriz de segunda derivada del sistema, , también conocida como matriz hessiana , que describe la curvatura del PES en $r$ , tiene todos los valores propios positivos (es definida positiva ). ${\begin{pmatrix}{\frac {\partial ^{2}E}{\partial r_{i}\,\partial r_{j}}}\end{pmatrix}}_{ij}$

Un caso especial de optimización de geometría es la búsqueda de la geometría de un estado de transición ; esto se analiza a continuación.

El modelo computacional que proporciona una $E$ $($ $r$ $)$ aproximada podría basarse en la mecánica cuántica (usando la teoría funcional de la densidad o métodos semiempíricos ), campos de fuerza o una combinación de ellos en el caso de QM/MM . Utilizando este modelo computacional y una estimación inicial (o ansatz ) de la geometría correcta, se sigue un procedimiento de optimización iterativo, por ejemplo:

Calcular la fuerza sobre cada átomo (es decir, $-\partial E /\partial r$ )
Si la fuerza es menor que un cierto umbral, terminar
de lo contrario, mueva los átomos mediante algún paso calculado $∆ r$ que se predice que reducirá la fuerza
repetir desde el principio

Aspectos prácticos de la optimización

Como se describió anteriormente, se puede utilizar algún método como la mecánica cuántica para calcular la energía, $E$ $($ $r$ $)$ , el gradiente del PES, es decir, la derivada de la energía con respecto a la posición de los átomos, $\partial$ $E$ $/\partial$ $r$ y la matriz de segunda derivada del sistema, $\partial\partial$ $E$ $/\partial$ $r$ $i$ $\partial$ $r$ $j$ , también conocida como matriz hessiana , que describe la curvatura del PES en $r$ .

Un algoritmo de optimización puede utilizar algunos o todos los $E$ $($ $r$ $)$ , $\partial$ $E$ $/\partial$ $r$ y $\partial\partial$ $E$ $/\partial$ $r$ $i$ $\partial$ $r$ $j$ para intentar minimizar las fuerzas y esto podría en teoría ser cualquier método como el descenso de gradiente, gradiente conjugado o el método de Newton, pero en la práctica, los algoritmos que utilizan el conocimiento de la curvatura PES, es decir, la matriz de Hesse, resultan ser superiores. Sin embargo, para la mayoría de los sistemas de interés práctico, puede ser prohibitivamente costoso calcular la matriz de la segunda derivada, y se estima a partir de valores sucesivos del gradiente, como es típico en una optimización Quasi-Newton .

La elección del sistema de coordenadas puede ser crucial para realizar una optimización exitosa. Las coordenadas cartesianas, por ejemplo, son redundantes ya que una molécula no lineal con $N$ átomos tiene $3$ $N$ $-6$ grados de libertad vibracional mientras que el conjunto de coordenadas cartesianas tiene $3$ $N$ dimensiones. Además, las coordenadas cartesianas están altamente correlacionadas, es decir, la matriz de Hesse tiene muchos términos no diagonales que no están cerca de cero. Esto puede conducir a problemas numéricos en la optimización, porque, por ejemplo, es difícil obtener una buena aproximación a la matriz de Hesse y calcularla con precisión es demasiado costoso computacionalmente. Sin embargo, en caso de que la energía se exprese con campos de fuerza estándar, se han desarrollado métodos computacionalmente eficientes ^[2] capaces de derivar analíticamente la matriz de Hesse en coordenadas cartesianas mientras se preserva una complejidad computacional del mismo orden que la de los cálculos de gradiente. Las coordenadas internas tienden a estar menos correlacionadas pero son más difíciles de configurar y puede ser difícil describir algunos sistemas, como aquellos con simetría o grandes fases condensadas. ^[3] Muchos paquetes de software de química computacional modernos contienen procedimientos automáticos para la generación automática de sistemas de coordenadas razonables para la optimización. ^[4]

Restricción del grado de libertad

Se pueden eliminar algunos grados de libertad de una optimización; por ejemplo, se pueden asignar valores fijos a las posiciones de los átomos o a las longitudes y ángulos de los enlaces. A veces, se los denomina grados de libertad congelados .

La figura 1 muestra una optimización geométrica de los átomos de un nanotubo de carbono en presencia de un campo electrostático externo. En esta optimización, los átomos de la izquierda tienen sus posiciones congeladas. Su interacción con los demás átomos del sistema se sigue calculando, pero se evita la alteración de la posición de los átomos durante la optimización.

Figura 1: Deflexiones electrostáticas de un nanotubo de carbono en un campo eléctrico .

Optimización del estado de transición

Las estructuras de estados de transición se pueden determinar buscando puntos de silla en el PES de las especies químicas de interés. ^[5] Un punto de silla de primer orden es una posición en el PES que corresponde a un mínimo en todas las direcciones excepto una; un punto de silla de segundo orden es un mínimo en todas las direcciones excepto dos, y así sucesivamente. Definido matemáticamente, un punto de silla de orden n se caracteriza por lo siguiente: $\partial E /\partial r = 0$ y la matriz hessiana, $\partial\partial E /\partial r i \partial r j$ , tiene exactamente n valores propios negativos.

Los algoritmos para localizar geometrías de estados de transición se dividen en dos categorías principales: métodos locales y métodos semiglobales. Los métodos locales son adecuados cuando el punto de partida de la optimización está muy cerca del estado de transición verdadero ( en breve se definirá " muy cerca" ) y los métodos semiglobales se aplican cuando se busca localizar el estado de transición con muy poco conocimiento a priori de su geometría. Algunos métodos, como el método Dimer (ver más abajo), se incluyen en ambas categorías.

Búsquedas locales

Una optimización local requiere una estimación inicial del estado de transición que sea muy cercana al estado de transición verdadero. Muy cercana significa que la estimación inicial debe tener una matriz hessiana correspondiente con un valor propio negativo, o bien, el valor propio negativo correspondiente a la coordenada de reacción debe ser mayor en magnitud que los otros valores propios negativos. Además, el vector propio con el valor propio más negativo debe corresponder a la coordenada de reacción, es decir, debe representar la transformación geométrica relacionada con el proceso cuyo estado de transición se busca.

Dados los requisitos previos anteriores, un algoritmo de optimización local puede entonces moverse "cuesta arriba" a lo largo del vector propio con el valor propio más negativo y "cuesta abajo" a lo largo de todos los demás grados de libertad, utilizando algo similar a un método cuasi-Newton.

Método del dímero

El método del dímero ^[6] se puede utilizar para encontrar posibles estados de transición sin conocer la estructura final o para refinar una buena estimación de una estructura de transición. El “dímero” está formado por dos imágenes muy cercanas entre sí en el PES. El método funciona moviendo el dímero hacia arriba desde la posición inicial mientras se gira para encontrar la dirección de curvatura más baja (en última instancia negativa).

Técnica de Relajación Activada (TRA)

La técnica de relajación de activación (ART) ^[7]^[8]^[9] también es un método abierto para encontrar nuevos estados de transición o refinar puntos de silla conocidos en el PES. El método sigue la dirección de la curvatura negativa más baja (calculada utilizando el algoritmo de Lanczos ) en el PES para alcanzar el punto de silla, relajándose en el hiperplano perpendicular entre cada "salto" (activación) en esta dirección.

Métodos de cadena de estados

Los métodos de cadena de estados ^[10] se pueden utilizar para encontrar la geometría aproximada del estado de transición en función de las geometrías del reactivo y del producto. La geometría aproximada generada puede servir entonces como punto de partida para el refinamiento mediante una búsqueda local, que se describió anteriormente.

Los métodos de cadena de estados utilizan una serie de vectores, es decir, puntos en el PES, que conectan el reactivo y el producto de la reacción de interés, $r$ $reactivo$ y $r$ $producto$ , discretizando así la vía de reacción. Muy comúnmente, estos puntos se conocen como perlas debido a una analogía de un conjunto de perlas conectadas por cuerdas o resortes, que conectan el reactivo y los productos. La serie de perlas a menudo se crea inicialmente interpolando entre $r$ $reactivo$ y $r$ $producto$ , por ejemplo, para una serie de $N$ $+ 1$ perlas, la perla $i$ podría estar dada por

$\mathbf {r} _{i}={\frac {i}{N}}\mathbf {r} _{\mathrm {producto} }+\left(1-{\frac {i}{N }}\right)\mathbf {r} _{\mathrm {reactante} }$

donde $i \in 0, 1, ..., N$ . Cada una de las perlas $r i$ tiene una energía, $E$ $($ $r$ $i$ $)$ , y fuerzas, $-\partial$ $E$ $/\partial$ $r$ $i$ y estas se tratan con un proceso de optimización restringida que busca obtener una representación lo más precisa posible de la vía de reacción. Para lograr esto, se deben aplicar restricciones de espaciado de modo que cada perla $r$ $i$ no se optimice simplemente para la geometría del reactivo y del producto.

A menudo, esta restricción se logra proyectando los componentes de la fuerza sobre cada perla $r i$ , o alternativamente el movimiento de cada perla durante la optimización, que son tangenciales a la trayectoria de reacción. Por ejemplo, si por conveniencia se define que $g i = \partial E /\partial r i$ , entonces el gradiente de energía en cada perla menos el componente del gradiente de energía que es tangencial a la trayectoria de reacción viene dado por

$\mathbf {g} _{i}^{\perp }=\mathbf {g} _{i}-\mathbf {\tau } _{i}(\mathbf {\tau } _{i}\cdot \mathbf {g} _{i})=\left(I-\mathbf {\tau } _{i}\mathbf {\tau } _{i}^{T}\right)\mathbf {g} _{i}$

donde $I$ es la matriz identidad y $τ i$ es un vector unitario que representa la tangente de la trayectoria de reacción en $r i$ . Al proyectar los componentes del gradiente de energía o el paso de optimización que son paralelos a la trayectoria de reacción, un algoritmo de optimización reduce significativamente la tendencia de cada una de las perlas a optimizarse directamente al mínimo.

Tránsito sincrónico

El método de cadena de estados más simple es el método de tránsito sincrónico lineal (LST). Funciona tomando puntos interpolados entre las geometrías del reactivo y del producto y eligiendo el que tiene la energía más alta para un refinamiento posterior mediante una búsqueda local. El método de tránsito sincrónico cuadrático (QST) extiende el LST al permitir una trayectoria de reacción parabólica, con optimización del punto de energía más alta ortogonalmente a la parábola.

Banda elástica con empujoncitos

En el método de banda elástica empujada (NEB) ^[11] , las perlas a lo largo de la vía de reacción tienen fuerzas de resorte simuladas además de las fuerzas químicas, $-\partial E /\partial r i$ , para hacer que el optimizador mantenga la restricción de espaciado. Específicamente, la fuerza $f i$ en cada punto i está dada por

$\mathbf {f} _{i}=\mathbf {f} _{i}^{\paralelo }-\mathbf {g} _{i}^{\perp }$

dónde

$\mathbf {f} _{i}^{\paralelo }=k\left[\left(\left(\mathbf {r} _{i+1}-\mathbf {r} _{i}\right)-\left(\mathbf {r} _{i}-\mathbf {r} _{i-1}\right)\right)\cdot \tau _{i}\right]\tau _{i}$

es la fuerza del resorte paralela a la trayectoria en cada punto $r i$ ( k es una constante del resorte y $τ i$ , como antes, es un vector unitario que representa la tangente de la trayectoria de reacción en $r i$ ).

En una implementación tradicional, el punto con la energía más alta se utiliza para el refinamiento posterior en una búsqueda local. Existen muchas variaciones del método NEB, ^[12] como la NEB de imagen ascendente, en la que el punto con la energía más alta se empuja hacia arriba durante el procedimiento de optimización para (con suerte) dar una geometría que sea aún más cercana a la del estado de transición. También ha habido extensiones ^[13] para incluir la regresión del proceso gaussiano para reducir el número de evaluaciones. Para sistemas con geometría no euclidiana (R^2), como los sistemas magnéticos, el método se modifica al enfoque de banda elástica empujada geodésica. ^[14]

Método de cadena

El método de cadenas ^[15]^[16]^[17] utiliza splines que conectan los puntos, $r i$ , para medir y aplicar restricciones de distancia entre los puntos y para calcular la tangente en cada punto. En cada paso de un procedimiento de optimización, los puntos pueden moverse de acuerdo con la fuerza que actúa sobre ellos perpendicularmente a la trayectoria y, luego, si la restricción de equidistancia entre los puntos ya no se satisface, los puntos pueden redistribuirse, utilizando la representación spline de la trayectoria para generar nuevos vectores con el espaciado requerido.

Las variaciones del método de cadenas incluyen el método de cadena creciente, ^[18] en el que la estimación de la ruta aumenta a partir de los puntos finales (es decir, el reactivo y los productos) a medida que avanza la optimización.

Comparación con otras técnicas

La optimización geométrica es fundamentalmente diferente de una simulación de dinámica molecular . Esta última simula el movimiento de las moléculas con respecto al tiempo, sujetas a la temperatura, fuerzas químicas, velocidades iniciales, movimiento browniano de un disolvente, etc., mediante la aplicación de las leyes de Newton del movimiento . Esto significa que las trayectorias de los átomos que se calculan tienen algún significado físico. La optimización geométrica, por el contrario, no produce una "trayectoria" con ningún significado físico; se ocupa de la minimización de las fuerzas que actúan sobre cada átomo de un conjunto de átomos, y la vía por la que se logra esto carece de significado. Diferentes algoritmos de optimización podrían dar el mismo resultado para la estructura de energía mínima, pero llegar a ella por una vía diferente.

Véase también

Gráfico compuesto con restricciones
Cortes de gráficos en visión artificial : aparatos para resolver problemas de visión artificial que pueden formularse en términos de minimización de energía
Principios de energía en la mecánica estructural

Referencias

^ "Referencia de entrada de la versión troncal de CP2K, Sección GEO_OPT, Palabra clave TYPE". CP2K . Consultado el 30 de abril de 2015 .
^ Chatzieleftheriou, S.; Adendorff, MR; Lagaros, ND (2016). "Elementos finitos de energía potencial generalizados para modelar nanoestructuras moleculares". J. Chem. Inf. Model . 56 (10): 1963–1978. doi :10.1021/acs.jcim.6b00356. PMID 27653992.
^ Peng, C.; Ayala, PY; Schlegel, HB (1996). "Uso de coordenadas internas redundantes para optimizar geometrías de equilibrio y estados de transición". Journal of Computational Chemistry . 17 (1): 49–56. doi :10.1002/(sici)1096-987x(19960115)17:1<49::aid-jcc5>3.3.co;2-#.
^ "Inicio". gaussian.com .
^ Frank Jensen (1999). Introducción a la química computacional . Inglaterra: John Wiley and Sons Ltd.
^ Graeme Henkelman; Hannes Jónsson (1999). "Un método de dímeros para encontrar puntos de silla en superficies de potencial de alta dimensión utilizando solo las primeras derivadas". J. Chem. Phys . 111 (15): 7010–7022. Bibcode :1999JChPh.111.7010H. doi :10.1063/1.480097.
^ GT Barkema; Normand Mousseau (1996). "Relajación basada en eventos de sistemas desordenados continuos". Phys. Rev. Lett . 77 (21): 4358–4361. arXiv : cond-mat/9607156 . Código Bibliográfico :1996PhRvL..77.4358B. doi :10.1103/PhysRevLett.77.4358. PMID 10062518. S2CID 27932059.
^ Rachid Malek; Normand Mousseau (2011). "Exploración optimizada del paisaje energético utilizando la técnica de activación-relajación basada en ab initio". Physical Review E . 135 (6): 7723–7728. arXiv : cond-mat/0006042 . Bibcode :2000PhRvE..62.7723M. doi :10.1103/PhysRevE.62.7723. PMID 11138044. S2CID 119453527.
^ Eduardo Machado-Charry; Laurent Karim Béland; Damien Caliste; Luigi Genovese; Thierry Deutsch; Normand Mousseau; Pascal Pochet (2011). "Exploración optimizada del paisaje energético utilizando la técnica de activación-relajación basada en ab initio". J. Chem. Phys . 62 (3): 034102–034112. Bibcode :2011JChPh.135c4102M. doi :10.1063/1.3609924. PMID 21786982.
^ Jensen, F. Introducción a la química computacional; Wiley: 2.ª edición; 2006
^ (a) G. Mills y H. Jónsson, Phys. Rev. Lett. 72, 1124 (1994) (b) Graeme Henkelman y Hannes Jónsson, Estimación de tangente mejorada en el método de banda elástica empujada para encontrar trayectorias de energía mínima y puntos de silla, J. Chem. Phys. 113, 9978 - 9985 (2000)
^ "Banda elástica con empujoncitos". UT Austin . Archivado desde el original el 3 de febrero de 2014.
^ Koistinen, Olli-Pekka; Dagbjartsdóttir, Freyja B.; Ásgeirsson, Vilhjálmur; Vehtari, Aki; Jónsson, Hannes (21 de octubre de 2017). "Los cálculos de bandas elásticas empujadas se aceleraron con la regresión del proceso gaussiano". La Revista de Física Química . 147 (15): 152720. arXiv : 1706.04606 . Código Bib :2017JChPh.147o2720K. doi : 10.1063/1.4986787. ISSN 0021-9606. PMID 29055305. S2CID 21822734.
^ Ivanov, AV; Dagbartsson, D; Tranchida, J; Uzdin, VM; Jónsson, H (12 de agosto de 2020). "Método de optimización eficiente para encontrar caminos de energía mínima de transiciones magnéticas". Journal of Physics: Condensed Matter . 32 (34): 345901. arXiv : 2001.10372 . Bibcode :2020JPCM...32H5901I. doi :10.1088/1361-648X/ab8b9c. ISSN 0953-8984. PMID 32316000. S2CID 210932577.
^ "Eventos raros, vías de transición y tasas de reacción".y "La página del método de cadena".
^ Weinan E, Weiqing Ren, Eric Vanden-Eijnden (2002). "Método de cuerdas para el estudio de eventos raros". Phys. Rev. B . 66 (5): 052301. arXiv : cond-mat/0205527 . Código Bibliográfico :2002PhRvB..66e2301E. doi :10.1103/PhysRevB.66.052301. S2CID 119326534.{{cite journal}}: CS1 maint: varios nombres: lista de autores ( enlace )
^ Amit Samanta; Weinan E (2010). "Método de cuerdas modificado para encontrar la ruta de energía mínima". arXiv : 1009.5612 [physics.comp-ph].
^ Baron Peters; Andreas Heyden; Alexis T. Bell; Arup Chakraborty (2004). "Un método de cuerdas crecientes para determinar estados de transición: Comparación con los métodos de cuerdas y bandas elásticas empujadas". J. Chem. Phys . 120 (17): 7877–7886. Bibcode :2004JChPh.120.7877P. doi :10.1063/1.1691018. PMID 15267702.

Enlaces externos

Recetas numéricas en Fortran 77

Referencias adicionales

Payne et al., "Técnicas de minimización iterativa para cálculos de energía total ab initio: dinámica molecular y gradientes conjugados", Reviews of Modern Physics 64 (4), pp. 1045–1097. (1992) (resumen)
Stich et al., "Minimización del gradiente conjugado del funcional energético: un nuevo método para el cálculo de la estructura electrónica", Physical Review B 39 (8), págs. 4997–5004, (1989)
Chadi, "Enfoque de minimización de energía para la geometría atómica de superficies de semiconductores", Physical Review Letters 41 (15), págs. 1062-1065 (1978)