Hispanopedia Wiki
Articulo de referencia

Restricción (química computacional)

En química computacional , un algoritmo de restricción es un método para satisfacer el movimiento newtoniano de un cuerpo rígido compuesto por puntos de masa. Se utiliza un algo...

Hispanopedia WikiContenido en espanolLectura gratuita

En química computacional , un algoritmo de restricción es un método para satisfacer el movimiento newtoniano de un cuerpo rígido compuesto por puntos de masa. Se utiliza un algoritmo de restricción para asegurar que se mantenga la distancia entre los puntos de masa. Los pasos generales involucrados son: (i) elegir nuevas coordenadas sin restricciones (coordenadas internas), (ii) introducir fuerzas de restricción explícitas, (iii) minimizar las fuerzas de restricción implícitamente mediante la técnica de multiplicadores de Lagrange o métodos de proyección.

Los algoritmos de restricción se aplican frecuentemente a simulaciones de dinámica molecular . Si bien estas simulaciones a veces se realizan utilizando coordenadas internas que satisfacen automáticamente las restricciones de longitud, ángulo y torsión de enlace, también pueden llevarse a cabo utilizando fuerzas de restricción explícitas o implícitas para estas tres restricciones. Sin embargo, las fuerzas de restricción explícitas generan ineficiencia, ya que requieren mayor potencia computacional para obtener una trayectoria de una longitud determinada. Por lo tanto, generalmente se prefieren los algoritmos de restricción que utilizan coordenadas internas y fuerzas implícitas.

Los algoritmos de restricción logran eficiencia computacional al ignorar el movimiento a lo largo de algunos grados de libertad. Por ejemplo, en la dinámica molecular atomística, normalmente se restringe la longitud de los enlaces covalentes con el hidrógeno; sin embargo, no se deben usar algoritmos de restricción si las vibraciones a lo largo de estos grados de libertad son importantes para el fenómeno que se estudia.

Formación matemática

El movimiento de un conjunto de N partículas se puede describir mediante un conjunto de ecuaciones diferenciales ordinarias de segundo orden, la segunda ley de Newton, que se puede escribir en forma matricial.

METRO⋅d2qdt2=F=−∂V∂q{\displaystyle \mathbf {M} \cdot {\frac {d^{2}\mathbf {q} }{dt^{2}}}=\mathbf {f} =-{\frac {\partial V}{\partial \mathbf {q} }}}

donde M es una matriz de masas y q es el vector de coordenadas generalizadas que describen las posiciones de las partículas. Por ejemplo, el vector q puede ser una matriz de coordenadas cartesianas 3N de las posiciones de las partículas r k , donde k va de 1 a N ; en ausencia de restricciones, M sería la matriz cuadrada diagonal 3N x 3N de las masas de las partículas. El vector f representa las fuerzas generalizadas y el escalar V ( q ) representa la energía potencial, ambas funciones de las coordenadas generalizadas q .

Si existen M restricciones, las coordenadas también deben satisfacer M ecuaciones algebraicas independientes del tiempo.

gramoj(q)=0{\displaystyle g_{j}(\mathbf {q} )=0}

donde el índice j va de 1 a M. Para simplificar, estas funciones g i se agrupan en un vector g de dimensión M que se muestra a continuación. La tarea consiste en resolver el conjunto combinado de ecuaciones diferenciales algebraicas (EDA), en lugar de solo las ecuaciones diferenciales ordinarias (EDO) de la segunda ley de Newton.

Este problema fue estudiado en detalle por Joseph Louis Lagrange , quien estableció la mayoría de los métodos para resolverlo. [ 1 ] El enfoque más simple consiste en definir nuevas coordenadas generalizadas sin restricciones; este enfoque elimina las ecuaciones algebraicas y reduce el problema nuevamente a la resolución de una ecuación diferencial ordinaria . Este enfoque se utiliza, por ejemplo, para describir el movimiento de un cuerpo rígido; la posición y orientación de un cuerpo rígido se pueden describir mediante seis coordenadas independientes y sin restricciones, en lugar de describir las posiciones de las partículas que lo componen y las restricciones entre ellas que mantienen sus distancias relativas. La desventaja de este enfoque es que las ecuaciones pueden volverse difíciles de manejar y complejas; por ejemplo, la matriz de masa M puede volverse no diagonal y depender de las coordenadas generalizadas.

Un segundo enfoque consiste en introducir fuerzas explícitas que mantengan la restricción; por ejemplo, se podrían introducir fuerzas elásticas intensas que impongan las distancias entre los puntos de masa dentro de un cuerpo "rígido". Las dos dificultades de este enfoque son que las restricciones no se satisfacen exactamente y que las fuerzas intensas pueden requerir pasos de tiempo muy cortos, lo que hace que las simulaciones sean computacionalmente ineficientes.

Un tercer enfoque consiste en utilizar un método como los multiplicadores de Lagrange o la proyección a la variedad de restricciones para determinar los ajustes de coordenadas necesarios para satisfacer las restricciones.

