Tasa de aceptación de Bennett

Algoritmo para estimar la diferencia de energía libre entre dos sistemas

El método de la relación de aceptación de Bennett ( BAR ) es un algoritmo para estimar la diferencia de energía libre entre dos sistemas (normalmente, los sistemas se simulan en la computadora). Fue sugerido por Charles H. Bennett en 1976. ^[1]

Preliminares

Tomemos un sistema en un determinado estado super (es decir, Gibbs). Al realizar un recorrido de Monte Carlo por Metropolis, es posible muestrear el paisaje de estados entre los que se mueve el sistema, utilizando la ecuación

${\displaystyle p({\text{State}}_{x}\rightarrow {\text{State}}_{y})=\min \left(e^{-\beta \,\Delta U},1\right)=M(\beta \,\Delta U)}$

donde Δ U =  U (Estado _y ) −  U (Estado _x ) es la diferencia en energía potencial, β = 1/ kT ( T es la temperatura en kelvins , mientras que k es la constante de Boltzmann ), y es la función Metropolis. Los estados resultantes se muestrean luego de acuerdo con la distribución de Boltzmann del superestado a la temperatura T. Alternativamente, si el sistema se simula dinámicamente en el conjunto canónico (también llamado conjunto NVT ), los estados resultantes a lo largo de la trayectoria simulada se distribuyen de la misma manera. El promedio a lo largo de la trayectoria (en cualquier formulación) se denota mediante corchetes angulares . ${\displaystyle M(x)\equiv \min(e^{-x},1)}$ ${\displaystyle \left\langle \cdots \right\rangle }$

Supongamos que se dan dos superestados de interés, A y B. Suponemos que tienen un espacio de configuración común, es decir, comparten todos sus microestados, pero las energías asociadas a estos (y, por lo tanto, las probabilidades) difieren debido a un cambio en algún parámetro (como la fuerza de una determinada interacción). La pregunta básica que debe abordarse es, entonces, ¿cómo se puede calcular el cambio de energía libre de Helmholtz (Δ F  =  F _B  −  F _A ) al moverse entre los dos superestados a partir del muestreo en ambos conjuntos? La parte de energía cinética en la energía libre es igual entre los estados, por lo que se puede ignorar. Además, la energía libre de Gibbs corresponde al conjunto NpT .

El caso general

Bennett muestra que para cada función f que satisface la condición (que es esencialmente la condición de equilibrio detallada ), y para cada compensación de energía C , se tiene la relación exacta ${\displaystyle f(x)/f(-x)\equiv e^{-x}}$

${\displaystyle e^{-\beta (\Delta F-C)}={\frac {\left\langle f\left(\beta (U_{\text{B}}-U_{\text{A}}-C)\right)\right\rangle _{\text{A}}}{\left\langle f\left(\beta (U_{\text{A}}-U_{\text{B}}+C)\right)\right\rangle _{\text{B}}}}}$

donde U _A y U _B son las energías potenciales de las mismas configuraciones, calculadas utilizando la función potencial A (cuando el sistema está en el superestado A) y la función potencial B (cuando el sistema está en el superestado B) respectivamente.

El caso básico

Sustituyendo f por la función Metropolis definida anteriormente (que satisface la condición de equilibrio detallada) y fijando C en cero, se obtiene

${\displaystyle e^{-\beta \Delta F}={\frac {\left\langle M\left(\beta (U_{\text{B}}-U_{\text{A}})\right)\right\rangle _{\text{A}}}{\left\langle M\left(\beta (U_{\text{A}}-U_{\text{B}})\right)\right\rangle _{\text{B}}}}}$

La ventaja de esta formulación (aparte de su simplicidad) es que se puede calcular sin realizar dos simulaciones, una en cada conjunto específico. De hecho, es posible definir un tipo adicional de movimiento de prueba de Metropolis de "cambio potencial" (que se realiza cada un número fijo de pasos), de modo que el muestreo único del conjunto "mixto" sea suficiente para el cálculo.

El caso más eficiente

Bennett explora qué expresión específica para Δ F es la más eficiente, en el sentido de producir el error estándar más pequeño para un tiempo de simulación determinado. Demuestra que la opción óptima es tomar

1. ${\displaystyle f(x)\equiv {\frac {1}{1+e^{x}}}}$, que es esencialmente la distribución de Fermi-Dirac (que de hecho satisface la condición de equilibrio detallada).
2. ${\displaystyle C\approx \Delta F}$Este valor, por supuesto, no se conoce (es exactamente lo que se intenta calcular), pero se puede elegir aproximadamente de manera autoconsistente.

Algunas suposiciones necesarias para la eficiencia son las siguientes:

1. Las densidades de los dos superestados (en su espacio de configuración común) deberían tener una gran superposición. De lo contrario, puede ser necesaria una cadena de superestados entre A y B, de modo que la superposición de cada dos superestados consecutivos sea adecuada.
2. El tamaño de la muestra debe ser grande. En particular, como los estados sucesivos están correlacionados, el tiempo de simulación debe ser mucho mayor que el tiempo de correlación.
3. El costo de simular ambos conjuntos debería ser aproximadamente igual y, de hecho, el sistema se muestrea aproximadamente por igual en ambos superestados. De lo contrario, se modifica la expresión óptima para C y el muestreo debería dedicar tiempos iguales (en lugar de un número igual de pasos de tiempo) a los dos conjuntos.

