[gmx-users] problem with linear molecules and suggested fix

[David Mobley]

All,

Earlier, I wrote about how I ran a large set of small molecules using
GAFF parameters and all of my small molecules with triple bonds
(nitriles, pentyne, propyne, etc) were crashing. This was not a
parameter problem, and not a timestep problem, and, perplexingly, only
seems to occur for these molecules when they are in vacuum or water,
but not in a binding site.

As it turns out, we (John Chodera, Chris Fennell and I) have now
traced the problem back to problems with the angles, and in
particular, how GROMACS handles angle bending for nearly linear
molecules. These molecules, since they have triple bonds, are
basically linear, of course.

Anyway, as it turns out, all of these molecules have one or several
angles with preferred angles of nearly pi (or 180, whichever units you
prefer). As it turns out, the angle bending code apparently can't
handle angles of exactly 180 degrees, which causes the problem. This
also explains why these molecules are OK when run in a binding site:
The binding site bends the angle so it isn't 180.

In gmxlib/bondfree.c, the relevant algorithm (in pseudocode) is
something like this (where
x_i, x_j, and x_k are atom vectors):

t_1 = x_i - x_j
t_2 = x_k - x_j
cos_theta = dot(t_1, t_2) / (|t_1| |t_2|)
theta = inverse_cosine(cos_theta)
V = (1/2) K (theta - theta_0)^2

Now, if theta_0 is anywhere near pi, this obviously isn't going to
work if inverse_cosine is defined on the domain [-pi,+pi).  Hence, the
problem we're seeing.

The solution is simple, as John points out to me: shift this
inverse_cosine domain to cover a range of [-pi,+pi) around theta_0.
Then only in the (extremely rare) case that a molecule is so strained
that it is pi away from its preferred value would there be
instabilities.

I will also submit a bugzilla.

Thank you,
David Mobley
UCSF
http://www.dillgroup.ucsf.edu/~dmobley