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 interconexión neta -La fuerza atómica sobre cada átomo es aceptablemente cercana a cero y la posición en la superficie de energía potencial (PES) es un punto estacionario (que se describe más adelante). El conjunto de átomos puede ser una sola molécula , un ion , una fase condensada , un estado de transición o incluso un conjunto de cualquiera de estos. El modelo computacional de enlace químico podría ser, por ejemplo, la mecánica cuántica.

Por ejemplo, al optimizar la geometría de una molécula de agua , se pretende obtener las longitudes de los enlaces hidrógeno-oxígeno y el ángulo del enlace hidrógeno-oxígeno-hidrógeno que minimicen las fuerzas que de otro modo juntarían o separarían los átomos.

La motivación para realizar una optimización de la geometría 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 estructura química , termodinámica , cinética química , espectroscopia y otros.

Normalmente, pero no siempre, el proceso busca encontrar la geometría de una disposición particular de los átomos que represente un mínimo de energía local o global. En lugar de buscar un mínimo de energía global, podría ser deseable optimizar a 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 las posiciones de los átomos. Este podría ser el conjunto de coordenadas cartesianas de los átomos o, cuando 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 diédricos.

Dado un conjunto de átomos y un vector, $r$ , que describe las posiciones de los átomos, se puede introducir el concepto de energía en función de las posiciones, $E$ $($ $r$ $)$ . La optimización de la geometría 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 segunda derivada del sistema, también conocida como matriz de Hesse , que describe la curvatura del PES en $r$ , tiene todos 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 la 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 suposició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 algún umbral, termine
de lo contrario, mueva los átomos en 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 segunda derivada del sistema, $\partial\partial$ $E$ $/\partial$ $r$ $i$ $\partial$ $r$ $j$ , también conocida como matriz de Hesse , que describe la curvatura del PES en $r$ .

Un algoritmo de optimización puede utilizar parte o la totalidad de $E$ $($ $r$ $)$ , $\partial$ $E$ $/\partial$ $r$ y $\partial\partial$ $E$ $/\partial$ $r$ $i$ $\partial$ $r$ $j$ para intentar minimizar las fuerzas y, en teoría, esto podría ser cualquier método como descenso de gradiente, conjugado gradiente 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 resultar 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 Cuasi-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 vibratorios , 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 desde el punto de vista computacional. Sin embargo, en el 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 preservando al mismo tiempo 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 resultar difícil describir algunos sistemas, como los que tienen 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 su optimización. ^[4]

Grado de restricción de libertad

Algunos grados de libertad se pueden eliminar de una optimización, por ejemplo, a las posiciones de los átomos o a las longitudes y ángulos de los enlaces se les pueden dar valores fijos. A veces se hace referencia a estos como grados de libertad congelados .

La Figura 1 muestra una optimización de la geometría de los átomos en 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 otros átomos del sistema aún se calcula, 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 los 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 silla de primer orden es una posición en el PES correspondiente a un mínimo en todas las direcciones excepto una; un punto silla de segundo orden es un mínimo en todas las direcciones excepto en 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 de Hesse, $\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 para la optimización está muy cerca del verdadero estado de transición ( muy cerca se definirá en breve) y los métodos semiglobales encuentran aplicación cuando se busca localizar el estado de transición con muy poco conocimiento a priori de él. su geometría. Algunos métodos, como el método Dimer (ver más abajo), se incluyen en ambas categorías.

Búsquedas locales

La denominada optimización local requiere una estimación inicial del estado de transición que esté muy cerca del estado de transición real. Muy cerca generalmente significa que la suposición inicial debe tener una matriz de Hesse correspondiente con un valor propio negativo, o 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 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 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 suposició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 cuesta arriba desde la posición inicial mientras lo gira para encontrar la dirección de curvatura más baja (en última instancia, negativa).

Técnica de Activación y Relajación (ART)

La Técnica de Activación y Relajació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

Se pueden utilizar métodos de cadena de estados ^{[10] 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, como 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 ruta de reacción. Muy comúnmente, estos puntos se denominan cuentas debido a una analogía de un conjunto de cuentas conectadas por hilos o resortes, que conectan el reactivo y los productos. La serie de cuentas a menudo se crea inicialmente interpolando entre $r$ $reactivo$ y $r$ $producto$ ; por ejemplo, para una serie de $N$ $+ 1$ cuentas, la cuenta $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 cuentas $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 restringido que busca obtener una representación lo más precisa posible de la ruta de reacción. Para lograr esto, se deben aplicar restricciones de espacio para que cada perla $ri$ $no$ se optimice simplemente según la geometría del reactivo y del producto.

A menudo, esta restricción se logra proyectando componentes de la fuerza en cada cordón $ri$ , o alternativamente el movimiento de cada cordón 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 cuenta menos el componente del gradiente de energía que es tangencial a la ruta 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 del camino de reacción en $r i$ . Al proyectar componentes del gradiente de energía o del paso de optimización que son paralelos a la ruta 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 síncrono lineal (LST). Opera tomando puntos interpolados entre las geometrías del reactivo y del producto y eligiendo el que tiene mayor energía para su posterior refinamiento mediante una búsqueda local. El método del tránsito síncrono cuadrático (QST) extiende el LST al permitir una trayectoria de reacción parabólica, con optimización del punto de mayor energía ortogonalmente a la parábola.

Banda elástica empujada

En el método de banda elástica empujada (NEB) ^[11] , las perlas a lo largo de la ruta de reacción han simulado fuerzas de resorte 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}^{\parallel }-\mathbf {g} _{i}^{\perp }$

