From d57e07180eb222aa0f8588cb0d4045c9eabca8bb Mon Sep 17 00:00:00 2001
From: Giacomo Fiorin <giacomo.fiorin@gmail.com>
Date: Thu, 12 Feb 2026 14:43:02 -0500
Subject: [PATCH 01/15] refactor: Use a dedicated class to store system
 boundaries

---
 devel-tools/compile-gromacs.sh             |   4 +-
 misc_interfaces/stubs/colvarproxy_stub.cpp |   7 +-
 src/colvarmodule.h                         |   1 +
 src/colvarproxy_system.cpp                 |  41 +----
 src/colvarproxy_system.h                   | 117 +-------------
 src/colvars_system.h                       | 172 +++++++++++++++++++++
 vmd/src/colvarproxy_vmd.C                  |   2 +-
 7 files changed, 189 insertions(+), 155 deletions(-)
 create mode 100644 src/colvars_system.h

diff --git a/devel-tools/compile-gromacs.sh b/devel-tools/compile-gromacs.sh
index e95ead5aa..396fd7d3e 100755
--- a/devel-tools/compile-gromacs.sh
+++ b/devel-tools/compile-gromacs.sh
@@ -45,8 +45,8 @@ compile_gromacs_target() {
     while [ $# -ge 1 ]; do
         if [ "${1,,}" = "debug" ]; then
             GMX_BUILD_TYPE=Debug
-            GMX_BUILD_OPTS+=(-DCMAKE_BUILD_TYPE=Debug -DCMAKE_VERBOSE_MAKEFILE=yes)
-            GMX_BUILD_OPTS+=(-DCOLVARS_DEBUG=ON)
+            GMX_BUILD_OPTS+=(-DCMAKE_VERBOSE_MAKEFILE=yes)
+            GMX_BUILD_OPTS+=(-DCMAKE_CXX_FLAGS="-DCOLVARS_DEBUG=1")
         else
             GMX_INSTALL_DIR=${1}
         fi
diff --git a/misc_interfaces/stubs/colvarproxy_stub.cpp b/misc_interfaces/stubs/colvarproxy_stub.cpp
index 643bc4f72..9469f1a2d 100644
--- a/misc_interfaces/stubs/colvarproxy_stub.cpp
+++ b/misc_interfaces/stubs/colvarproxy_stub.cpp
@@ -52,13 +52,14 @@ colvarproxy_stub::~colvarproxy_stub()
 
 int colvarproxy_stub::setup()
 {
-  boundaries_type = boundaries_non_periodic;
-  reset_pbc_lattice();
-  cvmodule->it = cvmodule->it_restart = 0;
+  // Set system boundaries as the default (non-periodic)
+  boundaries_.reset();
 
   if (cvmodule) {
+    cvmodule->it = cvmodule->it_restart = 0;
     return cvmodule->update_engine_parameters();
   }
+
   return COLVARS_OK;
 }
 
diff --git a/src/colvarmodule.h b/src/colvarmodule.h
index 16b0e0822..81e6c75ab 100644
--- a/src/colvarmodule.h
+++ b/src/colvarmodule.h
@@ -228,6 +228,7 @@ static inline real acos(real const &x)
   template <class T> class matrix2d;
   class quaternion;
   class rotation;
+  class system_boundary_conditions;
 
   class usage;
   class memory_stream;
diff --git a/src/colvarproxy_system.cpp b/src/colvarproxy_system.cpp
index 1daf31750..92c91f734 100644
--- a/src/colvarproxy_system.cpp
+++ b/src/colvarproxy_system.cpp
@@ -22,12 +22,10 @@ colvarproxy_system::colvarproxy_system()
   timestep_ = 1.0;
   target_temperature_ = 0.0;
   boltzmann_ = 0.001987191; // Default: kcal/mol/K
-  boundaries_type = boundaries_unsupported;
   total_force_requested = false;
   indirect_lambda_biasing_force = 0.0;
   cached_alch_lambda_changed = false;
   cached_alch_lambda = -1.0;
-  reset_pbc_lattice();
 }
 
 
@@ -90,47 +88,10 @@ bool colvarproxy_system::total_forces_same_step() const
 }
 
 
-void colvarproxy_system::update_pbc_lattice()
-{
-  // Periodicity is assumed in all directions
-
-  if (boundaries_type == boundaries_unsupported ||
-      boundaries_type == boundaries_non_periodic) {
-    cvm::error_static("Error: setting PBC lattice with unsupported boundaries.\n",
-               COLVARS_BUG_ERROR);
-    return;
-  }
-
-  {
-    cvm::rvector const v = cvm::rvector::outer(unit_cell_y, unit_cell_z);
-    reciprocal_cell_x = v/(v*unit_cell_x);
-  }
-  {
-    cvm::rvector const v = cvm::rvector::outer(unit_cell_z, unit_cell_x);
-    reciprocal_cell_y = v/(v*unit_cell_y);
-  }
-  {
-    cvm::rvector const v = cvm::rvector::outer(unit_cell_x, unit_cell_y);
-    reciprocal_cell_z = v/(v*unit_cell_z);
-  }
-}
-
-
-void colvarproxy_system::reset_pbc_lattice()
-{
-  unit_cell_x.reset();
-  unit_cell_y.reset();
-  unit_cell_z.reset();
-  reciprocal_cell_x.reset();
-  reciprocal_cell_y.reset();
-  reciprocal_cell_z.reset();
-}
-
-
 cvm::rvector colvarproxy_system::position_distance(cvm::atom_pos const &pos1,
                                                    cvm::atom_pos const &pos2) const
 {
-  return position_distance_internal(pos1, pos2);
+  return boundaries_.position_distance(pos1, pos2);
 }
 
 
diff --git a/src/colvarproxy_system.h b/src/colvarproxy_system.h
index 3746a5b66..4af4cc69a 100644
--- a/src/colvarproxy_system.h
+++ b/src/colvarproxy_system.h
@@ -10,6 +10,7 @@
 #ifndef COLVARPROXY_SYSTEM_H
 #define COLVARPROXY_SYSTEM_H
 
