Handle boundary conditions internally

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}")

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

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}
@@ -2012,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).
@@ -8223,6 +8212,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}
 
 

gromacs/src/applied_forces/colvars/colvarproxygromacs.cpp.patch

@@ -1,67 +1,64 @@
 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();
 -        errorCode |= colvars->setup_input();
 +        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;
      }
  }
- 
+

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,33 +46,78 @@ 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");
      }
 
      if (mpiComm.isMainRank())
      {
 -        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
 +    }
 
      // 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,11 +129,12 @@ 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)
  {
 -    colvars->write_state_buffer(stateToCheckpoint_.colvarStateFile_);
 +    cvmodule->write_state_buffer(stateToCheckpoint_.colvarStateFile_);
      stateToCheckpoint_.writeState(checkpointWriting.builder_, moduleName);
  }
+

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

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) +

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;
 }