Enable rtp angle/dihedral entries not connected by bonds

pdb2gmx has previously silently ignored any angles or
torsions not connected by bonds. This patch adds code
to mark with RTP/hackblock entries that have been
assigned in this search, and afterwards we add those
that have not yet been matched. In particular, this makes
it possible to have angles or torsions in RTP entries
even if those atoms are not connected by bonds.

Fixes #1276.

Change-Id: I61c3c3954ef8c4cd59956d88b73c6a69c1a66c65

diff --git a/src/gromacs/gmxpreprocess/gen_ad.c b/src/gromacs/gmxpreprocess/gen_ad.c
index 7e8bfad88c19e864e1ccfc00113f24c68a6b34fe..4fb527ac0a53a5e32e74ab57906f579fe854a08e 100644 (file)
--- a/src/gromacs/gmxpreprocess/gen_ad.c
+++ b/src/gromacs/gmxpreprocess/gen_ad.c
@@ -787,6 +787,7 @@ void gen_pad(t_nextnb *nnb, t_atoms *atoms, t_restp rtp[],
     t_param    *ang, *dih, *pai, *improper;
     t_rbondeds *hbang, *hbdih;
     char      **anm;
+    const char *p;
     int         res, minres, maxres;
     int         i, j, j1, k, k1, l, l1, m, n, i1, i2;
     int         ninc, maxang, maxdih, maxpai;
@@ -794,7 +795,6 @@ void gen_pad(t_nextnb *nnb, t_atoms *atoms, t_restp rtp[],
     int         nFound;
     gmx_bool    bFound, bExcl;
 
-
     /* These are the angles, dihedrals and pairs that we generate
      * from the bonds. The ones that are already there from the rtp file
      * will be retained.
@@ -817,6 +817,17 @@ void gen_pad(t_nextnb *nnb, t_atoms *atoms, t_restp rtp[],
     if (hb)
     {
         gen_excls(atoms, excls, hb, bAllowMissing);
+        /* mark all entries as not matched yet */
+        for (i = 0; i < atoms->nres; i++)
+        {
+            for (j = 0; j < ebtsNR; j++)
+            {
+                for (k = 0; k < hb[i].rb[j].nb; k++)
+                {
+                    hb[i].rb[j].b[k].match = FALSE;
+                }
+            }
+        }
     }
 
     /* Extract all i-j-k-l neighbours from nnb struct to generate all
@@ -877,6 +888,8 @@ void gen_pad(t_nextnb *nnb, t_atoms *atoms, t_restp rtp[],
                                         if (bFound)
                                         {
                                             set_p_string(&(ang[nang]), hbang->b[l].s);
+                                            /* Mark that we found a match for this entry */
+                                            hbang->b[l].match = TRUE;
                                         }
                                     }
                                 }
