How to calculate the nuclear energy derivatives in molecular mechanics?

Question

In molecular mechanics, the energy is often written as a sum of bonds, angles and torsion energies, and an electrostatic term, e.g.

$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}}]$

While this is pretty easy to evaluate, the nuclear energy derivatives (gradient) are more complex because they need to be written in the Cartesian coordinate system.

This reference:

First and second derivative matrix elements for the stretching, bending, and torsional energy, Kenneth J. Miller, Robert J. Hinde, and Janet Anderson, Journal of Computational Chemistry, vol 10, 63-76, 1989, https://www.onlinelibrary.wiley.com/doi/abs/10.1002/jcc.540100107

Outlines the mathematics for calculating the nuclear energy derivatives for the bonds, angles and torsions but it's rather complicated.

For example, I don't understand Table V.

Are there any open-source codes which implement this at a high-level?

I'd like to use this code to supplement the text and for educational purposes so the easier the code the better

Answer 1

You iterate through the bond, angle, torsion, etc. contributions and sum the gradients on each atom.

When we were implementing these for Open Babel, we found a very nice dissertation which gave some insight into MMFF94 gradients:

Dr. Andreas Moll BALLView : a molecular viewer and modeling tool

• Bonds - harmonic potential leads to, e.g.: $$ S_{i j} \frac{\hat{\mathbf{d}_{\mathrm{jj}}}}{\left|\mathbf{d}_{\mathrm{ij}}\right|} $$

• Angle Bend, e.g. $$ 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|} $$

• Torsion, e.g.; $$ 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|} $$

The code in Open Babel can be found starting here: forcefield.cpp

It's basically a bunch of cross-products and normalizations to get the force on each atom from a particular interaction.

OpenMM looks like they have their implementations here - fairly similar in approach.

Answer 2

Typically, force evaluation happens in Cartesian coordinates regardless of the approach (QM or MM). In QM codes the forces may then be projected into internal coordinates, which in turn is usually dependent on the geometry, see e.g. J. Chem. Phys. 110, 4986 (1999); even though this carries some cost, it is more than offset by the savings in fewer QM calculations. Also solids can be optimized in internal coordinates, see Chem. Phys. Lett. 335, 321 (2001). (AFAIK molecular dynamics codes don't use internal coordinates for geometry optimization, but I may be mistaken.)

The equations for the gradients and hessians for your force field are a bit complicated, but straightforward, since you're still in the Cartesian space. Table V is just a helper table of multiple-angle formulas for cosine:

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

and so on for $\cos 5\phi_{ijkl}$ and $\cos 6\phi_{ijkl}$. They say in the main text (in the section on torsional motion) that using these trigonometric identities avoids undefined values at $\phi_{ijkl}=0$ and $\phi_{ijkl}=\pi/2$ for the derivatives that get divided by $\sin \phi_{ijkl}$.

However, if you're hoping for high-level implementations - meaning simple and very intelligible code - my guess is that you're out of luck: since the evaluation of forces is a key bottleneck in MD codes, it is probably heavily optimized in all codes.

Your potential is quite simple, consisting of the harmonic bond and angle stretching, dihedral, and Lennard-Jones term, so it is probably available in almost any molecular dynamics code. (I think you are actually missing the Coulomb term that is typically included to model electrostatic interactions between non-bonded regions.)

GROMACS is one of the better-known open source molecular MD codes out there, and it is fast. It used to be written in plain C; I think it may have gotten some C++ later on. I haven't really looked at the source code in a decade...