Finalmente, existen diversos enfoques híbridos en los que diferentes conjuntos de restricciones se satisfacen mediante diferentes métodos, por ejemplo, coordenadas internas, fuerzas explícitas y soluciones de fuerza implícita.

Métodos de coordenadas internas

El método más sencillo para satisfacer las restricciones en la minimización de energía y la dinámica molecular consiste en representar el sistema mecánico mediante las llamadas coordenadas internas, que corresponden a grados de libertad independientes y sin restricciones del sistema. Por ejemplo, los ángulos diedros de una proteína constituyen un conjunto independiente de coordenadas que especifican las posiciones de todos los átomos sin necesidad de restricciones. La dificultad de estos métodos de coordenadas internas radica en dos aspectos: las ecuaciones de movimiento newtonianas se vuelven mucho más complejas y las coordenadas internas pueden ser difíciles de definir para sistemas cíclicos de restricciones, como en el plegamiento de anillos o cuando una proteína presenta un enlace disulfuro.

Los métodos originales para la minimización recursiva eficiente de la energía en coordenadas internas fueron desarrollados por Gō y colaboradores. [ 2 ] [ 3 ]

Los solucionadores recursivos eficientes con restricciones de coordenadas internas se extendieron a la dinámica molecular. [ 4 ] [ 5 ] Métodos análogos se aplicaron posteriormente a otros sistemas. [ 6 ] [ 7 ] [ 8 ]

métodos basados ​​en multiplicadores de Lagrange

Resolución de las restricciones de una molécula de agua rígida mediante multiplicadores de Lagrange : a) se obtienen las posiciones sin restricciones después de un paso de tiempo de simulación, b) se calculan los gradientes de cada restricción sobre cada partícula y c) se calculan los multiplicadores de Lagrange para cada gradiente de manera que se satisfagan las restricciones.

En la mayoría de las simulaciones de dinámica molecular que utilizan algoritmos de restricciones, las restricciones se imponen utilizando el método de los multiplicadores de Lagrange. Dado un conjunto de n restricciones lineales ( holonómicas ) en el tiempo t ,

σk(t):=‖incógnitakα(t)−incógnitakβ(t)‖2−dk2=0,k=1…norte{\displaystyle \sigma _{k}(t):=\|\mathbf {x} _{k\alpha }(t)-\mathbf {x} _{k\beta }(t)\|^{2}-d_{k}^{2}=0,\quad k=1\ldots n}

dóndeincógnitakα(t){\displaystyle \scriptstyle \mathbf {x} _{k\alpha }(t)}yincógnitakβ(t){\displaystyle \scriptstyle \mathbf {x} _{k\beta }(t)}son las posiciones de las dos partículas involucradas en la k -ésima restricción en el tiempo t ydk{\displaystyle d_{k}}es la distancia interpartícula prescrita.

Las fuerzas debidas a estas restricciones se añaden a las ecuaciones de movimiento, lo que resulta en, para cada una de las N partículas del sistema

∂2incógnitai(t)∂t2metroi=−∂∂incógnitai[V(incógnitai(t))−∑k=1norteλkσk(t)],i=1…norte.{\displaystyle {\frac {\partial ^{2}\mathbf {x} _{i}(t)}{\partial t^{2}}}m_{i}=-{\frac {\partial }{\partial \mathbf {x} _{i}}}\left[V(\mathbf {x} _{i}(t))-\sum _{k=1}^{n}\lambda _{k}\sigma _{k}(t)\right],\quad i=1\ldots N.}

Agregar las fuerzas de restricción no cambia la energía total, ya que el trabajo neto realizado por las fuerzas de restricción (tomado sobre el conjunto de partículas sobre las que actúan las restricciones) es cero. Tenga en cuenta que el signo enλk{\displaystyle \lambda _{k}}es arbitrario y algunas referencias [ 9 ] tienen un signo opuesto.

Al integrar ambos lados de la ecuación con respecto al tiempo, las coordenadas restringidas de las partículas en el tiempo,t+Δt{\displaystyle t+\Delta t}, se dan,

incógnitai(t+Δt)=incógnita^i(t+Δt)+∑k=1norteλk∂σk(t)∂incógnitai(Δt)2metroi−1,i=1…norte{\displaystyle \mathbf {x} _{i}(t+\Delta t)={\hat {\mathbf {x} }}_{i}(t+\Delta t)+\sum _{k=1}^{n}\lambda _{k}{\frac {\partial \sigma _{k}(t)}{\partial \mathbf {x} _{i}}}\left(\Delta t\right)^{2}m_{i}^{-1},\quad i=1\ldots N}

dóndeincógnita^i(t+Δt){\displaystyle {\hat {\mathbf {x} }}_{i}(t+\Delta t)}es la posición no restringida (o no corregida) de la i -ésima partícula después de integrar las ecuaciones de movimiento no restringidas.

