Bug #2862

Division by zero in restrained dihedrals

Added by David van der Spoel 10 days ago. Updated 4 days ago.

Target version:
Affected version - extra info:
Affected version:


The following code in restcbt.cpp, line 154 can lead to division by zero (variable sine_phi_sq is set to zero).

/*    It is possible that cosine_phi is slightly bigger than 1.0 due to round-off errors. */
if (sine_phi_sq < 0.0) {
sine_phi_sq = 0.0;
< snip >
/*      Computation of the prefactor - common term for all forces */
*prefactor_phi = -(k_torsion) * delta_cosine * norm_phi * term_phi_phi0 / (sine_phi_sq * sine_phi_sq);

Associated revisions

Revision dfa21149 (diff)
Added by David van der Spoel 4 days ago

Precision fix for rescbt code.

The compute_restangles function was not sufficiently
precise. In fact it return completely different
results in single and double precision. By making part
of the function work in double the issue is fixed.

New tests added.

Part of #2795
Fixes #2862

Change-Id: Iad7f2bc45be996ba3e16358aab838c5427b157b8


#1 Updated by Mark Abraham 10 days ago

Nico, should we improve something here?

#2 Updated by Berk Hess 10 days ago

Unrelated to the division by zero, but the factor (sine_theta_sq)^2 in the denominator of prefactor looks suspicious to me. I would expect a factor sin^2, not a factor sin^4.

#3 Updated by Berk Hess 10 days ago

Maybe the sin^4 factor could be correct, but that is very difficult to verify without the definitions of any of the variables in the code.

The potential diverges at theta=0, so you should simply never put in such a configuration!

#4 Updated by David van der Spoel 8 days ago

Here is an issue related to the same code, where single and double produce completely different results.

#5 Updated by David van der Spoel 8 days ago

For sure the code is not sufficiently robust. Maybe some variables need to be made double precision.

#6 Updated by Gerrit Code Review Bot 8 days ago

Gerrit received a related patchset '1' for Issue #2862.
Uploader: David van der Spoel ()
Change-Id: gromacs~master~Iad7f2bc45be996ba3e16358aab838c5427b157b8
Gerrit URL:

#7 Updated by David van der Spoel 4 days ago

  • Status changed from New to Resolved

Also available in: Atom PDF