Actually fix gmx helix segmentation fault

The residue indexing in hxprops was wrong because the residue indices
used for residue name lookup were wrong.

Added some documentation to clarify what this code does.

Fixes #2701

Change-Id: Ie6c002be6ff2139f6af27c7b891a180338ae5b20

diff --git a/docs/release-notes/2018/2018.4.rst b/docs/release-notes/2018/2018.4.rst
index 90fdddef108a740669788a996dd64acee795a1cd..79e90028e7e717c8b315d1dd2401d36c1e188ead 100644 (file)
--- a/docs/release-notes/2018/2018.4.rst
+++ b/docs/release-notes/2018/2018.4.rst
@@ -60,14 +60,14 @@ empty file name strings.
 
 :issue:`2653`
 
-Fix helix segmentation fault
+Fix gmx helix segmentation faults
 """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
-The .tpr file is now read correctly. 
+The .tpr file is now read correctly, and the helix analysis correctly
+handles selections that include proline residues.
 
 :issue:`2701`
 
-
 Fixes to improve portability
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

diff --git a/src/gromacs/gmxana/hxprops.cpp b/src/gromacs/gmxana/hxprops.cpp
index 9c6d7aea0a09b7951e7febf03d70b8ab29a10d41..bf8fbd6cdec55c7614a7d58369ae2d734e6bd71b 100644 (file)
--- a/src/gromacs/gmxana/hxprops.cpp
+++ b/src/gromacs/gmxana/hxprops.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.
@@ -339,7 +339,7 @@ t_bb *mkbbind(const char *fn, int *nres, int *nbb, int res0,
 #define NBB asize(bb_nm)
     t_bb               *bb;
     char               *grpname;
-    int                 ai, i, i0, i1, j, k, ri, rnr, gnx, r0, r1;
+    int                 ai, i, i0, i1, j, k, rnr, gnx, r0, r1;
 
     fprintf(stderr, "Please select a group containing the entire backbone\n");
     rd_index(fn, 1, &gnx, index, &grpname);
@@ -362,40 +362,45 @@ t_bb *mkbbind(const char *fn, int *nres, int *nbb, int res0,
     for (i = j = 0; (i < gnx); i++)
     {
         ai = (*index)[i];
-        ri = atom[ai].resind-r0;
-        if (std::strcmp(*(resinfo[ri].name), "PRO") == 0)
+        // Create an index into the residue index for the topology.
+        int resindex = atom[ai].resind;
+        // Create an index into the residues present in the selected
+        // index group.
+        int bbindex  = resindex -r0;
+        if (std::strcmp(*(resinfo[resindex].name), "PRO") == 0)
         {
+            // For PRO in a peptide, there is no H bound to backbone
+            // N, so use CD instead.
             if (std::strcmp(*(atomname[ai]), "CD") == 0)
             {
-                bb[ri].H = ai;
+                bb[bbindex].H = ai;
             }
         }
         for (k = 0; (k < NBB); k++)
         {
             if (std::strcmp(bb_nm[k], *(atomname[ai])) == 0)
             {
-                j++;
                 break;
             }
         }
         switch (k)
         {
             case 0:
-                bb[ri].N = ai;
+                bb[bbindex].N = ai;
                 break;
             case 1:
             case 5:
                 /* No attempt to address the case where some weird input has both H and HN atoms in the group */
-                bb[ri].H = ai;
+                bb[bbindex].H = ai;
                 break;
             case 2:
-                bb[ri].CA = ai;
+                bb[bbindex].CA = ai;
                 break;
             case 3:
-                bb[ri].C = ai;
+                bb[bbindex].C = ai;
                 break;
             case 4:
-                bb[ri].O = ai;
+                bb[bbindex].O = ai;
                 break;
             default:
                 break;
@@ -449,8 +454,8 @@ t_bb *mkbbind(const char *fn, int *nres, int *nbb, int res0,
     /* Set the labels */
     for (i = 0; (i < rnr); i++)
     {
-        ri = atom[bb[i].CA].resind;
-        sprintf(bb[i].label, "%s%d", *(resinfo[ri].name), ri+res0);
+        int resindex = atom[bb[i].CA].resind;
+        sprintf(bb[i].label, "%s%d", *(resinfo[resindex].name), resinfo[resindex].nr);
     }
 
     *nres = rnr;

diff --git a/src/gromacs/gmxana/hxprops.h b/src/gromacs/gmxana/hxprops.h
index 88de40811cbb4af95529ecf9f30fffc15212e7aa..25f2e492691c847e69e17de42c117a09abdbdf26 100644 (file)
--- a/src/gromacs/gmxana/hxprops.h
+++ b/src/gromacs/gmxana/hxprops.h
@@ -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, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,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.
@@ -58,16 +58,49 @@ extern "C"
 /* Canonical values of the helix phi/psi angles */
 
 
-typedef struct {
-    real     phi, psi, pprms2;
+/*! \internal \brief Struct containing properties of a residue in a protein backbone. */
+struct t_bb {
+    //! Protein backbone phi angle.
+    real     phi;
+    //! Protein backbone psi angle.
+    real     psi;
+    //! RMS distance of phi and psi angles from ideal helix
+    real     pprms2;
+    //! Estimated J-coupling value
     real     jcaha;
-    real     d3, d4, d5, rmsa;
+    //! Value of 3 turn helix?
+    real     d3;
+    //! Value of 4 turn helix?
+    real     d4;
+    //! Value of 5 turn?
+    real     d5;
+    //! Average of RMS for analysis.
+    real     rmsa;
+    //! If the structure is helical.
     gmx_bool bHelix;
+    //! Number of elliptical elements
     int      nhx;
-    int      nrms, resno;
-    int      Cprev, N, H, CA, C, O, Nnext;
+    //! Average RMS Deviation when atoms of this residue are fitted to ideal helix
+    int      nrms;
+    //! Residue index for output, relative to gmx_helix -r0 value
+    int      resno;
+    //! Index for previous carbon.
+    int      Cprev;
+    //! Index for backbone nitrogen.
+    int      N;
+    //! Index for backbone NH hydrogen.
+    int      H;
+    //! Index for alpha carbon.
+    int      CA;
+    //! Index for carbonyl carbon.
+    int      C;
+    //! Index for carbonyl oxygen.
+    int      O;
+    //! Index for next backbone nitrogen.
+    int      Nnext;
+    //! Name for this residue.
     char     label[32];
-} t_bb;
+};
 
 enum {
     efhRAD,  efhTWIST, efhRISE, efhLEN,
@@ -109,6 +142,17 @@ extern void av_hblen(FILE *fp3, FILE *fp3a,
 extern void av_phipsi(FILE *fphi, FILE *fpsi, FILE *fphi2, FILE *fpsi2,
                       real t, int nres, t_bb bb[]);
 
+/*! \brief Allocate and fill an array of information about residues in a protein backbone.
+ *
+ * The user is propted for an index group of protein residues (little
+ * error checking occurs). For the number of residues found in the
+ * selected group, nbb entries are made in the returned array.  Each
+ * entry contains the atom indices of the N, H, CA, C and O atoms (for
+ * PRO, H means CD), as well as the C of the previous residue and the
+ * N of the next (-1 if not found).
+ *
+ * In the output array, the first residue will be numbered starting
+ * from res0. */
 extern t_bb *mkbbind(const char *fn, int *nres, int *nbb, int res0,
                      int *nall, int **index,
                      char ***atomname, t_atom atom[],