Para satisfacer las restriccionesσk(t+Δt){\displaystyle \sigma _{k}(t+\Delta t)}En el siguiente paso de tiempo, los multiplicadores de Lagrange deben determinarse según la siguiente ecuación:

σk(t+Δt):=‖incógnitakα(t+Δt)−incógnitakβ(t+Δt)‖2−dk2=0.{\displaystyle \sigma _{k}(t+\Delta t):=\left\|\mathbf {x} _{k\alpha }(t+\Delta t)-\mathbf {x} _{k\beta }(t+\Delta t)\right\|^{2}-d_{k}^{2}=0.}

Esto implica resolver un sistema denorte{\displaystyle n}ecuaciones no lineales

σj(t+Δt):=‖incógnita^jα(t+Δt)−incógnita^jβ(t+Δt)+∑k=1norteλk(Δt)2[∂σk(t)∂incógnitajαmetrojα−1−∂σk(t)∂incógnitajβmetrojβ−1]‖2−dj2=0,j=1…norte{\displaystyle \sigma _{j}(t+\Delta t):=\left\|{\hat {\mathbf {x} }}_{j\alpha }(t+\Delta t)-{\hat {\mathbf {x} }}_{j\beta }(t+\Delta t)+\sum _{k=1}^{n}\lambda _{k}\left(\Delta t\right)^{2}\left[{\frac {\partial \sigma _{k}(t)}{\partial \mathbf {x} _{j\alpha }}}m_{j\alpha }^{-1}-{\frac {\partial \sigma _{k}(t)}{\partial \mathbf {x} _{j\beta }}}m_{j\beta }^{-1}\right]\right\|^{2}-d_{j}^{2}=0,\quad j=1\ldots n}

simultáneamente para elnorte{\displaystyle n}multiplicadores de Lagrange desconocidosλk{\displaystyle \lambda _{k}}.

Este sistema denorte{\displaystyle n}ecuaciones no lineales ennorte{\displaystyle n}Las incógnitas se resuelven comúnmente utilizando el método de Newton-Raphson, donde el vector soluciónλ_{\displaystyle {\underline {\lambda }}}se actualiza usando

λ_(l+1)←λ_(l)−Jσ−1σ_(t+Δt){\displaystyle {\underline {\lambda }}^{(l+1)}\leftarrow {\underline {\lambda }}^{(l)}-\mathbf {J} _{\sigma }^{-1}{\underline {\sigma }}(t+\Delta t)}

dóndeJσ{\displaystyle \mathbf {J} _ {\sigma }}es el jacobiano de las ecuaciones σ k :

J=(∂σ1∂λ1∂σ1∂λ2⋯∂σ1∂λnorte∂σ2∂λ1∂σ2∂λ2⋯∂σ2∂λnorte⋮⋮⋱⋮∂σnorte∂λ1∂σnorte∂λ2⋯∂σnorte∂λnorte).{\displaystyle \mathbf {J} =\left({\begin{array}{cccc}{\frac {\partial \sigma _{1}}{\partial \lambda _{1}}}&{\frac {\partial \sigma _{1}}{\partial \lambda _{2}}}&\cdots &{\frac {\partial \sigma _{1}}{\partial \lambda _{n}}}\\[5pt]{\frac {\partial \sigma _{2}}{\partial \lambda _{1}}}&{\frac {\partial \sigma _{2}}{\partial \lambda _{2}}}&\cdots &{\frac {\partial \sigma _{2}}{\partial \lambda _{n}}}\\[5pt]\vdots &\vdots &\ddots &\vdots \\[5pt]{\frac {\partial \sigma _{n}}{\partial \lambda _{1}}}&{\frac {\partial \sigma _{n}}{\partial \lambda _{2}}}&\cdots &{\frac {\partial \sigma _{n}}{\partial \lambda _{n}}}\end{array}}\right).}

Dado que no todas las partículas contribuyen a todas las restricciones,Jσ{\displaystyle \mathbf {J} _ {\sigma }}es una matriz de bloques y se puede resolver individualmente para cada bloque de la matriz. En otras palabras,Jσ{\displaystyle \mathbf {J} _ {\sigma }}se puede resolver individualmente para cada molécula.

En lugar de actualizar constantemente el vectorλ_{\displaystyle {\underline {\lambda }}}, la iteración puede comenzar conλ_(0)=0{\displaystyle {\underline {\lambda }}^{(0)}=\mathbf {0} }, lo que resulta en expresiones más simples paraσk(t){\displaystyle \sigma _{k}(t)}y∂σk(t)∂λj{\displaystyle {\frac {\partial \sigma _{k}(t)}{\partial \lambda _{j}}}}. En este caso

Jij=∂σj∂λi|λ=0=2[incógnita^jα−incógnita^jβ][∂σi∂incógnitajαmetrojα−1−∂σi∂incógnitajβmetrojβ−1].{\displaystyle J_{ij}=\left.{\frac {\partial \sigma _{j}}{\partial \lambda _{i}}}\right|_{\mathbf {\lambda } =0}=2\left[{\hat {x}}_{j\alpha }-{\hat {x}}_{j\beta }\right]\left[{\frac {\partial \sigma _{i}}{\partial x_{j\alpha }}}m_{j\alpha }^{-1}-{\frac {\partial \sigma _{i}}{\partial x_{j\beta }}}m_{j\beta }^{-1}\right].}