+#include "colvars_system.h"
 
 /// Methods for accessing the simulation system (PBCs, integrator, etc)
 class colvarproxy_system {
@@ -85,47 +86,17 @@ class colvarproxy_system {
   /// Pass restraint energy value for current timestep to MD engine
   virtual void add_energy(cvm::real energy);
 
-  /// Use the PBC functions from the Colvars library (as opposed to MD engine)
+  /// Account for system boundaries within the Colvars library (as opposed to using the MD engine)
   inline bool & use_internal_pbc() { return use_internal_pbc_; }
 
   /// Get the PBC-aware distance vector between two positions (using the MD engine's convention)
   virtual cvm::rvector position_distance(cvm::atom_pos const &pos1,
                                          cvm::atom_pos const &pos2) const;
 
-  /// Inline version of position_distance()
-  cvm::rvector position_distance_internal(cvm::atom_pos const &pos1,
-                                          cvm::atom_pos const &pos2) const;
-
-  /// Kernel used by position_distance_internal()
-  /// @param pos1 First position
-  /// @param pos2 Second position
-  /// @param a Unit cell vector
-  /// @param b Unit cell vector
-  /// @param c Unit cell vector
-  /// @param a_r Reciprocal cell vector
-  /// @param b_r Reciprocal cell vector
-  /// @param c_r Reciprocal cell vector
-  /// @param a_p Is this dimension periodic?
-  /// @param b_p Is this dimension periodic?
-  /// @param c_p Is this dimension periodic?
-  /// @return (pos2 - pos1) without periodicity, minimum-image difference otherwise
-  static cvm::rvector position_distance_kernel(cvm::atom_pos const &pos1,
-                                               cvm::atom_pos const &pos2,
-                                               cvm::rvector const &a,
-                                               cvm::rvector const &b,
-                                               cvm::rvector const &c,
-                                               cvm::rvector const &a_r,
-                                               cvm::rvector const &b_r,
-                                               cvm::rvector const &c_r,
-                                               bool a_p = false,
-                                               bool b_p = false,
-                                               bool c_p = false);
-
-  /// Recompute PBC reciprocal lattice (assumes XYZ periodicity)
-  void update_pbc_lattice();
-
-  /// Set the lattice vectors to zero
-  void reset_pbc_lattice();
+  /// Get the current system boundary conditions
+  inline cvm::system_boundary_conditions const &get_system_boundaries() const {
+    return boundaries_;
+  }
 
   /// \brief Tell the proxy whether total forces are needed (they may not
   /// always be available)
@@ -214,81 +185,9 @@ class colvarproxy_system {
   /// Use the PBC functions from the Colvars library (as opposed to MD engine)
   bool use_internal_pbc_ = false;
 
-  /// Type of boundary conditions defined for the current computation
-  ///
-  /// Orthogonal and triclinic cells are made available to objects.
-  /// For any other conditions (mixed periodicity, triclinic cells in LAMMPS)
-  /// minimum-image distances are computed by the host engine by default
-  enum Boundaries_type {
-    boundaries_non_periodic,
-    boundaries_pbc_ortho,
-    boundaries_pbc_triclinic,
-    boundaries_unsupported
-  };
-
-  /// Type of boundary conditions
-  Boundaries_type boundaries_type;
-
-  /// Bravais lattice vectors
-  cvm::rvector unit_cell_x, unit_cell_y, unit_cell_z;
-
-  /// Reciprocal lattice vectors
-  cvm::rvector reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z;
+  /// Current system boundary conditions
+  cvm::system_boundary_conditions boundaries_;
 };
 
 
-inline cvm::rvector colvarproxy_system::position_distance_internal(cvm::atom_pos const &pos1,
-                                                                   cvm::atom_pos const &pos2) const
-{
-  if (boundaries_type == boundaries_unsupported) {
-    cvm::error_static("Error: unsupported boundary conditions.\n", COLVARS_INPUT_ERROR);
-    return cvm::rvector({0.0, 0.0, 0.0});
-  }
-
-  if (boundaries_type == boundaries_non_periodic) {
-    return pos2 - pos1;
-  }
-
-  // Periodicity flags are hard-coded, because this is the only case supported so far other than
-  // the two above
-  return position_distance_kernel(pos1, pos2, unit_cell_x, unit_cell_y, unit_cell_z,
-                                  reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z,
-                                  true, true, true);
-}
-
-
-inline cvm::rvector colvarproxy_system::position_distance_kernel(cvm::atom_pos const &pos1,
-                                                                 cvm::atom_pos const &pos2,
-                                                                 cvm::rvector const &a,
-                                                                 cvm::rvector const &b,
-                                                                 cvm::rvector const &c,
-                                                                 cvm::rvector const &a_r,
-                                                                 cvm::rvector const &b_r,
-                                                                 cvm::rvector const &c_r,
-                                                                 bool a_p,
-                                                                 bool b_p,
-                                                                 bool c_p)
-{
-  cvm::rvector diff = (pos2 - pos1);
-
-  cvm::real const x_shift = std::floor(a_r * diff + 0.5);
-  cvm::real const y_shift = std::floor(b_r * diff + 0.5);
-  cvm::real const z_shift = std::floor(c_r * diff + 0.5);
-
-  if (a_p) {
-    diff.x -= x_shift * a.x + y_shift * b.x + z_shift * c.x;
-  }
-
-  if (b_p) {
-    diff.y -= x_shift * a.y + y_shift * b.y + z_shift * c.y;
-  }
-
-  if (c_p) {
-    diff.z -= x_shift * a.z + y_shift * b.z + z_shift * c.z;
-  }
-
-  return diff;
-}
-
-
 #endif
diff --git a/src/colvars_system.h b/src/colvars_system.h
new file mode 100644
index 000000000..f7b4d03b0
--- /dev/null
+++ b/src/colvars_system.h
@@ -0,0 +1,172 @@
+// -*- c++ -*-
+
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/Colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
+#ifndef COLVARS_SYSTEM_H
+#define COLVARS_SYSTEM_H
+
+#include "colvartypes.h"
+
+/// Class to store the system's boundary conditions
+class colvarmodule::system_boundary_conditions {
+public:
+
+  /// Type of boundary conditions defined for the current computation
+  enum class types {
+    non_periodic,   /// All three dimensions are non-periodic
+    mixed,          /// Some dimensions are periodic, but others are not
+    pbc_orthogonal, /// All three dimensions are periodic, lattice vectors are orthogonal
+    pbc_triclinic,  /// All three dimensions are periodic, lattice vectors are *not* orthogonal
+    unsupported     /// Unsupported boundary conditions
+  };
+
+  /// Type of boundary conditions in the current computation
+  inline COLVARS_HOST_DEVICE types type() const { return type_; }
+
+  /// Set the type of boundary explicitly
+  inline COLVARS_HOST_DEVICE void set_type(types t) { type_ = t; }
+
+  /// Compute the distance between two positions
+  cvm::rvector position_distance(cvm::atom_pos const &pos1, cvm::atom_pos const &pos2) const;
+
+  /// Reset to defaults (non-periodic)
+  inline COLVARS_HOST_DEVICE void reset() {
+    periodic_x = periodic_y = periodic_z = false;
+    type_ = types::non_periodic;
+    unit_cell_x.reset();
+    unit_cell_y.reset();
+    unit_cell_z.reset();
+    reciprocal_cell_x.reset();
+    reciprocal_cell_y.reset();
+    reciprocal_cell_z.reset();
+  }
+
+  /// Set from explicit boundary configuration
+  void set_boundaries(bool periodic_x_in, bool periodic_y_in, bool periodic_z_in,
+                      cvm::rvector const &A, cvm::rvector const &B, cvm::rvector const &C);
+
+protected:
+
+  /// Type of boundary conditions in the current computation
+  types type_ = types::non_periodic;
+
+  /// Bravais lattice vectors
+  cvm::rvector unit_cell_x, unit_cell_y, unit_cell_z;
+
+  /// Reciprocal lattice vectors
+  cvm::rvector reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z;
+
+  /// Periodic flags in each dimension
+  bool periodic_x = false, periodic_y = false, periodic_z = false;
+};
+
+
+/// Set from explicit boundary configuration
+inline COLVARS_HOST_DEVICE void
+cvm::system_boundary_conditions::set_boundaries(bool periodic_x_in, bool periodic_y_in,
+                                                bool periodic_z_in, cvm::rvector const &A,
+                                                cvm::rvector const &B, cvm::rvector const &C)
+{
+  constexpr double diagonal_tol2 = 1.0e-10;
+
+  periodic_x = periodic_x_in;
+  periodic_y = periodic_y_in;
+  periodic_z = periodic_z_in;
+
+  if ((periodic_x == periodic_y) && (periodic_x == periodic_z)) {
+    if (periodic_x) {
+      // Temporarily set as fully-periodic & orthogonal; will check below for triclinic
+      set_type(types::pbc_orthogonal);
+    } else {
+      set_type(types::non_periodic);
+      return;
+    }
+  } else {
+    set_type(types::mixed);
+  }
+
+  bool off_diagonal = false;
+
+  if (periodic_x) {
+    unit_cell_x = A;
+    if ((A.y * A.y) > diagonal_tol2 || (A.z * A.z) > diagonal_tol2) {
+      off_diagonal = true;
+    } else {
+      reciprocal_cell_x = unit_cell_x/unit_cell_x.norm2();
+    }
+  } else {
+    unit_cell_x = {0.0, 0.0, 0.0};
+  }
+
+  if (periodic_y) {
+    unit_cell_y = B;
+    if ((B.x * B.x) > diagonal_tol2 || (B.z * B.z) > diagonal_tol2) {
+      off_diagonal = true;
+    } else {
+      reciprocal_cell_y = unit_cell_y/unit_cell_y.norm2();
+    }
+  } else {
+    unit_cell_y = {0.0, 0.0, 0.0};
+  }
+
+  if (periodic_z) {
+    unit_cell_z = C;
+    if ((C.x * C.x) > diagonal_tol2 || (C.y * C.y) > diagonal_tol2) {
+      off_diagonal = true;
+    } else {
+      reciprocal_cell_z = unit_cell_z/unit_cell_z.norm2();
+    }
+  } else {
+    unit_cell_z = {0.0, 0.0, 0.0};
+  }
+
+  if (type() == types::pbc_orthogonal && off_diagonal) {
+    set_type(types::pbc_triclinic);
+  }
+
+  if (type() == types::pbc_triclinic) {
+    cvm::rvector const v_yz = cvm::rvector::outer(unit_cell_y, unit_cell_z);
+    reciprocal_cell_x = v_yz / (v_yz * unit_cell_x);
+    cvm::rvector const v_zx = cvm::rvector::outer(unit_cell_z, unit_cell_x);
+    reciprocal_cell_y = v_zx / (v_zx * unit_cell_y);
+    cvm::rvector const v_xy = cvm::rvector::outer(unit_cell_x, unit_cell_y);
+    reciprocal_cell_z = v_xy / (v_xy * unit_cell_z);
+  }
+}
+
+
+inline COLVARS_HOST_DEVICE
+cvm::rvector cvm::system_boundary_conditions::position_distance(cvm::atom_pos const &pos1,
+                                                                cvm::atom_pos const &pos2) const
+{
+  cvm::rvector diff = (pos2 - pos1);
+
+  if (type() == types::non_periodic || type() == types::unsupported) {
+    return diff;
+  }
+
+  cvm::real const x_shift = ::floor(reciprocal_cell_x * diff + 0.5);
+  cvm::real const y_shift = ::floor(reciprocal_cell_y * diff + 0.5);
+  cvm::real const z_shift = ::floor(reciprocal_cell_z * diff + 0.5);
+
+  if (periodic_x) {
+    diff.x -= x_shift * unit_cell_x.x + y_shift * unit_cell_y.x + z_shift * unit_cell_z.x;
+  }
+
+  if (periodic_y) {
+    diff.y -= x_shift * unit_cell_x.y + y_shift * unit_cell_y.y + z_shift * unit_cell_z.y;
+  }
+
+  if (periodic_z) {
+    diff.z -= x_shift * unit_cell_x.z + y_shift * unit_cell_y.z + z_shift * unit_cell_z.z;
+  }
+
+  return diff;
+}
+
+#endif
diff --git a/vmd/src/colvarproxy_vmd.C b/vmd/src/colvarproxy_vmd.C
index e23450066..5a0d4dd79 100644
--- a/vmd/src/colvarproxy_vmd.C
+++ b/vmd/src/colvarproxy_vmd.C
@@ -805,7 +805,7 @@ void colvarproxy_vmd::compute_voldata(VolumetricData const *voldata,
   cvm::atom_pos const origin(0.0, 0.0, 0.0);
   for (; i < atoms->size(); ++i) {
     // Wrap around the origin
-    cvm::rvector const wrapped_pos = position_distance(
+    cvm::rvector const wrapped_pos = boundaries_.position_distance(
       origin, cvm::atom_pos(atoms->pos_x(i), atoms->pos_y(i), atoms->pos_z(i)));
     coord[0] = internal_to_angstrom(wrapped_pos.x);
     coord[1] = internal_to_angstrom(wrapped_pos.y);

From 06a8479cc24e0ca885e010f60e759e86b68fd221 Mon Sep 17 00:00:00 2001
From: Giacomo Fiorin <giacomo.fiorin@gmail.com>
Date: Thu, 12 Feb 2026 15:29:45 -0500
Subject: [PATCH 02/15] refactor: Use boundary conditions object in coordNum

---
 src/colvarcomp_coordnums.cpp | 15 +++++++++------
 src/colvarcomp_coordnums.h   |  8 ++++----
 2 files changed, 13 insertions(+), 10 deletions(-)

diff --git a/src/colvarcomp_coordnums.cpp b/src/colvarcomp_coordnums.cpp
index 16dd3b22f..62ebcf4de 100644
--- a/src/colvarcomp_coordnums.cpp
+++ b/src/colvarcomp_coordnums.cpp
@@ -203,6 +203,8 @@ void inline colvar::coordnum::main_loop()
 
   bool *pairlist_elem = pairlist.get();
 
+  auto const &boundary_conditions = cvm::main()->proxy->get_system_boundaries();
+
   for (size_t i = 0; i < group1_num_coords; ++i) {
 
     cvm::real const x1 = use_group1_com ? group1_com.x : group1->pos_x(i);
@@ -233,12 +235,12 @@ void inline colvar::coordnum::main_loop()
                                                             x1, y1, z1, x2, y2, z2,
                                                             gx1, gy1, gz1, gx2, gy2, gz2,
                                                             tolerance, tolerance_l2_max,
-                                                            cvmodule) :
+                                                            boundary_conditions) :
          compute_pair_coordnum<flags>(inv_r0_vec, inv_r0sq_vec, en, ed,
                                       x1, y1, z1, x2, y2, z2,
                                       gx1, gy1, gz1, gx2, gy2, gz2,
                                       tolerance, tolerance_l2_max,
-                                      cvmodule) ) :
+                                      boundary_conditions) ) :
         0.0;
 
       if ((flags & ef_use_pairlist) && (flags & ef_rebuild_pairlist)) {
@@ -447,7 +449,7 @@ void colvar::h_bond::calc_value()
                                                         atom_groups[0]->grad_y(1),
                                                         atom_groups[0]->grad_z(1),
                                                         0.0, 1.0e20,
-                                                        cvmodule);
+                                                        cvm::main()->proxy->get_system_boundaries());
   // Skip the gradient
 }
 
@@ -479,7 +481,7 @@ void colvar::h_bond::calc_gradients()
                                          atom_groups[0]->grad_y(1),
                                          atom_groups[0]->grad_z(1),
                                          0.0, 1.0e20,
-                                         cvmodule);
+                                         cvm::main()->proxy->get_system_boundaries());
 }
 
 
@@ -493,6 +495,7 @@ template <int flags> inline void colvar::selfcoordnum::selfcoordnum_sequential_l
 {
   size_t const n = group1->size();
   bool *pairlist_elem = pairlist.get();
+  auto const &boundary_conditions = cvm::main()->proxy->get_system_boundaries();
 
   for (size_t i = 0; i < n - 1; i++) {
 
@@ -522,12 +525,12 @@ template <int flags> inline void colvar::selfcoordnum::selfcoordnum_sequential_l
                                                             x1, y1, z1, x2, y2, z2,
                                                             gx1, gy1, gz1, gx2, gy2, gz2,
                                                             tolerance, tolerance_l2_max,
-                                                            cvmodule) :
+                                                            boundary_conditions) :
          compute_pair_coordnum<flags>(inv_r0_vec, inv_r0sq_vec, en, ed,
                                       x1, y1, z1, x2, y2, z2,
                                       gx1, gy1, gz1, gx2, gy2, gz2,
                                       tolerance, tolerance_l2_max,
-                                      cvmodule) ) :
+                                      boundary_conditions) ) :
         0.0;
 
       if ((flags & ef_use_pairlist) && (flags & ef_rebuild_pairlist)) {
diff --git a/src/colvarcomp_coordnums.h b/src/colvarcomp_coordnums.h
index 650fa5b2d..e7dac7510 100644
--- a/src/colvarcomp_coordnums.h
+++ b/src/colvarcomp_coordnums.h
@@ -56,7 +56,7 @@ class colvar::coordnum : public colvar::cvc {
                                          cvm::real &g1x, cvm::real &g1y, cvm::real &g1z,
                                          cvm::real &g2x, cvm::real &g2y, cvm::real &g2z,
                                          cvm::real pairlist_tol, cvm::real pairlist_tol_l2_max,
-                                         colvarmodule *cvmodule);
+                                         cvm::system_boundary_conditions const &boundary_conditions);
 
   /// Workhorse function
   template <bool use_group1_com, bool use_group2_com, int flags> int compute_coordnum();
