Move QMMM charge setting to grompp

Also made sure that the const gmx_mtop_t does not get changed during
QMMM initialization.

Change-Id: I9bec4bbdb22ea6043019e887f8dbe44e7512ced5

diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp
index 94807e64580013046e663633ef77d703e8cea3fa..d95b9fdf32289c5a75bce580848ffda8196a581f 100644 (file)
--- a/src/gromacs/gmxpreprocess/grompp.cpp
+++ b/src/gromacs/gmxpreprocess/grompp.cpp
@@ -79,6 +79,7 @@
 #include "gromacs/mdlib/compute_io.h"
 #include "gromacs/mdlib/constr.h"
 #include "gromacs/mdlib/perf_est.h"
+#include "gromacs/mdlib/qmmm.h"
 #include "gromacs/mdlib/sim_util.h"
 #include "gromacs/mdrunutility/mdmodules.h"
 #include "gromacs/mdtypes/inputrec.h"
@@ -2132,7 +2133,7 @@ int gmx_grompp(int argc, char *argv[])
         pr_symtab(debug, 0, "After close", &sys.symtab);
     }
 
-    /* make exclusions between QM atoms */
+    /* make exclusions between QM atoms and remove charges if needed */
     if (ir->bQMMM)
     {
         if (ir->QMMMscheme == eQMMMschemenormal && ir->ns_type == ensSIMPLE)
@@ -2143,6 +2144,11 @@ int gmx_grompp(int argc, char *argv[])
         {
             generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_ORIGINAL);
         }
+        if (ir->QMMMscheme != eQMMMschemeoniom)
+        {
+            std::vector<int> qmmmAtoms = qmmmAtomIndices(*ir, sys);
+            removeQmmmAtomCharges(&sys, qmmmAtoms);
+        }
     }
 
     if (ir->eI == eiMimic)

diff --git a/src/gromacs/mdlib/qmmm.cpp b/src/gromacs/mdlib/qmmm.cpp
index d3c9013ec42c48aeb5fd46af0181818554988028..86460336e3e9dec01f4741000b16064b45586c41 100644 (file)
--- a/src/gromacs/mdlib/qmmm.cpp
+++ b/src/gromacs/mdlib/qmmm.cpp
@@ -351,6 +351,70 @@ t_QMMMrec *mk_QMMMrec()
 } /* mk_QMMMrec */
 #endif
 