entoncesλ{\displaystyle \lambda }se actualiza a

λj=−J−1[‖incógnita^jα(t+Δt)−incógnita^jβ(t+Δt)‖2−dj2].{\displaystyle \mathbf {\lambda } _{j}=-\mathbf {J} ^{-1}\left[\left\|{\hat {\mathbf {x} }}_{j\alpha }(t+\Delta t)-{\hat {\mathbf {x} }}_{j\beta }(t+\Delta t)\right\|^{2}-d_{j}^{2}\right].}

Después de cada iteración, las posiciones de las partículas no restringidas se actualizan utilizando

incógnita^i(t+Δt)←incógnita^i(t+Δt)+∑k=1norteλk∂σk∂incógnitai(Δt)2metroi−1.{\displaystyle {\hat {\mathbf {x} }}_{i}(t+\Delta t)\leftarrow {\hat {\mathbf {x} }}_{i}(t+\Delta t)+\sum _{k=1}^{n}\lambda _{k}{\frac {\partial \sigma _{k}}{\partial \mathbf {x} _{i}}}\left(\Delta t\right)^{2}m_{i}^{-1}.}

El vector se restablece entonces a

λ_=0.{\displaystyle {\underline {\lambda }}=\mathbf {0} .}

El procedimiento anterior se repite hasta que se encuentre la solución de las ecuaciones de restricción.σk(t+Δt){\displaystyle \sigma _{k}(t+\Delta t)}, converge a una tolerancia prescrita de un error numérico .

Si bien existen diversos algoritmos para calcular los multiplicadores de Lagrange, la diferencia radica únicamente en los métodos para resolver el sistema de ecuaciones. Para ello, se suelen utilizar los métodos cuasi-Newton .

El algoritmo SETTLE

El algoritmo SETTLE [ 9 ] resuelve analíticamente el sistema de ecuaciones no lineales paranorte=3{\displaystyle n=3}restricciones en tiempo constante. Aunque no se adapta a un mayor número de restricciones, se utiliza con mucha frecuencia para restringir moléculas de agua rígidas, que están presentes en casi todas las simulaciones biológicas y que normalmente se modelan utilizando tres restricciones (por ejemplo, los modelos de agua SPC/E y TIP3P ).

El algoritmo SHAKE

El algoritmo SHAKE se desarrolló inicialmente para satisfacer una restricción de geometría de enlace durante simulaciones de dinámica molecular. [ 10 ] Posteriormente, el método se generalizó para manejar cualquier restricción holonómica, como las necesarias para mantener ángulos de enlace constantes o rigidez molecular. [ 11 ]

En el algoritmo SHAKE, el sistema de ecuaciones de restricción no lineales se resuelve utilizando el método de Gauss-Seidel , que aproxima la solución del sistema de ecuaciones lineales utilizando el método de Newton-Raphson ;

λ_=−Jσ−1σ_.{\displaystyle {\underline {\lambda }}=-\mathbf {J} _{\sigma }^{-1}{\underline {\sigma }}.}

Esto equivale a suponer queJσ{\displaystyle \mathbf {J} _{\sigma }}es diagonalmente dominante y resolver elk{\displaystyle k}ecuación solo para lak{\displaystyle k}desconocido. En la práctica, calculamos

λk←σk(t)∂σk(t)/∂λk,incógnitakα←incógnitakα+λk∂σk(t)∂incógnitakα,incógnitakβ←incógnitakβ+λk∂σk(t)∂incógnitakβ,{\displaystyle {\begin{aligned}\lambda _{k}&\leftarrow {\frac {\sigma _{k}(t)}{\partial \sigma _{k}(t)/\partial \lambda _{k}}},\\[5pt]\mathbf {x} _{k\alpha }&\leftarrow \mathbf {x} _{k\alpha }+\lambda _{k}{\frac {\partial \sigma _{k}(t)}{\partial \mathbf {x} _{k\alpha }}},\\[5pt]\mathbf {x} _{k\beta }&\leftarrow \mathbf {x} _{k\beta }+\lambda _{k}{\frac {\partial \sigma _{k}(t)}{\partial \mathbf {x} _{k\beta }}},\end{aligned}}}

a pesar dek=1…norte{\displaystyle k=1\ldots n}iterativamente hasta que las ecuaciones de restricciónσk(t+Δt){\displaystyle \sigma _{k}(t+\Delta t)}se resuelven con una tolerancia determinada.

El coste de cálculo de cada iteración esO(norte){\displaystyle {\mathcal {O}}(n)}y las iteraciones mismas convergen linealmente.

Posteriormente se desarrolló una forma no iterativa de SHAKE. [ 12 ]