@@ -247,13 +247,13 @@ inline cvm::real colvar::coordnum::compute_pair_coordnum(cvm::rvector const &inv
                                                          cvm::real& g2z,
                                                          cvm::real pairlist_tol,
                                                          cvm::real pairlist_tol_l2_max,
-                                                         colvarmodule *cvmodule_in)
+                                                         cvm::system_boundary_conditions const &boundary_conditions)
 {
   const cvm::atom_pos pos1{a1x, a1y, a1z};
   const cvm::atom_pos pos2{a2x, a2y, a2z};
   cvm::rvector const diff = (flags & ef_use_internal_pbc)
-                                ? cvmodule_in->proxy->position_distance_internal(pos1, pos2)
-                                : cvmodule_in->proxy->position_distance(pos1, pos2);
+                                ? boundary_conditions.position_distance(pos1, pos2)
+                                : cvm::main()->proxy->position_distance(pos1, pos2);
 
   cvm::rvector const scal_diff(diff.x * inv_r0_vec.x,
                                diff.y * inv_r0_vec.y,

From 52c09e53024e97a086e923abdc7f2dd283daa14f Mon Sep 17 00:00:00 2001
From: Giacomo Fiorin <giacomo.fiorin@gmail.com>
Date: Fri, 13 Feb 2026 17:43:23 -0500
Subject: [PATCH 03/15] opt: Stop using the proxy directly to compute distances

---
 src/colvarcomp.cpp           |  2 ++
 src/colvarcomp.h             |  3 +++
 src/colvarcomp_angles.cpp    | 12 ++++++------
 src/colvarcomp_coordnums.cpp |  7 ++-----
 src/colvarcomp_distances.cpp | 36 ++++++++++++++++--------------------
 src/colvarmodule.cpp         |  7 -------
 src/colvarmodule.h           |  5 -----
 src/colvars_system.h         |  6 ++++++
 8 files changed, 35 insertions(+), 43 deletions(-)

diff --git a/src/colvarcomp.cpp b/src/colvarcomp.cpp
index 620206444..e0f6b69e7 100644
--- a/src/colvarcomp.cpp
+++ b/src/colvarcomp.cpp
@@ -456,6 +456,8 @@ int colvar::cvc::set_param(std::string const &param_name,
 
 void colvar::cvc::read_data()
 {
+  boundary_conditions = cvm::main()->proxy->get_system_boundaries();
+
   if (is_enabled(f_cvc_explicit_atom_groups)) {
     for (auto agi = atom_groups.begin(); agi != atom_groups.end(); agi++) {
       auto &atoms = *(*agi);
diff --git a/src/colvarcomp.h b/src/colvarcomp.h
index 4634ed13a..694580e2b 100644
--- a/src/colvarcomp.h
+++ b/src/colvarcomp.h
@@ -349,6 +349,9 @@ class colvar::cvc
 
   /// \brief CVC-specific default colvar width (default: not provided)
   cvm::real width = 0.0;
+
+  /// Boundary conditions for the system, copied from the engine when needed
+  cvm::system_boundary_conditions boundary_conditions;
 };
 
 
diff --git a/src/colvarcomp_angles.cpp b/src/colvarcomp_angles.cpp
index d6fc402c3..370a94993 100644
--- a/src/colvarcomp_angles.cpp
+++ b/src/colvarcomp_angles.cpp
@@ -64,11 +64,11 @@ void colvar::angle::calc_value()
   cvm::atom_pos const g3_pos = group3->center_of_mass();
 
   r21  = is_enabled(f_cvc_pbc_minimum_image) ?
-    cvmodule->position_distance(g2_pos, g1_pos) :
+    boundary_conditions.position_distance(g2_pos, g1_pos) :
     g1_pos - g2_pos;
   r21l = r21.norm();
   r23  = is_enabled(f_cvc_pbc_minimum_image) ?
-    cvmodule->position_distance(g2_pos, g3_pos) :
+    boundary_conditions.position_distance(g2_pos, g3_pos) :
     g3_pos - g2_pos;
   r23l = r23.norm();
 
@@ -159,7 +159,7 @@ void colvar::dipole_angle::calc_value()
   r21 = group1->dipole();
   r21l = r21.norm();
   r23  = is_enabled(f_cvc_pbc_minimum_image) ?
-    cvmodule->position_distance(g2_pos, g3_pos) :
+    boundary_conditions.position_distance(g2_pos, g3_pos) :
     g3_pos - g2_pos;
   r23l = r23.norm();
 
@@ -282,13 +282,13 @@ void colvar::dihedral::calc_value()
 
   // Usual sign convention: r12 = r2 - r1
   r12 = is_enabled(f_cvc_pbc_minimum_image) ?
-    cvmodule->position_distance(g1_pos, g2_pos) :
+    boundary_conditions.position_distance(g1_pos, g2_pos) :
     g2_pos - g1_pos;
   r23 = is_enabled(f_cvc_pbc_minimum_image) ?
-    cvmodule->position_distance(g2_pos, g3_pos) :
+    boundary_conditions.position_distance(g2_pos, g3_pos) :
     g3_pos - g2_pos;
   r34 = is_enabled(f_cvc_pbc_minimum_image) ?
-    cvmodule->position_distance(g3_pos, g4_pos) :
+    boundary_conditions.position_distance(g3_pos, g4_pos) :
     g4_pos - g3_pos;
 
   cvm::rvector const n1 = cvm::rvector::outer(r12, r23);
diff --git a/src/colvarcomp_coordnums.cpp b/src/colvarcomp_coordnums.cpp
index 62ebcf4de..b40418fd6 100644
--- a/src/colvarcomp_coordnums.cpp
+++ b/src/colvarcomp_coordnums.cpp
@@ -203,8 +203,6 @@ void inline colvar::coordnum::main_loop()
 
   bool *pairlist_elem = pairlist.get();
 
-  auto const &boundary_conditions = cvm::main()->proxy->get_system_boundaries();
-
   for (size_t i = 0; i < group1_num_coords; ++i) {
 
     cvm::real const x1 = use_group1_com ? group1_com.x : group1->pos_x(i);
@@ -449,7 +447,7 @@ void colvar::h_bond::calc_value()
                                                         atom_groups[0]->grad_y(1),
                                                         atom_groups[0]->grad_z(1),
                                                         0.0, 1.0e20,
-                                                        cvm::main()->proxy->get_system_boundaries());
+                                                        boundary_conditions);
   // Skip the gradient
 }
 
@@ -481,7 +479,7 @@ void colvar::h_bond::calc_gradients()
                                          atom_groups[0]->grad_y(1),
                                          atom_groups[0]->grad_z(1),
                                          0.0, 1.0e20,
-                                         cvm::main()->proxy->get_system_boundaries());
+                                         boundary_conditions);
 }
 
 
@@ -495,7 +493,6 @@ template <int flags> inline void colvar::selfcoordnum::selfcoordnum_sequential_l
 {
   size_t const n = group1->size();
   bool *pairlist_elem = pairlist.get();
-  auto const &boundary_conditions = cvm::main()->proxy->get_system_boundaries();
 
   for (size_t i = 0; i < n - 1; i++) {
 
diff --git a/src/colvarcomp_distances.cpp b/src/colvarcomp_distances.cpp
index fb6e55321..0c0c11631 100644
--- a/src/colvarcomp_distances.cpp
+++ b/src/colvarcomp_distances.cpp
@@ -54,8 +54,7 @@ void colvar::distance::calc_value()
   if (!is_enabled(f_cvc_pbc_minimum_image)) {
     dist_v = group2->center_of_mass() - group1->center_of_mass();
   } else {
-    dist_v = cvmodule->position_distance(group1->center_of_mass(),
-                                    group2->center_of_mass());
+    dist_v = boundary_conditions.position_distance(group1->center_of_mass(), group2->center_of_mass());
   }
   x.real_value = dist_v.norm();
 }
@@ -101,8 +100,8 @@ void colvar::distance_vec::calc_value()
   if (!is_enabled(f_cvc_pbc_minimum_image)) {
     x.rvector_value = group2->center_of_mass() - group1->center_of_mass();
   } else {
-    x.rvector_value = cvmodule->position_distance(group1->center_of_mass(),
-                                             group2->center_of_mass());
+    x.rvector_value =
+        boundary_conditions.position_distance(group1->center_of_mass(), group2->center_of_mass());
   }
 }
 
@@ -127,7 +126,7 @@ void colvar::distance_vec::apply_force(colvarvalue const &force)
 cvm::real colvar::distance_vec::dist2(colvarvalue const &x1, colvarvalue const &x2) const
 {
   if (is_enabled(f_cvc_pbc_minimum_image)) {
-    return (cvmodule->position_distance(x1.rvector_value, x2.rvector_value)).norm2();
+    return (boundary_conditions.position_distance(x1.rvector_value, x2.rvector_value)).norm2();
   }
   return (x2.rvector_value - x1.rvector_value).norm2();
 }
@@ -136,7 +135,7 @@ cvm::real colvar::distance_vec::dist2(colvarvalue const &x1, colvarvalue const &
 colvarvalue colvar::distance_vec::dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const
 {
   if (is_enabled(f_cvc_pbc_minimum_image)) {
-    return 2.0 * cvmodule->position_distance(x2.rvector_value, x1.rvector_value);
+    return 2.0 * boundary_conditions.position_distance(x2.rvector_value, x1.rvector_value);
   }
   return 2.0 * (x2.rvector_value - x1.rvector_value);
 }
@@ -203,7 +202,7 @@ void colvar::distance_z::calc_value()
     if (!is_enabled(f_cvc_pbc_minimum_image)) {
       dist_v = M - R1;
     } else {
-      dist_v = cvmodule->position_distance(R1, M);
+      dist_v = boundary_conditions.position_distance(R1, M);
     }
   } else {
     cvm::rvector const R2 = ref2->center_of_mass();
@@ -212,8 +211,8 @@ void colvar::distance_z::calc_value()
       dist_v = M - C;
       axis = R2 - R1;
     } else {
-      dist_v = cvmodule->position_distance(C, M);
-      axis = cvmodule->position_distance(R1, R2);
+      dist_v = boundary_conditions.position_distance(C, M);
+      axis = boundary_conditions.position_distance(R1, R2);
     }
     axis_norm = axis.norm();
     axis = axis.unit();
@@ -275,15 +274,13 @@ void colvar::distance_xy::calc_value()
   if (!is_enabled(f_cvc_pbc_minimum_image)) {
     dist_v = main->center_of_mass() - ref1->center_of_mass();
   } else {
-    dist_v = cvmodule->position_distance(ref1->center_of_mass(),
-                                    main->center_of_mass());
+    dist_v = boundary_conditions.position_distance(ref1->center_of_mass(), main->center_of_mass());
   }
   if (!fixed_axis) {
     if (!is_enabled(f_cvc_pbc_minimum_image)) {
       v12 = ref2->center_of_mass() - ref1->center_of_mass();
     } else {
-      v12 = cvmodule->position_distance(ref1->center_of_mass(),
-                                   ref2->center_of_mass());
+      v12 = boundary_conditions.position_distance(ref1->center_of_mass(), ref2->center_of_mass());
     }
     axis_norm = v12.norm();
     axis = v12.unit();
@@ -311,8 +308,7 @@ void colvar::distance_xy::calc_gradients()
     if (!is_enabled(f_cvc_pbc_minimum_image)) {
       v13 = main->center_of_mass() - ref1->center_of_mass();
     } else {
-      v13 = cvmodule->position_distance(ref1->center_of_mass(),
-                                   main->center_of_mass());
+      v13 = boundary_conditions.position_distance(ref1->center_of_mass(), main->center_of_mass());
     }
     A = (dist_v * axis) / axis_norm;
 
@@ -357,8 +353,8 @@ void colvar::distance_dir::calc_value()
   if (!is_enabled(f_cvc_pbc_minimum_image)) {
     dist_v = group2->center_of_mass() - group1->center_of_mass();
   } else {
-    dist_v = cvmodule->position_distance(group1->center_of_mass(),
-                                    group2->center_of_mass());
+    dist_v =
+        boundary_conditions.position_distance(group1->center_of_mass(), group2->center_of_mass());
   }
   x.rvector_value = dist_v.unit();
 }
@@ -465,7 +461,7 @@ void colvar::distance_inv::calc_value()
                                group2->pos_z(j));      \
       cvm::rvector dv;                                 \
       if (USE_PBC_MINIMUM_IMAGE) {                     \
-        dv = cvmodule->position_distance(pos1, pos2);       \
+        dv = boundary_conditions.position_distance(pos1, pos2);       \
       } else {                                         \
         dv = pos2 - pos1;                              \
       }                                                \
@@ -548,7 +544,7 @@ void colvar::distance_pairs::calc_value()
                                group2->pos_y(i2),                  \
                                group2->pos_z(i2));                 \
       const cvm::rvector dv = USE_PBC_MINIMUM_IMAGE ?              \
-                              cvmodule->position_distance(pos1, pos2) : \
+                              boundary_conditions.position_distance(pos1, pos2) : \
                               pos2 - pos1;                         \
       cvm::real const d = dv.norm();                               \
       x.vector1d_value[i1*group2->size() + i2] = d;                \
@@ -593,7 +589,7 @@ void colvar::distance_pairs::apply_force(colvarvalue const &force)
                                group2->pos_y(i2),                       \
                                group2->pos_z(i2));                      \
       const cvm::rvector dv = USE_PBC_MINIMUM_IMAGE ?                   \
-                              cvmodule->position_distance(pos1, pos2) :      \
+                              boundary_conditions.position_distance(pos1, pos2) :      \
                               pos2 - pos1;                              \
       cvm::real const d = dv.norm();                                    \
       x.vector1d_value[i1*group2->size() + i2] = d;                     \
