Restricción (química computacional)

Método para satisfacer el movimiento newtoniano de un cuerpo rígido que consta de puntos de masa

En química computacional , un algoritmo de restricción es un método para satisfacer el movimiento newtoniano de un cuerpo rígido que consta de puntos de masa. Un algoritmo de restricción se utiliza para garantizar que se mantenga la distancia entre los puntos de masa. Los pasos generales involucrados son: (i) elegir nuevas coordenadas no restringidas (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 a menudo a simulaciones de dinámica molecular . Aunque a veces dichas simulaciones se realizan utilizando coordenadas internas que satisfacen automáticamente las restricciones de longitud de enlace, ángulo de enlace y ángulo de torsión, también se pueden realizar simulaciones utilizando fuerzas de restricción explícitas o implícitas para estas tres restricciones. Sin embargo, las fuerzas de restricción explícitas dan lugar a ineficiencia; se requiere más potencia computacional para obtener una trayectoria de una longitud determinada. Por lo tanto, generalmente se prefieren las coordenadas internas y los solucionadores de restricciones de fuerza implícita.

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

Antecedentes matemáticos

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.

${\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 cartesiana de 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 de 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

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

donde el índice j va de 1 a M . Para abreviar, estas funciones g _i se agrupan en un vector g de dimensión M a continuación. La tarea consiste en resolver el conjunto combinado de ecuaciones diferenciales algebraicas (DAE), 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 expuso la mayoría de los métodos para resolverlo. ^[1] El enfoque más simple es definir nuevas coordenadas generalizadas que no tengan restricciones; este enfoque elimina las ecuaciones algebraicas y reduce el problema una vez más a la solució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 la 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. El inconveniente 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 trabajen para mantener la restricción; por ejemplo, se podrían introducir fuerzas de resorte fuertes 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 las fuerzas fuertes pueden requerir pasos de tiempo muy cortos, lo que hace que las simulaciones sean ineficientes desde el punto de vista computacional.

Un tercer enfoque es 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.

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

Métodos de coordenadas internas

El enfoque más simple para satisfacer las restricciones en la minimización de la energía y la dinámica molecular es representar el sistema mecánico en las llamadas coordenadas internas que corresponden a los grados de libertad independientes sin restricciones del sistema. Por ejemplo, los ángulos diedros de una proteína son un conjunto independiente de coordenadas que especifican las posiciones de todos los átomos sin requerir ninguna restricción. La dificultad de estos enfoques de coordenadas internas es doble: 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, por ejemplo, en el fruncimiento de anillos o cuando una proteína tiene un enlace disulfuro.

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

Los solucionadores de restricciones de coordenadas internas, recursivos y eficientes 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 utilizando multiplicadores de Lagrange : a) las posiciones sin restricciones se obtienen 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 modo que se satisfagan las restricciones.

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

${\displaystyle \sigma _{k}(t):=\|\mathbf {x} _{k\alpha }(t)-\mathbf {x} _{k\beta }(t)\|^{2}-d_{k}^{2}=0,\quad k=1\ldots n}$

donde y son las posiciones de las dos partículas involucradas en la restricción k en el tiempo t y es la distancia entre partículas prescrita. ${\displaystyle \scriptstyle \mathbf {x} _{k\alpha }(t)}$ ${\displaystyle \scriptstyle \mathbf {x} _{k\beta }(t)}$ ${\displaystyle d_{k}}$

Las fuerzas debidas a estas restricciones se suman en las ecuaciones de movimiento, lo que da como resultado, para cada una de las N partículas en el sistema

${\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.}$

La suma de las fuerzas de restricción no modifica 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. Nótese que el signo es arbitrario y algunas referencias ^[9] tienen un signo opuesto. ${\displaystyle \lambda _{k}}$

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

${\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}$

donde es la posición sin restricciones (o sin corregir) de la i -ésima partícula después de integrar las ecuaciones de movimiento sin restricciones. ${\displaystyle {\hat {\mathbf {x} }}_{i}(t+\Delta t)}$

Para satisfacer las restricciones en el siguiente paso de tiempo, los multiplicadores de Lagrange deben determinarse como la siguiente ecuación, ${\displaystyle \sigma _{k}(t+\Delta t)}$

${\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 de ecuaciones no lineales. ${\displaystyle n}$

${\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 los multiplicadores de Lagrange desconocidos . ${\displaystyle n}$ ${\displaystyle \lambda _{k}}$

Este sistema de ecuaciones no lineales con incógnitas se resuelve comúnmente utilizando el método de Newton-Raphson, donde el vector solución se actualiza utilizando ${\displaystyle n}$ ${\displaystyle n}$ ${\displaystyle {\underline {\lambda }}}$

${\displaystyle {\underline {\lambda }}^{(l+1)}\leftarrow {\underline {\lambda }}^{(l)}-\mathbf {J} _{\sigma }^{-1}{\underline {\sigma }}(t+\Delta t)}$

¿Dónde está el jacobiano de las ecuaciones σ _k : ${\displaystyle \mathbf {J} _{\sigma }}$

${\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, es una matriz de bloques y se puede resolver individualmente para obtener la unidad de bloques de la matriz. En otras palabras, se puede resolver individualmente para cada molécula. ${\displaystyle \mathbf {J} _{\sigma }}$ ${\displaystyle \mathbf {J} _{\sigma }}$

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

${\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].}$

luego se actualiza a ${\displaystyle \lambda }$

${\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 partículas sin restricciones se actualizan utilizando

${\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}.}$

A continuación, el vector se restablece a

${\displaystyle {\underline {\lambda }}=\mathbf {0} .}$

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

Aunque existen varios algoritmos para calcular los multiplicadores de Lagrange, estas diferencias dependen únicamente de los métodos para resolver el sistema de ecuaciones. Para estos métodos, se utilizan comúnmente los métodos cuasi-Newton .

El algoritmo SETTLE

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

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. ^[11] Luego, el método se generalizó para manejar cualquier restricción holonómica, como las requeridas para mantener ángulos de enlace constantes o rigidez molecular. ^[12]

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

${\displaystyle {\underline {\lambda }}=-\mathbf {J} _{\sigma }^{-1}{\underline {\sigma }}.}$

Esto equivale a suponer que es diagonalmente dominante y resolver la ecuación sólo para la incógnita. En la práctica, calculamos ${\displaystyle \mathbf {J} _{\sigma }}$ ${\displaystyle k}$ ${\displaystyle k}$

${\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}}}$

para todos iterativamente hasta que las ecuaciones de restricción se resuelvan a una tolerancia dada. ${\displaystyle k=1\ldots n}$ ${\displaystyle \sigma _{k}(t+\Delta t)}$

El coste de cálculo de cada iteración es , y las iteraciones en sí convergen linealmente. ${\displaystyle {\mathcal {O}}(n)}$

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

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

El algoritmo SHAKE original es capaz de restringir tanto moléculas rígidas como flexibles (por ejemplo, agua, benceno y bifenilo ) e introduce un error insignificante o una deriva de energía en una simulación de dinámica molecular. ^[14] 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 una precisión de computadora de 64 bits (una tolerancia relativa de ) en una simulación de dinámica molecular típica a una temperatura de 310 K, un modelo de agua de 3 sitios que tiene 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). ^[14] Por lo tanto, los requisitos de CPU del algoritmo SHAKE pueden llegar a ser significativos, particularmente si un modelo molecular tiene un alto grado de rigidez. ${\displaystyle \approx 10^{-16}}$

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 de uso general. ^[15] 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. ^[16]

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