Existen varias variantes del algoritmo SHAKE. Si bien difieren en la forma en que calculan o aplican las restricciones, estas se modelan mediante multiplicadores de Lagrange, los cuales se calculan utilizando el método de Gauss-Seidel.

El algoritmo SHAKE original es capaz de restringir moléculas tanto rígidas como flexibles (por ejemplo, agua, benceno y bifenilo ) e introduce un error o deriva de energía insignificante en una simulación de dinámica molecular. [ 13 ] Un problema con SHAKE es que el número de iteraciones necesarias para alcanzar un cierto nivel de convergencia aumenta a medida que la geometría molecular se vuelve más compleja. Para alcanzar la precisión de la computadora de 64 bits (una tolerancia relativa de≈10−16{\displaystyle \approx 10^{-16}}) en una simulación típica de dinámica molecular a una temperatura de 310 K, un modelo de agua de 3 sitios con 3 restricciones para mantener la geometría molecular requiere un promedio de 9 iteraciones (que es 3 por sitio por paso de tiempo). Un modelo de butano de 4 sitios con 5 restricciones necesita 17 iteraciones (22 por sitio), un modelo de benceno de 6 sitios con 12 restricciones necesita 36 iteraciones (72 por sitio), mientras que un modelo de bifenilo de 12 sitios con 29 restricciones requiere 92 iteraciones (229 por sitio por paso de tiempo). [ 13 ] Por lo tanto, los requisitos de CPU del algoritmo SHAKE pueden volverse significativos, particularmente si un modelo molecular tiene un alto grado de rigidez.

Una extensión posterior del método, QSHAKE ( Quaternion SHAKE), se desarrolló como una alternativa más rápida para moléculas compuestas de unidades rígidas, pero no es tan general. [ 14 ] Funciona satisfactoriamente para bucles rígidos como los sistemas de anillos aromáticos , pero QSHAKE falla para bucles flexibles, como cuando una proteína tiene un enlace disulfuro. [ 15 ]

Otras extensiones incluyen RATTLE, [ 16 ] WIGGLE, [ 17 ] y MSHAKE. [ 18 ]

Mientras que RATTLE funciona de la misma manera que SHAKE, [ 19 ] pero utilizando el esquema de integración temporal de Velocity Verlet , WIGGLE extiende SHAKE y RATTLE utilizando una estimación inicial para los multiplicadores de Lagrange.λk{\displaystyle \lambda _{k}}basado en las velocidades de las partículas. Cabe mencionar que MSHAKE calcula correcciones en las fuerzas de restricción , logrando una mejor convergencia.

Una modificación final del algoritmo SHAKE es el algoritmo P-SHAKE [ 20 ] que se aplica a moléculas muy rígidas o semirrígidas . P-SHAKE calcula y actualiza un precondicionador que se aplica a los gradientes de restricción antes de la iteración SHAKE, lo que provoca que el jacobianoJσ{\displaystyle \mathbf {J} _{\sigma }}para volverse diagonal o fuertemente diagonalmente dominante. Las restricciones desacopladas convergen mucho más rápido (cuadráticamente en lugar de linealmente) a costa deO(norte2){\displaystyle {\mathcal {O}}(n^{2})}.

El algoritmo M-SHAKE

El algoritmo M-SHAKE [ 21 ] resuelve el sistema de ecuaciones no lineales utilizando directamente el método de Newton . En cada iteración, el sistema de ecuaciones lineales

λ_=−Jσ−1σ_{\displaystyle {\underline {\lambda }}=-\mathbf {J} _{\sigma }^{-1}{\underline {\sigma }}}

se resuelve exactamente usando una descomposición LU . Cada iteración cuestaO(norte3){\displaystyle {\mathcal {O}}(n^{3})}operaciones, sin embargo, la solución converge cuadráticamente , requiriendo menos iteraciones que SHAKE.

Esta solución fue propuesta por primera vez en 1986 por Ciccotti y Ryckaert [ 11 ] bajo el título "el método matricial", pero difería en la solución del sistema lineal de ecuaciones. Ciccotti y Ryckaert sugieren invertir la matriz.Jσ{\displaystyle \mathbf {J} _{\sigma }}directamente, pero haciéndolo solo una vez, en la primera iteración. La primera iteración luego cuestaO(norte3){\displaystyle {\mathcal {O}}(n^{3})}operaciones, mientras que las siguientes iteraciones cuestan soloO(norte2){\displaystyle {\mathcal {O}}(n^{2})}operaciones (para la multiplicación matriz-vector). Sin embargo, esta mejora tiene un costo, ya que el jacobiano ya no se actualiza y la convergencia es solo lineal , aunque a un ritmo mucho más rápido que en el algoritmo SHAKE.

Barth et al. estudiaron varias variantes de este enfoque basadas en técnicas de matrices dispersas . [ 22 ]

Algoritmo SHAPE