diff --git a/src/colvarmodule.cpp b/src/colvarmodule.cpp
index 31b0ad6fe..781b11b34 100644
--- a/src/colvarmodule.cpp
+++ b/src/colvarmodule.cpp
@@ -2482,13 +2482,6 @@ void colvarmodule::request_total_force()
 }
 
 
-cvm::rvector colvarmodule::position_distance(cvm::atom_pos const &pos1,
-                                    cvm::atom_pos const &pos2)
-{
-  return proxy->position_distance(pos1, pos2);
-}
-
-
 cvm::real colvarmodule::rand_gaussian(void)
 {
   return proxy->rand_gaussian();
diff --git a/src/colvarmodule.h b/src/colvarmodule.h
index 81e6c75ab..498737779 100644
--- a/src/colvarmodule.h
+++ b/src/colvarmodule.h
@@ -821,11 +821,6 @@ static inline real acos(real const &x)
     return 5;
   }
 
-  /// \brief Get the distance between two atomic positions with pbcs handled
-  /// correctly
-  rvector position_distance(atom_pos const &pos1,
-                                   atom_pos const &pos2);
-
   /// \brief Names of .ndx files that have been loaded
   std::vector<std::string> index_file_names;
 
diff --git a/src/colvars_system.h b/src/colvars_system.h
index f7b4d03b0..7c23da343 100644
--- a/src/colvars_system.h
+++ b/src/colvars_system.h
@@ -25,6 +25,12 @@ class colvarmodule::system_boundary_conditions {
     unsupported     /// Unsupported boundary conditions
   };
 
+  /// Default constructor
+  inline COLVARS_HOST_DEVICE system_boundary_conditions() { reset(); }
+
+  /// Copy constructor
+  system_boundary_conditions(system_boundary_conditions const &) = default;
+
   /// Type of boundary conditions in the current computation
   inline COLVARS_HOST_DEVICE types type() const { return type_; }
 

From 4eef1b1c78e3fcd8044cd7b19d7c1363e6bfb031 Mon Sep 17 00:00:00 2001
From: Giacomo Fiorin <giacomo.fiorin@gmail.com>
Date: Mon, 16 Feb 2026 15:53:31 -0500
Subject: [PATCH 04/15] cleanup: Remove now redundant logic in CVC compute
 functions

---
 src/colvarcomp_angles.cpp    |  24 ++------
 src/colvarcomp_distances.cpp | 103 ++++++++---------------------------
 2 files changed, 28 insertions(+), 99 deletions(-)

diff --git a/src/colvarcomp_angles.cpp b/src/colvarcomp_angles.cpp
index 370a94993..2f32fc239 100644
--- a/src/colvarcomp_angles.cpp
+++ b/src/colvarcomp_angles.cpp
@@ -63,13 +63,9 @@ void colvar::angle::calc_value()
   cvm::atom_pos const g2_pos = group2->center_of_mass();
   cvm::atom_pos const g3_pos = group3->center_of_mass();
 
-  r21  = is_enabled(f_cvc_pbc_minimum_image) ?
-    boundary_conditions.position_distance(g2_pos, g1_pos) :
-    g1_pos - g2_pos;
+  r21 = boundary_conditions.position_distance(g2_pos, g1_pos);
   r21l = r21.norm();
-  r23  = is_enabled(f_cvc_pbc_minimum_image) ?
-    boundary_conditions.position_distance(g2_pos, g3_pos) :
-    g3_pos - g2_pos;
+  r23 = boundary_conditions.position_distance(g2_pos, g3_pos);
   r23l = r23.norm();
 
   cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
@@ -158,9 +154,7 @@ void colvar::dipole_angle::calc_value()
 
   r21 = group1->dipole();
   r21l = r21.norm();
-  r23  = is_enabled(f_cvc_pbc_minimum_image) ?
-    boundary_conditions.position_distance(g2_pos, g3_pos) :
-    g3_pos - g2_pos;
+  r23 = boundary_conditions.position_distance(g2_pos, g3_pos);
   r23l = r23.norm();
 
   cvm::real denom = (r21l*r23l);
@@ -281,15 +275,9 @@ void colvar::dihedral::calc_value()
   cvm::atom_pos const g4_pos = group4->center_of_mass();
 
   // Usual sign convention: r12 = r2 - r1
-  r12 = is_enabled(f_cvc_pbc_minimum_image) ?
-    boundary_conditions.position_distance(g1_pos, g2_pos) :
-    g2_pos - g1_pos;
-  r23 = is_enabled(f_cvc_pbc_minimum_image) ?
-    boundary_conditions.position_distance(g2_pos, g3_pos) :
-    g3_pos - g2_pos;
-  r34 = is_enabled(f_cvc_pbc_minimum_image) ?
-    boundary_conditions.position_distance(g3_pos, g4_pos) :
-    g4_pos - g3_pos;
+  r12 = boundary_conditions.position_distance(g1_pos, g2_pos);
+  r23 = boundary_conditions.position_distance(g2_pos, g3_pos);
+  r34 = boundary_conditions.position_distance(g3_pos, g4_pos);
 
   cvm::rvector const n1 = cvm::rvector::outer(r12, r23);
   cvm::rvector const n2 = cvm::rvector::outer(r23, r34);
diff --git a/src/colvarcomp_distances.cpp b/src/colvarcomp_distances.cpp
index 0c0c11631..2901ad1c2 100644
--- a/src/colvarcomp_distances.cpp
+++ b/src/colvarcomp_distances.cpp
@@ -51,11 +51,7 @@ int colvar::distance::init(std::string const &conf)
 
 void colvar::distance::calc_value()
 {
-  if (!is_enabled(f_cvc_pbc_minimum_image)) {
-    dist_v = group2->center_of_mass() - group1->center_of_mass();
-  } else {
-    dist_v = boundary_conditions.position_distance(group1->center_of_mass(), group2->center_of_mass());
-  }
+  dist_v = boundary_conditions.position_distance(group1->center_of_mass(), group2->center_of_mass());
   x.real_value = dist_v.norm();
 }
 
@@ -97,12 +93,8 @@ colvar::distance_vec::distance_vec()
 
 void colvar::distance_vec::calc_value()
 {
-  if (!is_enabled(f_cvc_pbc_minimum_image)) {
-    x.rvector_value = group2->center_of_mass() - group1->center_of_mass();
-  } else {
-    x.rvector_value =
-        boundary_conditions.position_distance(group1->center_of_mass(), group2->center_of_mass());
-  }
+  x.rvector_value =
+      boundary_conditions.position_distance(group1->center_of_mass(), group2->center_of_mass());
 }
 
 
@@ -125,19 +117,13 @@ void colvar::distance_vec::apply_force(colvarvalue const &force)
 
 cvm::real colvar::distance_vec::dist2(colvarvalue const &x1, colvarvalue const &x2) const
 {
-  if (is_enabled(f_cvc_pbc_minimum_image)) {
-    return (boundary_conditions.position_distance(x1.rvector_value, x2.rvector_value)).norm2();
-  }
-  return (x2.rvector_value - x1.rvector_value).norm2();
+  return (boundary_conditions.position_distance(x1.rvector_value, x2.rvector_value)).norm2();
 }
 
 
 colvarvalue colvar::distance_vec::dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const
 {
-  if (is_enabled(f_cvc_pbc_minimum_image)) {
-    return 2.0 * boundary_conditions.position_distance(x2.rvector_value, x1.rvector_value);
-  }
-  return 2.0 * (x2.rvector_value - x1.rvector_value);
+  return 2.0 * boundary_conditions.position_distance(x2.rvector_value, x1.rvector_value);
 }
 
 
@@ -199,21 +185,12 @@ void colvar::distance_z::calc_value()
   cvm::rvector const M = main->center_of_mass();
   cvm::rvector const R1 = ref1->center_of_mass();
   if (fixed_axis) {
-    if (!is_enabled(f_cvc_pbc_minimum_image)) {
-      dist_v = M - R1;
-    } else {
-      dist_v = boundary_conditions.position_distance(R1, M);
-    }
+    dist_v = boundary_conditions.position_distance(R1, M);
   } else {
     cvm::rvector const R2 = ref2->center_of_mass();
     cvm::rvector const C = 0.5 * (R1 + R2);
-    if (!is_enabled(f_cvc_pbc_minimum_image)) {
-      dist_v = M - C;
-      axis = R2 - R1;
-    } else {
-      dist_v = boundary_conditions.position_distance(C, M);
-      axis = boundary_conditions.position_distance(R1, R2);
-    }
+    dist_v = boundary_conditions.position_distance(C, M);
+    axis = boundary_conditions.position_distance(R1, R2);
     axis_norm = axis.norm();
     axis = axis.unit();
   }
@@ -271,17 +248,9 @@ colvar::distance_xy::distance_xy()
 
 void colvar::distance_xy::calc_value()
 {
-  if (!is_enabled(f_cvc_pbc_minimum_image)) {
-    dist_v = main->center_of_mass() - ref1->center_of_mass();
-  } else {
-    dist_v = boundary_conditions.position_distance(ref1->center_of_mass(), main->center_of_mass());
-  }
+  dist_v = boundary_conditions.position_distance(ref1->center_of_mass(), main->center_of_mass());
   if (!fixed_axis) {
-    if (!is_enabled(f_cvc_pbc_minimum_image)) {
-      v12 = ref2->center_of_mass() - ref1->center_of_mass();
-    } else {
-      v12 = boundary_conditions.position_distance(ref1->center_of_mass(), ref2->center_of_mass());
-    }
+    v12 = boundary_conditions.position_distance(ref1->center_of_mass(), ref2->center_of_mass());
     axis_norm = v12.norm();
     axis = v12.unit();
   }
@@ -305,11 +274,7 @@ void colvar::distance_xy::calc_gradients()
     ref1->set_weighted_gradient(-1.0 * x_inv * dist_v_ortho);
     main->set_weighted_gradient(       x_inv * dist_v_ortho);
   } else {
-    if (!is_enabled(f_cvc_pbc_minimum_image)) {
-      v13 = main->center_of_mass() - ref1->center_of_mass();
-    } else {
-      v13 = boundary_conditions.position_distance(ref1->center_of_mass(), main->center_of_mass());
-    }
+    v13 = boundary_conditions.position_distance(ref1->center_of_mass(), main->center_of_mass());
     A = (dist_v * axis) / axis_norm;
 
     ref1->set_weighted_gradient( (A - 1.0) * x_inv * dist_v_ortho);
@@ -350,12 +315,8 @@ colvar::distance_dir::distance_dir()
 
 void colvar::distance_dir::calc_value()
 {
-  if (!is_enabled(f_cvc_pbc_minimum_image)) {
-    dist_v = group2->center_of_mass() - group1->center_of_mass();
-  } else {
-    dist_v =
-        boundary_conditions.position_distance(group1->center_of_mass(), group2->center_of_mass());
-  }
+  dist_v =
+      boundary_conditions.position_distance(group1->center_of_mass(), group2->center_of_mass());
   x.rvector_value = dist_v.unit();
 }
 
@@ -448,7 +409,7 @@ int colvar::distance_inv::init(std::string const &conf)
 
 void colvar::distance_inv::calc_value()
 {
-#define CALL_KERNEL(USE_PBC_MINIMUM_IMAGE) do {        \
+#define CALL_KERNEL() do {                             \
   const int factor = -1*(exponent/2);                  \
   for (size_t i = 0; i < group1->size(); ++i) {        \
     const cvm::atom_pos pos1(group1->pos_x(i),         \
@@ -460,11 +421,7 @@ void colvar::distance_inv::calc_value()
                                group2->pos_y(j),       \
                                group2->pos_z(j));      \
       cvm::rvector dv;                                 \
-      if (USE_PBC_MINIMUM_IMAGE) {                     \
-        dv = boundary_conditions.position_distance(pos1, pos2);       \
-      } else {                                         \
-        dv = pos2 - pos1;                              \
-      }                                                \
+      dv = boundary_conditions.position_distance(pos1, pos2); \
       cvm::real const d2 = dv.norm2();                                      \
       cvm::real const dinv = cvm::integer_power(d2, factor);                \
       x.real_value += dinv;                                                 \
@@ -480,11 +437,7 @@ void colvar::distance_inv::calc_value()
   }                            \
 } while (0);
   x.real_value = 0.0;
-  if (!is_enabled(f_cvc_pbc_minimum_image)) {
-    CALL_KERNEL(false);
-  } else {
-    CALL_KERNEL(true);
-  }
+  CALL_KERNEL();
 
   x.real_value *= 1.0 / cvm::real(group1->size() * group2->size());
   x.real_value = cvm::pow(x.real_value, -1.0/cvm::real(exponent));
@@ -533,7 +486,7 @@ int colvar::distance_pairs::init(std::string const &conf)
 void colvar::distance_pairs::calc_value()
 {
   x.vector1d_value.resize(group1->size() * group2->size());
-#define CALL_KERNEL(USE_PBC_MINIMUM_IMAGE) do {                    \
+#define CALL_KERNEL() do {                                         \
   for (size_t i1 = 0; i1 < group1->size(); ++i1) {                 \
     const cvm::atom_pos pos1(group1->pos_x(i1),                    \
                              group1->pos_y(i1),                    \
@@ -543,9 +496,7 @@ void colvar::distance_pairs::calc_value()
       const cvm::atom_pos pos2(group2->pos_x(i2),                  \
                                group2->pos_y(i2),                  \
                                group2->pos_z(i2));                 \
-      const cvm::rvector dv = USE_PBC_MINIMUM_IMAGE ?              \
-                              boundary_conditions.position_distance(pos1, pos2) : \
-                              pos2 - pos1;                         \
+      const cvm::rvector dv = boundary_conditions.position_distance(pos1, pos2); \
       cvm::real const d = dv.norm();                               \
       x.vector1d_value[i1*group2->size() + i2] = d;                \
       const cvm::rvector g2 = dv.unit();                           \
@@ -559,11 +510,7 @@ void colvar::distance_pairs::calc_value()
     group1->grad_z(i1) += g1.z;*/                                    \
   }                                                                \
 } while (0);
-  if (!is_enabled(f_cvc_pbc_minimum_image)) {
-    CALL_KERNEL(false);
-  } else {
-    CALL_KERNEL(true);
-  }
+  CALL_KERNEL();
 #undef CALL_KERNEL
 }
 
@@ -576,7 +523,7 @@ void colvar::distance_pairs::calc_gradients()
 
 void colvar::distance_pairs::apply_force(colvarvalue const &force)
 {
-#define CALL_KERNEL(USE_PBC_MINIMUM_IMAGE) do {                         \
+#define CALL_KERNEL() do {                                              \
   auto group1_force_obj = group1->get_group_force_object();             \
   auto group2_force_obj = group2->get_group_force_object();             \
   for (size_t i1 = 0; i1 < group1->size(); i1++) {                      \
@@ -588,9 +535,7 @@ void colvar::distance_pairs::apply_force(colvarvalue const &force)
       const cvm::atom_pos pos2(group2->pos_x(i2),                       \
                                group2->pos_y(i2),                       \
                                group2->pos_z(i2));                      \
-      const cvm::rvector dv = USE_PBC_MINIMUM_IMAGE ?                   \
-                              boundary_conditions.position_distance(pos1, pos2) :      \
-                              pos2 - pos1;                              \
+      const cvm::rvector dv = boundary_conditions.position_distance(pos1, pos2); \
       cvm::real const d = dv.norm();                                    \
       x.vector1d_value[i1*group2->size() + i2] = d;                     \
       const cvm::rvector f2 = force[i1*group2->size() + i2] * dv.unit();\
@@ -600,11 +545,7 @@ void colvar::distance_pairs::apply_force(colvarvalue const &force)
     group1_force_obj.add_atom_force(i1, f1);                            \
   }                                                                     \
 } while (0);
-  if (!is_enabled(f_cvc_pbc_minimum_image)) {
-    CALL_KERNEL(false);
-  } else {
-    CALL_KERNEL(true);
-  }
+  CALL_KERNEL();
 #undef CALL_KERNEL
 }
 

