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

[Mark Abraham]

David Mobley wrote:
> 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.

Hmmmm, but arccos has domain [-1,1] and range [0,pi], and bondfree.c is 
calling acos() from the standard math library and that's the range it 
returns.

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

My brain is telling me you have a point, but that this isn't the fix. 
The value of theta returned by the above pseudocode is correct in all 
cases - any bond angle near +pi will give theta in the range 
[+pi-delta,+pi], and so V will return sensible values for theta_0 = +pi. 
Perhaps the problem lies in some assumption in the implementation of the 
force component calculation?

Mark