Remove unnecessary config.h includes
[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / gen_ad.c
index adee8b11b510468df32b6666aca806654ff4ecc4..c110bc7a3e4a31963416a849307e1a9256152a08 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014, 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.
  * the research papers on the package. Check out http://www.gromacs.org.
  */
 /* This file is completely threadsafe - keep it that way! */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
 
-#include <math.h>
 #include <ctype.h>
-#include "sysstuff.h"
-#include "macros.h"
-#include "smalloc.h"
-#include "string2.h"
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/cstringutil.h"
 #include "gromacs/fileio/confio.h"
-#include "vec.h"
-#include "pbc.h"
+#include "gromacs/math/vec.h"
 #include "toputil.h"
 #include "topio.h"
 #include "gpp_nextnb.h"
-#include "symtab.h"
-#include "macros.h"
-#include "gmx_fatal.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/utility/fatalerror.h"
 #include "pgutil.h"
 #include "resall.h"
 #include "gen_ad.h"
@@ -208,7 +206,7 @@ static void rm2par(t_param p[], int *np, peq eq)
             {
                 fprintf(debug,
                         "Something VERY strange is going on in rm2par (gen_ad.c)\n"
-                        "a[0] %u a[1] %u a[2] %u a[3] %u\n",
+                        "a[0] %d a[1] %d a[2] %d a[3] %d\n",
                         p[i].a[0], p[i].a[1], p[i].a[2], p[i].a[3]);
             }
             strcpy(p[i].s, "");
@@ -785,6 +783,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;
@@ -792,7 +791,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.
@@ -815,6 +813,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
@@ -875,6 +884,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;
                                         }
                                     }
                                 }
@@ -937,6 +948,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.
@@ -1025,6 +1038,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)
     {