Si bien RATTLE funciona de la misma manera que SHAKE, ^[20] pero utilizando el esquema de integración temporal de Velocity Verlet , WIGGLE extiende SHAKE y RATTLE al utilizar una estimación inicial para los multiplicadores de Lagrange basados ​​en las velocidades de las partículas. Vale la pena mencionar que MSHAKE calcula correcciones en las fuerzas de restricción , logrando una mejor convergencia. ${\displaystyle \lambda _{k}}$

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

El algoritmo M-SHAKE

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

${\displaystyle {\underline {\lambda }}=-\mathbf {J} _{\sigma }^{-1}{\underline {\sigma }}}$

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

Esta solución fue propuesta por primera vez en 1986 por Ciccotti y Ryckaert ^[12] 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 directamente, pero hacerlo solo una vez, en la primera iteración. La primera iteración entonces cuesta operaciones, mientras que las iteraciones siguientes cuestan solo operaciones (para la multiplicación matriz-vector). Sin embargo, esta mejora tiene un costo, ya que el jacobiano ya no se actualiza, la convergencia es solo lineal , aunque a una velocidad mucho más rápida que para el algoritmo SHAKE. ${\displaystyle \mathbf {J} _{\sigma }}$ ${\displaystyle {\mathcal {O}}(n^{3})}$ ${\displaystyle {\mathcal {O}}(n^{2})}$

Barth et al. estudiaron varias variantes de este enfoque basadas en técnicas de matriz dispersa . ^[23]

Algoritmo SHAPE

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

