[gmx-users] problem with linear molecules and suggested fix
Mark.Abraham at anu.edu.au
Wed Oct 17 10:09:18 CEST 2007
David Mobley wrote:
> 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
> 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
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?
More information about the gromacs.org_gmx-users