¿Cómo calcular los derivados de la energía nuclear en mecánica molecular?
En mecánica molecular, la energía a menudo se escribe como una suma de enlaces, ángulos y energías de torsión, y un término electrostático, p. Ej.
$V = \sum_{bonds} K_r (r-r_{eq})^2 +\sum_{angles}K_{\theta}(\theta -\theta_0)^2 +\sum_{dihedrals}V_n/2 [1+cos(n\phi -\gamma)] + \sum_{i<j}[\frac{A_{ij}}{R_{ij}^{12}} + \frac{B_{ij}}{R_{ij}^{6}} + \frac{q_1 q_2}{\epsilon R_{ij}}]$
Si bien esto es bastante fácil de evaluar, las derivadas de la energía nuclear (gradiente) son más complejas porque deben escribirse en el sistema de coordenadas cartesianas.
Esta referencia:
Elementos de matriz de primera y segunda derivada para la energía de estiramiento, flexión y torsión, Kenneth J. Miller, Robert J. Hinde y Janet Anderson, Journal of Computational Chemistry, vol 10, 63-76, 1989, https://www.onlinelibrary.wiley.com/doi/abs/10.1002/jcc.540100107
Describe las matemáticas para calcular las derivadas de energía nuclear para los enlaces, ángulos y torsiones, pero es bastante complicado.
Por ejemplo, no entiendo la Tabla V.
¿Existen códigos de fuente abierta que implementen esto a un alto nivel?
Me gustaría usar este código para complementar el texto y con fines educativos, de modo que cuanto más fácil sea el código, mejor.
Respuestas
Se itera a través de las contribuciones de enlace, ángulo, torsión, etc. y se suman los gradientes en cada átomo.
Cuando estábamos implementando estos para Open Babel, encontramos una disertación muy agradable que dio una idea de los gradientes de MMFF94:
Dr. Andreas Moll BALLView: un visor molecular y una herramienta de modelado
Enlaces: el potencial armónico conduce a, por ejemplo: $$ S_{i j} \frac{\hat{\mathbf{d}_{\mathrm{jj}}}}{\left|\mathbf{d}_{\mathrm{ij}}\right|} $$
Curva de ángulo, p. Ej. $$ B_{i j k} \frac{\hat{\mathbf{d}}_{\mathrm{ij}} \times \hat{\mathbf{d}}_{\mathrm{ki}} \times \hat{\mathbf{d}}_{\mathrm{ij}}}{\left|\mathbf{d}_{\mathrm{ij}}\right|} $$
Torsión, por ejemplo; $$ T_{i j k l} \frac{-\hat{\mathbf{d}}_{\mathrm{ij}} \times \hat{\mathrm{d}}_{\mathrm{jk}}}{\sin (\phi)^{2}\left|\mathbf{d}_{\mathrm{ij}}\right|} $$
El código en Open Babel se puede encontrar comenzando aquí: forcefield.cpp
Es básicamente un montón de productos cruzados y normalizaciones para obtener la fuerza en cada átomo de una interacción particular.
OpenMM parece que tienen sus implementaciones aquí , bastante similar en enfoque.
Normalmente, la evaluación de la fuerza se realiza en coordenadas cartesianas independientemente del enfoque (QM o MM). En los códigos QM, las fuerzas pueden luego proyectarse en coordenadas internas, que a su vez suelen depender de la geometría, véase, por ejemplo, J. Chem. Phys. 110, 4986 (1999) ; aunque esto conlleva algún costo, se compensa con creces por los ahorros en menos cálculos de QM. También los sólidos se pueden optimizar en coordenadas internas, ver Chem. Phys. Letón. 335, 321 (2001) . (Los códigos de dinámica molecular AFAIK no usan coordenadas internas para la optimización de la geometría, pero puedo estar equivocado).
Las ecuaciones para los gradientes y arpilleras para su campo de fuerza son un poco complicadas, pero sencillas, ya que todavía se encuentra en el espacio cartesiano. La Tabla V es solo una tabla auxiliar de fórmulas de múltiples ángulos para el coseno:
$\begin{aligned}\cos \phi_{ijkl} &= \cos \phi_{ijkl} \\ \cos 2\phi_{ijkl} &= -1 + 2 \cos^2 \phi_{ijkl} \\ \cos 3\phi_{ijkl} &= -3\cos \phi_{ijkl} + 4 \cos^3 \phi_{ijkl} \\ \cos 4\phi_{ijkl} &= 1 -8\cos^2 \phi_{ijkl} + 8 \cos^4 \phi_{ijkl} \end{aligned}$
y así sucesivamente para $\cos 5\phi_{ijkl}$ y $\cos 6\phi_{ijkl}$. Dicen en el texto principal (en la sección sobre movimiento de torsión) que el uso de estas identidades trigonométricas evita valores indefinidos en$\phi_{ijkl}=0$ y $\phi_{ijkl}=\pi/2$ para las derivadas que se dividen por $\sin \phi_{ijkl}$.
Sin embargo, si espera implementaciones de alto nivel, es decir, código simple y muy inteligible, supongo que no tiene suerte: dado que la evaluación de fuerzas es un cuello de botella clave en los códigos MD, probablemente esté muy optimizado en todos los códigos.
Su potencial es bastante simple, que consiste en el enlace armónico y el estiramiento del ángulo, diedro y el término de Lennard-Jones, por lo que probablemente esté disponible en casi cualquier código de dinámica molecular. (Creo que en realidad se está perdiendo el término de Coulomb que normalmente se incluye para modelar interacciones electrostáticas entre regiones no unidas).
GROMACS es uno de los códigos MD moleculares de código abierto más conocidos que existen, y es rápido. Solía estar escrito en C simple; Creo que pudo haber obtenido algo de C ++ más adelante. Realmente no he mirado el código fuente en una década ...