${\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 única diagonalización de la matriz 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 se proporciona con SHAKE iterativo completamente convergente, pero se ha descubierto que es más eficiente y más preciso que SHAKE cuando se aplica a sistemas que involucran tres o más centros. Amplía la capacidad de las restricciones como SHAKE a sistemas lineales con tres o más átomos, sistemas planos con cuatro o más átomos y a estructuras rígidas significativamente más grandes donde SHAKE es intratable. También permite vincular cuerpos rígidos con uno o dos centros comunes (por ejemplo, planos de péptidos) al resolver restricciones de cuerpos rígidos de manera iterativa de la misma manera básica en que SHAKE se usa para átomos que involucran 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, ^[25] y se basó en el método de 1986 de Edberg, Evans y Morriss (EEM), ^[26] y una modificación del mismo de Baranyai y Evans (BE). ^[27]

LINCS aplica multiplicadores de Lagrange a las fuerzas de restricción y resuelve los multiplicadores utilizando una expansión en serie para aproximar la inversa del jacobiano : ${\displaystyle \mathbf {J} _{\sigma }}$

${\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 solo para moléculas con baja conectividad.

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

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. ^{[28] [29] [30]} Este enfoque fue iniciado por Lagrange, ^[1] y da como resultado ecuaciones de Lagrange del tipo mixto . ^[31]

Véase también

• Dinámica molecular
• Software para modelado de mecánica molecular

Referencias y notas a pie de página

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". Revista de la Sociedad de Física de Japón . 52 (10): 3685–3690. Código Bibliográfico :1983JPSJ...52.3685N. doi :10.1143/JPSJ.52.3685.
3. Abe, H; Braun W; Noguti T; Gō N (1984). "Cálculo rápido de la primera y la segunda derivadas 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 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". Proteínas: estructura, función y genética . 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 multipolar celular para interacciones sin enlaces y el método del operador de masa inversa de Newton-Euler para dinámica de coordenadas internas". Proteínas: estructura, función y genética . 20 (3): 227–247. doi :10.1002/prot.340200304. PMID  7892172. S2CID  25753031.
8. Mazur, AK (1997). "Ecuaciones cuasi-hamiltonianas de movimiento para 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. Miyamoto, S; Kollman PA (1992). "SETTLE: Una versión analítica del algoritmo SHAKE y RATTLE para modelos de agua rígida". Revista de química computacional . 13 (8): 952–962. doi :10.1002/jcc.540130805. S2CID  122506495.
10. Miyamoto, S; Kollman PA (1992). "SETTLE: Una versión analítica del algoritmo SHAKE y RATTLE para modelos de agua rígida". Revista de química computacional . 13 (8): 952–962. doi :10.1002/jcc.540130805. S2CID  122506495.
11. 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". Revista de física computacional . 23 (3): 327–341. Bibcode :1977JCoPh..23..327R. CiteSeerX 10.1.1.399.6868 . doi :10.1016/0021-9991(77)90098-5. 
12. 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.
13. Yoneya, M; Berendsen HJC; Hirasawa K (1994). "Un método matricial no iterativo para simulaciones de dinámica molecular con restricciones". Simulaciones moleculares . 13 (6): 395–405. doi :10.1080/08927029408022001.
14. 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". Revista de física química . 152 (2): 024114_1–024114_15. doi :10.1063/1.5139708. PMID  31941339. S2CID  210333551.
15. Forester, TR; Smith W (1998). "SHAKE, Rattle, and Roll: Algoritmos de restricción eficientes para cuerpos rígidos vinculados". Journal of Computational Chemistry . 19 : 102–111. doi :10.1002/(SICI)1096-987X(19980115)19:1<102::AID-JCC9>3.0.CO;2-T.
16. McBride, C; Wilson MR; Howard JAK (1998). "Simulaciones de dinámica molecular de fases de cristal líquido utilizando potenciales atomísticos". Física molecular . 93 (6): 955–964. Código Bibliográfico :1998MolPh..93..955C. doi :10.1080/002689798168655.
17. 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. 
18. Lee, Sang-Ho; Kim Palmo; Samuel Krimm (2005). "WIGGLE: Un nuevo algoritmo de dinámica molecular restringida en coordenadas cartesianas". Revista de Física Computacional . 210 (1): 171–182. Código Bibliográfico :2005JCoPh.210..171L. doi :10.1016/j.jcp.2005.04.006.
19. 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.
20. 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.
21. Gonnet, Pedro (2007). "P-SHAKE: Un SHAKE cuadráticamente convergente en ". Journal of Computational Physics . 220 (2): 740–750. Bibcode :2007JCoPh.220..740G. doi :10.1016/j.jcp.2006.05.032. ${\displaystyle {\mathcal {O}}(n^{2})}$
22. 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". Revista de química computacional . 22 (5): 501–508. doi :10.1002/1096-987X(20010415)22:5<501::AID-JCC1021>3.0.CO;2-V. S2CID  6187100.
23. Barth, Eric; K. Kuczera; B. Leimkuhler; R. Skeel (1995). "Algoritmos para dinámica molecular restringida". Revista de química computacional . 16 (10): 1192–1209. doi :10.1002/jcc.540161003. S2CID  38109923.
24. 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. 
25. Hess, B; Bekker H; Berendsen HJC; Fraaije JGEM (1997). "LINCS: Un solucionador de restricciones lineales para simulaciones moleculares". Revista de química computacional . 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. 
26. 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.
27. Baranyai, A; Evans DJ (1990). "Nuevo algoritmo para simulación de dinámica molecular restringida de benceno líquido y naftaleno". Física molecular . 70 (1): 53–63. Código Bibliográfico :1990MolPh..70...53B. doi :10.1080/00268979000100841.
28. 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.
29. 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.
30. Rodríguez, G; Jain A; Kreutz-Delgado K (1991). "Un álgebra de operadores espaciales para el modelado y control de manipuladores". Revista internacional de investigación en robótica . 10 (4): 371–381. doi :10.1177/027836499101000406. hdl : 2060/19900020578 . S2CID  12166182.
31. Sommerfeld, Arnold (1952). Lecciones de física teórica, vol. I: Mecánica . Nueva York: Academic Press. ISBN 978-0-12-654670-5.