Fix entropy calculation in gmx anaeig

gmx anaeig reads eigenvectors produced by gmx covar and
can compute entropy according to Schlitter's formula or
based on the quasiharmonic method. If the number of
eigenvectors is not consistent with the number of atoms
the entropy calculation could use uninitialized variables.
Added a warning when this happens.

Fixes #2668

Change-Id: I4f265212ce0a7bf82e25e1aa6f9cbb544c45db3f

diff --git a/docs/release-notes/2018/2018.4.rst b/docs/release-notes/2018/2018.4.rst
index b18982d6fb07a13f92bbf98777943d429d45f3bc..2f83fe826e026e8b8102c1ea9a51e3e12374781a 100644 (file)
--- a/docs/release-notes/2018/2018.4.rst
+++ b/docs/release-notes/2018/2018.4.rst
@@ -79,6 +79,15 @@ handles selections that include proline residues.
 
 :issue:`2701`
 
+Fix bug in entropy calculation in gmx anaeig
+""""""""""""""""""""""""""""""""""""""""""""
+
+When gmx anaeig received an inconsistent number of atoms and
+eigenvectors (fewer eigenvectors than three times the number of
+atoms) the entropy calculations would use uninitialized values.
+
+:issue:`2668`
+
 Fixes to improve portability
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

diff --git a/src/gromacs/gmxana/gmx_anaeig.cpp b/src/gromacs/gmxana/gmx_anaeig.cpp
index bd390ebab75912d0b634e6190b41261d613f33fc..e894b17ffd7b02adedd0fa0fe56b6ca37aa1e2fc 100644 (file)
--- a/src/gromacs/gmxana/gmx_anaeig.cpp
+++ b/src/gromacs/gmxana/gmx_anaeig.cpp
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -1164,7 +1164,13 @@ int gmx_anaeig(int argc, char *argv[])
     read_eigenvectors(VecFile, &natoms, &bFit1,
                       &xref1, &bDMR1, &xav1, &bDMA1,
                       &nvec1, &eignr1, &eigvec1, &eigval1);
-    neig1 = DIM*natoms;
+    neig1 = std::min(nvec1, DIM*natoms);
+    if (nvec1 != DIM*natoms)
+    {
+        fprintf(stderr, "Warning: number of eigenvectors %d does not match three times\n"
+                "the number of atoms %d in %s. Using %d eigenvectors.\n\n",
+                nvec1, natoms, VecFile, neig1);
+    }
 
     /* Overwrite eigenvalues from separate files if the user provides them */
     if (EigFile != nullptr)