From 7583e455a592e51679cd0b22cc5684d9fcc359bc Mon Sep 17 00:00:00 2001
From: Berk Hess <hess@kth.se>
Date: Thu, 13 Aug 2020 15:13:10 +0200
Subject: [PATCH] Add new FEP perturbation check functions

Added functions to check whether the system has perturbed masses
and perturbed constraints. Simplified the perturbed atoms check.
---
 src/gromacs/nbnxm/nbnxm_setup.cpp  |  2 +-
 src/gromacs/topology/mtop_util.cpp | 54 ++++++++++++++++++++++++------
 src/gromacs/topology/mtop_util.h   | 14 ++++----
 3 files changed, 52 insertions(+), 18 deletions(-)

diff --git a/src/gromacs/nbnxm/nbnxm_setup.cpp b/src/gromacs/nbnxm/nbnxm_setup.cpp
index e7994d84ba..bcb0b7d926 100644
--- a/src/gromacs/nbnxm/nbnxm_setup.cpp
+++ b/src/gromacs/nbnxm/nbnxm_setup.cpp
@@ -392,7 +392,7 @@ std::unique_ptr<nonbonded_verlet_t> init_nb_verlet(const gmx::MDLogger& mdlog,
 
     const bool haveMultipleDomains = havePPDomainDecomposition(cr);
 
-    bool           bFEP_NonBonded = (fr->efep != efepNO) && haveFepPerturbedNBInteractions(mtop);
+    bool           bFEP_NonBonded = (fr->efep != efepNO) && haveFepPerturbedNBInteractions(*mtop);
     PairlistParams pairlistParams(kernelSetup.kernelType, bFEP_NonBonded, ir->rlist, haveMultipleDomains);
 
     setupDynamicPairlistPruning(mdlog, ir, mtop, box, fr->ic, &pairlistParams);
diff --git a/src/gromacs/topology/mtop_util.cpp b/src/gromacs/topology/mtop_util.cpp
index fb4a6ea7d6..789a35eced 100644
--- a/src/gromacs/topology/mtop_util.cpp
+++ b/src/gromacs/topology/mtop_util.cpp
@@ -1051,23 +1051,55 @@ void convertAtomsToMtop(t_symtab* symtab, char** name, t_atoms* atoms, gmx_mtop_
     mtop->finalize();
 }
 
-bool haveFepPerturbedNBInteractions(const gmx_mtop_t* mtop)
+bool haveFepPerturbedNBInteractions(const gmx_mtop_t& mtop)
 {
-    for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
+    for (const gmx_moltype_t& molt : mtop.moltype)
     {
-        const gmx_molblock_t& molb = mtop->molblock[mb];
-        const gmx_moltype_t&  molt = mtop->moltype[molb.type];
-        for (int m = 0; m < molb.nmol; m++)
+        for (int a = 0; a < molt.atoms.nr; a++)
         {
-            for (int a = 0; a < molt.atoms.nr; a++)
+            if (PERTURBED(molt.atoms.atom[a]))
             {
-                const t_atom& atom = molt.atoms.atom[a];
-                if (PERTURBED(atom))
-                {
-                    return true;
-                }
+                return true;
+            }
+        }
+    }
+
+    return false;
+}
+
+bool haveFepPerturbedMasses(const gmx_mtop_t& mtop)
+{
+    for (const gmx_moltype_t& molt : mtop.moltype)
+    {
+        for (int a = 0; a < molt.atoms.nr; a++)
+        {
+            const t_atom& atom = molt.atoms.atom[a];
+            if (atom.m != atom.mB)
+            {
+                return true;
             }
         }
     }
+
+    return false;
+}
+
+bool havePerturbedConstraints(const gmx_mtop_t& mtop)
+{
+    // This code assumes that all perturbed constraints parameters are actually used
+    const auto& ffparams = mtop.ffparams;
+
+    for (gmx::index i = 0; i < gmx::ssize(ffparams.functype); i++)
+    {
+        if (ffparams.functype[i] == F_CONSTR || ffparams.functype[i] == F_CONSTRNC)
+        {
+            const auto& iparams = ffparams.iparams[i];
+            if (iparams.constr.dA != iparams.constr.dB)
+            {
+                return true;
+            }
+        }
+    }
+
     return false;
 }
diff --git a/src/gromacs/topology/mtop_util.h b/src/gromacs/topology/mtop_util.h
index 9a027bac7d..7d0527e654 100644
--- a/src/gromacs/topology/mtop_util.h
+++ b/src/gromacs/topology/mtop_util.h
@@ -277,11 +277,13 @@ std::vector<int> get_atom_index(const gmx_mtop_t* mtop);
  */
 void convertAtomsToMtop(t_symtab* symtab, char** name, t_atoms* atoms, gmx_mtop_t* mtop);
 
-/*! \brief Checks if the non-bonded FEP should be performed in this run.
- *
- * \param[in]  mtop  Molecular topology.
- * \returns Whether FEP non-bonded is requested.
- */
-bool haveFepPerturbedNBInteractions(const gmx_mtop_t* mtop);
+//! Checks and returns whether non-bonded interactions are perturbed for free-energy calculations
+bool haveFepPerturbedNBInteractions(const gmx_mtop_t& mtop);
+
+//! Checks whether masses are perturbed for free-energy calculations
+bool haveFepPerturbedMasses(const gmx_mtop_t& mtop);
+
+//! Checks whether constraints are perturbed for free-energy calculations
+bool havePerturbedConstraints(const gmx_mtop_t& mtop);
 
 #endif
-- 
2.22.0