Dinámica de Langevin

En física, la dinámica de Langevin es un enfoque para el modelado matemático de la dinámica de los sistemas moleculares utilizando la ecuación de Langevin. Fue desarrollada originalmente por el físico francés Paul Langevin. El enfoque se caracteriza por el uso de modelos simplificados, al tiempo que tiene en cuenta los grados de libertad omitidos mediante el uso de ecuaciones diferenciales estocásticas. Las simulaciones de la dinámica de Langevin son un tipo de simulación de Monte Carlo.^[1]​

Los sistemas moleculares del mundo real se producen en el aire o en disolventes, en lugar de forma aislada, en el vacío. El choque de las moléculas del disolvente o del aire provoca fricción, y las colisiones ocasionales a alta velocidad perturban el sistema. La dinámica de Langevin intenta ampliar la dinámica molecular para tener en cuenta estos efectos. Además, la dinámica de Langevin permite controlar la temperatura como con un termostato, aproximándose así al conjunto canónico.

La dinámica de Langevin imita el aspecto viscoso de un disolvente. No modela completamente un disolvente implícito; concretamente, el modelo no tiene en cuenta el apantallamiento electrostático ni el efecto hidrofóbico. En el caso de disolventes más densos, las interacciones hidrodinámicas no se captan mediante la dinámica de Langevin.

Para un sistema de ${\displaystyle N}$ partículas con masas ${\displaystyle M}$, con coordenadas ${\displaystyle X=X(t)}$ que constituyen una variable aleatoria dependiente del tiempo, la ecuación de Langevin resultante es^[2]​^[3]​ $${\displaystyle M\,{\ddot {\mathbf {X} }}=-\mathbf {\nabla } U(\mathbf {X} )-\gamma \,M\,{\dot {\mathbf {X} }}+{\sqrt {2\,M\,\gamma \,k_{\rm {B}}T}}\,\mathbf {R} (t)\,,}$$ donde ${\displaystyle U(\mathbf {X} )}$ es el potencial de interacción entre partículas; ${\displaystyle \nabla }$ es el operador gradiente, de modo que ${\displaystyle -\mathbf {\nabla } U(\mathbf {X} )}$ es la fuerza calculada a partir de los potenciales de interacción entre partículas; el punto es una derivada temporal, de modo que ${\displaystyle {\dot {\mathbf {X} }}}$ es la velocidad y ${\displaystyle {\ddot {\mathbf {X} }}}$ es la aceleración; ${\displaystyle \gamma }$ es la constante de amortiguación (unidades de tiempo recíproco), también conocida como frecuencia de colisión; ${\displaystyle T}$ es la temperatura, ${\displaystyle k_{\rm {B}}}$ es la constante de Boltzmann; y ${\displaystyle \mathbf {R} (t)}$ es un estacionario correlacionado con delta proceso gaussiano con media cero, denominado ruido blanco gaussiano, que satisface $${\displaystyle \left\langle \mathbf {R} (t)\right\rangle =0}$$ $${\displaystyle \left\langle \mathbf {R} (t)\cdot \mathbf {R} (t')\right\rangle =\delta (t-t')}$$

Aquí, ${\displaystyle \delta }$ es la Delta de Dirac

Considerando la covarianza del movimiento browniano estándar o proceso de Wiener ${\displaystyle W_{t}}$, podemos encontrar que

$${\displaystyle \mathbb {E} (W_{t}W_{\tau })=\min(t,\tau )}$$

Definimos la matriz de covarianza de la derivada como $${\displaystyle \mathbb {E} ({\dot {W_{t}}}{\dot {W_{\tau }}})={\frac {\partial }{\partial t}}{\frac {\partial }{\partial \tau }}\mathbb {E} (W_{t}W_{\tau })={\frac {\partial }{\partial t}}{\frac {\partial }{\partial \tau }}\min(t,\tau )=\delta (t-\tau )}$$ Así, en el sentido de covarianza, podemos decir que $${\displaystyle {\rm {d}}W_{t}=\mathbf {R} (t){\rm {d}}t}$$

