Merge release-5-0 into master

[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / pdb2top.cpp

diff --git a/src/gromacs/gmxpreprocess/pdb2top.cpp b/src/gromacs/gmxpreprocess/pdb2top.cpp
index 7b08b863c42bc69caf7e78d21223ceba7401bb0b..89f059f7fd8c877edb443396eb590fa7e2035b17 100644 (file)
--- a/src/gromacs/gmxpreprocess/pdb2top.cpp
+++ b/src/gromacs/gmxpreprocess/pdb2top.cpp
@@ -734,22 +734,7 @@ void write_top(FILE *out, char *pr, char *molname,
     }
 }
 
-static atom_id search_res_atom(const char *type, int resind,
-                               t_atoms *atoms,
-                               const char *bondtype, gmx_bool bAllowMissing)
-{
-    int i;
-
-    for (i = 0; (i < atoms->nr); i++)
-    {
-        if (atoms->atom[i].resind == resind)
-        {
-            return search_atom(type, i, atoms, bondtype, bAllowMissing);
-        }
-    }
 
-    return NO_ATID;
-}
 
 static void do_ssbonds(t_params *ps, t_atoms *atoms,
                        int nssbonds, t_ssbond *ssbonds, gmx_bool bAllowMissing)
@@ -815,6 +800,7 @@ static void at2bonds(t_params *psb, t_hackblock *hb,
             {
                 dist2 = distance2(x[ai], x[aj]);
                 if (dist2 > long_bond_dist2)
+
                 {
                     fprintf(stderr, "Warning: Long Bond (%d-%d = %g nm)\n",
                             ai+1, aj+1, sqrt(dist2));
@@ -1425,10 +1411,11 @@ void match_atomnames_with_rtp(t_restp restp[], t_hackblock hb[],
 }
 
 #define NUM_CMAP_ATOMS 5
-static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuetype_t *rt)
+static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms)
 {
     int         residx, i, j, k;
     const char *ptr;
+    const char *pname;
     t_resinfo  *resinfo = atoms->resinfo;
     int         nres    = atoms->nres;
     gmx_bool    bAddCMAP;
@@ -1444,11 +1431,14 @@ static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuet
         ptr = "check";
     }
 
-    fprintf(stderr, "Making cmap torsions...");
+    fprintf(stderr, "Making cmap torsions...\n");
     i = 0;
-    /* End loop at nres-1, since the very last residue does not have a +N atom, and
-     * therefore we get a valgrind invalid 4 byte read error with atom am */
-    for (residx = 0; residx < nres-1; residx++)
+    /* Most cmap entries use the N atom from the next residue, so the last
+     * residue should not have its CMAP entry in that case, but for things like
+     * dipeptides we sometimes define a complete CMAP entry inside a residue,
+     * and in this case we need to process everything through the last residue.
+     */
+    for (residx = 0; residx < nres; residx++)
     {
         /* Add CMAP terms from the list of CMAP interactions */
         for (j = 0; j < restp[residx].rb[ebtsCMAP].nb; j++)
@@ -1459,7 +1449,21 @@ static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuet
              * from residues labelled as protein. */
             for (k = 0; k < NUM_CMAP_ATOMS && bAddCMAP; k++)
             {
-                cmap_atomid[k] = search_atom(restp[residx].rb[ebtsCMAP].b[j].a[k],
+                /* Assign the pointer to the name of the next reference atom.
+                 * This can use -/+ labels to refer to previous/next residue.
+                 */
+                pname = restp[residx].rb[ebtsCMAP].b[j].a[k];
+                /* Skip this CMAP entry if it refers to residues before the
+                 * first or after the last residue.
+                 */
+                if (((strchr(pname, '-') != NULL) && (residx == 0)) ||
+                    ((strchr(pname, '+') != NULL) && (residx == nres-1)))
+                {
+                    bAddCMAP = FALSE;
+                    break;
+                }
+
+                cmap_atomid[k] = search_atom(pname,
                                              i, atoms, ptr, TRUE);
                 bAddCMAP = bAddCMAP && (cmap_atomid[k] != NO_ATID);
                 if (!bAddCMAP)
@@ -1482,7 +1486,17 @@ static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuet
                     bAddCMAP = bAddCMAP &&
                         cmap_chainnum == resinfo[this_residue_index].chainnum;
                 }
-                bAddCMAP = bAddCMAP && gmx_residuetype_is_protein(rt, *(resinfo[this_residue_index].name));
+                /* Here we used to check that the residuetype was protein and
+                 * disable bAddCMAP if that was not the case. However, some
+                 * special residues (say, alanine dipeptides) might not adhere
+                 * to standard naming, and if we start calling them normal
+                 * protein residues the user will be bugged to select termini.
+                 *
+                 * Instead, I believe that the right course of action is to
+                 * keep the CMAP interaction if it is present in the RTP file
+                 * and we correctly identified all atoms (which is the case
+                 * if we got here).
+                 */
             }
 
             if (bAddCMAP)
@@ -1499,7 +1513,6 @@ static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuet
             }
         }
     }
-
     /* Start the next residue */
 }
 
@@ -1603,7 +1616,7 @@ void pdb2top(FILE *top_file, char *posre_fn, char *molname,
     /* Make CMAP */
     if (TRUE == bCmap)
     {
-        gen_cmap(&(plist[F_CMAP]), restp, atoms, rt);
+        gen_cmap(&(plist[F_CMAP]), restp, atoms);
         if (plist[F_CMAP].nr > 0)
         {
             fprintf(stderr, "There are %4d cmap torsion pairs\n",