From 3c6ba9c734d8b883eca6bc7919cdd97307798985 Mon Sep 17 00:00:00 2001
From: Berk Hess <hess@kth.se>
Date: Sat, 28 May 2016 09:22:37 +0200
Subject: [PATCH] Add grompp check for unbound atoms

grompp now print a note for atoms that are not connected by
a potential or constraint to any other atom in the same moleculetype,
since this often means the user made a mistake.

Refs #1958.

Change-Id: Iabb00563c76a9f7954f84d89d1c67d438f2c31ff
---
 src/gromacs/gmxpreprocess/grompp.cpp | 69 ++++++++++++++++++++++++++++++++++++
 1 file changed, 69 insertions(+)

diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp
index fdfa16c950..1da65c6e36 100644
--- a/src/gromacs/gmxpreprocess/grompp.cpp
+++ b/src/gromacs/gmxpreprocess/grompp.cpp
@@ -1354,6 +1354,73 @@ static real get_max_reference_temp(const t_inputrec *ir,
     return ref_t;
 }
 
+/* Checks if there are unbound atoms in moleculetype molt.
+ * Prints a note for each unbound atoms and a warning if any is present.
+ */
+static void checkForUnboundAtoms(const gmx_moltype_t *molt,
+                                 warninp_t            wi)
+{
+    const t_atoms *atoms = &molt->atoms;
+
+    if (atoms->nr == 1)
+    {
+        /* Only one atom, there can't be unbound atoms */
+        return;
+    }
+
+    std::vector<int> count(atoms->nr, 0);
+
+    for (int ftype = 0; ftype < F_NRE; ftype++)
+    {
+        if (((interaction_function[ftype].flags & IF_BOND) && ftype != F_CONNBONDS) ||
+            (interaction_function[ftype].flags & IF_CONSTRAINT) ||
+            ftype == F_SETTLE)
+        {
+            const t_ilist *il   = &molt->ilist[ftype];
+            int            nral = NRAL(ftype);
+
+            for (int i = 0; i < il->nr; i += 1 + nral)
+            {
+                for (int j = 0; j < nral; j++)
+                {
+                    count[il->iatoms[i + 1 + j]]++;
+                }
+            }
+        }
+    }
+
+    int numDanglingAtoms = 0;
+    for (int a = 0; a < atoms->nr; a++)
+    {
+        if (atoms->atom[a].ptype != eptVSite &&
+            count[a] == 0)
+        {
+            fprintf(stderr, "\nAtom %d '%s' in moleculetype '%s' is not bound by a potential or constraint to any other atom in the same moleculetype.\n",
+                    a + 1, *atoms->atomname[a], *molt->name);
+
+            numDanglingAtoms++;
+        }
+    }
+
+    if (numDanglingAtoms > 0)
+    {
+        char buf[STRLEN];
+        sprintf(buf, "In moleculetype '%s' %d atoms are not bound by a potential or constraint to any other atom in the same moleculetype. Although technically this might not cause issues in a simulation, this often means that the user forgot to add a bond/potential/constraint or put multiple molecules in the same moleculetype definition by mistake.",
+                *molt->name, numDanglingAtoms);
+        warning_note(wi, buf);
+    }
+}
+
+/* Checks all moleculetypes for unbound atoms */
+static void checkForUnboundAtoms(const gmx_mtop_t *mtop,
+                                 warninp_t         wi)
+{
+    for (int mt = 0; mt < mtop->nmoltype; mt++)
+    {
+        checkForUnboundAtoms(&mtop->moltype[mt], wi);
+    }
+}
+
 static void set_verlet_buffer(const gmx_mtop_t *mtop,
                               t_inputrec       *ir,
                               real              buffer_temp,
@@ -1810,6 +1877,8 @@ int gmx_grompp(int argc, char *argv[])
     /* check masses */
     check_mol(sys, wi);
 
+    checkForUnboundAtoms(sys, wi);
+
     for (i = 0; i < sys->nmoltype; i++)
     {
         check_cg_sizes(ftp2fn(efTOP, NFILE, fnm), &sys->moltype[i].cgs, wi);
-- 
2.11.4.GIT