Sin pérdida de generalidad, sea la masa ${\displaystyle M=1}$, ${\displaystyle \sigma ={\sqrt {M\gamma k_{\rm {B}}T}}}$, entonces la EDE original se convertirá en $${\displaystyle {\rm {d}}{\dot {\mathbf {X} }}=-\nabla U(\mathbf {X} ){\rm {d}}t-\gamma {\rm {d}}{\mathbf {X} }+{\sqrt {2}}\sigma {\rm {d}}\mathbf {W} (t)}$$

Si el objetivo principal es controlar la temperatura, se debe tener cuidado de utilizar una constante de amortiguación pequeña ${\displaystyle \gamma }$. A medida que ${\displaystyle \gamma }$ crece, la distancia se extiende desde el régimen inercial hasta el difusivo (Browniano). El límite de la dinámica de Langevin de la no inercia se describe comúnmente como Dinámica browniana. La dinámica browniana puede considerarse como una dinámica de Langevin sobredampeada, es decir, una dinámica de Langevin en la que no se produce una aceleración media. Bajo este límite, tenemos ${\displaystyle {\rm {d}}{\dot {X}}=0}$, por lo que la EDE original se convierte en

$${\displaystyle {\rm {d}}{\mathbf {X} }=-{\frac {1}{\gamma }}\nabla U(\mathbf {X} ){\rm {d}}t+{\frac {{\sqrt {2}}\sigma }{\gamma }}{\rm {d}}\mathbf {W} (t)}$$

La ecuación de Langevin traslacional se puede resolver utilizando varios métodos numéricos ^[4]​^[5]​ con diferencias en la sofisticación de las soluciones analíticas, los pasos de tiempo permitidos, la reversibilidad temporal (métodos simplécticos), en el límite de fricción cero, “'etc.”'

El termostato Langevin^[6]​ es un tipo de algoritmo termostato en dinámica molecular, que se utiliza para simular un conjunto canónico (NVT) a una temperatura deseada. Integra la siguiente ecuación de movimiento de Langevin:

${\displaystyle M{\ddot {\mathbf {X} }}=-\nabla U(\mathbf {X} )-\gamma {\dot {\mathbf {X} }}+{\sqrt {2\gamma k_{B}T}}{\textbf {R}}(t)}$

${\displaystyle -\nabla U(\mathbf {X} )}$ es el término de fuerza determinista; ${\displaystyle \gamma }$ es el coeficiente de fricción y ${\displaystyle \gamma {\dot {X}}}$ es el término de fricción o amortiguación; el último término es el término de fuerza aleatoria (${\displaystyle k_{B}}$: constante de Boltzmann, ${\displaystyle T}$: temperatura). Esta ecuación permite que el sistema se acople con un «baño de calor» imaginario: la energía cinética del sistema se disipa a partir del término de fricción/amortiguación y se gana a partir de la fuerza aleatoria/fluctuación; la fuerza del acoplamiento se controla mediante ${\displaystyle \gamma }$. Esta ecuación se puede simular con solucionadores de SDE, como el método de Euler-Maruyama, en el que el término de fuerza aleatoria se sustituye por un número aleatorio gaussiano en cada paso de integración (varianza ${\displaystyle \sigma ^{2}=2\gamma k_{B}T/\Delta t}$, ${\displaystyle \Delta t}$: paso de tiempo), o Método del salto de rana Langevin, etc. Este método también se conoce como integrador de Langevin.^[7]​

La ecuación de Langevin sobreamortiguada da

${\displaystyle {\rm {d}}\mathbf {x} _{t}=-{\frac {D}{k_{B}T}}\nabla _{\mathbf {x} }U(\mathbf {x} _{t}){\rm {d}}t+{\sqrt {2D}}{\rm {d}}W_{t}}$

Aquí, ${\displaystyle D=k_{B}T/\gamma }$ es el coeficiente de difusión de la relación de Einstein. Como se ha demostrado con la ecuación de Fokker-Planck, en condiciones adecuadas, la distribución estacionaria de ${\displaystyle \mathbf {x} _{t}}$ es la distribución de Boltzmann ${\displaystyle p(\mathbf {x} )\propto e^{-U(\mathbf {x} )/k_{B}T}}$.

