diff --git a/CHANGES/v2.9.md b/CHANGES/v2.9.md
index 2e69f926d8..414f2893e5 100644
--- a/CHANGES/v2.9.md
+++ b/CHANGES/v2.9.md
@@ -61,4 +61,5 @@ Changes from version 2.8 which are relevant for users:
 ## Version 2.9.2 (tba)
   - Patch for GROMACS 2023 updated to the latest version
   - new patch for GROMACS 2024 
+  - Small optimization of \ref COORDINATION and related collective variables (\issue{1096}).
 
diff --git a/src/colvar/CoordinationBase.cpp b/src/colvar/CoordinationBase.cpp
index 9a6cd66eee..e7ab51611c 100644
--- a/src/colvar/CoordinationBase.cpp
+++ b/src/colvar/CoordinationBase.cpp
@@ -149,17 +149,21 @@ void CoordinationBase::calculate()
   const unsigned nn=nl->size();
   if(nt*stride*10>nn) nt=1;
 
+  const unsigned elementsPerRank = std::ceil(double(nn)/stride);
+  const unsigned int start= rank*elementsPerRank;
+  const unsigned int end = ((start + elementsPerRank)< nn)?(start + elementsPerRank): nn;
+
   #pragma omp parallel num_threads(nt)
   {
     std::vector<Vector> omp_deriv(getPositions().size());
     Tensor omp_virial;
 
     #pragma omp for reduction(+:ncoord) nowait
-    for(unsigned int i=rank; i<nn; i+=stride) {
+    for(unsigned int i=start; i<end; ++i) {
 
       Vector distance;
-      unsigned i0=nl->getClosePair(i).first;
-      unsigned i1=nl->getClosePair(i).second;
+      const unsigned i0=nl->getClosePair(i).first;
+      const unsigned i1=nl->getClosePair(i).second;
 
       if(getAbsoluteIndex(i0)==getAbsoluteIndex(i1)) continue;
 
@@ -187,7 +191,8 @@ void CoordinationBase::calculate()
     }
     #pragma omp critical
     if(nt>1) {
-      for(unsigned i=0; i<getPositions().size(); i++) deriv[i]+=omp_deriv[i];
+      for(unsigned i=0; i<getPositions().size(); i++)
+        deriv[i]+=omp_deriv[i];
       virial+=omp_virial;
     }
   }