Actually fix gmx wham with angle geometries

A previous fix for converting the force constant in gmx wham
for angle and dihedral geometry applied an incorrect conversion
factor.

Fixes #2609 and #3094

Change-Id: I0aa9b848c7fa62306ccf15b327a650b1ad88e03f

diff --git a/docs/release-notes/2019/2019.3.rst b/docs/release-notes/2019/2019.3.rst
index c5e7b7cd83b10ce8dcb199380ffa199032ee9cec..80f501cb4da12cb1fc7f5be07f77511217b231d7 100644 (file)
--- a/docs/release-notes/2019/2019.3.rst
+++ b/docs/release-notes/2019/2019.3.rst
@@ -67,9 +67,11 @@ Fix gmx wham with angle geometries
 """"""""""""""""""""""""""""""""""
 
 gmx wham would mix up degree and radian units leading to no overlap
-or not-a-number output.
+or not-a-number output. **Note**: this fix is not correct, a correct
+fix is applied in the next patch release.
 
 :issue:`2609`
+:issue:`3094`
 
 Add some information for grompp error with wrong line endings
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

diff --git a/docs/release-notes/2019/2019.4.rst b/docs/release-notes/2019/2019.4.rst
index c9210825130f292e43bace29d26f573fdc474aff..da735b6132278a73979f7d3b35d4ad4971e2ce82 100644 (file)
--- a/docs/release-notes/2019/2019.4.rst
+++ b/docs/release-notes/2019/2019.4.rst
@@ -72,6 +72,15 @@ but not replaced by angle constraints.
 
 :issue:`3067`
 
+Fix gmx wham with angle and dihedral geometries
+"""""""""""""""""""""""""""""""""""""""""""""""
+
+gmx wham would apply an incorrect radian to degree unit conversion,
+leading to no overlap or not-a-number output.
+
+:issue:`2609`
+:issue:`3094`
+
 Fix bug in gmx xpm2ps
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 

diff --git a/src/gromacs/gmxana/gmx_wham.cpp b/src/gromacs/gmxana/gmx_wham.cpp
index 2bcc8e145b53ec18e2370a6931a93beeb3676410..5624c0ddb7788cbc39af7fd4c8ee45b9bb1a329c 100644 (file)
--- a/src/gromacs/gmxana/gmx_wham.cpp
+++ b/src/gromacs/gmxana/gmx_wham.cpp
@@ -2050,8 +2050,12 @@ static void read_tpr_header(const char *fn, t_UmbrellaHeader* header, t_Umbrella
         header->pcrd[i].pull_type     = ir->pull->coord[i].eType;
         header->pcrd[i].geometry      = ir->pull->coord[i].eGeom;
         header->pcrd[i].ngroup        = ir->pull->coord[i].ngroup;
-        /* Angle type coordinates are handled fully in degrees in gmx wham */
-        header->pcrd[i].k             = ir->pull->coord[i].k*pull_conversion_factor_internal2userinput(&ir->pull->coord[i]);
+        /* Angle type coordinates are handled fully in degrees in gmx wham.
+         * The pull "distance" appears with a power of -2 in the force constant,
+         * so we need to multiply with the internal units (radians for angle)
+         * to user units (degrees for an angle) with the same power.
+         */
+        header->pcrd[i].k             = ir->pull->coord[i].k/gmx::square(pull_conversion_factor_internal2userinput(&ir->pull->coord[i]));
         header->pcrd[i].init_dist     = ir->pull->coord[i].init;
 
         copy_ivec(ir->pull->coord[i].dim, header->pcrd[i].dim);