Dado que ${\displaystyle \nabla \log p(\mathbf {x} )=-\nabla U(\mathbf {x} )/k_{B}T}$, esta ecuación es equivalente a la siguiente forma:

${\displaystyle {\rm {d}}\mathbf {x} _{t}=\epsilon \nabla _{\mathbf {x} }\log p(\mathbf {x} _{t}){\rm {d}}t+{\sqrt {2\epsilon }}{\rm {d}}W_{t}}$

Y la distribución de ${\displaystyle \mathbf {x} _{t}(t\to \infty )}$ sigue ${\displaystyle p(\mathbf {x} )}$. En otras palabras, la dinámica de Langevin impulsa las partículas hacia una distribución estacionaria ${\displaystyle p(\mathbf {x} )}$ a lo largo de un flujo gradiente, debido al término ${\displaystyle \nabla \log p(\mathbf {x} )}$, al tiempo que permite algunas fluctuaciones aleatorias. Esto proporciona unos Métodos de Montecarlo basados en cadenas de Markov que se pueden utilizar para muestrear datos ${\displaystyle \mathbf {x} }$ de una distribución objetivo ${\displaystyle p(\mathbf {x} )}$, conocido como Langevin Monte Carlo.

En muchas aplicaciones, tenemos una distribución deseada ${\displaystyle p(\mathbf {x} )}$ de la que nos gustaría obtener una muestra ${\displaystyle \mathbf {x} }$, pero el muestreo directo puede resultar complicado o ineficaz. Langevin Monte Carlo ofrece otra forma de muestrear ${\displaystyle \mathbf {x} \sim p(\mathbf {x} )}$ mediante el muestreo de una cadena de Markov de acuerdo con la dinámica de Langevin, cuyo estado estacionario es ${\displaystyle p(\mathbf {x} )}$. El algoritmo de Langevin ajustado a Metropolis (MALA) es un ejemplo: dado un estado actual ${\displaystyle \mathbf {x} _{t}}$, el método MALA propone un nuevo estado ${\displaystyle {\tilde {x}}_{t+1}}$ utilizando la dinámica de Langevin anterior. A continuación, la propuesta se acepta o se rechaza en función del algoritmo de Metropolis-Hastings. La incorporación de la dinámica de Langevin en la elección de ${\displaystyle {\tilde {x}}_{t+1}}$ proporciona una mayor eficiencia computacional, ya que la dinámica impulsa las partículas hacia regiones con mayor probabilidad ${\displaystyle p(\mathbf {x} )}$ y, por lo tanto, es más probable que sean aceptadas. Más información en Algoritmo de Langevin ajustado por Metropolis.

La dinámica de Langevin es una de las bases de los modelos generativos basados en puntuación.^[8]​^[9]​ De la dinámica de Langevin (sobreamortiguada) dinámica de Langevin,

${\displaystyle {\rm {d}}\mathbf {x} _{t}=\epsilon \nabla _{\mathbf {x} }\log p(\mathbf {x} _{t}){\rm {d}}t+{\sqrt {2\epsilon }}{\rm {d}}W_{t}}$

Un modelo generativo tiene como objetivo generar muestras que sigan (una distribución de datos desconocida) ${\displaystyle p(\mathbf {x} )}$. Para lograrlo, un modelo basado en puntuaciones aprende una función de puntuación aproximada ${\displaystyle \mathbf {s} _{\theta }(\mathbf {x} )\approx \nabla _{\mathbf {x} }\log p(\mathbf {x} )}$ (un proceso denominado emparejamiento de puntuaciones). Con acceso a una función de puntuación, las muestras se generan mediante la siguiente iteración,^[8]​^[9]​

${\displaystyle \mathbf {x} _{i+1}\gets \mathbf {x} _{i}+\epsilon \nabla _{\mathbf {x} }\log p(\mathbf {x} _{i})+{\sqrt {2\epsilon }}\mathbf {z} _{i},\quad i=0,1,\cdots ,K}$

