Change grompp disres cost to linear

The cost of adding distance restraints (actually their parameters)
was quadratic in the number of restraints and has now been changed
to linear.

Fixes #3457

diff --git a/docs/release-notes/2021/major/performance.rst b/docs/release-notes/2021/major/performance.rst
index 0c895e763705df81a1e4c6295c126aac25c8819f..0931f3843f4aeeec7ef0a270264c21ed41f2f7e7 100644 (file)
--- a/docs/release-notes/2021/major/performance.rst
+++ b/docs/release-notes/2021/major/performance.rst
@@ -12,3 +12,11 @@ Extend supported use-cases for GPU version of update and constraints
 
 GPU version of update and constraints can now be used for FEP, except mass and constraints
 free-energy perturbation.
+       
+Reduce time spent in grompp with large numbers of distance restraints
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+The time `gmx grompp` spent processing distance restraint has been
+changed from quadratic in the number of restraints to linear.
+       
+:issue:`3457`

diff --git a/src/gromacs/gmxpreprocess/convparm.cpp b/src/gromacs/gmxpreprocess/convparm.cpp
index c9efb733cbe83608742057356d7cdd837eff6c5d..225086c52f6dd3aaad4e99bc7615b83666381425 100644 (file)
--- a/src/gromacs/gmxpreprocess/convparm.cpp
+++ b/src/gromacs/gmxpreprocess/convparm.cpp
@@ -451,7 +451,6 @@ static int enter_params(gmx_ffparams_t*           ffparams,
                         bool                      bAppend)
 {
     t_iparams newparam;
-    int       type;
     int       rc;
 
     if ((rc = assign_param(ftype, &newparam, forceparams, comb, reppow)) < 0)
@@ -462,21 +461,33 @@ static int enter_params(gmx_ffparams_t*           ffparams,
 
     if (!bAppend)
     {
-        for (type = start; (type < ffparams->numTypes()); type++)
+        if (ftype != F_DISRES)
         {
-            if (ffparams->functype[type] == ftype)
+            for (int type = start; type < ffparams->numTypes(); type++)
             {
-                if (memcmp(&newparam, &ffparams->iparams[type], static_cast<size_t>(sizeof(newparam))) == 0)
+                // Note that the first condition is always met by starting the loop at start
+                if (ffparams->functype[type] == ftype
+                    && memcmp(&newparam, &ffparams->iparams[type], static_cast<size_t>(sizeof(newparam))) == 0)
                 {
                     return type;
                 }
             }
         }
+        else
+        {
+            // Distance restraints should have unique labels and pairs with the same label
+            // should be consecutive, so we here we only need to check the last type in the list.
+            // This changes the complexity from quadratic to linear in the number of restraints.
+            const int type = ffparams->numTypes() - 1;
+            if (type >= 0 && ffparams->functype[type] == ftype
+                && memcmp(&newparam, &ffparams->iparams[type], static_cast<size_t>(sizeof(newparam))) == 0)
+            {
+                return type;
+            }
+        }
     }
-    else
-    {
-        type = ffparams->numTypes();
-    }
+
+    const int type = ffparams->numTypes();
 
     ffparams->iparams.push_back(newparam);
     ffparams->functype.push_back(ftype);