@@ -939,6 +952,8 @@ void gen_pad(t_nextnb *nnb, t_atoms *atoms, t_restp rtp[],
                                             if (bFound)
                                             {
                                                 set_p_string(&dih[ndih], hbdih->b[n].s);
+                                                /* Mark that we found a match for this entry */
+                                                hbdih->b[n].match = TRUE;
 
                                                 /* Set the last parameter to be able to see
                                                    if the dihedral was in the rtp list.
@@ -1027,6 +1042,106 @@ void gen_pad(t_nextnb *nnb, t_atoms *atoms, t_restp rtp[],
         }
     }
 
+    /* The above approach is great in that we double-check that e.g. an angle
+     * really corresponds to three atoms connected by bonds, but this is not
+     * generally true. Go through the angle and dihedral hackblocks to add
+     * entries that we have not yet marked as matched when going through bonds.
+     */
+    for (i = 0; i < atoms->nres; i++)
+    {
+        /* Add remaining angles from hackblock */
+        hbang = &hb[i].rb[ebtsANGLES];
+        for (j = 0; j < hbang->nb; j++)
+        {
+            if (hbang->b[j].match == TRUE)
+            {
+                /* We already used this entry, continue to the next */
+                continue;
+            }
+            /* Hm - entry not used, let's see if we can find all atoms */
+            if (nang == maxang)
+            {
+                maxang += ninc;
+                srenew(ang, maxang);
+            }
+            bFound = TRUE;
+            for (k = 0; k < 3 && bFound; k++)
+            {
+                p   = hbang->b[j].a[k];
+                res = i;
+                if (p[0] == '-')
+                {
+                    p++;
+                    res--;
+                }
+                else if (p[0] == '+')
+                {
+                    p++;
+                    res++;
+                }
+                ang[nang].a[k] = search_res_atom(p, res, atoms, "angle", TRUE);
+                bFound         = (ang[nang].a[k] != NO_ATID);
+            }
+            ang[nang].C0 = NOTSET;
+            ang[nang].C1 = NOTSET;
+
+            if (bFound)
+            {
+                set_p_string(&(ang[nang]), hbang->b[j].s);
+                hbang->b[j].match = TRUE;
+                /* Incrementing nang means we save this angle */
+                nang++;
+            }
+        }
+
+        /* Add remaining dihedrals from hackblock */
+        hbdih = &hb[i].rb[ebtsPDIHS];
+        for (j = 0; j < hbdih->nb; j++)
+        {
+            if (hbdih->b[j].match == TRUE)
+            {
+                /* We already used this entry, continue to the next */
+                continue;
+            }
+            /* Hm - entry not used, let's see if we can find all atoms */
+            if (ndih == maxdih)
+            {
+                maxdih += ninc;
+                srenew(dih, maxdih);
+            }
+            bFound = TRUE;
+            for (k = 0; k < 4 && bFound; k++)
+            {
+                p   = hbdih->b[j].a[k];
+                res = i;
+                if (p[0] == '-')
+                {
+                    p++;
+                    res--;
+                }
+                else if (p[0] == '+')
+                {
+                    p++;
+                    res++;
+                }
+                dih[ndih].a[k] = search_res_atom(p, res, atoms, "dihedral", TRUE);
+                bFound         = (dih[ndih].a[k] != NO_ATID);
+            }
+            for (m = 0; m < MAXFORCEPARAM; m++)
+            {
+                dih[ndih].c[m] = NOTSET;
+            }
+
+            if (bFound)
+            {
+                set_p_string(&(dih[ndih]), hbdih->b[j].s);
+                hbdih->b[j].match = TRUE;
+                /* Incrementing ndih means we save this dihedral */
+                ndih++;
+            }
+        }
+    }
+
     /* Sort angles with respect to j-i-k (middle atom first) */
     if (nang > 1)
     {

diff --git a/src/gromacs/gmxpreprocess/hackblock.c b/src/gromacs/gmxpreprocess/hackblock.c
index a5cdc7a0afd71dfef7304f7c44add61506eb9635..0e4e1c1cfa988a8723208fb7e5620eda7d03d1f0 100644 (file)
--- a/src/gromacs/gmxpreprocess/hackblock.c
+++ b/src/gromacs/gmxpreprocess/hackblock.c
@@ -44,6 +44,7 @@
 #include "gromacs/utility/smalloc.h"
 #include "vec.h"
 #include "macros.h"
+#include "names.h"
 
 /* these MUST correspond to the enum in hackblock.h */
 const char *btsNames[ebtsNR] = { "bonds", "angles", "dihedrals", "impropers", "exclusions", "cmap" };
@@ -176,7 +177,8 @@ static void copy_t_rbonded(t_rbonded *s, t_rbonded *d)
     {
         d->a[i] = safe_strdup(s->a[i]);
     }
-    d->s = safe_strdup(s->s);
+    d->s     = safe_strdup(s->s);
+    d->match = s->match;
 }
 
 static gmx_bool contains_char(t_rbonded *s, char c)
@@ -395,7 +397,7 @@ void dump_hb(FILE *out, int nres, t_hackblock hb[])
                     }
                     fprintf(out, " %s]", SS(hb[i].rb[j].b[k].s));
                 }