El algoritmo SHAPE [ 23 ] es un análogo multicéntrico de SHAKE para restringir cuerpos rígidos de tres o más centros. Al igual que SHAKE, se toma un paso sin restricciones y luego se corrige calculando y aplicando directamente la matriz de rotación del cuerpo rígido que satisface:

Lrígido(t+Δt2)=Lno rígido(t+Δt2){\displaystyle L^{\text{rigid}}\left(t+{\frac {\Delta t}{2}}\right)=L^{\text{nonrigid}}\left(t+{\frac {\Delta t}{2}}\right)}

Este enfoque implica una diagonalización de una matriz de 3×3 seguida de tres o cuatro iteraciones rápidas de Newton para determinar la matriz de rotación. SHAPE proporciona la misma trayectoria que SHAKE iterativo totalmente convergente, pero resulta más eficiente y preciso que SHAKE cuando se aplica a sistemas con tres o más centros. Extiende la capacidad de las restricciones tipo SHAKE a sistemas lineales con tres o más átomos, sistemas planares con cuatro o más átomos y a estructuras rígidas significativamente mayores donde SHAKE es intratable. También permite vincular cuerpos rígidos con uno o dos centros comunes (por ejemplo, planos peptídicos) resolviendo iterativamente las restricciones de cuerpo rígido de la misma manera básica en que se usa SHAKE para átomos con más de una restricción SHAKE.

Algoritmo LINCS

Un método de restricción alternativo, LINCS (Linear Constraint Solver), fue desarrollado en 1997 por Hess, Bekker, Berendsen y Fraaije, [ 24 ] y se basó en el método de 1986 de Edberg, Evans y Morriss (EEM), [ 25 ] y una modificación del mismo por Baranyai y Evans (BE). [ 26 ]

LINCS aplica multiplicadores de Lagrange a las fuerzas de restricción y resuelve los multiplicadores utilizando un desarrollo en serie para aproximar el inverso del jacobiano.Jσ{\displaystyle \mathbf {J} _{\sigma }}:

(I−Jσ)−1=I+Jσ+Jσ2+Jσ3+⋯{\displaystyle (\mathbf {I} -\mathbf {J} _{\sigma })^{-1}=\mathbf {I} +\mathbf {J} _{\sigma }+\mathbf {J} _{\sigma }^{2}+\mathbf {J} _{\sigma }^{3}+\cdots }

en cada paso de la iteración de Newton. Esta aproximación solo funciona para matrices con valores propios menores que 1, lo que hace que el algoritmo LINCS sea adecuado únicamente para moléculas con baja conectividad.

Se ha informado que LINCS es 3 a 4 veces más rápido que SHAKE. [ 24 ]

métodos híbridos

También se han introducido métodos híbridos en los que las restricciones se dividen en dos grupos; las restricciones del primer grupo se resuelven utilizando coordenadas internas, mientras que las del segundo grupo se resuelven utilizando fuerzas de restricción, por ejemplo, mediante un multiplicador de Lagrange o un método de proyección. [ 27 ] [ 28 ] [ 29 ] Este enfoque fue iniciado por Lagrange, [ 1 ] y da como resultado ecuaciones de Lagrange de tipo mixto . [ 30 ]

Véase también

