Sort all includes in src/gromacs
[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / toppush.c
index c823fb38674cdf7dbefa5bde32a5b31d241cd6b0..161940b72e0f2a73267baeab449016d12247b67d 100644 (file)
@@ -1,59 +1,60 @@
-/*  -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+/*
+ * This file is part of the GROMACS molecular simulation package.
  *
- *
- *                This source code is part of
- *
- *                 G   R   O   M   A   C   S
- *
- *          GROningen MAchine for Chemical Simulations
- *
- *                        VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
  * of the License, or (at your option) any later version.
  *
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
  *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  *
- * For more info, check our website at http://www.gromacs.org
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
  *
- * And Hey:
- * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
+ * 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 "toppush.h"
+
+#include <assert.h>
 #include <ctype.h>
 #include <math.h>
-
-#include "sysstuff.h"
-#include "smalloc.h"
-#include "macros.h"
-#include "string2.h"
-#include "names.h"
-#include "toputil.h"
-#include "toppush.h"
-#include "topdirs.h"
-#include "readir.h"
-#include "symtab.h"
-#include "gmx_fatal.h"
-#include "warninp.h"
-#include "gpp_atomtype.h"
-#include "gpp_bond_atomtype.h"
+#include <stdlib.h>
+
+#include "gromacs/gmxpreprocess/gpp_atomtype.h"
+#include "gromacs/gmxpreprocess/gpp_bond_atomtype.h"
+#include "gromacs/gmxpreprocess/readir.h"
+#include "gromacs/gmxpreprocess/topdirs.h"
+#include "gromacs/gmxpreprocess/toputil.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/warninp.h"
+#include "gromacs/topology/symtab.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
 
 void generate_nbparams(int comb, int ftype, t_params *plist, gpp_atomtype_t atype,
                        warninp_t wi)
@@ -93,7 +94,7 @@ void generate_nbparams(int comb, int ftype, t_params *plist, gpp_atomtype_t atyp
                     break;
 
                 case eCOMB_ARITHMETIC:
-                    /* c0 and c1 are epsilon and sigma */
+                    /* c0 and c1 are sigma and epsilon */
                     for (i = k = 0; (i < nr); i++)
                     {
                         for (j = 0; (j < nr); j++, k++)
@@ -102,14 +103,21 @@ void generate_nbparams(int comb, int ftype, t_params *plist, gpp_atomtype_t atyp
                             cj0                  = get_atomtype_nbparam(j, 0, atype);
                             ci1                  = get_atomtype_nbparam(i, 1, atype);
                             cj1                  = get_atomtype_nbparam(j, 1, atype);
-                            plist->param[k].c[0] = (ci0+cj0)*0.5;
+                            plist->param[k].c[0] = (fabs(ci0) + fabs(cj0))*0.5;
+                            /* Negative sigma signals that c6 should be set to zero later,
+                             * so we need to propagate that through the combination rules.
+                             */
+                            if (ci0 < 0 || cj0 < 0)
+                            {
+                                plist->param[k].c[0] *= -1;
+                            }
                             plist->param[k].c[1] = sqrt(ci1*cj1);
                         }
                     }
 
                     break;
                 case eCOMB_GEOM_SIG_EPS:
-                    /* c0 and c1 are epsilon and sigma */
+                    /* c0 and c1 are sigma and epsilon */
                     for (i = k = 0; (i < nr); i++)
                     {
                         for (j = 0; (j < nr); j++, k++)
@@ -118,7 +126,14 @@ void generate_nbparams(int comb, int ftype, t_params *plist, gpp_atomtype_t atyp
                             cj0                  = get_atomtype_nbparam(j, 0, atype);
                             ci1                  = get_atomtype_nbparam(i, 1, atype);
                             cj1                  = get_atomtype_nbparam(j, 1, atype);
-                            plist->param[k].c[0] = sqrt(ci0*cj0);
+                            plist->param[k].c[0] = sqrt(fabs(ci0*cj0));
+                            /* Negative sigma signals that c6 should be set to zero later,
+                             * so we need to propagate that through the combination rules.
+                             */
+                            if (ci0 < 0 || cj0 < 0)
+                            {
+                                plist->param[k].c[0] *= -1;
+                            }
                             plist->param[k].c[1] = sqrt(ci1*cj1);
                         }
                     }
@@ -950,6 +965,7 @@ void push_nbt(directive d, t_nbparam **nbt, gpp_atomtype_t atype,
         /* When the B topology parameters are not set,
          * copy them from topology A
          */
+        assert(nrfp <= 4);
         for (i = n; i < nrfp; i++)
         {
             c[i] = c[i-2];
@@ -1110,8 +1126,8 @@ push_cmaptype(directive d, t_params bt[], int nral, gpp_atomtype_t at,
     nrfp   = nrfpA+nrfpB;
 
     /* Allocate memory for the CMAP grid */
-    bt->ncmap += nrfp;
-    srenew(bt->cmap, bt->ncmap);
+    bt[F_CMAP].ncmap += nrfp;
+    srenew(bt[F_CMAP].cmap, bt[F_CMAP].ncmap);
 
     /* Read in CMAP parameters */
     sl = 0;
@@ -1123,7 +1139,7 @@ push_cmaptype(directive d, t_params bt[], int nral, gpp_atomtype_t at,
         }
         nn  = sscanf(line+start+sl, " %s ", s);
         sl += strlen(s);
-        bt->cmap[i+(bt->ncmap)-nrfp] = strtod(s, NULL);
+        bt[F_CMAP].cmap[i+(bt[F_CMAP].ncmap)-nrfp] = strtod(s, NULL);
 
         if (nn == 1)
         {
@@ -1141,7 +1157,7 @@ push_cmaptype(directive d, t_params bt[], int nral, gpp_atomtype_t at,
     {
         for (i = 0; i < ncmap; i++)
         {
-            bt->cmap[i+ncmap] = bt->cmap[i];
+            bt[F_CMAP].cmap[i+ncmap] = bt[F_CMAP].cmap[i];
         }
     }
     else
@@ -1164,12 +1180,12 @@ push_cmaptype(directive d, t_params bt[], int nral, gpp_atomtype_t at,
     /* Set grid spacing and the number of grids (we assume these numbers to be the same for all grids
      * so we can safely assign them each time
      */
-    bt->grid_spacing = nxcmap;     /* Or nycmap, they need to be equal */
-    bt->nc           = bt->nc + 1; /* Since we are incrementing here, we need to subtract later, see (*****) */
-    nct              = (nral+1) * bt->nc;
+    bt[F_CMAP].grid_spacing = nxcmap;            /* Or nycmap, they need to be equal */
+    bt[F_CMAP].nc           = bt[F_CMAP].nc + 1; /* Since we are incrementing here, we need to subtract later, see (*****) */
+    nct                     = (nral+1) * bt[F_CMAP].nc;
 
     /* Allocate memory for the cmap_types information */
-    srenew(bt->cmap_types, nct);
+    srenew(bt[F_CMAP].cmap_types, nct);
 
     for (i = 0; (i < nral); i++)
     {
@@ -1183,16 +1199,16 @@ push_cmaptype(directive d, t_params bt[], int nral, gpp_atomtype_t at,
         }
 
         /* Assign a grid number to each cmap_type */
-        bt->cmap_types[bt->nct++] = get_bond_atomtype_type(alc[i], bat);
+        bt[F_CMAP].cmap_types[bt[F_CMAP].nct++] = get_bond_atomtype_type(alc[i], bat);
     }
 
     /* Assign a type number to this cmap */
-    bt->cmap_types[bt->nct++] = bt->nc-1; /* Since we inremented earlier, we need to subtrac here, to get the types right (****) */
+    bt[F_CMAP].cmap_types[bt[F_CMAP].nct++] = bt[F_CMAP].nc-1; /* Since we inremented earlier, we need to subtrac here, to get the types right (****) */
 
     /* Check for the correct number of atoms (again) */
-    if (bt->nct != nct)
+    if (bt[F_CMAP].nct != nct)
     {
-        gmx_fatal(FARGS, "Incorrect number of atom types (%d) in cmap type %d\n", nct, bt->nc);
+        gmx_fatal(FARGS, "Incorrect number of atom types (%d) in cmap type %d\n", nct, bt[F_CMAP].nc);
     }
 
     /* Is this correct?? */
@@ -1530,7 +1546,7 @@ static gmx_bool default_cmap_params(t_params bondtype[],
     ct           = 0;
 
     /* Match the current cmap angle against the list of cmap_types */
-    for (i = 0; i < bondtype->nct && !bFound; i += 6)
+    for (i = 0; i < bondtype[F_CMAP].nct && !bFound; i += 6)
     {
         if (bB)
         {
@@ -1539,15 +1555,15 @@ static gmx_bool default_cmap_params(t_params bondtype[],
         else
         {
             if (
-                (get_atomtype_batype(at->atom[p->a[0]].type, atype) == bondtype->cmap_types[i])   &&
-                (get_atomtype_batype(at->atom[p->a[1]].type, atype) == bondtype->cmap_types[i+1]) &&
-                (get_atomtype_batype(at->atom[p->a[2]].type, atype) == bondtype->cmap_types[i+2]) &&
-                (get_atomtype_batype(at->atom[p->a[3]].type, atype) == bondtype->cmap_types[i+3]) &&
-                (get_atomtype_batype(at->atom[p->a[4]].type, atype) == bondtype->cmap_types[i+4]))
+                (get_atomtype_batype(at->atom[p->a[0]].type, atype) == bondtype[F_CMAP].cmap_types[i])   &&
+                (get_atomtype_batype(at->atom[p->a[1]].type, atype) == bondtype[F_CMAP].cmap_types[i+1]) &&
+                (get_atomtype_batype(at->atom[p->a[2]].type, atype) == bondtype[F_CMAP].cmap_types[i+2]) &&
+                (get_atomtype_batype(at->atom[p->a[3]].type, atype) == bondtype[F_CMAP].cmap_types[i+3]) &&
+                (get_atomtype_batype(at->atom[p->a[4]].type, atype) == bondtype[F_CMAP].cmap_types[i+4]))
             {
                 /* Found cmap torsion */
                 bFound       = TRUE;
-                ct           = bondtype->cmap_types[i+5];
+                ct           = bondtype[F_CMAP].cmap_types[i+5];
                 nparam_found = 1;
             }
         }
@@ -1821,7 +1837,8 @@ void push_bond(directive d, t_params bondtype[], t_params bond[],
         bFoundA = default_params(ftype, bondtype, at, atype, &param, FALSE, &param_defA, &nparam_defA);
         if (bFoundA)
         {
-            /* Copy the A-state and B-state default parameters */
+            /* Copy the A-state and B-state default parameters. */
+            assert(NRFPA(ftype)+NRFPB(ftype) <= MAXFORCEPARAM);
             for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++)
             {
                 param.c[j] = param_defA->c[j];
@@ -2513,14 +2530,14 @@ int add_atomtype_decoupled(t_symtab *symtab, gpp_atomtype_t at,
 static void convert_pairs_to_pairsQ(t_params *plist,
                                     real fudgeQQ, t_atoms *atoms)
 {
-    t_param *paramp1,*paramp2,*paramnew;
-    int      i,j,p1nr,p2nr,p2newnr;
+    t_param *paramp1, *paramp2, *paramnew;
+    int      i, j, p1nr, p2nr, p2newnr;
 
     /* Add the pair list to the pairQ list */
-    p1nr = plist[F_LJ14].nr;
-    p2nr = plist[F_LJC14_Q].nr;
+    p1nr    = plist[F_LJ14].nr;
+    p2nr    = plist[F_LJC14_Q].nr;
     p2newnr = p1nr + p2nr;
-    snew(paramnew,p2newnr);
+    snew(paramnew, p2newnr);
 
     paramp1             = plist[F_LJ14].param;
     paramp2             = plist[F_LJC14_Q].param;
@@ -2529,18 +2546,18 @@ static void convert_pairs_to_pairsQ(t_params *plist,
        it may be possible to just ADD the converted F_LJ14 array
        to the old F_LJC14_Q array, but since we have to create
        a new sized memory structure, better just to deep copy it all.
-    */
+     */
 
     for (i = 0; i < p2nr; i++)
     {
         /* Copy over parameters */
-        for (j=0;j<5;j++) /* entries are 0-4 for F_LJC14_Q */
+        for (j = 0; j < 5; j++) /* entries are 0-4 for F_LJC14_Q */
         {
             paramnew[i].c[j] = paramp2[i].c[j];
         }
 
         /* copy over atoms */
-        for (j=0;j<2;j++)
+        for (j = 0; j < 2; j++)
         {
             paramnew[i].a[j] = paramp2[i].a[j];
         }