Remove unnecessary config.h includes
[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / pdb2top.cpp
index bcbfc5857ffcbf836a9ab72201f788853d0cbf74..ec4f8927ce5d570c14cde224c8b774b9cb3e429c 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.
  * To help us fund GROMACS development, we humbly ask that you cite
  * the research papers on the package. Check out http://www.gromacs.org.
  */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
 
 #include <stdio.h>
 #include <math.h>
 #include <ctype.h>
 
-#include "vec.h"
-#include "smalloc.h"
-#include "macros.h"
-#include "symtab.h"
-#include "gromacs/fileio/futil.h"
-#include "gmx_fatal.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/utility/futil.h"
 #include "pdb2top.h"
 #include "gpp_nextnb.h"
 #include "topdirs.h"
 #include "pgutil.h"
 #include "resall.h"
 #include "topio.h"
-#include "string2.h"
-#include "physics.h"
+#include "gromacs/utility/cstringutil.h"
 #include "gromacs/fileio/pdbio.h"
 #include "gen_ad.h"
 #include "gromacs/fileio/filenm.h"
-#include "index.h"
 #include "gen_vsite.h"
 #include "add_par.h"
 #include "toputil.h"
 #include "fflibutil.h"
-#include "strdb.h"
-#include "copyrite.h"
+#include "gromacs/legacyheaders/copyrite.h"
 
+#include "gromacs/fileio/strdb.h"
+#include "gromacs/topology/residuetypes.h"
+#include "gromacs/topology/symtab.h"
 #include "gromacs/utility/exceptions.h"
-#include "gromacs/utility/programinfo.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/programcontext.h"
+#include "gromacs/utility/smalloc.h"
 
 /* this must correspond to enum in pdb2top.h */
 const char *hh[ehisNR]   = { "HISD", "HISE", "HISH", "HIS1" };
@@ -175,7 +172,7 @@ choose_ff(const char *ffsel,
         ptr = strrchr(ffdirs[i], '/');
         if (ptr == NULL)
         {
-            ffs[i]     = strdup(ffdirs[i]);
+            ffs[i]     = gmx_strdup(ffdirs[i]);
             ffs_dir[i] = low_gmxlibfn(ffdirs[i], FALSE, FALSE);
             if (ffs_dir[i] == NULL)
             {
@@ -184,8 +181,8 @@ choose_ff(const char *ffsel,
         }
         else
         {
-            ffs[i]     = strdup(ptr+1);
-            ffs_dir[i] = strdup(ffdirs[i]);
+            ffs[i]     = gmx_strdup(ptr+1);
+            ffs_dir[i] = gmx_strdup(ffdirs[i]);
         }
         ffs_dir[i][strlen(ffs_dir[i])-strlen(ffs[i])-1] = '\0';
         /* Remove the extension from the ffdir name */
@@ -250,14 +247,14 @@ choose_ff(const char *ffsel,
             if (gmx_fexist(buf))
             {
                 /* We don't use fflib_open, because we don't want printf's */
-                fp = ffopen(buf, "r");
+                fp = gmx_ffopen(buf, "r");
                 snew(desc[i], STRLEN);
                 get_a_line(fp, desc[i], STRLEN);
-                ffclose(fp);
+                gmx_ffclose(fp);
             }
             else
             {
-                desc[i] = strdup(ffs[i]);
+                desc[i] = gmx_strdup(ffs[i]);
             }
         }
         /* Order force fields from the same dir alphabetically
@@ -374,7 +371,7 @@ void choose_watermodel(const char *wmsel, const char *ffdir,
     }
     else if (strcmp(wmsel, "select") != 0)
     {
-        *watermodel = strdup(wmsel);
+        *watermodel = gmx_strdup(wmsel);
 
         return;
     }
@@ -410,7 +407,7 @@ void choose_watermodel(const char *wmsel, const char *ffdir,
             sfree(model[nwm]);
         }
     }
-    ffclose(fp);
+    gmx_ffclose(fp);
     fprintf(stderr, "%2d: %s\n", nwm+1, "None");
 
     sel = -1;
@@ -432,7 +429,7 @@ void choose_watermodel(const char *wmsel, const char *ffdir,
     }
     else
     {
-        *watermodel = strdup(model[sel]);
+        *watermodel = gmx_strdup(model[sel]);
     }
 
     for (i = 0; i < nwm; i++)
@@ -442,8 +439,8 @@ void choose_watermodel(const char *wmsel, const char *ffdir,
     sfree(model);
 }
 
-static int name2type(t_atoms *at, int **cgnr, gpp_atomtype_t atype,
-                     t_restp restp[], gmx_residuetype_t rt)
+static int name2type(t_atoms *at, int **cgnr,
+                     t_restp restp[], gmx_residuetype_t *rt)
 {
     int         i, j, prevresind, resind, i0, prevcg, cg, curcg;
     char       *name;
@@ -492,9 +489,8 @@ static int name2type(t_atoms *at, int **cgnr, gpp_atomtype_t atype,
             j                = search_jtype(&restp[resind], name, bNterm);
             at->atom[i].type = restp[resind].atom[j].type;
             at->atom[i].q    = restp[resind].atom[j].q;
-            at->atom[i].m    = get_atomtype_massA(restp[resind].atom[j].type,
-                                                  atype);
-            cg = restp[resind].cgnr[j];
+            at->atom[i].m    = restp[resind].atom[j].m;
+            cg               = restp[resind].cgnr[j];
             /* A charge group number -1 signals a separate charge group
              * for this atom.
              */
@@ -560,8 +556,7 @@ void print_top_comment(FILE       *out,
         gmx::BinaryInformationSettings settings;
         settings.generatedByHeader(true);
         settings.linePrefix(";\t");
-        gmx::printBinaryInformation(out, gmx::ProgramInfo::getInstance(),
-                                    settings);
+        gmx::printBinaryInformation(out, gmx::getProgramContext(), settings);
     }
     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 
@@ -737,22 +732,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)
@@ -818,6 +798,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));
@@ -1019,7 +1000,7 @@ void add_atom_to_restp(t_restp *restp, int at_start, const t_hack *hack)
         }
         snew( restp->atomname[at_start+1+k], 1);
         restp->atom     [at_start+1+k] = *hack->atom;