+std::vector<int> qmmmAtomIndices(const t_inputrec &ir, const gmx_mtop_t &mtop)
+{
+    const int           numQmmmGroups = ir.opts.ngQM;
+    const gmx_groups_t &groups        = mtop.groups;
+    std::vector<int>    qmmmAtoms;
+    for (int i = 0; i < numQmmmGroups; i++)
+    {
+        SystemAtomIterator aloop(mtop);
+        while (aloop.nextAtom())
+        {
+            int index = aloop.globalAtomNumber();
+            if (getGroupType(groups, egcQMMM, index) == i)
+            {
+                qmmmAtoms.push_back(index);
+            }
+        }
+        if (ir.QMMMscheme == eQMMMschemeoniom)
+        {
+            /* I assume that users specify the QM groups from small to
+             * big(ger) in the mdp file
+             */
+            gmx_mtop_ilistloop_all_t iloop = gmx_mtop_ilistloop_all_init(&mtop);
+            int nral1                      = 1 + NRAL(F_VSITE2);
+            int atomOffset                 = 0;
+            while (const InteractionLists *ilists = gmx_mtop_ilistloop_all_next(iloop, &atomOffset))
+            {
+                const InteractionList &ilist = (*ilists)[F_VSITE2];
+                for (int j = 0; j < ilist.size(); j += nral1)
+                {
+                    const int vsite = atomOffset + ilist.iatoms[j  ]; /* the vsite         */
+                    const int ai    = atomOffset + ilist.iatoms[j+1]; /* constructing atom */
+                    const int aj    = atomOffset + ilist.iatoms[j+2]; /* constructing atom */
+                    if (getGroupType(groups, egcQMMM, vsite) == getGroupType(groups, egcQMMM, ai)
+                        &&
+                        getGroupType(groups, egcQMMM, vsite) == getGroupType(groups, egcQMMM, aj))
+                    {
+                        /* this dummy link atom needs to be removed from qmmmAtoms
+                         * before making the QMrec of this layer!
+                         */
+                        qmmmAtoms.erase(std::remove_if(qmmmAtoms.begin(),
+                                                       qmmmAtoms.end(),
+                                                       [&vsite](int atom){return atom == vsite; }),
+                                        qmmmAtoms.end());
+                    }
+                }
+            }
+        }
+    }
+    return qmmmAtoms;
+}
+
+void removeQmmmAtomCharges(gmx_mtop_t *mtop, gmx::ArrayRef<const int> qmmmAtoms)
+{
+    int molb = 0;
+    for (int i = 0; i < qmmmAtoms.size(); i++)
+    {
+        int     indexInMolecule;
+        mtopGetMolblockIndex(mtop, qmmmAtoms[i], &molb, nullptr, &indexInMolecule);
+        t_atom *atom = &mtop->moltype[mtop->molblock[molb].type].atoms.atom[indexInMolecule];
+        atom->q  = 0.0;
+        atom->qB = 0.0;
+    }
+}
+
 void init_QMMMrec(const t_commrec  *cr,
                   const gmx_mtop_t *mtop,
                   const t_inputrec *ir,
@@ -362,11 +426,8 @@ void init_QMMMrec(const t_commrec  *cr,
      * simply contains true/false for QM and MM (the other) atoms.
      */
 
-    int                     *qm_arr = nullptr, vsite, ai, aj;
-    int                      qm_max = 0, qm_nr = 0, jmax;
     t_QMMMrec               *qr;
     t_MMrec                 *mm;
-    int                      a_offset;
 
     if (!GMX_QMMM)
     {
@@ -408,44 +469,28 @@ void init_QMMMrec(const t_commrec  *cr,
 
     /* small problem if there is only QM.... so no MM */
 
-    jmax = ir->opts.ngQM;
+    int numQmmmGroups = ir->opts.ngQM;
 
     if (qr->QMMMscheme == eQMMMschemeoniom)
     {
-        qr->nrQMlayers = jmax;
+        qr->nrQMlayers = numQmmmGroups;
     }
     else
     {
         qr->nrQMlayers = 1;
     }
 
-    const gmx_groups_t &groups = mtop->groups;
-
-    /* there are jmax groups of QM atoms. In case of multiple QM groups
+    /* there are numQmmmGroups groups of QM atoms. In case of multiple QM groups
      * I assume that the users wants to do ONIOM. However, maybe it
      * should also be possible to define more than one QM subsystem with
      * independent neighbourlists. I have to think about
      * that.. 11-11-2003
      */
-    snew(qr->qm, jmax);
-    for (int j = 0; j < jmax; j++)
+    std::vector<int> qmmmAtoms = qmmmAtomIndices(*ir, *mtop);
+    snew(qr->qm, numQmmmGroups);
+    for (int i = 0; i < numQmmmGroups; i++)
     {
         /* new layer */
-        SystemAtomIterator aloop(*mtop);
-        while (aloop.nextAtom())
-        {
-            int i = aloop.globalAtomNumber();
-            if (qm_nr >= qm_max)
-            {
-                qm_max += 1000;
-                srenew(qm_arr, qm_max);
-            }
-            if (getGroupType(groups, egcQMMM, i) == j)
-            {
-                /* hack for tip4p */
-                qm_arr[qm_nr++] = i;
-            }
-        }
         if (qr->QMMMscheme == eQMMMschemeoniom)
         {
             /* add the atoms to the bQMMM array
@@ -454,50 +499,10 @@ void init_QMMMrec(const t_commrec  *cr,
             /* I assume that users specify the QM groups from small to
              * big(ger) in the mdp file
              */
-            qr->qm[j] = mk_QMrec();
-            /* we need to throw out link atoms that in the previous layer
-             * existed to separate this QMlayer from the previous
-             * QMlayer. We use the iatoms array in the idef for that
-             * purpose. If all atoms defining the current Link Atom (Dummy2)
-             * are part of the current QM layer it needs to be removed from
-             * qm_arr[].  */
-
-            gmx_mtop_ilistloop_all_t iloop = gmx_mtop_ilistloop_all_init(mtop);
-            int nral1 = 1 + NRAL(F_VSITE2);
-            while (const InteractionLists *ilists = gmx_mtop_ilistloop_all_next(iloop, &a_offset))
-            {
-                const InteractionList &ilist = (*ilists)[F_VSITE2];
-                for (int i = 0; i < ilist.size(); i += nral1)
-                {
-                    vsite = a_offset + ilist.iatoms[i  ]; /* the vsite         */
-                    ai    = a_offset + ilist.iatoms[i+1]; /* constructing atom */
-                    aj    = a_offset + ilist.iatoms[i+2]; /* constructing atom */
-                    if (getGroupType(groups, egcQMMM, vsite) == getGroupType(groups, egcQMMM, ai)
-                        &&
-                        getGroupType(groups, egcQMMM, vsite) == getGroupType(groups, egcQMMM, aj))
-                    {
-                        /* this dummy link atom needs to be removed from the qm_arr
-                         * before making the QMrec of this layer!
-                         */
-                        for (int j = 0; j < qm_nr; j++)
-                        {
-                            if (qm_arr[j] == vsite)
-                            {
-                                /* drop the element */
-                                for (int l = i; l < qm_nr; l++)
-                                {
-                                    qm_arr[l] = qm_arr[l+1];
-                                }
-                                qm_nr--;
-                            }
-                        }
-                    }
-                }
-            }
-
+            qr->qm[i] = mk_QMrec();
             /* store QM atoms in this layer in the QMrec and initialise layer
              */
-            init_QMrec(j, qr->qm[j], qm_nr, qm_arr, mtop, ir);
+            init_QMrec(i, qr->qm[i], qmmmAtoms.size(), qmmmAtoms.data(), mtop, ir);
         }
     }
     if (qr->QMMMscheme != eQMMMschemeoniom)
@@ -510,19 +515,10 @@ void init_QMMMrec(const t_commrec  *cr,
          * TODO: Consider doing this in grompp instead.
          */
 
-        int molb = 0;
-        for (int k = 0; k < qm_nr; k++)
-        {
-            int     indexInMolecule;
-            mtopGetMolblockIndex(mtop, qm_arr[k], &molb, nullptr, &indexInMolecule);
-            t_atom *atom = &mtop->moltype[mtop->molblock[molb].type].atoms.atom[indexInMolecule];
-            atom->q  = 0.0;
-            atom->qB = 0.0;
-        }
         qr->qm[0] = mk_QMrec();
         /* store QM atoms in the QMrec and initialise
          */
-        init_QMrec(0, qr->qm[0], qm_nr, qm_arr, mtop, ir);
+        init_QMrec(0, qr->qm[0], qmmmAtoms.size(), qmmmAtoms.data(), mtop, ir);
 
         /* MM rec creation */
         mm               = mk_MMrec();

diff --git a/src/gromacs/mdlib/qmmm.h b/src/gromacs/mdlib/qmmm.h
index c8fb094c32bb99595166cb3c7a0d014676528888..78fb3b66b1151089eb3d69db39263881eba1824e 100644 (file)
--- a/src/gromacs/mdlib/qmmm.h
+++ b/src/gromacs/mdlib/qmmm.h
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2008, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2017,2018,2019, 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.
@@ -39,8 +39,11 @@
 
 #include "config.h"
 
+#include <vector>
+
 #include "gromacs/math/vectypes.h"
 #include "gromacs/mdlib/tgroup.h"
+#include "gromacs/utility/arrayref.h"
 
 #define GMX_QMMM (GMX_QMMM_MOPAC || GMX_QMMM_GAMESS || GMX_QMMM_GAUSSIAN || GMX_QMMM_ORCA)
 
@@ -146,4 +149,21 @@ real calculate_QMMM(const t_commrec  *cr,
  * called by system().
  */
 
+/*! \brief
+ * Return vector of atom indices for atoms in the QMMM region.
+ *
+ * \param[in] mtop Topology to use for populating array.
+ * \param[in] ir   Inputrec used in simulation.
+ * \returns Vector of atoms.
+ */
+std::vector<int> qmmmAtomIndices(const t_inputrec &ir, const gmx_mtop_t &mtop);
+
+/*! \brief
+ * Remove charges from QMMM atoms.
+ *
+ * \param[in] mtop Topology used for removing atoms.
+ * \param[in] qmmmAtoms ArrayRef to vector conatining qmmm atom indices.
+ */
+void removeQmmmAtomCharges(gmx_mtop_t *mtop, gmx::ArrayRef<const int> qmmmAtoms);
+
 #endif