From fb40f7b6f95b847ebac9dce004ae1129493b444c Mon Sep 17 00:00:00 2001
From: Giacomo Fiorin <giacomo.fiorin@gmail.com>
Date: Thu, 12 Feb 2026 17:58:34 -0500
Subject: [PATCH 05/15] refactor: Get unit cell from MD engines (also fix
 LAMMPS and support GROMACS)

The LAMMPS interface was not importing tilt factors before this fix
---
 .../colvars/colvarproxygromacs.cpp.patch      | 45 +++++++------
 .../colvars/colvarsforceprovider.cpp.patch    | 64 ++++++++++++++++---
 lammps/src/COLVARS/colvarproxy_lammps.cpp     | 24 ++-----
 namd/src/colvarproxy_namd.C                   | 27 ++------
 vmd/src/colvarproxy_vmd.C                     | 40 ++++--------
 5 files changed, 98 insertions(+), 102 deletions(-)

diff --git a/gromacs/src/applied_forces/colvars/colvarproxygromacs.cpp.patch b/gromacs/src/applied_forces/colvars/colvarproxygromacs.cpp.patch
index ed58908ca..ac936cdf9 100644
--- a/gromacs/src/applied_forces/colvars/colvarproxygromacs.cpp.patch
+++ b/gromacs/src/applied_forces/colvars/colvarproxygromacs.cpp.patch
@@ -1,43 +1,40 @@
 diff --git a/src/gromacs/applied_forces/colvars/colvarproxygromacs.cpp b/src/gromacs/applied_forces/colvars/colvarproxygromacs.cpp
