Merge release-5-0 into master
[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / toppush.c
index 18787878f9629204a9587f036ac1a145d713862c..8272b7787ecf3f2dad111040e9d2ca76dca6b884 100644 (file)
  * 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 <stdlib.h>
 
-#include "macros.h"
-#include "names.h"
-#include "toputil.h"
-#include "toppush.h"
-#include "topdirs.h"
-#include "readir.h"
-#include "symtab.h"
-#include "warninp.h"
-#include "gpp_atomtype.h"
-#include "gpp_bond_atomtype.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"
@@ -96,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++)
@@ -105,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++)
@@ -121,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);
                         }
                     }
@@ -587,9 +599,11 @@ static void push_bondtype(t_params     *       bt,
                     {
                         sprintf(errbuf, "Overriding %s parameters.%s",
                                 interaction_function[ftype].longname,
-                                (ftype == F_PDIHS) ? "\nUse dihedraltype 4 to allow several multiplicity terms." : "");
+                                (ftype == F_PDIHS) ?
+                                "\nUse dihedraltype 9 to allow several multiplicity terms. Only consecutive lines are combined. Non-consective lines overwrite each other."
+                                : "");
                         warning(wi, errbuf);
-                        fprintf(stderr, "  old:");
+                        fprintf(stderr, "  old:                                         ");
                         for (j = 0; (j < nrfp); j++)
                         {
                             fprintf(stderr, " %g", bt->param[i].c[j]);
@@ -1114,8 +1128,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;
@@ -1127,7 +1141,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)
         {
@@ -1145,7 +1159,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
@@ -1168,12 +1182,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++)
     {
@@ -1187,16 +1201,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?? */
@@ -1534,7 +1548,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)
         {
@@ -1543,15 +1557,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;
             }
         }