Tasa de aceptación de Bennett en varios estados

El índice de aceptación de Bennett multiestado ( MBAR ) es una generalización del índice de aceptación de Bennett que calcula las energías libres (relativas) de varios multiestados. Básicamente, se reduce al método BAR cuando solo intervienen dos superestados.

Relación con otros métodos

El método de la teoría de perturbaciones

Este método, también llamado perturbación de energía libre (o FEP), implica el muestreo únicamente del estado A. Requiere que todas las configuraciones de alta probabilidad del superestado B estén contenidas en las configuraciones de alta probabilidad del superestado A, lo que es un requisito mucho más estricto que la condición de superposición mencionada anteriormente.

El resultado exacto (orden infinito)

${\displaystyle e^{-\beta \Delta F}=\left\langle e^{-\beta (U_{\text{B}}-U_{\text{A}})}\right\rangle _{\text{A}}}$

o

${\displaystyle \Delta F=-kT\cdot \log \left\langle e^{\beta (U_{\text{A}}-U_{\text{B}})}\right\rangle _{\text{A}}}$

Este resultado exacto se puede obtener a partir del método BAR general, utilizando (por ejemplo) la función Metropolis, en el límite . De hecho, en ese caso, el denominador de la expresión del caso general anterior tiende a 1, mientras que el numerador tiende a . Sin embargo, una derivación directa a partir de las definiciones es más sencilla. ${\displaystyle C\rightarrow -\infty }$ ${\displaystyle e^{\beta C}\left\langle e^{-\beta (U_{\text{B}}-U_{\text{A}})}\right\rangle _{\text{A}}}$

El resultado de segundo orden (aproximado)

Suponiendo que Taylor expande la segunda expresión exacta de la teoría de perturbación al segundo orden, se obtiene la aproximación ${\displaystyle U_{\text{B}}-U_{\text{A}}\ll kT}$

${\displaystyle \Delta F\approx \left\langle U_{\text{B}}-U_{\text{A}}\right\rangle _{\text{A}}-{\frac {\beta }{2}}\left(\left\langle (U_{\text{B}}-U_{\text{A}})^{2}\right\rangle _{\text{A}}-\left(\left\langle (U_{\text{B}}-U_{\text{A}})\right\rangle _{\text{A}}\right)^{2}\right)}$

Nótese que el primer término es el valor esperado de la diferencia de energía, mientras que el segundo es esencialmente su varianza.

Las desigualdades de primer orden

Utilizando la convexidad de la función logarítmica que aparece en el resultado del análisis de perturbación exacta, junto con la desigualdad de Jensen , se obtiene una desigualdad en el nivel lineal; combinado con el resultado análogo para el conjunto B se obtiene la siguiente versión de la desigualdad de Gibbs-Bogoliubov :

${\displaystyle \langle U_{\text{B}}-U_{\text{A}}\rangle _{\text{B}}\leq \Delta F\leq \langle U_{\text{B}}-U_{\text{A}}\rangle _{\text{A}}}$

Nótese que la desigualdad concuerda con el signo negativo del coeficiente del término de varianza (positivo) en el resultado de segundo orden.

El método de integración termodinámica

Escribiendo la energía potencial como dependiente de un parámetro continuo, ${\displaystyle U_{\text{A}}=U(\lambda =0),U_{\text{B}}=U(\lambda =1),}$

se tiene el resultado exacto Esto se puede verificar directamente a partir de las definiciones o verse a partir del límite de las desigualdades de Gibbs-Bogoliubov anteriores cuando . por lo tanto, podemos escribir ${\displaystyle {\frac {\partial F(\lambda )}{\partial \lambda }}=\left\langle {\frac {\partial U(\lambda )}{\partial \lambda }}\right\rangle _{\lambda }}$ ${\displaystyle {\text{A}}=\lambda ^{+},{\text{B}}=\lambda ^{-}}$

${\displaystyle \Delta F=\int _{0}^{1}\left\langle {\frac {\partial U(\lambda )}{\partial \lambda }}\right\rangle \,d\lambda }$

que es el resultado de la integración termodinámica (o TI). Se puede aproximar dividiendo el rango entre los estados A y B en muchos valores de λ en los que se estima el valor esperado y realizando una integración numérica.

Implementación

El método de la relación de aceptación de Bennett se implementa en sistemas de dinámica molecular modernos , como Gromacs . El código basado en Python para MBAR y BAR está disponible para descargar en [2].

Véase también

• Templado paralelo

Referencias

1. Charles H. Bennett (1976) Estimación eficiente de las diferencias de energía libre a partir de datos de Monte Carlo. Journal of Computational Physics 22: 245–268 [1]

Enlaces externos

• Tasa de aceptación de Bennett de AlchemistryWiki.
• Índice de aceptación de Bennett multiestatal de AlchemistryWiki.
• Método de análisis de histograma ponderado (MBAR es el caso sin agrupar) de AlchemistryWiki.