con ${\displaystyle \mathbf {z} _{i}\sim N(0,1)}$. A medida que ${\displaystyle \epsilon \to 0}$ y ${\displaystyle K\to \infty }$, los ${\displaystyle \mathbf {x} _{K}}$ generados convergen en la distribución objetivo ${\displaystyle p(\mathbf {x} )}$. Los modelos basados en puntuaciones utilizan ${\displaystyle \mathbf {s} _{\theta }(\mathbf {x} )\approx \nabla _{\mathbf {x} }\log p(\mathbf {x} )}$ como aproximación.^[8]​

Como ecuación diferencial estocástica (SDE), la ecuación dinámica de Langevin tiene su correspondiente ecuación diferencial parcial (PDE), la ecuación de Klein-Kramers, una ecuación de Fokker-Planck especial que rige la distribución de probabilidad de las partículas en el espacio de fase. La ecuación original de la “”'dinámica de Langevin'“” puede reformularse como las siguientes EDP de primer orden: $${\displaystyle {\rm {d}}\mathbf {X} =\mathbf {P} {\rm {d}}t}$$ $${\displaystyle {\rm {d}}\mathbf {P} =-\gamma \mathbf {P} {\rm {d}}t-\nabla U(\mathbf {X} ){\rm {d}}t+{\sqrt {2}}\sigma {\rm {d}}\mathbf {W} (t)}$$ Ahora consideremos los siguientes casos y su ley de ${\displaystyle (\mathbf {X} ,\mathbf {P} )}$:

1.${\displaystyle \mathbf {{\rm {d}}{X}} =\mathbf {P} {\rm {d}}t,\mathbf {{\rm {d}}{P}} =-\gamma \mathbf {P} {\rm {d}}t-\nabla U(\mathbf {X} ){\rm {d}}t+{\sqrt {2}}\sigma {\rm {d}}\mathbf {W} (t)}$ con ${\displaystyle (\mathbf {X} _{0},\mathbf {P} _{0})\sim \rho _{0}}$

2. ${\displaystyle {\frac {\partial \rho }{\partial t}}=-\mathbf {P} \nabla _{\mathbf {X} }\rho +\nabla _{\mathbf {P} }(\gamma \mathbf {P} \rho +\nabla _{\mathbf {X} }U(\mathbf {X} )\rho )+\nabla _{\mathbf {P} }^{2}(\sigma _{T}^{2}\rho )}$ con ${\displaystyle \rho (t=0,\mathbf {X} ,\mathbf {P} )=\rho _{0}}$

Consideremos una función general de momento y posición

$${\displaystyle \Psi _{t}=\Psi (\mathbf {X} ,\mathbf {P} )}$$

El valor esperado de la función será

