Merge branch release-5-0
[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / hackblock.c
index c30a69851b82c6b85476b52c3bab218c1dfd0571..98148c5a1be280b2daa8ddbad6d6e96ae7d4fbd4 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) 2011,2014, by the GROMACS development team, led by
+ * Copyright (c) 2011,2014,2015, 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.
@@ -44,6 +44,7 @@
 #include "gromacs/legacyheaders/names.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/smalloc.h"
 
 /* these MUST correspond to the enum in hackblock.h */
@@ -198,13 +199,14 @@ static gmx_bool contains_char(t_rbonded *s, char c)
     return bRet;
 }
 
-gmx_bool rbonded_atoms_exist_in_list(t_rbonded *b, t_rbonded blist[], int nlist, int natoms)
+int
+rbonded_find_atoms_in_list(t_rbonded *b, t_rbonded blist[], int nlist, int natoms)
 {
     int      i, k;
-    gmx_bool matchFound = FALSE;
+    int      foundPos = -1;
     gmx_bool atomsMatch;
 
-    for (i = 0; i < nlist && !matchFound; i++)
+    for (i = 0; i < nlist && foundPos < 0; i++)
     {
         atomsMatch = TRUE;
         for (k = 0; k < natoms && atomsMatch; k++)
@@ -220,30 +222,79 @@ gmx_bool rbonded_atoms_exist_in_list(t_rbonded *b, t_rbonded blist[], int nlist,
                 atomsMatch = atomsMatch && !strcmp(b->a[k], blist[i].a[natoms-1-k]);
             }
         }
-        matchFound = atomsMatch;
+        if (atomsMatch)
+        {
+            foundPos = i;
+            /* If all the atoms AND all the parameters match, it is likely that
+             * the user made a copy-and-paste mistake (since it would be much cheaper
+             * to just bump the force constant 2x if you really want it twice).
+             * Since we only have the unparsed string here we can only detect
+             * EXACT matches (including identical whitespace).
+             */
+            if (!strcmp(b->s, blist[i].s))
+            {
+                gmx_warning("Duplicate line found in or between hackblock and rtp entries");
+            }
+        }
     }
-    return matchFound;
+    return foundPos;
 }
 
 gmx_bool merge_t_bondeds(t_rbondeds s[], t_rbondeds d[], gmx_bool bMin, gmx_bool bPlus)
 {
     int      i, j;
     gmx_bool bBondsRemoved;
+    int      nbHackblockStart;
+    int      index;
 
     bBondsRemoved = FALSE;
     for (i = 0; i < ebtsNR; i++)
     {
         if (s[i].nb > 0)
         {
+            /* Record how many bonds we have in the destination when we start.
+             *
+             * If an entry is present in the hackblock (destination), we will
+             * not add the one from the main rtp, since the point is for hackblocks
+             * to overwrite it. However, if there is no hackblock entry we do
+             * allow multiple main rtp entries since some forcefield insist on that.
+             *
+             * We accomplish this by checking the position we find an entry in,
+             * rather than merely checking whether it exists at all.
+             * If that index is larger than the original (hackblock) destination
+             * size, it was added from the main rtp, and then we will allow more
+             * such entries. In contrast, if the entry found has a lower index
+             * it is a hackblock entry meant to override the main rtp, and then
+             * we don't add the main rtp one.
+             */
+            nbHackblockStart = d[i].nb;
+
             /* make space */
             srenew(d[i].b, d[i].nb + s[i].nb);
             for (j = 0; j < s[i].nb; j++)
             {
                 /* Check if this bonded string already exists before adding.
-                 * We are merging from the main rtp to the hackblocks, so this
-                 * will mean the hackblocks overwrite the man rtp, as intended.
+                 * We are merging from the main RTP to the hackblocks, so this
+                 * will mean the hackblocks overwrite the man RTP, as intended.
+                 */
+                index = rbonded_find_atoms_in_list(&s[i].b[j], d[i].b, d[i].nb, btsNiatoms[i]);
+                /* - If we did not find this interaction at all, the index will be -1,
+                 *   and then we should definitely add it to the merged hackblock and rtp.
+                 *
+                 * Alternatively, if it was found, index will be >=0.
+                 * - In case this index is lower than the original number of entries,
+                 *   it is already present as a *hackblock* entry, and those should
+                 *   always override whatever we have listed in the RTP. Thus, we
+                 *   should just keep that one and not add anything from the RTP.
+                 * - Finally, if it was found, but with an index higher than
+                 *   the original number of entries, it comes from the RTP rather
+                 *   than hackblock, and then we must have added it ourselves
+                 *   in a previous iteration. In that case it is a matter of
+                 *   several entries for the same sequence of atoms, and we allow
+                 *   that in the RTP. In this case we should simply copy all of
+                 *   them, including this one.
                  */
-                if (!rbonded_atoms_exist_in_list(&s[i].b[j], d[i].b, d[i].nb, btsNiatoms[i]))
+                if (index < 0 || index >= nbHackblockStart)
                 {
                     if (!(bMin && contains_char(&s[i].b[j], '-'))
                         && !(bPlus && contains_char(&s[i].b[j], '+')))
@@ -256,6 +307,12 @@ gmx_bool merge_t_bondeds(t_rbondeds s[], t_rbondeds d[], gmx_bool bMin, gmx_bool
                         bBondsRemoved = TRUE;
                     }
                 }
+                else
+                {
+                    /* This is the common case where a hackblock entry simply
+                     * overrides the RTP, so we cannot warn here.
+                     */
+                }
             }
         }
     }