Referencias y notas al pie

  1. ^ Lagrange , GL (1788). Mécanique analytique .
  2. ↑ Noguti T, Toshiyuki; Gō N (1983). "Un método de cálculo rápido de una matriz de segunda derivada de energía conformacional para moléculas grandes". Journal of the Physical Society of Japan . 52 (10): 3685– 3690. Bibcode : 1983JPSJ...52.3685N . doi : 10.1143/JPSJ.52.3685 .
  3. ↑ Abe, H; Braun W; Noguti T; Gō N (1984). "Cálculo rápido de las derivadas 1.ª y 2.ª de la energía conformacional con respecto a los ángulos diedros para proteínas: ecuaciones recurrentes generales". Computers and Chemistry . 8 (4): 239– 247. doi : 10.1016/0097-8485(84)85015-9 .
  4. ↑ Bae, DS; Haug EJ (1988). "Una formulación recursiva para la dinámica de sistemas mecánicos restringidos: Parte I. Sistemas de lazo abierto". Mecánica de estructuras y máquinas . 15 (3): 359– 382. doi : 10.1080/08905458708905124 .
  5. ↑ Jain, A; Vaidehi N; Rodriguez G (1993). "Un algoritmo recursivo rápido para la simulación de dinámica molecular". Journal of Computational Physics . 106 (2): 258– 268. Bibcode : 1993JCoPh.106..258J . doi : 10.1006/jcph.1993.1106 .
  6. ↑ Rice, LM; Brünger AT (1994). "Dinámica del ángulo de torsión: el muestreo conformacional variable reducido mejora el refinamiento de la estructura cristalográfica". Proteins: Structure, Function, and Genetics . 19 (4): 277– 290. doi : 10.1002/prot.340190403 . PMID 7984624 . S2CID 25080482 .  
  7. ↑ Mathiowetz, AM; Jain A; Karasawa N; Goddard III, WA (1994). "Simulaciones de proteínas utilizando técnicas adecuadas para sistemas muy grandes: el método de multipolos celulares para interacciones no enlazantes y el método del operador de masa inversa de Newton-Euler para la dinámica de coordenadas internas". Proteins: Structure, Function, and Genetics . 20 (3): 227– 247. doi : 10.1002/prot.340200304 . PMID 7892172 . S2CID 25753031 .  
  8. ↑ Mazur, AK (1997). "Ecuaciones cuasi-hamiltonianas de movimiento para la dinámica molecular de coordenadas internas de polímeros". Journal of Computational Chemistry . 18 (11): 1354– 1364. arXiv : physics/9703019 . doi : 10.1002/(SICI)1096-987X(199708)18:11 < 1354::AID-JCC3 > 3.0.CO ; 2-K .
  9. 1 2 Miyamoto, S; Kollman PA (1992). "SETTLE: Una versión analítica del algoritmo SHAKE and RATTLE para modelos de agua rígida". Journal of Computational Chemistry . 13 (8): 952– 962. Bibcode : 1992JCoCh..13..952M . doi : 10.1002/jcc.540130805 . S2CID 122506495 . 
  10. ↑ Ryckaert, JP; Ciccotti G; Berendsen HJC (1977). "Integración numérica de las ecuaciones cartesianas de movimiento de un sistema con restricciones: dinámica molecular de n -alcanos". Journal of Computational Physics . 23 (3): 327– 341. Bibcode : 1977JCoPh..23..327R . CiteSeerX 10.1.1.399.6868 . doi : 10.1016/0021-9991(77)90098-5 . 
  11. 1 2 Ciccotti, G.; JP Ryckaert (1986). "Simulación de dinámica molecular de moléculas rígidas". Computer Physics Reports . 4 (6): 345– 392. Bibcode : 1986CoPhR...4..346C . doi : 10.1016/0167-7977(86)90022-5 .
  12. ↑ Yoneya, M; Berendsen HJC; Hirasawa K (1994). "Un método matricial no iterativo para simulaciones de dinámica molecular con restricciones". Molecular Simulations . 13 (6): 395– 405. doi : 10.1080/08927029408022001 .
  13. 1 2 Hammonds, KD; Heyes DM (2020). "Hamiltoniano de sombra en simulaciones de dinámica molecular NVE clásicas: un camino hacia la estabilidad a largo plazo". Journal of Chemical Physics . 152 (2): 024114_1–024114_15. doi : 10.1063/1.5139708 . PMID 31941339 . S2CID 210333551 .  
  14. ↑ Forester, TR; Smith W (1998). "SHAKE, Rattle, and Roll: Efficient Constraint Algorithms for Linked Rigid Bodies". Journal of Computational Chemistry . 19 : 102– 111. doi : 10.1002/(SICI)1096-987X(19980115)19:1 < 102::AID-JCC9 > 3.0.CO ; 2-T .
  15. ↑ McBride, C; Wilson MR; Howard JAK (1998). "Simulaciones de dinámica molecular de fases de cristal líquido utilizando potenciales atomísticos". Molecular Physics . 93 (6): 955– 964. Bibcode : 1998MolPh..93..955C . doi : 10.1080/002689798168655 .
  16. ↑ Andersen, Hans C. (1983). "RATTLE: Una versión de "velocidad" del algoritmo SHAKE para cálculos de dinámica molecular". Journal of Computational Physics . 52 (1): 24– 34. Bibcode : 1983JCoPh..52...24A . CiteSeerX 10.1.1.459.5668 . doi : 10.1016/0021-9991(83)90014-1 . 
  17. ↑ Lee, Sang-Ho; Kim Palmo; Samuel Krimm (2005). "WIGGLE: Un nuevo algoritmo de dinámica molecular con restricciones en coordenadas cartesianas". Journal of Computational Physics . 210 (1): 171– 182. Bibcode : 2005JCoPh.210..171L . doi : 10.1016/j.jcp.2005.04.006 .
  18. ↑ Lambrakos, SG; JP Boris; ES Oran; I. Chandrasekhar; M. Nagumo (1989). "Un algoritmo SHAKE modificado para mantener enlaces rígidos en simulaciones de dinámica molecular de moléculas grandes". Journal of Computational Physics . 85 (2): 473– 486. Bibcode : 1989JCoPh..85..473L . doi : 10.1016/0021-9991(89)90160-5 .
  19. ↑ Leimkuhler, Benedict; Robert Skeel (1994). "Integradores numéricos simplécticos en sistemas hamiltonianos restringidos". Journal of Computational Physics . 112 (1): 117– 125. Bibcode : 1994JCoPh.112..117L . doi : 10.1006/jcph.1994.1085 .
  20. ↑ Gonnet, Pedro (2007). "P-SHAKE: Un SHAKE cuadráticamente convergente enO(norte2){\displaystyle {\mathcal {O}}(n^{2})}". Journal of Computational Physics . 220 (2): 740– 750. Bibcode : 2007JCoPh.220..740G . doi : 10.1016/j.jcp.2006.05.032 .
  21. ↑ Kräutler, Vincent; WF van Gunsteren; PH Hünenberger (2001). "Un algoritmo SHAKE rápido para resolver ecuaciones de restricción de distancia para moléculas pequeñas en simulaciones de dinámica molecular". Journal of Computational Chemistry . 22 (5): 501– 508. doi : 10.1002/1096-987X(20010415)22:5 < 501::AID-JCC1021 > 3.0.CO ; 2-V . S2CID 6187100 . 
  22. ↑ Barth, Eric; K. Kuczera; B. Leimkuhler; R. Skeel (1995). "Algoritmos para dinámica molecular restringida". Journal of Computational Chemistry . 16 (10): 1192– 1209. Bibcode : 1995JCoCh..16.1192B . doi : 10.1002/jcc.540161003 . S2CID 38109923 . 
  23. ↑ Tao, Peng; Xiongwu Wu; Bernard R. Brooks (2012). "Mantener estructuras rígidas en simulaciones de dinámica molecular cartesiana basadas en Verlet" . The Journal of Chemical Physics . 137 (13): 134110. Bibcode : 2012JChPh.137m4110T . doi : 10.1063/1.4756796 . PMC 3477181. PMID 23039588 .  
  24. 1 2 Hess, B; Bekker H; Berendsen HJC; Fraaije JGEM (1997). "LINCS: Un solucionador de restricciones lineales para simulaciones moleculares". Journal of Computational Chemistry . 18 (12): 1463– 1472. CiteSeerX 10.1.1.48.2727 . doi : 10.1002/(SICI)1096-987X(199709)18:12 < 1463::AID-JCC4 > 3.0.CO ; 2-H . 
  25. ↑ Edberg, R; Evans DJ; Morriss GP (1986). "Simulaciones de dinámica molecular restringida de alcanos líquidos con un nuevo algoritmo". Journal of Chemical Physics . 84 (12): 6933– 6939. Bibcode : 1986JChPh..84.6933E . doi : 10.1063/1.450613 .
  26. ↑ Baranyai, A; Evans DJ (1990). "Nuevo algoritmo para la simulación de dinámica molecular restringida de benceno y naftaleno líquidos". Molecular Physics . 70 (1): 53– 63. Bibcode : 1990MolPh..70...53B . doi : 10.1080/00268979000100841 .
  27. ↑ Mazur, AK (1999). "Integración simpléctica de la dinámica de cuerpos rígidos de cadena cerrada con ecuaciones de movimiento de coordenadas internas". Journal of Chemical Physics . 111 (4): 1407– 1414. Bibcode : 1999JChPh.111.1407M . doi : 10.1063/1.479399 .
  28. ↑ Bae, DS; Haug EJ (1988). "Una formulación recursiva para la dinámica de sistemas mecánicos restringidos: Parte II. Sistemas de bucle cerrado". Mecánica de estructuras y máquinas . 15 (4): 481– 506. doi : 10.1080/08905458708905130 .
  29. ↑ Rodríguez, G; Jain A; Kreutz-Delgado K (1991). "Un álgebra de operadores espaciales para el modelado y control de manipuladores". The International Journal of Robotics Research . 10 (4): 371– 381. doi : 10.1177/027836499101000406 . hdl : 2060/19900020578 . S2CID 12166182 . 
  30. ↑ Sommerfeld, Arnold (1952). Lecciones de física teórica, vol. I: mecánica . Nueva York: Academic Press. ISBN 978-0-12-654670-5.{{cite book}}: Incompatibilidad de ISBN/Fecha ( ayuda )