$${\displaystyle \mathbb {E} [\Psi _{t}]=\int \rho (t,\mathbf {X} ,\mathbf {P} )\Psi (\mathbf {X} ,\mathbf {P} ){\rm {d}}\mathbf {P} {\rm {d}}\mathbf {X} }$$ Derivando con respecto al tiempo ${\displaystyle t}$ y aplicando la fórmula de Itô, obtenemos $${\displaystyle \mathbb {E} [{\frac {\rm {d}}{{\rm {d}}t}}\Psi (\mathbf {X} ,\mathbf {P} )]=\mathbb {E} [\nabla _{\mathbf {X} }\Psi {\frac {{\rm {d}}\mathbf {X} }{{\rm {d}}t}}+\nabla _{\mathbf {P} }\Psi {\frac {{\rm {d}}\mathbf {P} }{{\rm {d}}t}}+\sigma _{T}^{2}\nabla _{\mathbf {P} }^{2}\Psi {\frac {1}{{\rm {d}}t}}({\rm {d}}\mathbf {W} (t))^{2}]}$$ que se puede simplificar a $${\displaystyle \int ({\frac {\partial }{\partial t}}\rho )\Psi (\mathbf {X} ,\mathbf {P} ){\rm {d}}\mathbf {X} {\rm {d}}\mathbf {P} =\mathbb {E} [(\nabla _{\mathbf {X} }\Psi )\mathbf {P} +\nabla _{\mathbf {P} }\Psi (-\gamma \mathbf {P} -\nabla _{\mathbf {X} }U(\mathbf {X} ))+\sigma _{T}^{2}\nabla _{\mathbf {P} }^{2}\Psi ]}$$ Integración por partes en el lado derecho, debido a la densidad nula para el momento o la velocidad infinitos, tenemos $${\displaystyle ({\frac {\partial }{\partial t}}\rho )\Psi (\mathbf {X} ,\mathbf {P} ){\rm {d}}\mathbf {X} {\rm {d}}\mathbf {P} =\int (-\mathbf {P} \nabla _{\mathbf {X} }\rho +\nabla _{\mathbf {P} }(\gamma \mathbf {P} \rho +\nabla _{\mathbf {X} }U(\mathbf {X} )\rho )+\nabla _{\mathbf {P} }^{2}(\sigma _{T}^{2}\rho ))\Psi (\mathbf {X} ,\mathbf {P} ){\rm {d}}\mathbf {X} {\rm {d}}\mathbf {P} }$$ Esta ecuación se cumple para cualquier $\Psi$, por lo que requerimos que la densidad satisfaga $${\displaystyle {\frac {\partial \rho }{\partial t}}=-\mathbf {P} \nabla _{\mathbf {X} }\rho +\nabla _{\mathbf {P} }(\gamma \mathbf {P} \rho +\nabla _{\mathbf {X} }U(\mathbf {X} )\rho )+\nabla _{\mathbf {P} }^{2}(\sigma _{T}^{2}\rho )}$$ Esta ecuación se denomina ecuación de Klein-Kramers, una versión especial de la ecuación de Fokker-Planck. Es una ecuación diferencial parcial que describe la evolución de la densidad de probabilidad del sistema en el espacio de fases.

Para el límite sobredampeado, tenemos ${\displaystyle {\rm {d}}\mathbf {P} =0}$, por lo que la evolución del sistema puede reducirse al subespacio de posición. Siguiendo una lógica similar, podemos demostrar que la ecuación diferencial estocástica para la posición, $${\displaystyle {\rm {d}}\mathbf {X} =-{\frac {1}{\gamma }}\nabla U(\mathbf {X} ){\rm {d}}t+{\sqrt {2}}{\frac {\sigma }{\gamma }}\mathbf {R} (t){\rm {d}}t}$$

corresponde a la ecuación de Fokker-Planck para la densidad de probabilidad $${\displaystyle {\frac {\partial \rho (t,\mathbf {X} )}{\partial t}}=\nabla _{\mathbf {X} }({\frac {1}{\gamma }}\nabla _{\mathbf {X} }U(\mathbf {X} )\rho (t,\mathbf {X} ))+\Delta _{\mathbf {X} }({\frac {\sigma ^{2}}{\gamma ^{2}}}\rho (t,\mathbf {X} ))}$$

Consideremos la dinámica de Langevin de una partícula libre (es decir, ${\displaystyle U(\mathbf {X} )=0}$), entonces la ecuación para el momento será

$${\displaystyle {\rm {d}}\mathbf {P} =-{\frac {1}{\gamma }}\mathbf {P} {\rm {d}}t+{\frac {{\sqrt {2}}\sigma }{\gamma }}{\rm {d}}\mathbf {W} _{t}}$$

la solución analítica de esta SDE es

