Merge branch release-5-1

[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / toppush.cpp

diff --git a/src/gromacs/gmxpreprocess/toppush.cpp b/src/gromacs/gmxpreprocess/toppush.cpp
index 9bc3a53aecc10af74ad4615915f824cf682d81dc..3763c2e55a7319e94931b14046c164f97a992a19 100644 (file)
--- a/src/gromacs/gmxpreprocess/toppush.cpp
+++ b/src/gromacs/gmxpreprocess/toppush.cpp
@@ -1576,13 +1576,37 @@ static gmx_bool default_cmap_params(t_params bondtype[],
     return bFound;
 }
 
+/* Returns the number of exact atom type matches, i.e. non wild-card matches,
+ * returns -1 when there are no matches at all.
+ */
+static int natom_match(t_param *pi,
+                       int type_i, int type_j, int type_k, int type_l,
+                       const gpp_atomtype_t atype)
+{
+    if ((pi->ai() == -1 || get_atomtype_batype(type_i, atype) == pi->ai()) &&
+        (pi->aj() == -1 || get_atomtype_batype(type_j, atype) == pi->aj()) &&
+        (pi->ak() == -1 || get_atomtype_batype(type_k, atype) == pi->ak()) &&
+        (pi->al() == -1 || get_atomtype_batype(type_l, atype) == pi->al()))
+    {
+        return
+            (pi->ai() == -1 ? 0 : 1) +
+            (pi->aj() == -1 ? 0 : 1) +
+            (pi->ak() == -1 ? 0 : 1) +
+            (pi->al() == -1 ? 0 : 1);
+    }
+    else
+    {
+        return -1;
+    }
+}
+
 static gmx_bool default_params(int ftype, t_params bt[],
                                t_atoms *at, gpp_atomtype_t atype,
                                t_param *p, gmx_bool bB,
                                t_param **param_def,
                                int *nparam_def)
 {
-    int          i, j, nparam_found;
+    int          nparam_found;
     gmx_bool     bFound, bSame;
     t_param     *pi    = NULL;
     t_param     *pj    = NULL;
@@ -1597,55 +1621,54 @@ static gmx_bool default_params(int ftype, t_params bt[],
     }
 
 
-    /* We allow wildcards now. The first type (with or without wildcards) that
-     * fits is used, so you should probably put the wildcarded bondtypes
-     * at the end of each section.
-     */
     bFound       = FALSE;
     nparam_found = 0;
-    /* OPLS uses 1000s of dihedraltypes, so in order to speed up the scanning we have a
-     * special case for this. Check for B state outside loop to speed it up.
-     */
     if (ftype == F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS || ftype == F_PIDIHS)
     {
-        if (bB)
+        int nmatch_max = -1;
+        int i          = -1;
+        int t;
+
+        /* For dihedrals we allow wildcards. We choose the first type
+         * that has the most real matches, i.e. non-wildcard matches.
+         */
+        for (t = 0; ((t < nr) && nmatch_max < 4); t++)
         {
-            for (i = 0; ((i < nr) && !bFound); i++)
+            int      nmatch;
+            t_param *pt;
+
+            pt = &(bt[ftype].param[t]);
+            if (bB)
             {
-                pi     = &(bt[ftype].param[i]);
-                bFound =
-                    (
-                        ((pi->ai() == -1) || (get_atomtype_batype(at->atom[p->ai()].typeB, atype) == pi->ai())) &&
-                        ((pi->aj() == -1) || (get_atomtype_batype(at->atom[p->aj()].typeB, atype) == pi->aj())) &&
-                        ((pi->ak() == -1) || (get_atomtype_batype(at->atom[p->ak()].typeB, atype) == pi->ak())) &&
-                        ((pi->al() == -1) || (get_atomtype_batype(at->atom[p->al()].typeB, atype) == pi->al()))
-                    );
+                nmatch = natom_match(pt, at->atom[p->ai()].typeB, at->atom[p->aj()].typeB, at->atom[p->ak()].typeB, at->atom[p->al()].typeB, atype);
             }
-        }
-        else
-        {
-            /* State A */
-            for (i = 0; ((i < nr) && !bFound); i++)
+            else
+            {
+                nmatch = natom_match(pt, at->atom[p->ai()].type, at->atom[p->aj()].type, at->atom[p->ak()].type, at->atom[p->al()].type, atype);
+            }
+            if (nmatch > nmatch_max)
             {
-                pi     = &(bt[ftype].param[i]);
-                bFound =
-                    (
-                        ((pi->ai() == -1) || (get_atomtype_batype(at->atom[p->ai()].type, atype) == pi->ai())) &&
-                        ((pi->aj() == -1) || (get_atomtype_batype(at->atom[p->aj()].type, atype) == pi->aj())) &&
-                        ((pi->ak() == -1) || (get_atomtype_batype(at->atom[p->ak()].type, atype) == pi->ak())) &&
-                        ((pi->al() == -1) || (get_atomtype_batype(at->atom[p->al()].type, atype) == pi->al()))
-                    );
+                nmatch_max = nmatch;
+                i          = t;
+                bFound     = TRUE;
             }
         }
-        /* Find additional matches for this dihedral - necessary for ftype==9 which is used e.g. for charmm.
-         * The rules in that case is that additional matches HAVE to be on adjacent lines!
-         */
+
         if (bFound == TRUE)
         {
+            int j;
+
+            pi    = &(bt[ftype].param[i]);
             nparam_found++;
+
+            /* Find additional matches for this dihedral - necessary
+             * for ftype==9.
+             * The rule in that case is that additional matches
+             * HAVE to be on adjacent lines!
+             */
             bSame = TRUE;
             /* Continue from current i value */
-            for (j = i+1; j < nr && bSame; j += 2)
+            for (j = i + 2; j < nr && bSame; j += 2)
             {
                 pj    = &(bt[ftype].param[j]);
                 bSame = (pi->ai() == pj->ai() && pi->aj() == pj->aj() && pi->ak() == pj->ak() && pi->al() == pj->al());
@@ -1659,6 +1682,8 @@ static gmx_bool default_params(int ftype, t_params bt[],
     }
     else   /* Not a dihedral */
     {
+        int i, j;
+
         for (i = 0; ((i < nr) && !bFound); i++)
         {
             pi = &(bt[ftype].param[i]);