# How to calculate the nuclear energy derivatives in molecular mechanics?

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

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

• There is no separation between electrons and nucleus in Molecular Mechanics. They are both "fuse" in balls. Also, why the need to get the equations in Cartesian coordinates? The torsion (as all of other therms) are well represented in internal or generalized coordinates.
– Camps
Aug 21 '20 at 18:21
• What is "high-level"?
– Camps
Aug 21 '20 at 18:23
• 1) the reason you need the gradient in Cartesian coordinates is because internal coordinates are redundant and most optimizers or molecular dynamics work in Cartesian coordinates. 2) For a definition of high-level see high-level programming language. So for example, I'd consider python a high-level and C++ a little lower because you need to define memory and stuff like that. Aug 21 '20 at 18:26
• @CodyAldaz which text are you supplementing? You mean just to understand the Miller et al paper better? Aug 21 '20 at 20:06
• @SusiLehtola The Miller paper which is referenced since I'm having a hard time doing all the calculus. But I guess it would be a good supplement for any molecular mechanics work Aug 21 '20 at 20:08

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.

• +1 Wow! Two excellent answers. I should have thought to check OpenBabel, that code is very informative. Thanks Aug 21 '20 at 22:14
• @CodyAldaz - that file needs to be split, but that part is readable. Aug 21 '20 at 22:15
• @TristanMaxson - thanks - should be fixed now - accidentally added an extra '/' Aug 22 '20 at 17:04

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...

• +1 nice information! I was missing the electrostatic term. I figured here might not be any simple implementations. But I'm still hopeful since it is an excellent exercise, and relatively easier than SCF. Aug 21 '20 at 20:01