-        *restp->atomname[at_start+1+k] = strdup(buf);
+        *restp->atomname[at_start+1+k] = gmx_strdup(buf);
         if (hack->cgnr != NOTSET)
         {
             restp->cgnr   [at_start+1+k] = hack->cgnr;
@@ -1193,7 +1174,7 @@ void get_hackblocks_rtp(t_hackblock **hb, t_restp **restp,
                             }
                             snew( (*restp)[i].atomname[l], 1);
                             (*restp)[i].atom[l]      =       *(*hb)[i].hack[j].atom;
-                            *(*restp)[i].atomname[l] = strdup((*hb)[i].hack[j].nname);
+                            *(*restp)[i].atomname[l] = gmx_strdup((*hb)[i].hack[j].nname);
                             if ( (*hb)[i].hack[j].cgnr != NOTSET)
                             {
                                 (*restp)[i].cgnr   [l] = (*hb)[i].hack[j].cgnr;
@@ -1350,7 +1331,7 @@ static gmx_bool match_atomnames_with_rtp_atom(t_atoms *pdba, rvec *x, int atind,
             }
             /* Rename the atom in pdba */
             snew(pdba->atomname[atind], 1);
-            *pdba->atomname[atind] = strdup(newnm);
+            *pdba->atomname[atind] = gmx_strdup(newnm);
         }
         else if (hbr->hack[j].oname != NULL && hbr->hack[j].nname == NULL &&
                  gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
@@ -1428,10 +1409,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;
@@ -1447,11 +1429,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++)
@@ -1462,7 +1447,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)
@@ -1485,7 +1484,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)
@@ -1502,7 +1511,6 @@ static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuet
             }
         }
     }
-
     /* Start the next residue */
 }
 
@@ -1542,7 +1550,7 @@ void pdb2top(FILE *top_file, char *posre_fn, char *molname,
     int              *vsite_type;
     int               i, nmissat;
     int               bts[ebtsNR];
-    gmx_residuetype_t rt;
+    gmx_residuetype_t*rt;
 
     init_plist(plist);
     gmx_residuetype_init(&rt);
@@ -1563,7 +1571,7 @@ void pdb2top(FILE *top_file, char *posre_fn, char *molname,
                atoms, nssbonds, ssbonds,
                bAllowMissing);
 
-    nmissat = name2type(atoms, &cgnr, atype, restp, rt);
+    nmissat = name2type(atoms, &cgnr, restp, rt);
     if (nmissat)
     {
         if (bAllowMissing)
@@ -1606,7 +1614,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",