-                fprintf(out, "\n");
+                fprintf(out, " Entry matched: %s\n", yesno_names[hb[i].rb[j].b[k].match]);
             }
         }
         fprintf(out, "\n");

diff --git a/src/gromacs/gmxpreprocess/hackblock.h b/src/gromacs/gmxpreprocess/hackblock.h
index 8d9fee3c1705e5932b2c97733b19f245651c29f1..4bbbb7d9e5e8690efbe07c9cad3d2147c71c9357 100644 (file)
--- a/src/gromacs/gmxpreprocess/hackblock.h
+++ b/src/gromacs/gmxpreprocess/hackblock.h
@@ -65,6 +65,7 @@ typedef struct {
     char  *s;              /* optional define string which gets copied from
                               .rtp/.tdb to .top and will be parsed by cpp
                               during grompp */
+    gmx_bool match;        /* boolean to mark that the entry has been found */
 } t_rbonded;
 
 typedef struct {

diff --git a/src/gromacs/gmxpreprocess/pdb2top.cpp b/src/gromacs/gmxpreprocess/pdb2top.cpp
index 9729f84d51b3bb4a11a4a1ecdad76c29ac8de382..85eb9c4c9b569a1224681eba5444b0b16518ed16 100644 (file)
--- a/src/gromacs/gmxpreprocess/pdb2top.cpp
+++ b/src/gromacs/gmxpreprocess/pdb2top.cpp
@@ -735,22 +735,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)
@@ -816,6 +801,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));

diff --git a/src/gromacs/gmxpreprocess/pgutil.c b/src/gromacs/gmxpreprocess/pgutil.c
index 80088875f820b240cfd02beeeeb03f4e497f8959..76ea6c83b40e11ab2b066acbe068565f30525113 100644 (file)
--- a/src/gromacs/gmxpreprocess/pgutil.c
+++ b/src/gromacs/gmxpreprocess/pgutil.c
@@ -148,6 +148,25 @@ atom_id search_atom(const char *type, int start,
     return NO_ATID;
 }
 
+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;
+}
+
+
 void set_at(t_atom *at, real m, real q, int type, int resind)
 {
     at->m      = m;

diff --git a/src/gromacs/gmxpreprocess/pgutil.h b/src/gromacs/gmxpreprocess/pgutil.h
index 3631ade8cb808721ca9db1e2c767a96c54aaaac8..8a5a622fc4a09b6366e00317f0144082c0b7ba38 100644 (file)
--- a/src/gromacs/gmxpreprocess/pgutil.h
+++ b/src/gromacs/gmxpreprocess/pgutil.h
@@ -45,18 +45,28 @@ extern "C"
 {
 #endif
 
-atom_id search_atom(const char *type, int start,
-                    t_atoms *atoms,
-                    const char *bondtype, gmx_bool bAllowMissing);
 /* Search an atom in array of pointers to strings, starting from start
  * if type starts with '-' then searches backwards from start.
  * bondtype is only used for printing the error/warning string,
  * when bondtype="check" no error/warning is issued.
  * When bAllowMissing=FALSE an fatal error is issued, otherwise a warning.
  */
+atom_id search_atom(const char *type, int start,
+                    t_atoms *atoms,
+                    const char *bondtype, gmx_bool bAllowMissing);
+
+/* Similar to search_atom, but this routine searches for the specified
+ * atom in residue resind.
+ */
+atom_id
+search_res_atom(const char *type, int resind,
+                t_atoms *atoms,
+                const char *bondtype, gmx_bool bAllowMissing);
+
 
 void set_at(t_atom *at, real m, real q, int type, int resind);
 
+
 #ifdef __cplusplus
 }
 #endif