Lecturas relacionadas
6 articulos
Química computacionalQuímica computacionalLa química computacional es una rama de la química que utiliza simulaciones por computadora para ayudar a resolver problemas químicos. [2] Utiliza ...{\displaystyle {\mathcal {L}}(x,\lambda )\equiv f(x)+\langle \lambda ,g(x)\rangle }Multiplicador de LagrangeEn optimización matemática , el método de los multiplicadores de Lagrange es una estrategia para encontrar los máximos y mínimos locales de una fun...Dinámica molecularDinámica molecularLa dinámica molecular ( MD ) es un método de simulación por computadora para analizar los movimientos físicos de átomos y moléculas . Se permite qu...Vector euclidianoUn vector que apunta de A a B En matemáticas , física e ingeniería , un vector euclidiano o simplemente un vector (a veces llamado vector geométric...{\displaystyle {\textbf {F}}={\frac {d\mathbf {p} }{dt}}}Coordenadas generalizadasEn mecánica analítica , las coordenadas generalizadas son un conjunto de parámetros utilizados para representar el estado de un sistema en un espac...Matriz cuadradaMatriz cuadradaEn matemáticas , una matriz cuadrada es una matriz con el mismo número de filas y columnas. Una matriz n por n se conoce como matriz cuadrada de or...