dónde

$\mathbf {f} _{i}^{\parallel }=k\left[\left(\left(\mathbf {r} _{i+1}-\mathbf {r} _{i}\ derecha)-\izquierda(\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 trayectoria de reacción tangente en $r i$ ).

En una implementación tradicional, el punto con mayor energía se utiliza para su posterior refinamiento en una búsqueda local. Hay muchas variaciones del método NEB, ^[12] como la imagen de ascenso NEB, en la que el punto con mayor energía se empuja hacia arriba durante el procedimiento de optimización para (con suerte) dar una geometría aún más cercana a la de el 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 con empuje geodésico. ^[14]

método de cadena

El método de cuerda ^[15]^[16]^[17] utiliza splines que conectan los puntos, $r i$ , para medir y aplicar restricciones de distancia entre los puntos y 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 cumple, los puntos pueden redistribuirse utilizando la spline. representación del camino para generar nuevos vectores con el espaciado requerido.

Las variaciones del método de cadena incluyen el método de cadena creciente, ^[18] en el que la estimación de la ruta aumenta desde 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 de la geometría es fundamentalmente diferente de una simulación de dinámica molecular . Este último 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 del movimiento de Newton . Esto significa que las trayectorias de los átomos que se calculan tienen algún significado físico. La optimización de la geometría, 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 en una colección de átomos, y el camino a través del cual 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 él por un camino diferente.

Ver también

Gráfico compuesto de restricciones
Recortes de gráficos en visión por computadora : aparato para resolver problemas de visión por computadora que se pueden formular en términos de minimización de energía
Principios energéticos en mecánica estructural.

Referencias

^ "Referencia de entrada de la troncal de la versión CP2K, sección GEO_OPT, TIPO de palabra clave". CP2K . Consultado el 30 de abril de 2015 .
^ Chatzieleftheriou, S.; Adendorff, señor; Lagaros, Dakota del Norte (2016). "Elementos finitos de energía potencial generalizada para el modelado de nanoestructuras moleculares". J. química. inf. Modelo . 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". Revista de Química Computacional . 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ímero para encontrar puntos de silla en superficies potenciales de alta dimensión utilizando sólo primeras derivadas". J. química. Física . 111 (15): 7010–7022. Código Bib :1999JChPh.111.7010H. doi : 10.1063/1.480097.
^ GT Barkema; Normand Mousseau (1996). "Relajación basada en eventos de sistemas continuos desordenados". Física. 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 mediante la técnica de activación-relajación basada ab initio". Revisión física E. 135 (6): 7723–7728. arXiv : cond-mat/0006042 . Código bibliográfico : 2000PhRvE..62.7723M. doi : 10.1103/PhysRevE.62.7723. PMID 11138044. S2CID 119453527.
^ Eduardo Machado-Charry; Laurent Karim Béland; Damián Caliste; Luigi Genovese; Thierry Deutsch; Normand Mousseau; Pascal Pochet (2011). "Exploración optimizada del paisaje energético mediante la técnica de activación-relajación basada ab initio". J. química. Física . 62 (3): 034102–034112. Código Bib :2011JChPh.135c4102M. doi :10.1063/1.3609924. PMID 21786982.
^ Jensen, F. Introducción a la química computacional; Wiley: 2 ed.; 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 caminos de energía mínima y puntos de silla, J. Chem. Física. 113, 9978 - 9985 (2000)
^ "Banda elástica empujada". 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". Revista de Física: Materia Condensada . 32 (34): 345901. arXiv : 2001.10372 . Código Bib : 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". Física. Rev. B. 66 (5): 052301. arXiv : cond-mat/0205527 . Código Bib : 2002PhRvB..66e2301E. doi : 10.1103/PhysRevB.66.052301. S2CID 119326534.{{cite journal}}: Mantenimiento CS1: varios nombres: lista de autores ( enlace )
^ Amit Samanta; Weinan E (2010). "Método de cadena modificado para encontrar la ruta de energía mínima". arXiv : 1009.5612 [física.comp-ph].
^ Barón Peters; Andreas Heyden; Alexis T. Bell; Arup Chakraborty (2004). "Un método de cuerda en crecimiento para determinar estados de transición: comparación con los métodos de cuerda y banda elástica empujada". J. química. Física . 120 (17): 7877–7886. Código Bib : 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 ab initio de energía total: dinámica molecular y gradientes conjugados", Reviews of Modern Physics 64 (4), págs. (1992) (resumen)
Stich et al., "Minimización del gradiente conjugado de la energía funcional: 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 semiconductoras", Physical Review Letters 41 (15), págs. 1062-1065 (1978)