$${\displaystyle \mathbf {P} =\mathbf {P} _{0}e^{-t/\gamma }+{\frac {{\sqrt {2}}\sigma }{\gamma }}\int _{0}^{t}{\rm {e}}^{-(t-t')/\gamma }{\rm {d}}\mathbf {W} _{t}'}$$

por lo que el valor medio del segundo momento del momento será (aquí aplicamos la isometría de Itô) $${\displaystyle \mathbb {E} (\mathbf {P} ^{2})=\mathbf {P} _{0}^{2}{\rm {e}}^{-2t/\gamma }+{\frac {\sigma ^{2}}{\gamma }}(1-{\rm {e}}^{-2t/\gamma }){\overset {t\to \infty }{\to }}{\frac {\sigma ^{2}}{\gamma }}}$$ Es decir, el comportamiento límite cuando el tiempo se aproxima al infinito positivo, la fluctuación del momento de este sistema está relacionada con la disipación de energía (parámetro del término fricción ${\displaystyle \gamma }$) de este sistema. Combinando este resultado con el teorema de equipartición, que relaciona el valor medio de la energía cinética de las partículas con la temperatura

$${\displaystyle \langle v^{2}\rangle =k_{\rm {B}}T}$$

podemos determinar el valor de la varianza ${\displaystyle \sigma }$ en aplicaciones como el termostato de Langevin.

$${\displaystyle \sigma ^{2}/\gamma =k_{B}T\to \sigma ={\sqrt {k_{\rm {B}}T\gamma }}}$$ Esto es coherente con la definición original que asume ${\displaystyle M=1}$.

La formulación de la integral de camino proviene de la mecánica cuántica. Pero para una EDE de Langevin también podemos inducir una integral de camino correspondiente. Considerando la siguiente ecuación de Langevin sobredampeada bajo, donde sin pérdida de generalidad tomamos ${\displaystyle \gamma =\sigma =1}$,

$${\displaystyle {\rm {d}}{X}=-\nabla U({X}){\rm {d}}t+{\sqrt {2}}{\rm {d}}W_{t}}$$

Discretizamos y definimos ${\displaystyle t_{n}=n\Delta t}$, obtenemos

$${\displaystyle {X}_{n+1}-{X}_{n}+\nabla U({X})\Delta t={\sqrt {2}}(W_{t_{n}}-W_{t_{n-1}})\sim {\mathcal {N}}(0,2{\sqrt {\Delta t}})}$$ Por lo tanto, la probabilidad de propagación será $${\displaystyle P({X}_{n+1}|{X}_{n})=\int {\rm {d}}\xi {\frac {1}{2{\sqrt {\pi \Delta t}}}}{\rm {e}}^{-{\frac {\xi ^{2}}{4\Delta t}}}\delta ({X}_{n+1}-{X}_{n}+\nabla U({X})\Delta t-\xi )}$$ Aplicando la transformada de Fourier de la función delta, obtendremos $${\displaystyle P=\int {\frac {{\rm {d}}k}{2\pi }}{\rm {e}}^{{\rm {i}}k({X}_{n+1}-{X}_{n}+\nabla U({X})\Delta t)}\int {\rm {d}}\xi {\frac {1}{2{\sqrt {\pi \Delta t}}}}{\rm {e}}^{-{\frac {\xi ^{2}}{4\Delta t}}}{\rm {e}}^{-{\rm {i}}k\xi }}$$ La segunda parte es una integral gaussiana, que da como resultado

$${\displaystyle P=\int {\frac {{\rm {d}}k}{2\pi }}{\rm {e}}^{{\rm {i}}k({X}_{n+1}-{X}_{n}+\nabla U({X})\Delta t)}{\rm {e}}^{-k^{2}\Delta t}}$$ Ahora consideremos la probabilidad desde el valor inicial ${\displaystyle X_{0}}$ hasta el valor final ${\displaystyle X_{n}}$.

$${\displaystyle P(\mathbf {X} _{n}|\mathbf {X} _{0})=\int {\frac {1}{2\pi }}\prod _{i}^{N-1}{\rm {d}}k_{i}{\rm {e}}^{({\rm {i}}k_{i}({\dot {X}}+\nabla U(X))-k_{i}^{2})\Delta t}}$$ tomando el límite de ${\displaystyle \Delta t\to 0}$, obtendremos $${\displaystyle P(\mathbf {X} _{n}|\mathbf {X} _{0})=\int {\mathcal {D}}[k]{\rm {e}}^{\int _{0}^{t_{n}}({\rm {i}}k({\dot {X}}+\nabla U(X))-k^{2}){\rm {d}}t}}$$

• Mecánica hamiltoniana
• Mecánica estadística
• Solvación implícita
• Ecuación diferencial estocástica
• Ecuación de Langevin
• Langevin Monte Carlo
• Ecuación de Klein-Kramers

• Langevin Dynamics (LD) Simulation