-index 76b350c611..37e5b53647 100644
+index 76b350c611..9272e9b365 100644
 --- a/src/gromacs/applied_forces/colvars/colvarproxygromacs.cpp
 +++ b/src/gromacs/applied_forces/colvars/colvarproxygromacs.cpp
 @@ -79,9 +79,6 @@ ColvarProxyGromacs::ColvarProxyGromacs(const std::string& colvarsConfigString,
      // Retrieve masses and charges from input file
      updated_masses_ = updated_charges_ = true;
- 
+
 -    // User-scripted forces are not available in GROMACS
 -    have_scripts = false;
 -
      angstrom_value_ = 0.1;
- 
+
      // From Gnu units
-index 76b350c611..e7cbe0dc28 100644
---- a/src/gromacs/applied_forces/colvars/colvarproxygromacs.cpp
-+++ b/src/gromacs/applied_forces/colvars/colvarproxygromacs.cpp
-@@ -117,23 +117,23 @@ ColvarProxyGromacs::ColvarProxyGromacs(const std::string& colvarsConfigString,
+@@ -117,23 +114,23 @@ ColvarProxyGromacs::ColvarProxyGromacs(const std::string& colvarsConfigString,
              input_streams_[inputName] = new std::istringstream(content);
          }
- 
+
 -        colvars = new colvarmodule(this);
 -        cvm::log(cvm::line_marker);
 -        cvm::log("Start colvars Initialization.");
 +        cvmodule = new colvarmodule(this);
 +        cvmodule->log(cvmodule->line_marker);
 +        cvmodule->log("Start colvars Initialization.");
- 
- 
+
+
 -        colvars->cite_feature("GROMACS engine");
 -        colvars->cite_feature("Colvars-GROMACS interface");
 +        cvmodule->cite_feature("GROMACS engine");
 +        cvmodule->cite_feature("Colvars-GROMACS interface");
- 
+
          if (cvm::debug())
          {
 -            cvm::log("Initializing the colvars proxy object.\n");
 +            cvmodule->log("Initializing the colvars proxy object.\n");
          }
- 
+
          int errorCode = colvarproxy::setup();
 -        errorCode |= colvars->read_config_string(colvarsConfigString);
 -        errorCode |= colvars->update_engine_parameters();
@@ -45,23 +42,23 @@ index 76b350c611..e7cbe0dc28 100644
 +        errorCode |= cvmodule->read_config_string(colvarsConfigString);
 +        errorCode |= cvmodule->update_engine_parameters();
 +        errorCode |= cvmodule->setup_input();
- 
+
          if (errorCode != COLVARS_OK)
          {
-@@ -141,10 +141,10 @@ ColvarProxyGromacs::ColvarProxyGromacs(const std::string& colvarsConfigString,
+@@ -141,10 +138,10 @@ ColvarProxyGromacs::ColvarProxyGromacs(const std::string& colvarsConfigString,
          }
- 
+
          // Citation Reporter
 -        cvm::log(std::string("\n") + colvars->feature_report(0) + std::string("\n"));
 +        cvmodule->log(std::string("\n") + cvmodule->feature_report(0) + std::string("\n"));
- 
+
          // TODO get initial step number from MDModules
 -        // colvars->set_initial_step(static_cast<cvm::step_number>(0L));
 +        // cvmodule->set_initial_step(static_cast<cvm::step_number>(0L));
      }
  }
- 
-@@ -196,7 +196,7 @@ int ColvarProxyGromacs::set_unit_system(std::string const& unitsIn, bool /*colva
+
+@@ -196,7 +193,7 @@ int ColvarProxyGromacs::set_unit_system(std::string const& unitsIn, bool /*colva
  {
      if (unitsIn != "gromacs")
      {
@@ -70,17 +67,19 @@ index 76b350c611..e7cbe0dc28 100644
                  "Specified unit system \"" + unitsIn
                  + "\" is unsupported in Gromacs. Supported units are \"gromacs\" (nm, kJ/mol).\n");
          return COLVARS_ERROR;
-@@ -218,7 +218,7 @@ int ColvarProxyGromacs::check_atom_id(int atomNumber)
+@@ -218,8 +215,8 @@ int ColvarProxyGromacs::check_atom_id(int atomNumber)
      }
      if ((aid < 0) || (aid >= gmxAtoms_.nr))
      {
 -        cvm::error("Error: invalid atom number specified, " + cvm::to_str(atomNumber) + "\n",
+-                   COLVARS_INPUT_ERROR);
 +        cvmodule->error("Error: invalid atom number specified, " + cvm::to_str(atomNumber) + "\n",
-                    COLVARS_INPUT_ERROR);
++                        COLVARS_INPUT_ERROR);
          return COLVARS_INPUT_ERROR;
      }
-@@ -273,10 +273,10 @@ void ColvarProxyGromacs::updateAtomProperties(int index)
- 
+
+@@ -273,10 +270,10 @@ void ColvarProxyGromacs::updateAtomProperties(int index)
+
  ColvarProxyGromacs::~ColvarProxyGromacs()
  {
 -    if (colvars != nullptr)
@@ -92,4 +91,4 @@ index 76b350c611..e7cbe0dc28 100644
 +        cvmodule = nullptr;
      }
  }
- 
+
diff --git a/gromacs/src/applied_forces/colvars/colvarsforceprovider.cpp.patch b/gromacs/src/applied_forces/colvars/colvarsforceprovider.cpp.patch
index 85ddb97fe..10194e024 100644
--- a/gromacs/src/applied_forces/colvars/colvarsforceprovider.cpp.patch
+++ b/gromacs/src/applied_forces/colvars/colvarsforceprovider.cpp.patch
@@ -1,5 +1,5 @@
 diff --git a/src/gromacs/applied_forces/colvars/colvarsforceprovider.cpp b/src/gromacs/applied_forces/colvars/colvarsforceprovider.cpp
-index 5860eed790..8a3232131c 100644
+index 5860eed790..1eee50f51a 100644
 --- a/src/gromacs/applied_forces/colvars/colvarsforceprovider.cpp
 +++ b/src/gromacs/applied_forces/colvars/colvarsforceprovider.cpp
 @@ -201,7 +201,7 @@ ColvarsForceProvider::ColvarsForceProvider(const std::string& colvarsConfigStrin
@@ -46,7 +46,7 @@ index 5860eed790..8a3232131c 100644
 +        cvmodule->log("positions = " + cvm::to_str(atoms_positions) + "\n");
 +        cvmodule->log("total_forces = " + cvm::to_str(atoms_total_forces) + "\n");
 +        cvmodule->log("atoms_new_colvar_forces = " + cvm::to_str(atoms_new_colvar_forces) + "\n");
-+        cvmodule->log(cvmodule->line_marker);
++        cvmodule->log(cvm::line_marker);
          log("Done initializing the colvars proxy object.\n");
      }
 
@@ -54,17 +54,62 @@ index 5860eed790..8a3232131c 100644
      {
 -        cvm::log(cvm::line_marker);
 -        cvm::log("End colvars Initialization.\n\n");
-+        cvmodule->log(cvmodule->line_marker);
++        cvmodule->log(cvm::line_marker);
 +        cvmodule->log("End colvars Initialization.\n\n");
      }
  }
 
-@@ -347,8 +347,10 @@ void ColvarsForceProvider::calculateForces(const ForceProviderInput& forceProvid
-     const rvec* xPointer = &(x.data()->as_vec());
-     const auto& box      = forceProviderInput.box_;
+@@ -336,19 +336,52 @@ ColvarsForceProvider::~ColvarsForceProvider()
+ void ColvarsForceProvider::calculateForces(const ForceProviderInput& forceProviderInput,
+                                            ForceProviderOutput*      forceProviderOutput)
+ {
++    auto const& box = forceProviderInput.box_;
 
--    colvars->it = forceProviderInput.step_;
+-    // Construct t_pbc struct
+-    set_pbc(&gmxPbc_, pbcType_, forceProviderInput.box_);
++    // Set PBC struct for use in ColvarProxyGromacs::position_distance()
++    set_pbc(&gmxPbc_, pbcType_, box);
++
++    // Share PBCs with the library as well
++    switch (pbcType_)
++    {
++        case PbcType::Unset:
++            GMX_RELEASE_ASSERT(false, "PBC type not set when being used by ColvarsForceProvider.");
++            break;
++        case PbcType::XY:
++            boundaries_.set_boundaries(true,
++                                       true,
++                                       false,
++                                       cvm::rvector{ box[0][0], box[0][1], box[0][2] },
++                                       cvm::rvector{ box[1][0], box[1][1], box[1][2] },
++                                       cvm::rvector{ 0.0, 0.0, 0.0 });
++            break;
++        case PbcType::Xyz:
++            boundaries_.set_boundaries(true,
++                                       true,
++                                       true,
++                                       cvm::rvector{ box[0][0], box[0][1], box[0][2] },
++                                       cvm::rvector{ box[1][0], box[1][1], box[1][2] },
++                                       cvm::rvector{ box[2][0], box[2][1], box[2][2] });
++            break;
++        case PbcType::Screw:
++            boundaries_.reset();
++            boundaries_.set_type(cvm::system_boundary_conditions::types::unsupported);
++            break;
++        case PbcType::No: boundaries_.reset(); break;
++        default:
++            GMX_RELEASE_ASSERT(false, "Invalid pbcType in ColvarsForceProvider::calculateForces");
++    }
+
+     const MpiComm& mpiComm = forceProviderInput.mpiComm_;
+     // Local atom coords
+     const gmx::ArrayRef<const gmx::RVec> x = forceProviderInput.x_;
+     // Local atom coords (coerced into into old gmx type)
+     const rvec* xPointer = &(x.data()->as_vec());
+-    const auto& box      = forceProviderInput.box_;
 -
+-    colvars->it = forceProviderInput.step_;
+
 +    if (cvmodule) // Only one thread has a pointer to the module
 +    {
 +        cvmodule->it = forceProviderInput.step_; // 'it' is not a static class member anymore
@@ -72,7 +117,7 @@ index 5860eed790..8a3232131c 100644
 
      // Eventually there needs to be an interface to update local data upon neighbor search
      // We could check if by chance all atoms are in one node, and skip communication
-@@ -389,9 +391,9 @@ void ColvarsForceProvider::calculateForces(const ForceProviderInput& forceProvid
+@@ -389,9 +422,9 @@ void ColvarsForceProvider::calculateForces(const ForceProviderInput& forceProvid
 
          biasEnergy = 0.0;
          // Call the collective variable module to fill atoms_new_colvar_forces
@@ -84,7 +129,7 @@ index 5860eed790..8a3232131c 100644
          }
 
          // Copy the forces to C array for broadcasting
-@@ -462,7 +464,7 @@ void ColvarsForceProvider::add_energy(cvm::real energy)
+@@ -462,7 +495,7 @@ void ColvarsForceProvider::add_energy(cvm::real energy)
  void ColvarsForceProvider::writeCheckpointData(MDModulesWriteCheckpointData checkpointWriting,
                                                 std::string_view             moduleName)
  {
@@ -92,3 +137,4 @@ index 5860eed790..8a3232131c 100644
 +    cvmodule->write_state_buffer(stateToCheckpoint_.colvarStateFile_);
      stateToCheckpoint_.writeState(checkpointWriting.builder_, moduleName);
  }
+
diff --git a/lammps/src/COLVARS/colvarproxy_lammps.cpp b/lammps/src/COLVARS/colvarproxy_lammps.cpp
index 8cc40de7e..9c00d6bc1 100644
--- a/lammps/src/COLVARS/colvarproxy_lammps.cpp
+++ b/lammps/src/COLVARS/colvarproxy_lammps.cpp
@@ -133,24 +133,12 @@ double colvarproxy_lammps::compute()
   }
   previous_step = _lmp->update->ntimestep;
 
-  unit_cell_x.set(_lmp->domain->xprd, 0.0, 0.0);
-  unit_cell_y.set(0.0, _lmp->domain->yprd, 0.0);
-  unit_cell_z.set(0.0, 0.0, _lmp->domain->zprd);
-
-  if (_lmp->domain->xperiodic == 0 && _lmp->domain->yperiodic == 0 &&
-      _lmp->domain->zperiodic == 0) {
-    boundaries_type = boundaries_non_periodic;
-    reset_pbc_lattice();
-  } else if ((_lmp->domain->nonperiodic == 0) &&
-             (_lmp->domain->dimension == 3) &&
-             (_lmp->domain->triclinic == 0)) {
-    // Orthogonal unit cell
-    boundaries_type = boundaries_pbc_ortho;
-    colvarproxy_system::update_pbc_lattice();
-    // It is safer to let LAMMPS deal with high-tilt triclinic boxes
-  } else {
-    boundaries_type = boundaries_unsupported;
-  }
+  boundaries_.set_boundaries(_lmp->domain->xperiodic,
+                             _lmp->domain->yperiodic,
+                             _lmp->domain->zperiodic,
+                             cvm::rvector{_lmp->domain->xprd, 0.0, 0.0},
+                             cvm::rvector{_lmp->domain->xy, _lmp->domain->yprd, 0.0},
+                             cvm::rvector{_lmp->domain->xz, _lmp->domain->yz, _lmp->domain->zprd});
 
   if (cvmodule->debug()) {
     cvmodule->log(std::string(cvm::line_marker) +
diff --git a/namd/src/colvarproxy_namd.C b/namd/src/colvarproxy_namd.C
index cdda407ea..584f47eae 100644
--- a/namd/src/colvarproxy_namd.C
+++ b/namd/src/colvarproxy_namd.C
@@ -388,29 +388,10 @@ void colvarproxy_namd::calculate()
   if (accelMDOn) update_accelMD_info();
 
   auto *lattice = globalmaster->get_lattice();
-
-  {
-    Vector const a = lattice->a();
-    Vector const b = lattice->b();
-    Vector const c = lattice->c();
-    unit_cell_x.set(a.x, a.y, a.z);
-    unit_cell_y.set(b.x, b.y, b.z);
-    unit_cell_z.set(c.x, c.y, c.z);
-  }
-
-  if (!lattice->a_p() && !lattice->b_p() && !lattice->c_p()) {
-    boundaries_type = boundaries_non_periodic;
-    reset_pbc_lattice();
-  } else if (lattice->a_p() && lattice->b_p() && lattice->c_p()) {
-    if (lattice->orthogonal()) {
-      boundaries_type = boundaries_pbc_ortho;
-    } else {
-      boundaries_type = boundaries_pbc_triclinic;
-    }
-    colvarproxy_system::update_pbc_lattice();
-  } else {
-    boundaries_type = boundaries_unsupported;
-  }
+  boundaries_.set_boundaries(lattice->a_p(), lattice->b_p(), lattice->c_p(),
+                             cvm::rvector{lattice->a().x, lattice->a().y, lattice->a().z},
+                             cvm::rvector{lattice->b().x, lattice->b().y, lattice->b().z},
+                             cvm::rvector{lattice->c().x, lattice->c().y, lattice->c().z});
 
   if (cvm::debug()) {
     cvmodule->log(std::string(cvm::line_marker)+
diff --git a/vmd/src/colvarproxy_vmd.C b/vmd/src/colvarproxy_vmd.C
index 5a0d4dd79..87cf05c99 100644
--- a/vmd/src/colvarproxy_vmd.C
+++ b/vmd/src/colvarproxy_vmd.C
@@ -227,36 +227,18 @@ int colvarproxy_vmd::update_input()
                                        angstrom_to_internal(vmdpos[atoms_ids[i]*3+2]));
   }
 
-
   Timestep const *ts = vmdmol->get_frame(vmdmol_frame);
-  {
-    // Get lattice vectors
-    float A[3];
-    float B[3];
-    float C[3];
-    ts->get_transform_vectors(A, B, C);
-    unit_cell_x.set(angstrom_to_internal(A[0]), angstrom_to_internal(A[1]), angstrom_to_internal(A[2]));
-    unit_cell_y.set(angstrom_to_internal(B[0]), angstrom_to_internal(B[1]), angstrom_to_internal(B[2]));
-    unit_cell_z.set(angstrom_to_internal(C[0]), angstrom_to_internal(C[1]), angstrom_to_internal(C[2]));
-  }
-
-  if (ts->a_length + ts->b_length + ts->c_length < 1.0e-6) {
-    boundaries_type = boundaries_non_periodic;
-    reset_pbc_lattice();
-  } else if ((ts->a_length > 1.0e-6) &&
-             (ts->b_length > 1.0e-6) &&
-             (ts->c_length > 1.0e-6)) {
-    if (((ts->alpha-90.0)*(ts->alpha-90.0)) +
-        ((ts->beta-90.0)*(ts->beta-90.0)) +
-        ((ts->gamma-90.0)*(ts->gamma-90.0)) < 1.0e-6) {
-      boundaries_type = boundaries_pbc_ortho;
-    } else {
-      boundaries_type = boundaries_pbc_triclinic;
-    }
-    colvarproxy_system::update_pbc_lattice();
-  } else {
-    boundaries_type = boundaries_unsupported;
-  }
+
+  float A[3];
+  float B[3];
+  float C[3];
+  ts->get_transform_vectors(A, B, C);
+  boundaries_.set_boundaries((ts->a_length > 0.0),
+                             (ts->b_length > 0.0),
+                             (ts->c_length > 0.0),
+                             cvm::rvector{A[0], A[1], A[2]},
+                             cvm::rvector{B[0], B[1], B[2]},
+                             cvm::rvector{C[0], C[1], C[2]});
 
   return error_code;
 }

From e5766c247d32fbc941c211ea3d4d871217f4c0ea Mon Sep 17 00:00:00 2001
From: Giacomo Fiorin <giacomo.fiorin@gmail.com>
Date: Mon, 16 Feb 2026 13:53:00 -0500
Subject: [PATCH 06/15] build: Add script to create patches to GROMACS
 interface files

---
 gromacs/update_patches_from_git.sh | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)
 create mode 100644 gromacs/update_patches_from_git.sh

diff --git a/gromacs/update_patches_from_git.sh b/gromacs/update_patches_from_git.sh
new file mode 100644
index 000000000..484e97f39
--- /dev/null
+++ b/gromacs/update_patches_from_git.sh
@@ -0,0 +1,23 @@
+#!/bin/bash
+
+# Run this script with a modified GROMACS worktree as argument to import the patches
+
+GROMACS_SOURCE=${1}
+if [ -n "${GROMACS_SOURCE}" ] ; then
+    GROMACS_SOURCE=$(git -C ${GROMACS_SOURCE} rev-parse --show-toplevel)
+fi
+
+if [ $? != 0 ] || [ ! -s ${GROMACS_SOURCE}/src/gromacs/version.h.cmakein ] ; then
+    echo >&2
+    echo "Usage: $0 <path_to_GROMACS_tree>" >& 2
+    exit 1
+fi
+
+COLVARS_SOURCE=$(dirname $(realpath $0))
+COLVARS_SOURCE=$(git -C ${COLVARS_SOURCE} rev-parse --show-toplevel)
+
+modified_files=($(git -C ${GROMACS_SOURCE} ls-files -m src/gromacs/applied_forces/colvars))
+
+for file in ${modified_files[@]} ; do
+    git -C ${GROMACS_SOURCE} diff ${file} > ${COLVARS_SOURCE}/gromacs/src/${file#src\/gromacs\/}.patch
+done

From 220a361554382b0ba59fa25e78c640319bdc50b4 Mon Sep 17 00:00:00 2001
From: Giacomo Fiorin <giacomo.fiorin@gmail.com>
Date: Mon, 16 Feb 2026 15:24:38 -0500
Subject: [PATCH 07/15] doc: Document handling of boundary conditions

---
 doc/colvars-refman-main.tex | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/doc/colvars-refman-main.tex b/doc/colvars-refman-main.tex
index 9194535e8..946939bfa 100644
--- a/doc/colvars-refman-main.tex
+++ b/doc/colvars-refman-main.tex
@@ -1495,15 +1495,16 @@
 
 \cvsubsec{Treatment of periodic boundary conditions}{sec:cvc_pbcs}
 
-In all colvar components described below, the following rules apply concerning periodic boundary
-conditions (PBCs):
+When copying over atomic coordinates from \MDENGINE, Colvars also copies its boundary conditions and unit cell vectors.  Periodic or non-periodic boundary conditions are supported for each direction independently, thus allowing for mixed-periodicity conditions if \MDENGINE{} allows them.  When using triclinic unit cells, tilt factors up to 1/2 are supported by looping over one additional layer of adjacent periodic images around the simulation cell to find the minimum image.  The use of tilt factors beyond 1/2 is \emph{not recommended}, and may lead to slower performance\cvlammpsonly{ (LAMMPS will also print warnings to that sense)}.
+
+For all colvar components implemented by Colvars internal code, the following rules apply:
 \begin{enumerate}
 \item Distance vectors between two coordinates
   $\mathbf{d}_{i,j} = \left(\mathbf{x}_1 - \mathbf{x}_2\right)$, are calculated following the
   minimum-image convention by default, unless \refkey{forceNoPBC}{colvar|cvc|forceNoPBC} is enabled.
   ($\mathbf{x}_1$ and $\mathbf{x}_2$ may be either individual atomic coordinates, or centers of mass
   of two groups.)
-\item For all other functions of individual atomic coordinates,
+\item For all other functions of multiple atomic coordinates,
   $f\left(\mathbf{x}_1, \mathbf{x}_2, \ldots\right)$, it is assumed that all atoms that are part of
   the same group are in the same periodic unit cell (see \ref{sec:colvar_atom_groups_wrapping}).
 \end{enumerate}
@@ -8223,6 +8224,9 @@
   The legacy keyword \texttt{disableForces} for atom groups is now deprecated and will be discontinued in a future release.
   Atom groups now have an automated way to save computation if forces are not used, and enabling this option otherwise would lead to incorrect behavior.
 
+\item \textbf{Colvars version 2026-XX-XX or later.}\\
+  Periodic boundary conditions are handled internally, rather than using code by \MDENGINE{} to compute the minimum-image convention (see \ref{sec:cvc_pbcs}).
+
 \end{itemize}
 
 

From 05cc6db20976d1d2ff68d41d4696e93b550f83ad Mon Sep 17 00:00:00 2001
From: Giacomo Fiorin <giacomo.fiorin@gmail.com>
Date: Mon, 16 Feb 2026 15:23:58 -0500
Subject: [PATCH 08/15] feat: Support forceNoPBC for all variables

---
 src/colvarcomp.cpp           | 8 +++++++-
 src/colvarcomp_coordnums.cpp | 4 ----
 2 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/src/colvarcomp.cpp b/src/colvarcomp.cpp
index e0f6b69e7..75edadfc8 100644
--- a/src/colvarcomp.cpp
+++ b/src/colvarcomp.cpp
@@ -456,7 +456,13 @@ int colvar::cvc::set_param(std::string const &param_name,
 
 void colvar::cvc::read_data()
 {
-  boundary_conditions = cvm::main()->proxy->get_system_boundaries();
+  if (!is_enabled(f_cvc_pbc_minimum_image)) {
+    // Copy boundary conditions from the proxy
+    boundary_conditions = cvm::main()->proxy->get_system_boundaries();
+  } else {
+    // Set as non-periodic boundary conditions (default) for this CVC
+    boundary_conditions.reset();
+  }
 
   if (is_enabled(f_cvc_explicit_atom_groups)) {
     for (auto agi = atom_groups.begin(); agi != atom_groups.end(); agi++) {
diff --git a/src/colvarcomp_coordnums.cpp b/src/colvarcomp_coordnums.cpp
index b40418fd6..53a677b8d 100644
--- a/src/colvarcomp_coordnums.cpp
+++ b/src/colvarcomp_coordnums.cpp
@@ -135,10 +135,6 @@ int colvar::coordnum::init(std::string const &conf)
                              COLVARS_INPUT_ERROR);
   }
 
-  if (!is_enabled(f_cvc_pbc_minimum_image)) {
-    cvmodule->log("Warning: only minimum-image distances are used by this variable.\n");
-  }
-
   if (function_type() != "groupCoord") {
     // All coordNum variables may benefit from a pairlist, except groupCoord
     get_keyval(conf, "tolerance", tolerance, tolerance);

From 88b2d0e23f926a3a3c63448407f41a2693c678e9 Mon Sep 17 00:00:00 2001
From: Giacomo Fiorin <giacomo.fiorin@gmail.com>
Date: Mon, 16 Feb 2026 11:18:03 -0500
Subject: [PATCH 09/15] opt: Compute shifts regardless of periodic flags (they
 are zero without PBCs)

---
 src/colvars_system.h | 23 +++++++++--------------
 1 file changed, 9 insertions(+), 14 deletions(-)

diff --git a/src/colvars_system.h b/src/colvars_system.h
index 7c23da343..d93a4c14c 100644
--- a/src/colvars_system.h
+++ b/src/colvars_system.h
@@ -106,7 +106,8 @@ cvm::system_boundary_conditions::set_boundaries(bool periodic_x_in, bool periodi
       reciprocal_cell_x = unit_cell_x/unit_cell_x.norm2();
     }
   } else {
-    unit_cell_x = {0.0, 0.0, 0.0};
+    unit_cell_x.reset();
+    reciprocal_cell_x.reset();
   }
 
   if (periodic_y) {
@@ -117,7 +118,8 @@ cvm::system_boundary_conditions::set_boundaries(bool periodic_x_in, bool periodi
       reciprocal_cell_y = unit_cell_y/unit_cell_y.norm2();
     }
   } else {
-    unit_cell_y = {0.0, 0.0, 0.0};
+    unit_cell_y.reset();
+    reciprocal_cell_y.reset();
   }
 
   if (periodic_z) {
@@ -128,7 +130,8 @@ cvm::system_boundary_conditions::set_boundaries(bool periodic_x_in, bool periodi
       reciprocal_cell_z = unit_cell_z/unit_cell_z.norm2();
     }
   } else {
-    unit_cell_z = {0.0, 0.0, 0.0};
+    unit_cell_z.reset();
+    reciprocal_cell_z.reset();
   }
 
   if (type() == types::pbc_orthogonal && off_diagonal) {
@@ -160,17 +163,9 @@ cvm::rvector cvm::system_boundary_conditions::position_distance(cvm::atom_pos co
   cvm::real const y_shift = ::floor(reciprocal_cell_y * diff + 0.5);
   cvm::real const z_shift = ::floor(reciprocal_cell_z * diff + 0.5);
 
-  if (periodic_x) {
-    diff.x -= x_shift * unit_cell_x.x + y_shift * unit_cell_y.x + z_shift * unit_cell_z.x;
-  }
-
-  if (periodic_y) {
-    diff.y -= x_shift * unit_cell_x.y + y_shift * unit_cell_y.y + z_shift * unit_cell_z.y;
-  }
-
-  if (periodic_z) {
-    diff.z -= x_shift * unit_cell_x.z + y_shift * unit_cell_y.z + z_shift * unit_cell_z.z;
-  }
+  diff.x -= x_shift * unit_cell_x.x + y_shift * unit_cell_y.x + z_shift * unit_cell_z.x;
+  diff.y -= x_shift * unit_cell_x.y + y_shift * unit_cell_y.y + z_shift * unit_cell_z.y;
+  diff.z -= x_shift * unit_cell_x.z + y_shift * unit_cell_y.z + z_shift * unit_cell_z.z;
 
   return diff;
 }

From de312209fa876d4ede077e440479e97a24375c25 Mon Sep 17 00:00:00 2001
From: Giacomo Fiorin <giacomo.fiorin@gmail.com>
Date: Mon, 16 Feb 2026 11:05:07 -0500
Subject: [PATCH 10/15] opt: Remove now obsoleted internalPBC flag

---
 doc/colvars-refman-main.tex  | 12 ---------
 src/colvarcomp_coordnums.cpp | 47 ++++++++++++------------------------
 src/colvarcomp_coordnums.h   |  8 +-----
 3 files changed, 16 insertions(+), 51 deletions(-)

diff --git a/doc/colvars-refman-main.tex b/doc/colvars-refman-main.tex
index 946939bfa..2b6f45c46 100644
--- a/doc/colvars-refman-main.tex
+++ b/doc/colvars-refman-main.tex
@@ -2013,18 +2013,6 @@
     100}{This controls the pairlist feature, dictating how many steps are taken between regenerating pair  lists if the tolerance is greater than 0.
   }
 
-\item %
-  \labelkey{colvar|coordNum|useInternalPBC}%
-  \keydef{useInternalPBC}{%
-    \texttt{coordNum}}{%
-    Use the PBC functions from the Colvars library (as opposed to MD engine)}{%
-    boolean}{%
-    engine-dependent}{%
-    If this keyword is set to \texttt{yes}, Colvars will use internal code to compute mininum-image distances, instead of calling the function used by \MDENGINE{}: enabling this flag may speed up computation.
-    \cvvmdonly{In VMD, the default value of this flag is true.} 
-    \cvlammpsonly{Please note that, unlike the LAMMPS function, the current internal code is not robust against large tilt factors: use this feature only for unit cells near an orthogonal shape.}
-    }
-
 \end{cvcoptions}
 
 This component returns a dimensionless number, which ranges from 0 (all interatomic distances are much larger than the cutoff) to $N_{\mathtt{group1}} \times N_{\mathtt{group2}}$ (all possible pairwise distances are much smaller than the cutoff).
diff --git a/src/colvarcomp_coordnums.cpp b/src/colvarcomp_coordnums.cpp
index 53a677b8d..7fbc13bd4 100644
--- a/src/colvarcomp_coordnums.cpp
+++ b/src/colvarcomp_coordnums.cpp
@@ -21,7 +21,6 @@ colvar::coordnum::coordnum()
   x.type(colvarvalue::type_scalar);
   cvm::real const r0 = cvmodule->proxy->angstrom_to_internal(4.0);
   update_cutoffs({r0, r0, r0});
-  b_use_internal_pbc = cvm::main()->proxy->use_internal_pbc();
   // Boundaries will be set later, when the number of pairs is known
 }
 
@@ -48,8 +47,6 @@ int colvar::coordnum::init(std::string const &conf)
 {
   int error_code = cvc::init(conf);
 
-  get_keyval(conf, "useInternalPBC", b_use_internal_pbc, b_use_internal_pbc);
-
   group1 = parse_group(conf, "group1");
 
   if (!group1) {
@@ -223,19 +220,12 @@ void inline colvar::coordnum::main_loop()
           ((flags & ef_use_pairlist) && (*pairlist_elem || (flags & ef_rebuild_pairlist))) ||
           !(flags & ef_use_pairlist);
 
-      cvm::real const partial = within ?
-        (b_use_internal_pbc ?
-         compute_pair_coordnum<flags | ef_use_internal_pbc>(inv_r0_vec, inv_r0sq_vec, en, ed,
-                                                            x1, y1, z1, x2, y2, z2,
-                                                            gx1, gy1, gz1, gx2, gy2, gz2,
-                                                            tolerance, tolerance_l2_max,
-                                                            boundary_conditions) :
-         compute_pair_coordnum<flags>(inv_r0_vec, inv_r0sq_vec, en, ed,
-                                      x1, y1, z1, x2, y2, z2,
-                                      gx1, gy1, gz1, gx2, gy2, gz2,
-                                      tolerance, tolerance_l2_max,
-                                      boundary_conditions) ) :
-        0.0;
+      cvm::real const partial =
+          within ? compute_pair_coordnum<flags>(inv_r0_vec, inv_r0sq_vec, en, ed,
+                                                x1, y1, z1, x2, y2, z2,
+                                                gx1, gy1, gz1, gx2, gy2, gz2,
+                                                tolerance, tolerance_l2_max, boundary_conditions)
+                 : 0.0;
 
       if ((flags & ef_use_pairlist) && (flags & ef_rebuild_pairlist)) {
         *pairlist_elem = partial > 0.0 ? true : false;
@@ -509,22 +499,15 @@ template <int flags> inline void colvar::selfcoordnum::selfcoordnum_sequential_l
       cvm::real &gz2 = group1->grad_z(j);
 
       bool const within =
-        ((flags & ef_use_pairlist) && (*pairlist_elem || (flags & ef_rebuild_pairlist))) ||
-        !(flags & ef_use_pairlist);
-
-      cvm::real const partial = within ?
-        (b_use_internal_pbc ?
-         compute_pair_coordnum<flags | ef_use_internal_pbc>(inv_r0_vec, inv_r0sq_vec, en, ed,
-                                                            x1, y1, z1, x2, y2, z2,
-                                                            gx1, gy1, gz1, gx2, gy2, gz2,
-                                                            tolerance, tolerance_l2_max,
-                                                            boundary_conditions) :
-         compute_pair_coordnum<flags>(inv_r0_vec, inv_r0sq_vec, en, ed,
-                                      x1, y1, z1, x2, y2, z2,
-                                      gx1, gy1, gz1, gx2, gy2, gz2,
-                                      tolerance, tolerance_l2_max,
-                                      boundary_conditions) ) :
-        0.0;
+          ((flags & ef_use_pairlist) && (*pairlist_elem || (flags & ef_rebuild_pairlist))) ||
+          !(flags & ef_use_pairlist);
+
+      cvm::real const partial =
+          within ? compute_pair_coordnum<flags>(inv_r0_vec, inv_r0sq_vec, en, ed,
+                                                x1, y1, z1, x2, y2, z2,
+                                                gx1, gy1, gz1, gx2, gy2, gz2,
+                                                tolerance, tolerance_l2_max, boundary_conditions)
+                 : 0.0;
 
       if ((flags & ef_use_pairlist) && (flags & ef_rebuild_pairlist)) {
         *pairlist_elem = partial > 0.0 ? true : false;
diff --git a/src/colvarcomp_coordnums.h b/src/colvarcomp_coordnums.h
index e7dac7510..d6e4b0857 100644
--- a/src/colvarcomp_coordnums.h
+++ b/src/colvarcomp_coordnums.h
@@ -32,7 +32,6 @@ class colvar::coordnum : public colvar::cvc {
   enum {
     ef_null = 0,
     ef_gradients = 1,
-    ef_use_internal_pbc = (1 << 8),
     ef_use_pairlist = (1 << 9),
     ef_rebuild_pairlist = (1 << 10)
   };
@@ -96,9 +95,6 @@ class colvar::coordnum : public colvar::cvc {
   /// If true, group2 will be treated as a single atom
   bool b_group2_center_only = false;
 
-  /// Use the PBC functions from the Colvars library (as opposed to MD engine)
-  bool b_use_internal_pbc = false;
-
   /// Tolerance for the pair list
   cvm::real tolerance = 0.0;
 
@@ -251,10 +247,8 @@ inline cvm::real colvar::coordnum::compute_pair_coordnum(cvm::rvector const &inv
 {
   const cvm::atom_pos pos1{a1x, a1y, a1z};
   const cvm::atom_pos pos2{a2x, a2y, a2z};
-  cvm::rvector const diff = (flags & ef_use_internal_pbc)
-                                ? boundary_conditions.position_distance(pos1, pos2)
-                                : cvm::main()->proxy->position_distance(pos1, pos2);
 
+  cvm::rvector const diff = boundary_conditions.position_distance(pos1, pos2);
   cvm::rvector const scal_diff(diff.x * inv_r0_vec.x,
                                diff.y * inv_r0_vec.y,
                                diff.z * inv_r0_vec.z);

From 2b205ba029230676ba2bac3f1f8c2eb63d6282d6 Mon Sep 17 00:00:00 2001
From: Giacomo Fiorin <giacomo.fiorin@gmail.com>
Date: Sat, 14 Feb 2026 17:06:59 -0500
Subject: [PATCH 11/15] fix: Silence warning regarding shadowing of class
 member

---
 src/colvarcomp_coordnums.h | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/colvarcomp_coordnums.h b/src/colvarcomp_coordnums.h
index d6e4b0857..0bca40969 100644
--- a/src/colvarcomp_coordnums.h
+++ b/src/colvarcomp_coordnums.h
@@ -55,7 +55,7 @@ class colvar::coordnum : public colvar::cvc {
                                          cvm::real &g1x, cvm::real &g1y, cvm::real &g1z,
                                          cvm::real &g2x, cvm::real &g2y, cvm::real &g2z,
                                          cvm::real pairlist_tol, cvm::real pairlist_tol_l2_max,
-                                         cvm::system_boundary_conditions const &boundary_conditions);
+                                         cvm::system_boundary_conditions const &bc);
 
   /// Workhorse function
   template <bool use_group1_com, bool use_group2_com, int flags> int compute_coordnum();
@@ -243,12 +243,12 @@ inline cvm::real colvar::coordnum::compute_pair_coordnum(cvm::rvector const &inv
                                                          cvm::real& g2z,
                                                          cvm::real pairlist_tol,
                                                          cvm::real pairlist_tol_l2_max,
-                                                         cvm::system_boundary_conditions const &boundary_conditions)
+                                                         cvm::system_boundary_conditions const &bc)
 {
   const cvm::atom_pos pos1{a1x, a1y, a1z};
   const cvm::atom_pos pos2{a2x, a2y, a2z};
 
-  cvm::rvector const diff = boundary_conditions.position_distance(pos1, pos2);
+  cvm::rvector const diff = bc.position_distance(pos1, pos2);
   cvm::rvector const scal_diff(diff.x * inv_r0_vec.x,
                                diff.y * inv_r0_vec.y,
                                diff.z * inv_r0_vec.z);

From 3bc036afcad1d5252e9f0acab63490f864313bc9 Mon Sep 17 00:00:00 2001
From: Giacomo Fiorin <giacomo.fiorin@gmail.com>
Date: Fri, 13 Feb 2026 16:37:00 -0500
Subject: [PATCH 12/15] chore: Update clang-format config to keep return type
 with function

---
 src/.clang-format | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/src/.clang-format b/src/.clang-format
index 66e78c8b5..e3f508c76 100644
--- a/src/.clang-format
+++ b/src/.clang-format
@@ -1,14 +1,16 @@
 ---
-Language: Cpp
-ColumnLimit: 100
-UseTab: Never
-TabWidth: 2
-IndentWidth: 2
-ObjCBlockIndentWidth: 2
-PenaltyBreakAssignment: 2
 AccessModifierOffset: -2
+BinPackParameters: false
 BreakBeforeBraces: WebKit
+ColumnLimit: 100
 EmptyLineAfterAccessModifier: Leave
-NamespaceIndentation: None
+IndentWidth: 2
+Language: Cpp
 MaxEmptyLinesToKeep: 2
+NamespaceIndentation: None
+ObjCBlockIndentWidth: 2
+PenaltyBreakAssignment: 2
+PenaltyReturnTypeOnItsOwnLine: 10000
+TabWidth: 2
+UseTab: Never
 ...

From 5d5444eb3d6436b8858153d3ca10acf145d52bcf Mon Sep 17 00:00:00 2001
From: Giacomo Fiorin <giacomo.fiorin@gmail.com>
Date: Fri, 13 Feb 2026 16:00:11 -0500
Subject: [PATCH 13/15] fix: Handle mininum-image for triclinic cells with
 moderate tilt in internal code

---
 src/colvars_system.h | 37 +++++++++++++++++++++++++++++++++++++
 1 file changed, 37 insertions(+)

diff --git a/src/colvars_system.h b/src/colvars_system.h
index d93a4c14c..6728f8164 100644
--- a/src/colvars_system.h
+++ b/src/colvars_system.h
@@ -40,6 +40,9 @@ class colvarmodule::system_boundary_conditions {
   /// Compute the distance between two positions
   cvm::rvector position_distance(cvm::atom_pos const &pos1, cvm::atom_pos const &pos2) const;
 
+  /// Compute a shift vector that accounts for tilt factors up to 0.5
+  cvm::rvector get_triclinic_shift(cvm::rvector const &diff) const;
+
   /// Reset to defaults (non-periodic)
   inline COLVARS_HOST_DEVICE void reset() {
     periodic_x = periodic_y = periodic_z = false;
@@ -167,7 +170,41 @@ cvm::rvector cvm::system_boundary_conditions::position_distance(cvm::atom_pos co
   diff.y -= x_shift * unit_cell_x.y + y_shift * unit_cell_y.y + z_shift * unit_cell_z.y;
   diff.z -= x_shift * unit_cell_x.z + y_shift * unit_cell_y.z + z_shift * unit_cell_z.z;
 
+  if (type() != types::pbc_orthogonal) {
+    // Matches both "mixed" and "pbc_triclinic", because reciprocal cell vectors are not used
+    diff += get_triclinic_shift(diff);
+  }
+
   return diff;
 }
 
+
+inline COLVARS_HOST_DEVICE cvm::rvector
+cvm::system_boundary_conditions::get_triclinic_shift(cvm::rvector const &diff) const
+{
+  cvm::real min_dist2 = diff.norm2();
+  cvm::rvector result{0.0, 0.0, 0.0};
+
+  int const nx = periodic_x ? 1 : 0;
+  int const ny = periodic_y ? 1 : 0;
+  int const nz = periodic_z ? 1 : 0;
+
+  // Loop over neighboring cells to find a shorter distance
+  for (int ix = -nx; ix <= nx; ix++) {
+    for (int iy = -ny; iy <= ny; iy++) {
+      for (int iz = -nz; iz <= nz; iz++) {
+        cvm::rvector const shift = ix * unit_cell_x + iy * unit_cell_y + iz * unit_cell_z;
+        cvm::real const this_dist2 = (diff + shift).norm2();
+        if (this_dist2 < min_dist2) {
+          result = shift;
+          min_dist2 = this_dist2;
+        }
+      }
+    }
+  }
+
+  return result;
+}
+
+
 #endif

From edb7b3baa728556b6256ad62e8a55bab2f72dc7f Mon Sep 17 00:00:00 2001
From: Giacomo Fiorin <giacomo.fiorin@gmail.com>
Date: Sun, 15 Feb 2026 20:31:18 -0500
Subject: [PATCH 14/15] chore: Make cvm::debug() a constexpr

(Leftover from a previous version of this branch but useful in general)
---
 src/colvarmodule.h | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/colvarmodule.h b/src/colvarmodule.h
index 498737779..135cf73e0 100644
--- a/src/colvarmodule.h
+++ b/src/colvarmodule.h
@@ -374,9 +374,13 @@ static inline real acos(real const &x)
   std::vector<colvarbias *> *biases_active();
 
   /// \brief Whether debug output should be enabled (compile-time option)
-  static inline bool debug()
+  static inline bool constexpr debug()
   {
+#if (defined(__HIP_DEVICE_COMPILE__)) || (defined(__CUDA_ARCH__))
+    return false;
+#else
     return COLVARS_DEBUG;
+#endif
   }
 
   /// How many objects (variables and biases) are configured yet?

From 815dfd4a962c71429c21083700cce610a3cd32e6 Mon Sep 17 00:00:00 2001
From: Giacomo Fiorin <giacomo.fiorin@gmail.com>
Date: Tue, 16 Jun 2026 17:49:47 -0400
Subject: [PATCH 15/15] chore: Add script useful for cherry-picking

[skip CI]
---
 devel-tools/compare_branch_commits.sh | 32 +++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)
 create mode 100644 devel-tools/compare_branch_commits.sh

diff --git a/devel-tools/compare_branch_commits.sh b/devel-tools/compare_branch_commits.sh
new file mode 100644
index 000000000..6477a4664
--- /dev/null
+++ b/devel-tools/compare_branch_commits.sh
@@ -0,0 +1,32 @@
+#!/usr/bin/env bash
+
+# Iterates over the commits in the first argument (ideally: "master...branch1") and prints out
+# those whose subject lines are not found in the second argument (ideally: "master...branch2")
+
+source_branch="$1"
+target_branch="$2"
+
+mapfile -t target_branch_messages < <(git log --format=%s "${target_branch}")
+mapfile -t target_branch_hashes < <(git log --format=%h "${target_branch}")
+
+
+declare -A target_branch_commits
+for i in "${!target_branch_hashes[@]}"; do
+    target_branch_commits["${target_branch_hashes[$i]}"]="${target_branch_messages[$i]}"
+done
+
+while IFS= read -r msg; do
+    if [ -z "$msg" ] ; then
+        continue
+    fi
+    found=0
+    for target_branch_hash in "${!target_branch_commits[@]}"; do
+        if [ "${target_branch_commits[${target_branch_hash}]}" == "${msg}" ] ; then
+            found=1
+            break
+        fi
+    done
+    if [ ${found} == 0 ] ; then
+        git log --format="%h %s" "${source_branch}" | grep "${msg}"
+    fi
+done < <(git log --format=%s "${source_branch}")