From b62dd060d011bfc38ed369f424757d74ac76a4c3 Mon Sep 17 00:00:00 2001
From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com>
Date: Thu, 28 Aug 2025 11:17:11 +0200
Subject: [PATCH 1/9] preparing multicolvarTemplate for an alternative PTM

---
 src/colvar/MultiColvarTemplate.h | 49 ++++++++++++++++----------------
 src/core/ParallelTaskManager.h   |  5 ++++
 2 files changed, 30 insertions(+), 24 deletions(-)

diff --git a/src/colvar/MultiColvarTemplate.h b/src/colvar/MultiColvarTemplate.h
index b8db6885fc..82e5ea4b7b 100644
--- a/src/colvar/MultiColvarTemplate.h
+++ b/src/colvar/MultiColvarTemplate.h
@@ -47,11 +47,12 @@ struct MultiColvarInput {
   }
 };
 
-template <class T>
+template <class T, typename myPTM=defaultPTM>
 class MultiColvarTemplate : public ActionWithVector {
 public:
   using input_type = MultiColvarInput;
-  using PTM = ParallelTaskManager<MultiColvarTemplate<T>>;
+  using mytype=MultiColvarTemplate<T,myPTM>;
+  using PTM =typename myPTM::template PTM<mytype>;
   constexpr static size_t virialSize = 9;
 private:
 /// The parallel task manager
@@ -89,8 +90,8 @@ class MultiColvarTemplate : public ActionWithVector {
                                ForceIndexHolder force_indices );
 };
 
-template <class T>
-void MultiColvarTemplate<T>::registerKeywords(Keywords& keys ) {
+template <class T, typename myPTM>
+void MultiColvarTemplate<T,myPTM>::registerKeywords(Keywords& keys ) {
   T::registerKeywords( keys );
   PTM::registerKeywords( keys );
   keys.addInputKeyword("optional","MASK","vector","the label for a sparse vector that should be used to determine which elements of the vector should be computed");
@@ -105,8 +106,8 @@ void MultiColvarTemplate<T>::registerKeywords(Keywords& keys ) {
   }
 }
 
-template <class T>
-MultiColvarTemplate<T>::MultiColvarTemplate(const ActionOptions&ao):
+template <class T, typename myPTM>
+MultiColvarTemplate<T,myPTM>::MultiColvarTemplate(const ActionOptions&ao):
   Action(ao),
   ActionWithVector(ao),
   taskmanager(this),
@@ -170,40 +171,40 @@ MultiColvarTemplate<T>::MultiColvarTemplate(const ActionOptions&ao):
   taskmanager.setActionInput( MultiColvarInput{ usepbc, mode, natoms_per_task });
 }
 
-template <class T>
-unsigned MultiColvarTemplate<T>::getNumberOfDerivatives() {
+template <class T, typename myPTM>
+unsigned MultiColvarTemplate<T,myPTM>::getNumberOfDerivatives() {
   return 3*getNumberOfAtoms()+9;
 }
 
-template <class T>
-void MultiColvarTemplate<T>::calculate() {
+template <class T, typename myPTM>
+void MultiColvarTemplate<T,myPTM>::calculate() {
   if( wholemolecules ) {
     makeWhole();
   }
   taskmanager.runAllTasks();
 }
 
-template <class T>
-void MultiColvarTemplate<T>::applyNonZeroRankForces( std::vector<double>& outforces ) {
+template <class T, typename myPTM>
+void MultiColvarTemplate<T,myPTM>::applyNonZeroRankForces( std::vector<double>& outforces ) {
   taskmanager.applyForces( outforces );
 }
 
-template <class T>
-void MultiColvarTemplate<T>::addValueWithDerivatives( const std::vector<std::size_t>& shape ) {
+template <class T, typename myPTM>
+void MultiColvarTemplate<T,myPTM>::addValueWithDerivatives( const std::vector<std::size_t>& shape ) {
   std::vector<std::size_t> s(1);
   s[0]=getNumberOfAtoms() / natoms_per_task;
   addValue( s );
 }
 
-template <class T>
-void MultiColvarTemplate<T>::addComponentWithDerivatives( const std::string& name, const std::vector<std::size_t>& shape ) {
+template <class T, typename myPTM>
+void MultiColvarTemplate<T,myPTM>::addComponentWithDerivatives( const std::string& name, const std::vector<std::size_t>& shape ) {
   std::vector<std::size_t> s(1);
   s[0]=getNumberOfAtoms() / natoms_per_task;
   addComponent( name, s );
 }
 
-template <class T>
-void MultiColvarTemplate<T>::getInputData( std::vector<double>& inputdata ) const {
+template <class T, typename myPTM>
+void MultiColvarTemplate<T,myPTM>::getInputData( std::vector<double>& inputdata ) const {
   std::size_t ntasks = getConstPntrToComponent(0)->getNumberOfStoredValues();
   if( inputdata.size()!=5*natoms_per_task*ntasks ) {
     inputdata.resize( 5*natoms_per_task*ntasks );
@@ -231,8 +232,8 @@ void MultiColvarTemplate<T>::getInputData( std::vector<double>& inputdata ) cons
   }
 }
 
-template <class T>
-void MultiColvarTemplate<T>::performTask( std::size_t task_index,
+template <class T, typename myPTM>
+void MultiColvarTemplate<T,myPTM>::performTask( std::size_t task_index,
     const MultiColvarInput& actiondata,
     ParallelActionsInput& input,
     ParallelActionsOutput& output ) {
@@ -293,14 +294,14 @@ void MultiColvarTemplate<T>::performTask( std::size_t task_index,
                   cvout );
 }
 
-template <class T>
-int MultiColvarTemplate<T>::getNumberOfValuesPerTask( std::size_t task_index,
+template <class T, typename myPTM>
+int MultiColvarTemplate<T,myPTM>::getNumberOfValuesPerTask( std::size_t task_index,
     const MultiColvarInput& actiondata ) {
   return 1;
 }
 
-template <class T>
-void MultiColvarTemplate<T>::getForceIndices( std::size_t task_index,
+template <class T, typename myPTM>
+void MultiColvarTemplate<T,myPTM>::getForceIndices( std::size_t task_index,
     std::size_t colno,
     std::size_t ntotal_force,
     const MultiColvarInput& actiondata,
diff --git a/src/core/ParallelTaskManager.h b/src/core/ParallelTaskManager.h
index b40bfaf5c9..dad0391801 100644
--- a/src/core/ParallelTaskManager.h
+++ b/src/core/ParallelTaskManager.h
@@ -442,6 +442,11 @@ class ParallelTaskManager {
                                         View<double> forces );
 };
 
+struct defaultPTM {
+  template <typename ACC>
+  using PTM=ParallelTaskManager<ACC>;
+};
+
 template <class T>
 void ParallelTaskManager<T>::registerKeywords( Keywords& keys ) {
   keys.addFlag("USEGPU",false,"run this calculation on the GPU");

From 3ef4328a33cccef32141ba9a4cc75981d3268a1a Mon Sep 17 00:00:00 2001
From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com>
Date: Thu, 28 Aug 2025 16:20:21 +0200
Subject: [PATCH 2/9] preparing distance for the new plugin

---
 src/colvar/Distance.cpp        | 333 +-----------------------------
 src/colvar/Distance.h          | 359 +++++++++++++++++++++++++++++++++
 src/core/ParallelTaskManager.h |   2 +-
 3 files changed, 364 insertions(+), 330 deletions(-)
 create mode 100644 src/colvar/Distance.h

diff --git a/src/colvar/Distance.cpp b/src/colvar/Distance.cpp
index 5dad010d39..b321a6973a 100644
--- a/src/colvar/Distance.cpp
+++ b/src/colvar/Distance.cpp
@@ -19,344 +19,19 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-//this is temporary:
-#ifdef __PLUMED_HAS_OPENACC
-#define __PLUMED_USE_OPENACC 1
-#endif //__PLUMED_HAS_OPENACC
-#include "Colvar.h"
+
+#include "Distance.h"
 #include "ColvarShortcut.h"
-#include "MultiColvarTemplate.h"
 #include "core/ActionRegister.h"
-#include "tools/Pbc.h"
 
 namespace PLMD {
 namespace colvar {
 
-//+PLUMEDOC COLVAR DISTANCE
-/*
-Calculate the distance/s between pairs of atoms.
-
-The following example illustrates how this action can be used to calculate and print the distance between atom 1
-and atom 2.
-
-```plumed
-d: DISTANCE ATOMS=1,2
-PRINT ARG=d FILE=colvar
-```
-
-By default the distance is computed in a way that takes periodic
-boundary conditions in account. This behavior can be changed by using the NOPBC flag.
-Furthermore, if you wish to calculate the vector connecting a pair of atoms you can use the
-`COMPONENTS` flag as shown below:
-
-```plumed
-d: DISTANCE ATOMS=1,2 COMPONENTS
-PRINT ARG=d.x,d.y,d.z FILE=colvar
-```
-
-Alternatively, you can calculate the components projected on the lattice vector by using the `SCALED_COMPONENTS`
-flag as shown below;
-
-```plumed
-d: DISTANCE ATOMS=1,2 SCALED_COMPONENTS
-PRINT ARG=d.a,d.b,d.c FILE=colvar
-```
-
-The advantage of using `SCALED_COMPONENTS` over `COMPONENTS` is that the a, b and c variables
-that are calculated when `SCALED_COMPONENTS` is employed have the proper periodicity. This feature is useful
-if you wish to study the motion of a molecule across a membrane.
-
-You can also use this command to calculate multiple indistinguishable distances or vectors with a single
-line of PLUMED input.  For example, the following input calculates and outputs the distances between four
-pairs of atoms:
-
-```plumed
-d: DISTANCE ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 ATOMS4=7,8
-PRINT ARG=d FILE=colvar
-```
-
-By a similar token, the following input outputs three four dimensional vectors that contain the x, y and z
-components of the vectors connecting the four atoms:
-
-```plumed
-d: DISTANCE COMPONENTS ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 ATOMS4=7,8
-PRINT ARG=d.x,d.y,d.z FILE=colvar
-```
-
-You can also replace COMPONENTS with SCALED_COMPONENTS in the above input and obtain the projects of these vectors
-on the lattice vectors.
-
-## Managing periodic boundary conditions
-
-When using the DISTANCE command to calculate the end-to-end distance for a large polymer you need to ensure that you
-are managing PBCs correctly.  This problems that can occur with these calculations are explained at length in the
-early parts of the document that is referenced in the bibliography. Notice, however, that the input
-provides an example of an input that could be used to compute the end-to-end distance for a polymer
-of 100 atoms and keeps it at a value around 5.
-
-```plumed
-WHOLEMOLECULES ENTITY0=1-100
-e2e: DISTANCE ATOMS=1,100 NOPBC
-RESTRAINT ARG=e2e KAPPA=1 AT=5
-```
-
-Notice that NOPBC is used here so as to ensure that the distance is calculated correctely even if the end-to-end distance is larger than half the simulation
-Also notice that, since many MD codes break molecules across cell boundary, it might be necessary to
-use the [WHOLEMOLECULES](WHOLEMOLECULES.md) keyword (also notice that it should be _before_ distance). The list of atoms provided to [WHOLEMOLECULES](WHOLEMOLECULES.md)
-here contains all the atoms between 1 and 100. Strictly speaking, this
-is not necessary. If you know for sure that atoms with difference in
-the index say equal to 10 are _not_ going to be farther than half cell
-you can e.g. use
-
-```plumed
-WHOLEMOLECULES ENTITY0=1,10,20,30,40,50,60,70,80,90,100
-e2e: DISTANCE ATOMS=1,100 NOPBC
-RESTRAINT ARG=e2e KAPPA=1 AT=5
-```
-
-Just be sure that the ordered list provide to [WHOLEMOLECULES](WHOLEMOLECULES.md) has the following
-properties:
-
-- Consecutive atoms should be closer than half-cell throughout the entire simulation.
-- Atoms required later for the distance (e.g. 1 and 100) should be included in the list
-
-The following example shows how to take periodicity into account when computing the z-component of a distance
-
-```plumed
-# this is a center of mass of a large group
-c: COM ATOMS=1-100
-# this is the distance between atom 101 and the group
-d: DISTANCE ATOMS=c,101 COMPONENTS
-# this makes a new variable, dd, equal to d and periodic, with domain -10,10
-# this is the right choise if e.g. the cell is orthorombic and its size in
-# z direction is 20.
-dz: COMBINE ARG=d.z PERIODIC=-10,10
-# metadynamics on dd
-METAD ARG=dz SIGMA=0.1 HEIGHT=0.1 PACE=200
-```
-
-You can use the same input even if the DISTANCE command is calculating the vectors connecting multiple pairs of atoms.
-However, using SCALED_COMPONENTS ensures this problem does not arise because these variables are always periodic
-with domain (-0.5,+0.5).
-
-*/
-//+ENDPLUMEDOC
-
-class Distance : public Colvar {
-  bool components;
-  bool scaled_components;
-  bool pbc;
-
-  std::vector<double> value;
-  std::vector<double> derivs;
-public:
-  static void registerKeywords( Keywords& keys );
-  explicit Distance(const ActionOptions&);
-  static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
-  static unsigned getModeAndSetupValues( ActionWithValue* av );
-// active methods:
-  void calculate() override;
-  static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout );
-};
-
 typedef ColvarShortcut<Distance> DistanceShortcut;
 PLUMED_REGISTER_ACTION(DistanceShortcut,"DISTANCE")
 PLUMED_REGISTER_ACTION(Distance,"DISTANCE_SCALAR")
 typedef MultiColvarTemplate<Distance> DistanceMulti;
 PLUMED_REGISTER_ACTION(DistanceMulti,"DISTANCE_VECTOR")
 
-void Distance::registerKeywords( Keywords& keys ) {
-  Colvar::registerKeywords( keys );
-  keys.setDisplayName("DISTANCE");
-  constexpr auto scalarOrVector = Keywords::componentType::scalar | Keywords::componentType::vector;
-  keys.add("atoms","ATOMS","the pair of atom that we are calculating the distance between");
-  keys.addFlag("COMPONENTS",false,"calculate the x, y and z components of the distance separately and store them as label.x, label.y and label.z");
-  keys.addFlag("SCALED_COMPONENTS",false,"calculate the a, b and c scaled components of the distance separately and store them as label.a, label.b and label.c");
-  keys.addOutputComponent("x","COMPONENTS",scalarOrVector,"the x-component of the vector connecting the two atoms");
-  keys.addOutputComponent("y","COMPONENTS",scalarOrVector,"the y-component of the vector connecting the two atoms");
-  keys.addOutputComponent("z","COMPONENTS",scalarOrVector,"the z-component of the vector connecting the two atoms");
-  keys.addOutputComponent("a","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the first lattice vector of the vector connecting the two atoms");
-  keys.addOutputComponent("b","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the second lattice vector of the vector connecting the two atoms");
-  keys.addOutputComponent("c","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the third lattice vector of the vector connecting the two atoms");
-  keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
-  keys.setValueDescription(scalarOrVector,"the DISTANCE between this pair of atoms");
-  keys.addDOI("10.1007/978-1-4939-9608-7_21");
-  keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
-}
-
-Distance::Distance(const ActionOptions&ao):
-  PLUMED_COLVAR_INIT(ao),
-  components(false),
-  scaled_components(false),
-  pbc(true),
-  value(1) {
-  std::vector<AtomNumber> atoms;
-  parseAtomList(-1,atoms,this);
-  if(atoms.size()!=2) {
-    error("Number of specified atoms should be 2");
-  }
-
-  bool nopbc=!pbc;
-  parseFlag("NOPBC",nopbc);
-  pbc=!nopbc;
-
-  if(pbc) {
-    log.printf("  using periodic boundary conditions\n");
-  } else {
-    log.printf("  without periodic boundary conditions\n");
-  }
-
-  unsigned mode = getModeAndSetupValues( this );
-  if(mode==1) {
-    components=true;
-  } else if(mode==2) {
-    scaled_components=true;
-  }
-  if( components || scaled_components ) {
-    value.resize(3);
-  }
-  requestAtoms(atoms);
-}
-
-void Distance::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
-  aa->parseAtomList("ATOMS",num,t);
-  if( t.size()==2 ) {
-    aa->log.printf("  between atoms %d %d\n",t[0].serial(),t[1].serial());
-  }
-}
-
-unsigned Distance::getModeAndSetupValues( ActionWithValue* av ) {
-  bool c;
-  av->parseFlag("COMPONENTS",c);
-  bool sc;
-  av->parseFlag("SCALED_COMPONENTS",sc);
-  if( c && sc ) {
-    av->error("COMPONENTS and SCALED_COMPONENTS are not compatible");
-  }
-
-  if(c) {
-    av->addComponentWithDerivatives("x");
-    av->componentIsNotPeriodic("x");
-    av->addComponentWithDerivatives("y");
-    av->componentIsNotPeriodic("y");
-    av->addComponentWithDerivatives("z");
-    av->componentIsNotPeriodic("z");
-    av->log<<"  WARNING: components will not have the proper periodicity - see manual\n";
-    return 1;
-  } else if(sc) {
-    av->addComponentWithDerivatives("a");
-    av->componentIsPeriodic("a","-0.5","+0.5");
-    av->addComponentWithDerivatives("b");
-    av->componentIsPeriodic("b","-0.5","+0.5");
-    av->addComponentWithDerivatives("c");
-    av->componentIsPeriodic("c","-0.5","+0.5");
-    return 2;
-  }
-  av->addValueWithDerivatives();
-  av->setNotPeriodic();
-  return 0;
-}
-
-// calculator
-void Distance::calculate() {
-
-  if(pbc) {
-    makeWhole();
-  }
-
-  if( components ) {
-    ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) );
-    calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), cvout );
-    Value* valuex=getPntrToComponent("x");
-    Value* valuey=getPntrToComponent("y");
-    Value* valuez=getPntrToComponent("z");
-
-    for(unsigned i=0; i<2; ++i) {
-      setAtomsDerivatives(valuex,i,cvout.getAtomDerivatives(0,i) );
-    }
-    setBoxDerivatives(valuex,cvout.virial[0]);
-    valuex->set(value[0]);
-
-    for(unsigned i=0; i<2; ++i) {
-      setAtomsDerivatives(valuey,i,cvout.getAtomDerivatives(1,i) );
-    }
-    setBoxDerivatives(valuey,cvout.virial[1]);
-    valuey->set(value[1]);
-
-    for(unsigned i=0; i<2; ++i) {
-      setAtomsDerivatives(valuez,i,cvout.getAtomDerivatives(2,i) );
-    }
-    setBoxDerivatives(valuez,cvout.virial[2]);
-    valuez->set(value[2]);
-  } else if( scaled_components ) {
-    ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) );
-    calculateCV( ColvarInput::createColvarInput( 2, getPositions(), this ), cvout );
-
-    Value* valuea=getPntrToComponent("a");
-    Value* valueb=getPntrToComponent("b");
-    Value* valuec=getPntrToComponent("c");
-    for(unsigned i=0; i<2; ++i) {
-      setAtomsDerivatives(valuea,i,cvout.getAtomDerivatives(0,i) );
-    }
-    valuea->set(value[0]);
-    for(unsigned i=0; i<2; ++i) {
-      setAtomsDerivatives(valueb,i,cvout.getAtomDerivatives(1,i) );
-    }
-    valueb->set(value[1]);
-    for(unsigned i=0; i<2; ++i) {
-      setAtomsDerivatives(valuec,i,cvout.getAtomDerivatives(2,i) );
-    }
-    valuec->set(value[2]);
-  } else  {
-    ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) );
-    calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout );
-    for(unsigned i=0; i<2; ++i) {
-      setAtomsDerivatives(i,cvout.getAtomDerivatives(0,i) );
-    }
-    setBoxDerivatives(cvout.virial[0]);
-    setValue           (value[0]);
-  }
-}
-
-void Distance::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
-  Vector distance=delta(cvin.pos[0],cvin.pos[1]);
-
-  if(cvin.mode==1) {
-    cvout.derivs[0][0] = Vector(-1,0,0);
-    cvout.derivs[0][1] = Vector(+1,0,0);
-    cvout.values[0] = distance[0];
-
-    cvout.derivs[1][0] = Vector(0,-1,0);
-    cvout.derivs[1][1] = Vector(0,+1,0);
-    cvout.values[1] = distance[1];
-
-    cvout.derivs[2][0] = Vector(0,0,-1);
-    cvout.derivs[2][1] = Vector(0,0,+1);
-    cvout.values[2] = distance[2];
-    ColvarInput::setBoxDerivativesNoPbc( cvin, cvout );
-  } else if(cvin.mode==2) {
-    Vector d=cvin.pbc.realToScaled(distance);
-    cvout.derivs[0][0] = matmul(cvin.pbc.getInvBox(),Vector(-1,0,0));
-    cvout.derivs[0][1] = matmul(cvin.pbc.getInvBox(),Vector(+1,0,0));
-    cvout.values[0] = Tools::pbc(d[0]);
-    cvout.derivs[1][0] = matmul(cvin.pbc.getInvBox(),Vector(0,-1,0));
-    cvout.derivs[1][1] = matmul(cvin.pbc.getInvBox(),Vector(0,+1,0));
-    cvout.values[1] = Tools::pbc(d[1]);
-    cvout.derivs[2][0] = matmul(cvin.pbc.getInvBox(),Vector(0,0,-1));
-    cvout.derivs[2][1] = matmul(cvin.pbc.getInvBox(),Vector(0,0,+1));
-    cvout.values[2] = Tools::pbc(d[2]);
-  } else {
-    const double value=distance.modulo();
-    const double invvalue=1.0/value;
-    cvout.derivs[0][0] = -invvalue*distance;
-    cvout.derivs[0][1] = invvalue*distance;
-    ColvarInput::setBoxDerivativesNoPbc( cvin, cvout );
-    cvout.values[0] = value;
-  }
-}
-
-}
-}
-
-
-
+} //namespace colvar
+} //namespace PLMD
diff --git a/src/colvar/Distance.h b/src/colvar/Distance.h
new file mode 100644
index 0000000000..2c0128cb46
--- /dev/null
+++ b/src/colvar/Distance.h
@@ -0,0 +1,359 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2011-2023 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_colvar_Distance_h
+#define __PLUMED_colvar_Distance_h
+#include "Colvar.h"
+#include "MultiColvarTemplate.h"
+#include "tools/Pbc.h"
+
+namespace PLMD {
+namespace colvar {
+
+//+PLUMEDOC COLVAR DISTANCE
+/*
+Calculate the distance/s between pairs of atoms.
+
+The following example illustrates how this action can be used to calculate and print the distance between atom 1
+and atom 2.
+
+```plumed
+d: DISTANCE ATOMS=1,2
+PRINT ARG=d FILE=colvar
+```
+
+By default the distance is computed in a way that takes periodic
+boundary conditions in account. This behavior can be changed by using the NOPBC flag.
+Furthermore, if you wish to calculate the vector connecting a pair of atoms you can use the
+`COMPONENTS` flag as shown below:
+
+```plumed
+d: DISTANCE ATOMS=1,2 COMPONENTS
+PRINT ARG=d.x,d.y,d.z FILE=colvar
+```
+
+Alternatively, you can calculate the components projected on the lattice vector by using the `SCALED_COMPONENTS`
+flag as shown below;
+
+```plumed
+d: DISTANCE ATOMS=1,2 SCALED_COMPONENTS
+PRINT ARG=d.a,d.b,d.c FILE=colvar
+```
+
+The advantage of using `SCALED_COMPONENTS` over `COMPONENTS` is that the a, b and c variables
+that are calculated when `SCALED_COMPONENTS` is employed have the proper periodicity. This feature is useful
+if you wish to study the motion of a molecule across a membrane.
+
+You can also use this command to calculate multiple indistinguishable distances or vectors with a single
+line of PLUMED input.  For example, the following input calculates and outputs the distances between four
+pairs of atoms:
+
+```plumed
+d: DISTANCE ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 ATOMS4=7,8
+PRINT ARG=d FILE=colvar
+```
+
+By a similar token, the following input outputs three four dimensional vectors that contain the x, y and z
+components of the vectors connecting the four atoms:
+
+```plumed
+d: DISTANCE COMPONENTS ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 ATOMS4=7,8
+PRINT ARG=d.x,d.y,d.z FILE=colvar
+```
+
+You can also replace COMPONENTS with SCALED_COMPONENTS in the above input and obtain the projects of these vectors
+on the lattice vectors.
+
+## Managing periodic boundary conditions
+
+When using the DISTANCE command to calculate the end-to-end distance for a large polymer you need to ensure that you
+are managing PBCs correctly.  This problems that can occur with these calculations are explained at length in the
+early parts of the document that is referenced in the bibliography. Notice, however, that the input
+provides an example of an input that could be used to compute the end-to-end distance for a polymer
+of 100 atoms and keeps it at a value around 5.
+
+```plumed
+WHOLEMOLECULES ENTITY0=1-100
+e2e: DISTANCE ATOMS=1,100 NOPBC
+RESTRAINT ARG=e2e KAPPA=1 AT=5
+```
+
+Notice that NOPBC is used here so as to ensure that the distance is calculated correctely even if the end-to-end distance is larger than half the simulation
+Also notice that, since many MD codes break molecules across cell boundary, it might be necessary to
+use the [WHOLEMOLECULES](WHOLEMOLECULES.md) keyword (also notice that it should be _before_ distance). The list of atoms provided to [WHOLEMOLECULES](WHOLEMOLECULES.md)
+here contains all the atoms between 1 and 100. Strictly speaking, this
+is not necessary. If you know for sure that atoms with difference in
+the index say equal to 10 are _not_ going to be farther than half cell
+you can e.g. use
+
+```plumed
+WHOLEMOLECULES ENTITY0=1,10,20,30,40,50,60,70,80,90,100
+e2e: DISTANCE ATOMS=1,100 NOPBC
+RESTRAINT ARG=e2e KAPPA=1 AT=5
+```
+
+Just be sure that the ordered list provide to [WHOLEMOLECULES](WHOLEMOLECULES.md) has the following
+properties:
+
+- Consecutive atoms should be closer than half-cell throughout the entire simulation.
+- Atoms required later for the distance (e.g. 1 and 100) should be included in the list
+
+The following example shows how to take periodicity into account when computing the z-component of a distance
+
+```plumed
+# this is a center of mass of a large group
+c: COM ATOMS=1-100
+# this is the distance between atom 101 and the group
+d: DISTANCE ATOMS=c,101 COMPONENTS
+# this makes a new variable, dd, equal to d and periodic, with domain -10,10
+# this is the right choise if e.g. the cell is orthorombic and its size in
+# z direction is 20.
+dz: COMBINE ARG=d.z PERIODIC=-10,10
+# metadynamics on dd
+METAD ARG=dz SIGMA=0.1 HEIGHT=0.1 PACE=200
+```
+
+You can use the same input even if the DISTANCE command is calculating the vectors connecting multiple pairs of atoms.
+However, using SCALED_COMPONENTS ensures this problem does not arise because these variables are always periodic
+with domain (-0.5,+0.5).
+
+*/
+//+ENDPLUMEDOC
+
+class Distance : public Colvar {
+  bool components;
+  bool scaled_components;
+  bool pbc;
+
+  std::vector<double> value;
+  std::vector<double> derivs;
+public:
+  static void registerKeywords( Keywords& keys );
+  explicit Distance(const ActionOptions&);
+  static void parseAtomList( const int& num,
+                             std::vector<AtomNumber>& t,
+                             ActionAtomistic* aa );
+  static unsigned getModeAndSetupValues( ActionWithValue* av );
+// active methods:
+  void calculate() override;
+  static void calculateCV( const ColvarInput& cvin,
+                           ColvarOutput& cvout );
+};
+
+void Distance::registerKeywords( Keywords& keys ) {
+  Colvar::registerKeywords( keys );
+  keys.setDisplayName("DISTANCE");
+  constexpr auto scalarOrVector = Keywords::componentType::scalar
+                                  | Keywords::componentType::vector;
+  keys.add("atoms","ATOMS","the pair of atom that we are calculating the distance between");
+  keys.addFlag("COMPONENTS",false,"calculate the x, y and z components of the distance separately and store them as label.x, label.y and label.z");
+  keys.addFlag("SCALED_COMPONENTS",false,"calculate the a, b and c scaled components of the distance separately and store them as label.a, label.b and label.c");
+  keys.addOutputComponent("x","COMPONENTS",scalarOrVector,"the x-component of the vector connecting the two atoms");
+  keys.addOutputComponent("y","COMPONENTS",scalarOrVector,"the y-component of the vector connecting the two atoms");
+  keys.addOutputComponent("z","COMPONENTS",scalarOrVector,"the z-component of the vector connecting the two atoms");
+  keys.addOutputComponent("a","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the first lattice vector of the vector connecting the two atoms");
+  keys.addOutputComponent("b","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the second lattice vector of the vector connecting the two atoms");
+  keys.addOutputComponent("c","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the third lattice vector of the vector connecting the two atoms");
+  keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
+  keys.setValueDescription(scalarOrVector,"the DISTANCE between this pair of atoms");
+  keys.addDOI("10.1007/978-1-4939-9608-7_21");
+  keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
+}
+
+Distance::Distance(const ActionOptions&ao):
+  PLUMED_COLVAR_INIT(ao),
+  components(false),
+  scaled_components(false),
+  pbc(true),
+  value(1) {
+  std::vector<AtomNumber> atoms;
+  parseAtomList(-1,atoms,this);
+  if(atoms.size()!=2) {
+    error("Number of specified atoms should be 2");
+  }
+
+  bool nopbc=!pbc;
+  parseFlag("NOPBC",nopbc);
+  pbc=!nopbc;
+
+  if(pbc) {
+    log.printf("  using periodic boundary conditions\n");
+  } else {
+    log.printf("  without periodic boundary conditions\n");
+  }
+
+  unsigned mode = getModeAndSetupValues( this );
+  if(mode==1) {
+    components=true;
+  } else if(mode==2) {
+    scaled_components=true;
+  }
+  if( components || scaled_components ) {
+    value.resize(3);
+  }
+  requestAtoms(atoms);
+}
+
+void Distance::parseAtomList( const int& num,
+                              std::vector<AtomNumber>& t,
+                              ActionAtomistic* aa ) {
+  aa->parseAtomList("ATOMS",num,t);
+  if( t.size()==2 ) {
+    aa->log.printf("  between atoms %d %d\n",t[0].serial(),t[1].serial());
+  }
+}
+
+unsigned Distance::getModeAndSetupValues( ActionWithValue* av ) {
+  bool c;
+  av->parseFlag("COMPONENTS",c);
+  bool sc;
+  av->parseFlag("SCALED_COMPONENTS",sc);
+  if( c && sc ) {
+    av->error("COMPONENTS and SCALED_COMPONENTS are not compatible");
+  }
+
+  if(c) {
+    av->addComponentWithDerivatives("x");
+    av->componentIsNotPeriodic("x");
+    av->addComponentWithDerivatives("y");
+    av->componentIsNotPeriodic("y");
+    av->addComponentWithDerivatives("z");
+    av->componentIsNotPeriodic("z");
+    av->log<<"  WARNING: components will not have the proper periodicity - see manual\n";
+    return 1;
+  } else if(sc) {
+    av->addComponentWithDerivatives("a");
+    av->componentIsPeriodic("a","-0.5","+0.5");
+    av->addComponentWithDerivatives("b");
+    av->componentIsPeriodic("b","-0.5","+0.5");
+    av->addComponentWithDerivatives("c");
+    av->componentIsPeriodic("c","-0.5","+0.5");
+    return 2;
+  }
+  av->addValueWithDerivatives();
+  av->setNotPeriodic();
+  return 0;
+}
+
+// calculator
+void Distance::calculate() {
+
+  if(pbc) {
+    makeWhole();
+  }
+
+  if( components ) {
+    ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) );
+    calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), cvout );
+    Value* valuex=getPntrToComponent("x");
+    Value* valuey=getPntrToComponent("y");
+    Value* valuez=getPntrToComponent("z");
+
+    for(unsigned i=0; i<2; ++i) {
+      setAtomsDerivatives(valuex,i,cvout.getAtomDerivatives(0,i) );
+    }
+    setBoxDerivatives(valuex,cvout.virial[0]);
+    valuex->set(value[0]);
+
+    for(unsigned i=0; i<2; ++i) {
+      setAtomsDerivatives(valuey,i,cvout.getAtomDerivatives(1,i) );
+    }
+    setBoxDerivatives(valuey,cvout.virial[1]);
+    valuey->set(value[1]);
+
+    for(unsigned i=0; i<2; ++i) {
+      setAtomsDerivatives(valuez,i,cvout.getAtomDerivatives(2,i) );
+    }
+    setBoxDerivatives(valuez,cvout.virial[2]);
+    valuez->set(value[2]);
+  } else if( scaled_components ) {
+    ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) );
+    calculateCV( ColvarInput::createColvarInput( 2, getPositions(), this ), cvout );
+
+    Value* valuea=getPntrToComponent("a");
+    Value* valueb=getPntrToComponent("b");
+    Value* valuec=getPntrToComponent("c");
+    for(unsigned i=0; i<2; ++i) {
+      setAtomsDerivatives(valuea,i,cvout.getAtomDerivatives(0,i) );
+    }
+    valuea->set(value[0]);
+    for(unsigned i=0; i<2; ++i) {
+      setAtomsDerivatives(valueb,i,cvout.getAtomDerivatives(1,i) );
+    }
+    valueb->set(value[1]);
+    for(unsigned i=0; i<2; ++i) {
+      setAtomsDerivatives(valuec,i,cvout.getAtomDerivatives(2,i) );
+    }
+    valuec->set(value[2]);
+  } else  {
+    ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) );
+    calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout );
+    for(unsigned i=0; i<2; ++i) {
+      setAtomsDerivatives(i,cvout.getAtomDerivatives(0,i) );
+    }
+    setBoxDerivatives(cvout.virial[0]);
+    setValue           (value[0]);
+  }
+}
+
+void Distance::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
+  Vector distance=delta(cvin.pos[0],cvin.pos[1]);
+
+  if(cvin.mode==1) {
+    cvout.derivs[0][0] = Vector(-1,0,0);
+    cvout.derivs[0][1] = Vector(+1,0,0);
+    cvout.values[0] = distance[0];
+
+    cvout.derivs[1][0] = Vector(0,-1,0);
+    cvout.derivs[1][1] = Vector(0,+1,0);
+    cvout.values[1] = distance[1];
+
+    cvout.derivs[2][0] = Vector(0,0,-1);
+    cvout.derivs[2][1] = Vector(0,0,+1);
+    cvout.values[2] = distance[2];
+    ColvarInput::setBoxDerivativesNoPbc( cvin, cvout );
+  } else if(cvin.mode==2) {
+    Vector d=cvin.pbc.realToScaled(distance);
+    cvout.derivs[0][0] = matmul(cvin.pbc.getInvBox(),Vector(-1,0,0));
+    cvout.derivs[0][1] = matmul(cvin.pbc.getInvBox(),Vector(+1,0,0));
+    cvout.values[0] = Tools::pbc(d[0]);
+    cvout.derivs[1][0] = matmul(cvin.pbc.getInvBox(),Vector(0,-1,0));
+    cvout.derivs[1][1] = matmul(cvin.pbc.getInvBox(),Vector(0,+1,0));
+    cvout.values[1] = Tools::pbc(d[1]);
+    cvout.derivs[2][0] = matmul(cvin.pbc.getInvBox(),Vector(0,0,-1));
+    cvout.derivs[2][1] = matmul(cvin.pbc.getInvBox(),Vector(0,0,+1));
+    cvout.values[2] = Tools::pbc(d[2]);
+  } else {
+    const double value=distance.modulo();
+    const double invvalue=1.0/value;
+    cvout.derivs[0][0] = -invvalue*distance;
+    cvout.derivs[0][1] = invvalue*distance;
+    ColvarInput::setBoxDerivativesNoPbc( cvin, cvout );
+    cvout.values[0] = value;
+  }
+}
+
+} //namespace colvar
+} //namespace PLMD
+
+#endif //__PLUMED_colvar_Distance_h
+
+
diff --git a/src/core/ParallelTaskManager.h b/src/core/ParallelTaskManager.h
index dad0391801..d203e44acc 100644
--- a/src/core/ParallelTaskManager.h
+++ b/src/core/ParallelTaskManager.h
@@ -381,7 +381,7 @@ class ParallelTaskManager {
 //  static constexpr bool has_custom_gather=PTMUtils::has_gatherForces_custom<T>;
 //  static constexpr bool has_GPU_gather=PTMUtils::has_gatherForces_GPU<T>;
 //  static constexpr size_t virialSize = PTMUtils::virialSize<T>;
-private:
+protected:
 /// The underlying action for which we are managing parallel tasks
   ActionWithVector* action;
 /// The MPI communicator

From 77b4ac2f2ef35a4beb07cd4a81ef380170904463 Mon Sep 17 00:00:00 2001
From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com>
Date: Thu, 28 Aug 2025 16:23:25 +0200
Subject: [PATCH 3/9] compile, do not link

---
 plugins/openaccPTM/ACCParallelTaskManager.h | 345 ++++++++++++++++++++
 plugins/openaccPTM/Distance.cpp             |  34 ++
 plugins/openaccPTM/Makefile                 |  47 +++
 plugins/openaccPTM/configure.sh             |  30 ++
 plugins/openaccPTM/plumed                   |   1 +
 5 files changed, 457 insertions(+)
 create mode 100644 plugins/openaccPTM/ACCParallelTaskManager.h
 create mode 100644 plugins/openaccPTM/Distance.cpp
 create mode 100644 plugins/openaccPTM/Makefile
 create mode 100755 plugins/openaccPTM/configure.sh
 create mode 120000 plugins/openaccPTM/plumed

diff --git a/plugins/openaccPTM/ACCParallelTaskManager.h b/plugins/openaccPTM/ACCParallelTaskManager.h
new file mode 100644
index 0000000000..cc1480d7a3
--- /dev/null
+++ b/plugins/openaccPTM/ACCParallelTaskManager.h
@@ -0,0 +1,345 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2017-2023 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_acc_ParallelTaskManager_h
+#define __PLUMED_acc_ParallelTaskManager_h
+
+#include "plumed/core/ActionWithVector.h"
+#include "plumed/core/ActionWithMatrix.h"
+#include "plumed/core/ParallelTaskManager.h"
+#include "plumed/tools/Communicator.h"
+#include "plumed/tools/OpenMP.h"
+#include "plumed/tools/View.h"
+#include "plumed/tools/View2D.h"
+
+#include "plumed/tools/ColvarOutput.h"
+#include "plumed/tools/OpenACC.h"
+
+namespace PLMD {
+
+template <class T>
+class AccParallelTaskManager : public ParallelTaskManager<T> {
+public:
+  static void registerKeywords( Keywords& keys ) {
+    ParallelTaskManager<T>::registerKeywords(keys);
+  }
+  AccParallelTaskManager(ActionWithVector* av):
+    PLMD::ParallelTaskManager<T>(av) {}
+/// This runs all the tasks
+  void runAllTasks();
+/// Apply the forces on the parallel object
+  void applyForces( std::vector<double>& forcesForApply );
+};
+
+struct ACCPTM {
+  template <typename T>
+  using PTM=AccParallelTaskManager<T>;
+};
+
+//use the __PLUMED_USE_OPENACC_TASKSMINE macro to debug the ptm ins a single file
+//so that compiling witha a small modification will be faster (the ptm is included nearly everywhere)
+#ifndef __PLUMED_USE_OPENACC_TASKSMINE
+template <class T>
+void runAllTasksACC(typename T::input_type actiondata,
+                    ParallelActionsInput myinput,
+                    std::vector<double>& value_stash,
+                    const std::vector<unsigned> & partialTaskList,
+                    const unsigned nactive_tasks,
+                    const std::size_t nderivatives_per_task,
+                    const std::size_t workspace_size
+                   ) {
+  auto myinput_acc = OpenACC::fromToDataHelper(myinput);
+  auto actiondata_acc = OpenACC::fromToDataHelper(actiondata);
+
+  //template type is deduced
+  OpenACC::memoryManager vs{value_stash};
+  auto value_stash_data = vs.devicePtr();
+
+  OpenACC::memoryManager ptl{partialTaskList};
+  auto partialTaskList_data = ptl.devicePtr();
+
+  OpenACC::memoryManager<double> buff{workspace_size*nactive_tasks};
+
+  auto buffer = buff.devicePtr();
+  OpenACC::memoryManager<double> dev(nderivatives_per_task*nactive_tasks);
+  auto derivatives = dev.devicePtr();
+#pragma acc parallel loop present(myinput, actiondata) \
+                           copyin(nactive_tasks, \
+                                 nderivatives_per_task, \
+                                 workspace_size)\
+                        deviceptr(derivatives, \
+                                  partialTaskList_data, \
+                                  value_stash_data, \
+                                  buffer) \
+                          default(none)
+  for(unsigned i=0; i<nactive_tasks; ++i) {
+    std::size_t task_index = partialTaskList_data[i];
+    std::size_t val_pos = task_index*myinput.nscalars;
+    auto myout = ParallelActionsOutput::create (myinput.nscalars,
+                 value_stash_data+val_pos,
+                 nderivatives_per_task,
+                 derivatives+nderivatives_per_task*i,
+                 workspace_size,
+                 (workspace_size>0)?
+                 buffer+workspace_size*i
+                 :nullptr  );
+    // Calculate the stuff in the loop for this action
+    T::performTask( task_index, actiondata, myinput, myout );
+  }
+  vs.copyFromDevice(value_stash.data());
+}
+#else
+template <class T>
+void runAllTasksACC(typename T::input_type actiondata,
+                    ParallelActionsInput myinput,
+                    std::vector<double>& value_stash,
+                    const std::vector<unsigned> & partialTaskList,
+                    const unsigned nactive_tasks,
+                    const std::size_t nderivatives_per_task,
+                    const std::size_t workspace_size
+                   ) ;
+#endif //__PLUMED_USE_OPENACC_TASKSMINE
+
+template <class T>
+void AccParallelTaskManager<T>::runAllTasks() {
+  // Get the list of active tasks
+  std::vector<unsigned> & partialTaskList( this->action->getListOfActiveTasks( this->action ) );
+  unsigned nactive_tasks=partialTaskList.size();
+  // Get all the input data so we can broadcast it to the GPU
+  this->myinput.noderiv = true;
+  this->action->getInputData( this->input_buffer );
+  this->myinput.dataSize = this->input_buffer.size();
+  this->myinput.inputdata = this->input_buffer.data();
+  // Transfer all the bookeeping information about the arguments
+  this->argumentsMap.setupArguments( this->action );
+  this->myinput.setupArguments( this->argumentsMap );
+  // Reset the values at the start of the task loop
+  std::size_t totalvals=this->getValueStashSize();
+  if( this->value_stash.size()!=totalvals ) {
+    this->value_stash.resize(totalvals);
+  }
+  std::fill (this->value_stash.begin(),this->value_stash.end(), 0.0);
+  if (this->comm.Get_rank()== 0) {// no multigpu shenanigans until this works
+    runAllTasksACC<T>(
+      this->actiondata,
+      this->myinput,
+      this->value_stash,
+      partialTaskList,
+      nactive_tasks,
+      this->nderivatives_per_task,
+      this->workspace_size
+    );
+  }
+  this->comm.Bcast( this->value_stash.data(), this->value_stash.size(), 0);
+
+  // MPI Gather everything
+  if( !this->action->runInSerial() ) {
+    this->comm.Sum( this->value_stash );
+  }
+
+// And transfer the value to the output values
+  this->action->transferStashToValues( this->value_stash );
+}
+
+//use the __PLUMED_USE_OPENACC_FORCESMINE macro to debug the ptm ins a single file
+//so that compiling witha a small modification will be faster (the ptm is included nearly everywhere)
+#ifndef __PLUMED_USE_OPENACC_FORCESMINE
+template <class T>
+void applyForcesWithACC(PLMD::View<double> forcesForApply,
+                        typename T::input_type actiondata,
+                        ParallelActionsInput myinput,
+                        const std::vector<double>& value_stash,
+                        const std::vector<unsigned> & partialTaskList,
+                        const unsigned nactive_tasks,
+                        const std::size_t nderivatives_per_task,
+                        const std::size_t workspace_size
+                       ) {
+  auto myinput_acc = OpenACC::fromToDataHelper(myinput);
+  auto actiondata_acc = OpenACC::fromToDataHelper(actiondata);
+
+  //template type is deduced
+  OpenACC::memoryManager vs{value_stash};
+  auto value_stash_data = vs.devicePtr();
+
+  OpenACC::memoryManager ptl{partialTaskList};
+  auto partialTaskList_data = ptl.devicePtr();
+
+  OpenACC::memoryManager ffa {forcesForApply};
+  auto forcesForApply_data = ffa.devicePtr();
+  const auto forcesForApply_size = ffa.size();
+  const auto nind_per_scalar = ForceIndexHolder::indexesPerScalar(myinput);
+  //nscalars is >=ncomponents (see setupParallelTaskManager )
+  const auto nind_per_task = nind_per_scalar*myinput.nscalars;
+
+  OpenACC::memoryManager<double> dev{nderivatives_per_task*nactive_tasks};
+  auto derivatives = dev.devicePtr();
+  OpenACC::memoryManager<std::size_t> ind{nind_per_task*nactive_tasks};
+  auto indices = ind.devicePtr();
+  OpenACC::memoryManager<double> vtmp{myinput.sizeOfFakeVals()*nactive_tasks};
+  auto valstmp = vtmp.devicePtr();
+  OpenACC::memoryManager<double> buff{workspace_size*nactive_tasks};
+  auto buffer = buff.devicePtr();
+
+#define forces_indicesArg(taskID,scalarID) ForceIndexHolder::create(myinput, \
+                          indices + taskID*nind_per_task + scalarID*nind_per_scalar)
+#define derivativeDrift(taskID,scalarID)  taskID*nderivatives_per_task \
+                           + scalarID*myinput.ncomponents*myinput.nderivatives_per_scalar
+#define stashDrift(taskID,scalarID) taskID*myinput.nscalars \
+                           + scalarID*myinput.ncomponents
+
+#pragma acc data present(myinput,actiondata) \
+                  copyin(nactive_tasks, \
+                         forcesForApply_size, \
+                         nderivatives_per_task, nind_per_task,nind_per_scalar, \
+                         workspace_size) \
+               deviceptr(derivatives, \
+                         indices, \
+                         value_stash_data, \
+                         partialTaskList_data, \
+                         forcesForApply_data, \
+                         valstmp, \
+                         buffer) \
+                 default(none)
+  {
+#pragma acc parallel loop
+    for(unsigned t=0; t<nactive_tasks; ++t) {
+      std::size_t task_index = partialTaskList_data[t];
+      auto myout = ParallelActionsOutput::create( myinput.nscalars,
+                   valstmp+myinput.nscalars*t,
+                   nderivatives_per_task,
+                   derivatives+nderivatives_per_task*t,
+                   workspace_size,
+                   (workspace_size>0)?buffer+workspace_size*t:nullptr);
+      // Calculate the stuff in the loop for this action
+      T::performTask( task_index, actiondata, myinput, myout );
+      // If this is a matrix this returns a number that isn't one as we have to loop over the columns
+      const std::size_t nvpt = T::getNumberOfValuesPerTask( task_index, actiondata );
+#pragma acc loop seq
+      for(unsigned vID=0; vID<nvpt; ++vID) {
+        auto force_indices = forces_indicesArg(t,vID);
+        // Create a force index holder
+        // Get the indices for forces
+        T::getForceIndices( task_index,
+                            vID,
+                            forcesForApply_size,
+                            actiondata,
+                            myinput,
+                            force_indices );
+
+        // Create a force input object
+        auto finput = ForceInput::create ( myinput.nscalars,
+                                           value_stash_data + stashDrift(task_index,vID),
+                                           myinput.nderivatives_per_scalar,
+                                           derivatives + derivativeDrift(t,vID));
+
+        // Gather forces that can be gathered locally
+        ParallelTaskManager<T>::gatherThreadSafeForces( myinput,
+            force_indices,
+            finput,
+            View<double>(forcesForApply_data,
+                         forcesForApply_size));
+      }
+    }
+
+#pragma acc parallel loop
+    for(unsigned v=myinput.threadunsafe_forces_start; v<forcesForApply_size; ++v) {
+      double tmp = 0.0;
+#pragma acc loop reduction(+:tmp)
+      for(unsigned t=0; t<nactive_tasks; ++t) {
+        const std::size_t task_index = partialTaskList_data[t];
+        const std::size_t nvpt = T::getNumberOfValuesPerTask( task_index, actiondata );
+        for(unsigned vID=0; vID<nvpt; ++vID) {
+          auto force_indices = forces_indicesArg(t,vID);
+
+          auto fdata = ForceInput::create( myinput.nscalars,
+                                           value_stash_data + stashDrift(task_index,vID),
+                                           myinput.nderivatives_per_scalar,
+                                           derivatives + derivativeDrift(t,vID));
+          for(unsigned i=0; i<myinput.ncomponents; ++i) {
+            const double ff = fdata.force[i];
+            for(unsigned d=force_indices.threadsafe_derivatives_end[i];
+                d<force_indices.tot_indices[i]; ++d) {
+              if( force_indices.indices[i][d]==v ) {
+                tmp += ff*fdata.deriv[i][d];
+                // break;
+              }
+            }
+          }
+        }
+      }
+      forcesForApply_data[v] = tmp;
+    }
+  }
+#undef forces_indicesArg
+#undef derivativeDrift
+#undef stashDrift
+  ffa.copyFromDevice(forcesForApply.data());
+}
+#else
+template <class T>
+void applyForcesWithACC(PLMD::View<double> forcesForApply,
+                        typename T::input_type actiondata,
+                        ParallelActionsInput myinput,
+                        const std::vector<double>& value_stash,
+                        const std::vector<unsigned> & partialTaskList,
+                        const unsigned nactive_tasks,
+                        const std::size_t nderivatives_per_task,
+                        const std::size_t workspace_size
+                       );
+#endif //__PLUMED_USE_OPENACC_FORCESMINE
+
+template <class T>
+void AccParallelTaskManager<T>::applyForces( std::vector<double>& forcesForApply ) {
+  // Get the list of active tasks
+  std::vector<unsigned> & partialTaskList= this->action->getListOfActiveTasks( this->action ) ;
+  unsigned nactive_tasks=partialTaskList.size();
+  // Clear force buffer
+  forcesForApply.assign( forcesForApply.size(), 0.0 );
+  //TODO: check if std::fill is faster (i get conflicting answers on the net)
+  //std::fill (forcesForApply.begin(),forcesForApply.end(), 0.0);
+  // Get all the input data so we can broadcast it to the GPU
+  this->myinput.noderiv = false;
+  // Retrieve the forces from the values
+  this->action->transferForcesToStash( this->value_stash );
+
+  std::fill (this->omp_forces[0].begin(),this->omp_forces[0].end(), 0.0);
+  if (this->comm.Get_rank() == 0) {
+    applyForcesWithACC<T>(
+      PLMD::View<double> { forcesForApply.data(), forcesForApply.size() },
+      this->actiondata,
+      this->myinput,
+      this->value_stash,
+      partialTaskList,
+      nactive_tasks,
+      this->nderivatives_per_task,
+      this->workspace_size
+    );
+  }
+  // MPI Gather everything (this must be extended to the gpu thing, after makning it mpi-aware)
+  if( !this->action->runInSerial() ) {
+    this->comm.Sum( forcesForApply );
+  }
+}
+
+} // namespace PLMD
+#endif // __PLUMED_acc_ParallelTaskManager_h
+
+
diff --git a/plugins/openaccPTM/Distance.cpp b/plugins/openaccPTM/Distance.cpp
new file mode 100644
index 0000000000..53b546db38
--- /dev/null
+++ b/plugins/openaccPTM/Distance.cpp
@@ -0,0 +1,34 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2011-2023 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "plumed/core/ActionRegister.h"
+#include "plumed/colvar/MultiColvarTemplate.h"
+#include "plumed/colvar/Distance.h"
+
+#include "ACCParallelTaskManager.h"
+
+template <typename T>
+using MCT= PLMD::colvar::MultiColvarTemplate<T,PLMD::ACCPTM>;
+
+typedef MCT<PLMD::colvar::Distance> DistanceMultiAcc;
+
+PLUMED_REGISTER_ACTION(DistanceMultiAcc,"DISTANCE_VECTOR_OPENACC")
+
diff --git a/plugins/openaccPTM/Makefile b/plugins/openaccPTM/Makefile
new file mode 100644
index 0000000000..b13e8bb8c2
--- /dev/null
+++ b/plugins/openaccPTM/Makefile
@@ -0,0 +1,47 @@
+include ./Makefile.conf
+
+#Dependency tracking based on https://make.mad-scientist.net/papers/advanced-auto-dependency-generation/#tldr
+#this assumes gcc
+DEPDIR := .deps
+DEPFLAGS = -MT $@ -MMD -MP -MF $(DEPDIR)/$*.d
+#tested with nvcc with :"Build cuda_11.7.r11.7/compiler.31442593_0"
+#-dc adds relocatable device code 
+#-dlto Perform link-time optimization of device code. 
+OBJS = Distance.o
+
+ADDCPPFLAGS=$(PLUMED_INCLUDE) -D__PLUMED_HAS_OPENACC=1 -acc -acc="gpu"
+ADDCLDFLAGS=-acc -static-nvidia
+
+ifeq ($(SOEXT),dylib)
+  SONAME_OPTION:=-Xlinker -install_name
+else
+  SONAME_OPTION:=-Xlinker -soname
+endif
+
+PLUGIN=plumedOpenACC.$(SOEXT)
+
+NVCXX ?= nvc++
+NVLDSHARED = $(LDSHARED:$(CXX)=$(NVCXX))
+all: $(PLUGIN)
+
+%.o: %.cpp $(DEPDIR)/%.d | $(DEPDIR)
+	@echo Compiling object $@
+	$(NVCXX) -c $(DEPFLAGS) $(CPPFLAGS) $(ADDCPPFLAGS) $(CXXFLAGS) $< -o $@
+
+$(DEPDIR): ; @mkdir -p $@
+
+DEPFILES := $(OBJS:%.o=$(DEPDIR)/%.d)
+$(DEPFILES):
+include $(wildcard $(DEPFILES))
+
+$(PLUGIN): $(OBJS)
+	@echo Linking $@
+	$(NVLDSHARED)  $(ADDCLDFLAGS) $(SONAME_OPTION),"$(notdir $@)" $(DYNAMIC_LIBS) $(PLUMED_KERNEL) -o $@ $^
+
+clean:
+	@rm -fv $(OBJS) $(PLUGIN)
+	@rm -rf .deps
+
+check: all
+	$(MAKE) -C regtest
+	$(MAKE) -C regtest checkfail
diff --git a/plugins/openaccPTM/configure.sh b/plugins/openaccPTM/configure.sh
new file mode 100755
index 0000000000..f6651959b1
--- /dev/null
+++ b/plugins/openaccPTM/configure.sh
@@ -0,0 +1,30 @@
+#! /usr/bin/env bash
+
+if [[ -z $PLUMED_KERNEL ]]; then
+    echo "$(basename $0) can work only if \"PLUMED_KERNEL\" is defined"
+    echo "either via module load or sourceme.sh"
+fi
+
+{
+    plumed config makefile_conf
+    echo "PLUMED_INCLUDE=-I$(plumed info --include-dir)"
+    echo "PLUMED_KERNEL=-L${PLUMED_KERNEL}"
+} >Make.tmp
+
+#pendantic adds a unuseful FOR EACH line with
+#"" warning: style of line directive is a GCC extension"
+{
+    #    grep CXXFLAGS Make.tmp \
+    #        | sed -e 's/-f/-Xcompiler -f/g' \
+    #            -e 's/-pedantic//g' \
+    #            -e 's/-W/-Xcompiler -W/g'
+    #    grep -eDYNAMIC_LIBS -eLDFLAGS Make.tmp \
+    #        | sed -e 's/-rdynamic/-Xcompiler -rdynamic/g' \
+    #            -e 's/-Wl,/-Xlinker /g' \
+    #            -e 's/-f/-Xcompiler -f/g'
+    #    #prints the rest of the file
+    #    grep -eDYNAMIC_LIBS -eLDFLAGS -eCXXFLAGS Make.tmp -v
+    sed -e 's/-fno-gnu-unique//g' Make.tmp
+} >Makefile.conf
+
+rm Make.tmp
diff --git a/plugins/openaccPTM/plumed b/plugins/openaccPTM/plumed
new file mode 120000
index 0000000000..929cb3dc9b
--- /dev/null
+++ b/plugins/openaccPTM/plumed
@@ -0,0 +1 @@
+../../src
\ No newline at end of file

From 13ca1dd41d3ed0e34bae251b0e45537a62ee228f Mon Sep 17 00:00:00 2001
From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com>
Date: Thu, 28 Aug 2025 16:59:55 +0200
Subject: [PATCH 4/9] not sure if it is a good idea, but the link problem was
 workarounded

---
 plugins/openaccPTM/Makefile | 14 +++++++++++++-
 1 file changed, 13 insertions(+), 1 deletion(-)

diff --git a/plugins/openaccPTM/Makefile b/plugins/openaccPTM/Makefile
index b13e8bb8c2..f3aaf1442d 100644
--- a/plugins/openaccPTM/Makefile
+++ b/plugins/openaccPTM/Makefile
@@ -7,7 +7,7 @@ DEPFLAGS = -MT $@ -MMD -MP -MF $(DEPDIR)/$*.d
 #tested with nvcc with :"Build cuda_11.7.r11.7/compiler.31442593_0"
 #-dc adds relocatable device code 
 #-dlto Perform link-time optimization of device code. 
-OBJS = Distance.o
+OBJS = Distance.o extra/colvar/ColvarInput.o extra/tools/Pbc.o
 
 ADDCPPFLAGS=$(PLUMED_INCLUDE) -D__PLUMED_HAS_OPENACC=1 -acc -acc="gpu"
 ADDCLDFLAGS=-acc -static-nvidia
@@ -28,7 +28,19 @@ all: $(PLUGIN)
 	@echo Compiling object $@
 	$(NVCXX) -c $(DEPFLAGS) $(CPPFLAGS) $(ADDCPPFLAGS) $(CXXFLAGS) $< -o $@
 
+extra/colvar/ColvarInput.o: plumed/colvar/ColvarInput.cpp $(DEPDIR)/extra/colvar/ColvarInput.d | $(DEPDIR)/extra/colvar extra/colvar
+	@echo Compiling object $@
+	$(NVCXX) -c $(DEPFLAGS) $(CPPFLAGS) $(ADDCPPFLAGS) $(CXXFLAGS) $< -o $@
+
+extra/tools/Pbc.o: plumed/tools/Pbc.cpp $(DEPDIR)/extra/tools/Pbc.d | $(DEPDIR)/extra/tools extra/tools
+	@echo Compiling object $@
+	$(NVCXX) -c $(DEPFLAGS) $(CPPFLAGS) $(ADDCPPFLAGS) $(CXXFLAGS) $< -o $@
+
 $(DEPDIR): ; @mkdir -p $@
+$(DEPDIR)/extra/colvar: ; @mkdir -p $@
+$(DEPDIR)/extra/tools: ; @mkdir -p $@
+extra/colvar: ; mkdir -p $@
+extra/tools: ; mkdir -p $@
 
 DEPFILES := $(OBJS:%.o=$(DEPDIR)/%.d)
 $(DEPFILES):

From b533e3ee1b49a5f0aa58b4c7dac8198f50d2e01b Mon Sep 17 00:00:00 2001
From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com>
Date: Fri, 29 Aug 2025 09:36:39 +0200
Subject: [PATCH 5/9] moving back the manual in the Distance.cpp

---
 src/colvar/Distance.cpp | 110 ++++++++++++++++++++++++++++++++++++++++
 src/colvar/Distance.h   | 110 ----------------------------------------
 2 files changed, 110 insertions(+), 110 deletions(-)

diff --git a/src/colvar/Distance.cpp b/src/colvar/Distance.cpp
index b321a6973a..92cd4933f5 100644
--- a/src/colvar/Distance.cpp
+++ b/src/colvar/Distance.cpp
@@ -27,6 +27,116 @@
 namespace PLMD {
 namespace colvar {
 
+//+PLUMEDOC COLVAR DISTANCE
+/*
+Calculate the distance/s between pairs of atoms.
+
+The following example illustrates how this action can be used to calculate and print the distance between atom 1
+and atom 2.
+
+```plumed
+d: DISTANCE ATOMS=1,2
+PRINT ARG=d FILE=colvar
+```
+
+By default the distance is computed in a way that takes periodic
+boundary conditions in account. This behavior can be changed by using the NOPBC flag.
+Furthermore, if you wish to calculate the vector connecting a pair of atoms you can use the
+`COMPONENTS` flag as shown below:
+
+```plumed
+d: DISTANCE ATOMS=1,2 COMPONENTS
+PRINT ARG=d.x,d.y,d.z FILE=colvar
+```
+
+Alternatively, you can calculate the components projected on the lattice vector by using the `SCALED_COMPONENTS`
+flag as shown below;
+
+```plumed
+d: DISTANCE ATOMS=1,2 SCALED_COMPONENTS
+PRINT ARG=d.a,d.b,d.c FILE=colvar
+```
+
+The advantage of using `SCALED_COMPONENTS` over `COMPONENTS` is that the a, b and c variables
+that are calculated when `SCALED_COMPONENTS` is employed have the proper periodicity. This feature is useful
+if you wish to study the motion of a molecule across a membrane.
+
+You can also use this command to calculate multiple indistinguishable distances or vectors with a single
+line of PLUMED input.  For example, the following input calculates and outputs the distances between four
+pairs of atoms:
+
+```plumed
+d: DISTANCE ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 ATOMS4=7,8
+PRINT ARG=d FILE=colvar
+```
+
+By a similar token, the following input outputs three four dimensional vectors that contain the x, y and z
+components of the vectors connecting the four atoms:
+
+```plumed
+d: DISTANCE COMPONENTS ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 ATOMS4=7,8
+PRINT ARG=d.x,d.y,d.z FILE=colvar
+```
+
+You can also replace COMPONENTS with SCALED_COMPONENTS in the above input and obtain the projects of these vectors
+on the lattice vectors.
+
+## Managing periodic boundary conditions
+
+When using the DISTANCE command to calculate the end-to-end distance for a large polymer you need to ensure that you
+are managing PBCs correctly.  This problems that can occur with these calculations are explained at length in the
+early parts of the document that is referenced in the bibliography. Notice, however, that the input
+provides an example of an input that could be used to compute the end-to-end distance for a polymer
+of 100 atoms and keeps it at a value around 5.
+
+```plumed
+WHOLEMOLECULES ENTITY0=1-100
+e2e: DISTANCE ATOMS=1,100 NOPBC
+RESTRAINT ARG=e2e KAPPA=1 AT=5
+```
+
+Notice that NOPBC is used here so as to ensure that the distance is calculated correctely even if the end-to-end distance is larger than half the simulation
+Also notice that, since many MD codes break molecules across cell boundary, it might be necessary to
+use the [WHOLEMOLECULES](WHOLEMOLECULES.md) keyword (also notice that it should be _before_ distance). The list of atoms provided to [WHOLEMOLECULES](WHOLEMOLECULES.md)
+here contains all the atoms between 1 and 100. Strictly speaking, this
+is not necessary. If you know for sure that atoms with difference in
+the index say equal to 10 are _not_ going to be farther than half cell
+you can e.g. use
+
+```plumed
+WHOLEMOLECULES ENTITY0=1,10,20,30,40,50,60,70,80,90,100
+e2e: DISTANCE ATOMS=1,100 NOPBC
+RESTRAINT ARG=e2e KAPPA=1 AT=5
+```
+
+Just be sure that the ordered list provide to [WHOLEMOLECULES](WHOLEMOLECULES.md) has the following
+properties:
+
+- Consecutive atoms should be closer than half-cell throughout the entire simulation.
+- Atoms required later for the distance (e.g. 1 and 100) should be included in the list
+
+The following example shows how to take periodicity into account when computing the z-component of a distance
+
+```plumed
+# this is a center of mass of a large group
+c: COM ATOMS=1-100
+# this is the distance between atom 101 and the group
+d: DISTANCE ATOMS=c,101 COMPONENTS
+# this makes a new variable, dd, equal to d and periodic, with domain -10,10
+# this is the right choise if e.g. the cell is orthorombic and its size in
+# z direction is 20.
+dz: COMBINE ARG=d.z PERIODIC=-10,10
+# metadynamics on dd
+METAD ARG=dz SIGMA=0.1 HEIGHT=0.1 PACE=200
+```
+
+You can use the same input even if the DISTANCE command is calculating the vectors connecting multiple pairs of atoms.
+However, using SCALED_COMPONENTS ensures this problem does not arise because these variables are always periodic
+with domain (-0.5,+0.5).
+
+*/
+//+ENDPLUMEDOC
+
 typedef ColvarShortcut<Distance> DistanceShortcut;
 PLUMED_REGISTER_ACTION(DistanceShortcut,"DISTANCE")
 PLUMED_REGISTER_ACTION(Distance,"DISTANCE_SCALAR")
diff --git a/src/colvar/Distance.h b/src/colvar/Distance.h
index 2c0128cb46..cd4eaaef87 100644
--- a/src/colvar/Distance.h
+++ b/src/colvar/Distance.h
@@ -28,116 +28,6 @@
 namespace PLMD {
 namespace colvar {
 
-//+PLUMEDOC COLVAR DISTANCE
-/*
-Calculate the distance/s between pairs of atoms.
-
-The following example illustrates how this action can be used to calculate and print the distance between atom 1
-and atom 2.
-
-```plumed
-d: DISTANCE ATOMS=1,2
-PRINT ARG=d FILE=colvar
-```
-
-By default the distance is computed in a way that takes periodic
-boundary conditions in account. This behavior can be changed by using the NOPBC flag.
-Furthermore, if you wish to calculate the vector connecting a pair of atoms you can use the
-`COMPONENTS` flag as shown below:
-
-```plumed
-d: DISTANCE ATOMS=1,2 COMPONENTS
-PRINT ARG=d.x,d.y,d.z FILE=colvar
-```
-
-Alternatively, you can calculate the components projected on the lattice vector by using the `SCALED_COMPONENTS`
-flag as shown below;
-
-```plumed
-d: DISTANCE ATOMS=1,2 SCALED_COMPONENTS
-PRINT ARG=d.a,d.b,d.c FILE=colvar
-```
-
-The advantage of using `SCALED_COMPONENTS` over `COMPONENTS` is that the a, b and c variables
-that are calculated when `SCALED_COMPONENTS` is employed have the proper periodicity. This feature is useful
-if you wish to study the motion of a molecule across a membrane.
-
-You can also use this command to calculate multiple indistinguishable distances or vectors with a single
-line of PLUMED input.  For example, the following input calculates and outputs the distances between four
-pairs of atoms:
-
-```plumed
-d: DISTANCE ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 ATOMS4=7,8
-PRINT ARG=d FILE=colvar
-```
-
-By a similar token, the following input outputs three four dimensional vectors that contain the x, y and z
-components of the vectors connecting the four atoms:
-
-```plumed
-d: DISTANCE COMPONENTS ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 ATOMS4=7,8
-PRINT ARG=d.x,d.y,d.z FILE=colvar
-```
-
-You can also replace COMPONENTS with SCALED_COMPONENTS in the above input and obtain the projects of these vectors
-on the lattice vectors.
-
-## Managing periodic boundary conditions
-
-When using the DISTANCE command to calculate the end-to-end distance for a large polymer you need to ensure that you
-are managing PBCs correctly.  This problems that can occur with these calculations are explained at length in the
-early parts of the document that is referenced in the bibliography. Notice, however, that the input
-provides an example of an input that could be used to compute the end-to-end distance for a polymer
-of 100 atoms and keeps it at a value around 5.
-
-```plumed
-WHOLEMOLECULES ENTITY0=1-100
-e2e: DISTANCE ATOMS=1,100 NOPBC
-RESTRAINT ARG=e2e KAPPA=1 AT=5
-```
-
-Notice that NOPBC is used here so as to ensure that the distance is calculated correctely even if the end-to-end distance is larger than half the simulation
-Also notice that, since many MD codes break molecules across cell boundary, it might be necessary to
-use the [WHOLEMOLECULES](WHOLEMOLECULES.md) keyword (also notice that it should be _before_ distance). The list of atoms provided to [WHOLEMOLECULES](WHOLEMOLECULES.md)
-here contains all the atoms between 1 and 100. Strictly speaking, this
-is not necessary. If you know for sure that atoms with difference in
-the index say equal to 10 are _not_ going to be farther than half cell
-you can e.g. use
-
-```plumed
-WHOLEMOLECULES ENTITY0=1,10,20,30,40,50,60,70,80,90,100
-e2e: DISTANCE ATOMS=1,100 NOPBC
-RESTRAINT ARG=e2e KAPPA=1 AT=5
-```
-
-Just be sure that the ordered list provide to [WHOLEMOLECULES](WHOLEMOLECULES.md) has the following
-properties:
-
-- Consecutive atoms should be closer than half-cell throughout the entire simulation.
-- Atoms required later for the distance (e.g. 1 and 100) should be included in the list
-
-The following example shows how to take periodicity into account when computing the z-component of a distance
-
-```plumed
-# this is a center of mass of a large group
-c: COM ATOMS=1-100
-# this is the distance between atom 101 and the group
-d: DISTANCE ATOMS=c,101 COMPONENTS
-# this makes a new variable, dd, equal to d and periodic, with domain -10,10
-# this is the right choise if e.g. the cell is orthorombic and its size in
-# z direction is 20.
-dz: COMBINE ARG=d.z PERIODIC=-10,10
-# metadynamics on dd
-METAD ARG=dz SIGMA=0.1 HEIGHT=0.1 PACE=200
-```
-
-You can use the same input even if the DISTANCE command is calculating the vectors connecting multiple pairs of atoms.
-However, using SCALED_COMPONENTS ensures this problem does not arise because these variables are always periodic
-with domain (-0.5,+0.5).
-
-*/
-//+ENDPLUMEDOC
-
 class Distance : public Colvar {
   bool components;
   bool scaled_components;

From 933d393a6f5e721bb1ea523d86b9b988ce7e6625 Mon Sep 17 00:00:00 2001
From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com>
Date: Fri, 29 Aug 2025 14:29:31 +0200
Subject: [PATCH 6/9] deactivating a check in ActionShortcut::readInputLine

---
 src/core/ActionShortcut.cpp | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/src/core/ActionShortcut.cpp b/src/core/ActionShortcut.cpp
index ab03e9d01e..b694cd509a 100644
--- a/src/core/ActionShortcut.cpp
+++ b/src/core/ActionShortcut.cpp
@@ -149,7 +149,12 @@ void ActionShortcut::readInputLine( const std::string& input, bool saveline ) {
       Keywords thiskeys;
       actionRegister().getKeywords( av->getName(), thiskeys );
       if( thiskeys.getDisplayName()!=getName() ) {
-        error("mismatch between display name of hidden action " + thiskeys.getDisplayName() + " and shortcut that creates it " + getName() );
+        //NOTE: this was an error
+        // I deactivated it because it clashes with the gpu-openacc-by-LOAD does not load
+        // at moment of writing this I do not know if there is a way to to the check by passing the
+        // dlloader object to this function.
+        // bypassing this check makes my contraption work
+        warning("mismatch between display name of hidden action " + thiskeys.getDisplayName() + " and shortcut that creates it " + getName() );
       }
     }
   } else {

From 63a90c5f05c590a5e139d204ad0691cfa057fb4b Mon Sep 17 00:00:00 2001
From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com>
Date: Fri, 29 Aug 2025 14:37:01 +0200
Subject: [PATCH 7/9] setting up USEGPU for distance_vector

---
 plugins/openaccPTM/.gitignore   |  2 ++
 plugins/openaccPTM/Distance.cpp |  2 +-
 src/colvar/ColvarShortcut.h     | 10 +++++++++-
 3 files changed, 12 insertions(+), 2 deletions(-)
 create mode 100644 plugins/openaccPTM/.gitignore

diff --git a/plugins/openaccPTM/.gitignore b/plugins/openaccPTM/.gitignore
new file mode 100644
index 0000000000..92dd77603a
--- /dev/null
+++ b/plugins/openaccPTM/.gitignore
@@ -0,0 +1,2 @@
+Makefile.conf
+
diff --git a/plugins/openaccPTM/Distance.cpp b/plugins/openaccPTM/Distance.cpp
index 53b546db38..0c8d49bea2 100644
--- a/plugins/openaccPTM/Distance.cpp
+++ b/plugins/openaccPTM/Distance.cpp
@@ -30,5 +30,5 @@ using MCT= PLMD::colvar::MultiColvarTemplate<T,PLMD::ACCPTM>;
 
 typedef MCT<PLMD::colvar::Distance> DistanceMultiAcc;
 
-PLUMED_REGISTER_ACTION(DistanceMultiAcc,"DISTANCE_VECTOR_OPENACC")
+PLUMED_REGISTER_ACTION(DistanceMultiAcc,"DISTANCE_VECTORACC")
 
diff --git a/src/colvar/ColvarShortcut.h b/src/colvar/ColvarShortcut.h
index 0bbc799d33..8cd8918445 100644
--- a/src/colvar/ColvarShortcut.h
+++ b/src/colvar/ColvarShortcut.h
@@ -45,6 +45,10 @@ void ColvarShortcut<T>::registerKeywords(Keywords& keys ) {
   }
   keys.addActionNameSuffix("_SCALAR");
   keys.addActionNameSuffix("_VECTOR");
+  //GPU related settings
+  keys.addFlag("USEGPU",false,"run this calculation on the GPU");
+  keys.addLinkInDocForFlag("USEGPU","gpu.md");
+  keys.addActionNameSuffix("_VECTORACC");
 }
 
 template <class T>
@@ -64,8 +68,12 @@ ColvarShortcut<T>::ColvarShortcut(const ActionOptions&ao):
     if( keywords.style( key, "atoms" ) ) {
       std::string inpt;
       parseNumbered( key, 1, inpt );
+      bool usegpuFLAG=false;
+      parseFlag("USEGPU",usegpuFLAG);
       if( inpt.length()>0 ) {
-        readInputLine( getShortcutLabel() + ": " + getName() + "_VECTOR " + key + "1=" + inpt + " " + convertInputLineToString() );
+        readInputLine( getShortcutLabel() + ": "
+                       + getName()  + "_VECTOR" + (usegpuFLAG ? "ACC ":" ")
+                       + key + "1=" + inpt + " " + convertInputLineToString() );
         scalar=false;
         break;
       }

From 93d5dfdc7ee6be92f5fded0e5f380429d4d3297c Mon Sep 17 00:00:00 2001
From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com>
Date: Fri, 29 Aug 2025 16:19:47 +0200
Subject: [PATCH 8/9] moved the multicolvar with openacc in a plugin

---
 plugins/openaccPTM/.gitignore              |   2 +
 plugins/openaccPTM/Angle.cpp               |  30 +++
 plugins/openaccPTM/DihedralCorrelation.cpp |  30 +++
 plugins/openaccPTM/Dipole.cpp              |  30 +++
 plugins/openaccPTM/Distance.cpp            |   6 +-
 plugins/openaccPTM/Makefile                |  10 +-
 plugins/openaccPTM/Plane.cpp               |  30 +++
 plugins/openaccPTM/Position.cpp            |  30 +++
 plugins/openaccPTM/Torsion.cpp             |  30 +++
 src/colvar/Angle.cpp                       | 105 +----------
 src/colvar/Angle.h                         | 125 +++++++++++++
 src/colvar/DihedralCorrelation.cpp         | 128 +------------
 src/colvar/DihedralCorrelation.h           | 154 ++++++++++++++++
 src/colvar/Dipole.cpp                      | 159 +---------------
 src/colvar/Dipole.h                        | 184 +++++++++++++++++++
 src/colvar/Plane.cpp                       | 121 +-----------
 src/colvar/Plane.h                         | 146 +++++++++++++++
 src/colvar/Position.cpp                    | 158 +---------------
 src/colvar/Position.h                      | 186 +++++++++++++++++++
 src/colvar/Torsion.cpp                     | 180 +-----------------
 src/colvar/Torsion.h                       | 204 +++++++++++++++++++++
 21 files changed, 1199 insertions(+), 849 deletions(-)
 create mode 100644 plugins/openaccPTM/Angle.cpp
 create mode 100644 plugins/openaccPTM/DihedralCorrelation.cpp
 create mode 100644 plugins/openaccPTM/Dipole.cpp
 create mode 100644 plugins/openaccPTM/Plane.cpp
 create mode 100644 plugins/openaccPTM/Position.cpp
 create mode 100644 plugins/openaccPTM/Torsion.cpp
 create mode 100644 src/colvar/Angle.h
 create mode 100644 src/colvar/DihedralCorrelation.h
 create mode 100644 src/colvar/Dipole.h
 create mode 100644 src/colvar/Plane.h
 create mode 100644 src/colvar/Position.h
 create mode 100644 src/colvar/Torsion.h

diff --git a/plugins/openaccPTM/.gitignore b/plugins/openaccPTM/.gitignore
index 92dd77603a..f11b086c1d 100644
--- a/plugins/openaccPTM/.gitignore
+++ b/plugins/openaccPTM/.gitignore
@@ -1,2 +1,4 @@
 Makefile.conf
+extra/
+.deps/
 
diff --git a/plugins/openaccPTM/Angle.cpp b/plugins/openaccPTM/Angle.cpp
new file mode 100644
index 0000000000..4779e6b121
--- /dev/null
+++ b/plugins/openaccPTM/Angle.cpp
@@ -0,0 +1,30 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2011-2023 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "plumed/core/ActionRegister.h"
+#include "plumed/colvar/MultiColvarTemplate.h"
+#include "plumed/colvar/Angle.h"
+
+#include "ACCParallelTaskManager.h"
+
+typedef PLMD::colvar::MultiColvarTemplate<PLMD::colvar::Angle,PLMD::ACCPTM> AngleMultiAcc;
+
+PLUMED_REGISTER_ACTION(AngleMultiAcc,"ANGLE_VECTORACC")
diff --git a/plugins/openaccPTM/DihedralCorrelation.cpp b/plugins/openaccPTM/DihedralCorrelation.cpp
new file mode 100644
index 0000000000..f77585c304
--- /dev/null
+++ b/plugins/openaccPTM/DihedralCorrelation.cpp
@@ -0,0 +1,30 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2011-2023 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "plumed/core/ActionRegister.h"
+#include "plumed/colvar/MultiColvarTemplate.h"
+#include "plumed/colvar/DihedralCorrelation.h"
+
+#include "ACCParallelTaskManager.h"
+
+typedef PLMD::colvar::MultiColvarTemplate<PLMD::colvar::DihedralCorrelation,PLMD::ACCPTM> DihedralCorrelationMultiAcc;
+
+PLUMED_REGISTER_ACTION(DihedralCorrelationMultiAcc,"DIHEDRALCORRELATION_VECTORACC")
diff --git a/plugins/openaccPTM/Dipole.cpp b/plugins/openaccPTM/Dipole.cpp
new file mode 100644
index 0000000000..af5801e1a0
--- /dev/null
+++ b/plugins/openaccPTM/Dipole.cpp
@@ -0,0 +1,30 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2011-2023 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "plumed/core/ActionRegister.h"
+#include "plumed/colvar/MultiColvarTemplate.h"
+#include "plumed/colvar/Dipole.h"
+
+#include "ACCParallelTaskManager.h"
+
+typedef PLMD::colvar::MultiColvarTemplate<PLMD::colvar::Dipole,PLMD::ACCPTM> DipoleMultiAcc;
+
+PLUMED_REGISTER_ACTION(DipoleMultiAcc,"DIPOLE_VECTORACC")
diff --git a/plugins/openaccPTM/Distance.cpp b/plugins/openaccPTM/Distance.cpp
index 0c8d49bea2..03d4c21dfc 100644
--- a/plugins/openaccPTM/Distance.cpp
+++ b/plugins/openaccPTM/Distance.cpp
@@ -25,10 +25,6 @@
 
 #include "ACCParallelTaskManager.h"
 
-template <typename T>
-using MCT= PLMD::colvar::MultiColvarTemplate<T,PLMD::ACCPTM>;
-
-typedef MCT<PLMD::colvar::Distance> DistanceMultiAcc;
+typedef PLMD::colvar::MultiColvarTemplate<PLMD::colvar::Distance,PLMD::ACCPTM> DistanceMultiAcc;
 
 PLUMED_REGISTER_ACTION(DistanceMultiAcc,"DISTANCE_VECTORACC")
-
diff --git a/plugins/openaccPTM/Makefile b/plugins/openaccPTM/Makefile
index f3aaf1442d..8d2e9391da 100644
--- a/plugins/openaccPTM/Makefile
+++ b/plugins/openaccPTM/Makefile
@@ -7,7 +7,7 @@ DEPFLAGS = -MT $@ -MMD -MP -MF $(DEPDIR)/$*.d
 #tested with nvcc with :"Build cuda_11.7.r11.7/compiler.31442593_0"
 #-dc adds relocatable device code 
 #-dlto Perform link-time optimization of device code. 
-OBJS = Distance.o extra/colvar/ColvarInput.o extra/tools/Pbc.o
+OBJS = Angle.o DihedralCorrelation.o Dipole.o Distance.o Plane.o Position.o Torsion.o extra/colvar/ColvarInput.o extra/tools/Pbc.o extra/tools/Torsion.o extra/tools/Angle.o
 
 ADDCPPFLAGS=$(PLUMED_INCLUDE) -D__PLUMED_HAS_OPENACC=1 -acc -acc="gpu"
 ADDCLDFLAGS=-acc -static-nvidia
@@ -36,6 +36,14 @@ extra/tools/Pbc.o: plumed/tools/Pbc.cpp $(DEPDIR)/extra/tools/Pbc.d | $(DEPDIR)/
 	@echo Compiling object $@
 	$(NVCXX) -c $(DEPFLAGS) $(CPPFLAGS) $(ADDCPPFLAGS) $(CXXFLAGS) $< -o $@
 
+extra/tools/Torsion.o: plumed/tools/Torsion.cpp $(DEPDIR)/extra/tools/Torsion.d | $(DEPDIR)/extra/tools extra/tools
+	@echo Compiling object $@
+	$(NVCXX) -c $(DEPFLAGS) $(CPPFLAGS) $(ADDCPPFLAGS) $(CXXFLAGS) $< -o $@
+
+extra/tools/Angle.o: plumed/tools/Angle.cpp $(DEPDIR)/extra/tools/Angle.d | $(DEPDIR)/extra/tools extra/tools
+	@echo Compiling object $@
+	$(NVCXX) -c $(DEPFLAGS) $(CPPFLAGS) $(ADDCPPFLAGS) $(CXXFLAGS) $< -o $@
+
 $(DEPDIR): ; @mkdir -p $@
 $(DEPDIR)/extra/colvar: ; @mkdir -p $@
 $(DEPDIR)/extra/tools: ; @mkdir -p $@
diff --git a/plugins/openaccPTM/Plane.cpp b/plugins/openaccPTM/Plane.cpp
new file mode 100644
index 0000000000..a771f4980a
--- /dev/null
+++ b/plugins/openaccPTM/Plane.cpp
@@ -0,0 +1,30 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2011-2023 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "plumed/core/ActionRegister.h"
+#include "plumed/colvar/MultiColvarTemplate.h"
+#include "plumed/colvar/Plane.h"
+
+#include "ACCParallelTaskManager.h"
+
+typedef PLMD::colvar::MultiColvarTemplate<PLMD::colvar::Plane,PLMD::ACCPTM> PlaneMultiAcc;
+
+PLUMED_REGISTER_ACTION(PlaneMultiAcc,"PLANE_VECTORACC")
diff --git a/plugins/openaccPTM/Position.cpp b/plugins/openaccPTM/Position.cpp
new file mode 100644
index 0000000000..41bd8e6e44
--- /dev/null
+++ b/plugins/openaccPTM/Position.cpp
@@ -0,0 +1,30 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2011-2023 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "plumed/core/ActionRegister.h"
+#include "plumed/colvar/MultiColvarTemplate.h"
+#include "plumed/colvar/Position.h"
+
+#include "ACCParallelTaskManager.h"
+
+typedef PLMD::colvar::MultiColvarTemplate<PLMD::colvar::Position,PLMD::ACCPTM> PositionMultiAcc;
+
+PLUMED_REGISTER_ACTION(PositionMultiAcc,"POSITION_VECTORACC")
diff --git a/plugins/openaccPTM/Torsion.cpp b/plugins/openaccPTM/Torsion.cpp
new file mode 100644
index 0000000000..709680057c
--- /dev/null
+++ b/plugins/openaccPTM/Torsion.cpp
@@ -0,0 +1,30 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2011-2023 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "plumed/core/ActionRegister.h"
+#include "plumed/colvar/MultiColvarTemplate.h"
+#include "plumed/colvar/Torsion.h"
+
+#include "ACCParallelTaskManager.h"
+
+typedef PLMD::colvar::MultiColvarTemplate<PLMD::colvar::Torsion,PLMD::ACCPTM> TorsionMultiAcc;
+
+PLUMED_REGISTER_ACTION(TorsionMultiAcc,"TORSION_VECTORACC")
diff --git a/src/colvar/Angle.cpp b/src/colvar/Angle.cpp
index b0e0866725..473cd5addd 100644
--- a/src/colvar/Angle.cpp
+++ b/src/colvar/Angle.cpp
@@ -19,10 +19,7 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#ifdef __PLUMED_HAS_OPENACC
-#define __PLUMED_USE_OPENACC 1
-#endif //__PLUMED_HAS_OPENACC
-#include "Colvar.h"
+#include "Angle.h"
 #include "ColvarShortcut.h"
 #include "MultiColvarTemplate.h"
 #include "core/ActionRegister.h"
@@ -100,108 +97,12 @@ PRINT ARG=a4,a5,a6 FILE=colvar
 */
 //+ENDPLUMEDOC
 
-class Angle : public Colvar {
-  bool pbc;
-  std::vector<double> value;
-  std::vector<double> derivs;
-public:
-  explicit Angle(const ActionOptions&);
-// active methods:
-  void calculate() override;
-  static void registerKeywords( Keywords& keys );
-  static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
-  static unsigned getModeAndSetupValues( ActionWithValue* av );
-  static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout );
-};
-
 typedef ColvarShortcut<Angle> AngleShortcut;
 PLUMED_REGISTER_ACTION(AngleShortcut,"ANGLE")
 PLUMED_REGISTER_ACTION(Angle,"ANGLE_SCALAR")
 typedef MultiColvarTemplate<Angle> AngleMulti;
 PLUMED_REGISTER_ACTION(AngleMulti,"ANGLE_VECTOR")
 
-void Angle::registerKeywords( Keywords& keys ) {
-  Colvar::registerKeywords(keys);
-  keys.setDisplayName("ANGLE");
-  keys.add("atoms","ATOMS","the list of atoms involved in this collective variable (either 3 or 4 atoms)");
-  keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
-  keys.setValueDescription("scalar/vector","the ANGLE involving these atoms");
-  keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
-}
-
-void Angle::parseAtomList( const int& num, std::vector<AtomNumber>& atoms, ActionAtomistic* aa ) {
-  aa->parseAtomList("ATOMS",num,atoms);
-  if(atoms.size()==3) {
-    aa->log.printf("  between atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial());
-    atoms.resize(4);
-    atoms[3]=atoms[2];
-    atoms[2]=atoms[1];
-  } else if(atoms.size()==4) {
-    aa->log.printf("  between lines %d-%d and %d-%d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial());
-  } else if( num<0 || atoms.size()>0 ) {
-    aa->error("Number of specified atoms should be either 3 or 4");
-  }
-}
-
-unsigned Angle::getModeAndSetupValues( ActionWithValue* av ) {
-  av->addValueWithDerivatives();
-  av->setNotPeriodic();
-  return 0;
-}
-
-Angle::Angle(const ActionOptions&ao):
-  PLUMED_COLVAR_INIT(ao),
-  pbc(true),
-  value(1) {
-  std::vector<AtomNumber> atoms;
-  parseAtomList( -1, atoms, this );
-  bool nopbc=!pbc;
-  parseFlag("NOPBC",nopbc);
-  pbc=!nopbc;
-
-  if(pbc) {
-    log.printf("  using periodic boundary conditions\n");
-  } else {
-    log.printf("  without periodic boundary conditions\n");
-  }
-
-  addValueWithDerivatives();
-  setNotPeriodic();
-  requestAtoms(atoms);
-  checkRead();
-}
-
-// calculator
-void Angle::calculate() {
-
-  if(pbc) {
-    makeWhole();
-  }
-  ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this);
-  calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout );
-  setValue( value[0] );
-  for(unsigned i=0; i<getPositions().size(); ++i) {
-    setAtomsDerivatives( i, cvout.getAtomDerivatives(0,i) );
-  }
-  setBoxDerivatives( cvout.virial[0] );
-}
-
-void Angle::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
-  Vector dij,dik;
-  dij=delta(cvin.pos[2],cvin.pos[3]);
-  dik=delta(cvin.pos[1],cvin.pos[0]);
-  Vector ddij,ddik;
-  PLMD::Angle a;
-  cvout.values[0]=a.compute(dij,dik,ddij,ddik);
-  cvout.derivs[0][0]=ddik;
-  cvout.derivs[0][1]=-ddik;
-  cvout.derivs[0][2]=-ddij;
-  cvout.derivs[0][3]=ddij;
-  ColvarInput::setBoxDerivativesNoPbc( cvin, cvout );
-}
-
-}
-}
-
-
+} //namespace colvar
+} //namespace PLMD
 
diff --git a/src/colvar/Angle.h b/src/colvar/Angle.h
new file mode 100644
index 0000000000..6d9d0981be
--- /dev/null
+++ b/src/colvar/Angle.h
@@ -0,0 +1,125 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2011-2023 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "Colvar.h"
+#include "ColvarInput.h"
+#include "tools/Angle.h"
+
+namespace PLMD {
+namespace colvar {
+
+class Angle : public Colvar {
+  bool pbc;
+  std::vector<double> value;
+  std::vector<double> derivs;
+public:
+  explicit Angle(const ActionOptions&);
+// active methods:
+  void calculate() override;
+  static void registerKeywords( Keywords& keys );
+  static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
+  static unsigned getModeAndSetupValues( ActionWithValue* av );
+  static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout );
+};
+
+void Angle::registerKeywords( Keywords& keys ) {
+  Colvar::registerKeywords(keys);
+  keys.setDisplayName("ANGLE");
+  keys.add("atoms","ATOMS","the list of atoms involved in this collective variable (either 3 or 4 atoms)");
+  keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
+  keys.setValueDescription("scalar/vector","the ANGLE involving these atoms");
+  keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
+}
+
+void Angle::parseAtomList( const int& num, std::vector<AtomNumber>& atoms, ActionAtomistic* aa ) {
+  aa->parseAtomList("ATOMS",num,atoms);
+  if(atoms.size()==3) {
+    aa->log.printf("  between atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial());
+    atoms.resize(4);
+    atoms[3]=atoms[2];
+    atoms[2]=atoms[1];
+  } else if(atoms.size()==4) {
+    aa->log.printf("  between lines %d-%d and %d-%d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial());
+  } else if( num<0 || atoms.size()>0 ) {
+    aa->error("Number of specified atoms should be either 3 or 4");
+  }
+}
+
+unsigned Angle::getModeAndSetupValues( ActionWithValue* av ) {
+  av->addValueWithDerivatives();
+  av->setNotPeriodic();
+  return 0;
+}
+
+Angle::Angle(const ActionOptions&ao):
+  PLUMED_COLVAR_INIT(ao),
+  pbc(true),
+  value(1) {
+  std::vector<AtomNumber> atoms;
+  parseAtomList( -1, atoms, this );
+  bool nopbc=!pbc;
+  parseFlag("NOPBC",nopbc);
+  pbc=!nopbc;
+
+  if(pbc) {
+    log.printf("  using periodic boundary conditions\n");
+  } else {
+    log.printf("  without periodic boundary conditions\n");
+  }
+
+  addValueWithDerivatives();
+  setNotPeriodic();
+  requestAtoms(atoms);
+  checkRead();
+}
+
+// calculator
+void Angle::calculate() {
+
+  if(pbc) {
+    makeWhole();
+  }
+  ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this);
+  calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout );
+  setValue( value[0] );
+  for(unsigned i=0; i<getPositions().size(); ++i) {
+    setAtomsDerivatives( i, cvout.getAtomDerivatives(0,i) );
+  }
+  setBoxDerivatives( cvout.virial[0] );
+}
+
+void Angle::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
+  Vector dij,dik;
+  dij=delta(cvin.pos[2],cvin.pos[3]);
+  dik=delta(cvin.pos[1],cvin.pos[0]);
+  Vector ddij,ddik;
+  PLMD::Angle a;
+  cvout.values[0]=a.compute(dij,dik,ddij,ddik);
+  cvout.derivs[0][0]=ddik;
+  cvout.derivs[0][1]=-ddik;
+  cvout.derivs[0][2]=-ddij;
+  cvout.derivs[0][3]=ddij;
+  ColvarInput::setBoxDerivativesNoPbc( cvin, cvout );
+}
+
+} //namespace colvar
+} //namespace PLMD
+
diff --git a/src/colvar/DihedralCorrelation.cpp b/src/colvar/DihedralCorrelation.cpp
index b5a93ff599..65cb5598ab 100644
--- a/src/colvar/DihedralCorrelation.cpp
+++ b/src/colvar/DihedralCorrelation.cpp
@@ -19,18 +19,12 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#ifdef __PLUMED_HAS_OPENACC
-#define __PLUMED_USE_OPENACC 1
-#endif //__PLUMED_HAS_OPENACC
-#include "Colvar.h"
+#include "DihedralCorrelation.h"
 #include "ColvarShortcut.h"
 #include "MultiColvarTemplate.h"
 #include "tools/Torsion.h"
 #include "core/ActionRegister.h"
 
-#include <string>
-#include <cmath>
-
 namespace PLMD {
 namespace colvar {
 
@@ -80,19 +74,6 @@ PRINT ARG=d FILE=colvar
 */
 //+ENDPLUMEDOC
 
-class DihedralCorrelation : public Colvar {
-private:
-  bool pbc;
-  std::vector<double> value;
-  std::vector<double> derivs;
-public:
-  static void registerKeywords( Keywords& keys );
-  explicit DihedralCorrelation(const ActionOptions&);
-  static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
-  static unsigned getModeAndSetupValues( ActionWithValue* av );
-  void calculate() override;
-  static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout );
-};
 
 typedef ColvarShortcut<DihedralCorrelation> DihedralCorrelationShortcut;
 PLUMED_REGISTER_ACTION(DihedralCorrelationShortcut,"DIHEDRAL_CORRELATION")
@@ -100,112 +81,5 @@ PLUMED_REGISTER_ACTION(DihedralCorrelation,"DIHEDRAL_CORRELATION_SCALAR")
 typedef MultiColvarTemplate<DihedralCorrelation> DihedralCorrelationMulti;
 PLUMED_REGISTER_ACTION(DihedralCorrelationMulti,"DIHEDRAL_CORRELATION_VECTOR")
 
-void DihedralCorrelation::registerKeywords( Keywords& keys ) {
-  Colvar::registerKeywords( keys );
-  keys.setDisplayName("DIHEDRAL_CORRELATION");
-  keys.add("atoms","ATOMS","the set of 8 atoms that are being used to calculate this quantity");
-  keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
-  keys.setValueDescription("scalar/vector","the DIHEDRAL_CORRELATION for these atoms");
-  keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
-}
-
-DihedralCorrelation::DihedralCorrelation(const ActionOptions&ao):
-  PLUMED_COLVAR_INIT(ao),
-  pbc(true),
-  value(1) {
-  std::vector<AtomNumber> atoms;
-  parseAtomList(-1,atoms,this);
-  if( atoms.size()!=8 ) {
-    error("Number of specified atoms should be 8");
-  }
-
-  bool nopbc=!pbc;
-  parseFlag("NOPBC",nopbc);
-  pbc=!nopbc;
-
-  if(pbc) {
-    log.printf("  using periodic boundary conditions\n");
-  } else {
-    log.printf("  without periodic boundary conditions\n");
-  }
-  addValueWithDerivatives();
-  setNotPeriodic();
-}
-
-void DihedralCorrelation::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
-  aa->parseAtomList("ATOMS",num,t);
-  if( num<0 && t.size()!=8 ) {
-    aa->error("Number of specified atoms should be 8");
-  }
-  if( t.size()==8 ) {
-    aa->log.printf("  correlation between dihedral angle for atoms %d %d %d %d and atoms %d %d %d %d\n",
-                   t[0].serial(),t[1].serial(),t[2].serial(),t[3].serial(),t[4].serial(),t[5].serial(),t[6].serial(),t[7].serial());
-  }
-}
-
-unsigned DihedralCorrelation::getModeAndSetupValues( ActionWithValue* av ) {
-  av->addValueWithDerivatives();
-  av->setNotPeriodic();
-  return 0;
-}
-
-void DihedralCorrelation::calculate() {
-
-  if(pbc) {
-    makeWhole();
-  }
-  ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this);
-  calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout );
-  setValue( value[0] );
-  for(unsigned i=0; i<getPositions().size(); ++i) {
-    setAtomsDerivatives( i, cvout.getAtomDerivatives(0, i) );
-  }
-  setBoxDerivatives( cvout.virial[0] );
-}
-
-void DihedralCorrelation::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
-  const Vector d10=delta(cvin.pos[1],cvin.pos[0]);
-  const Vector d11=delta(cvin.pos[2],cvin.pos[1]);
-  const Vector d12=delta(cvin.pos[3],cvin.pos[2]);
-
-  Vector dd10,dd11,dd12;
-  PLMD::Torsion t1;
-  const double phi1  = t1.compute(d10,d11,d12,dd10,dd11,dd12);
-
-  const Vector d20=delta(cvin.pos[5],cvin.pos[4]);
-  const Vector d21=delta(cvin.pos[6],cvin.pos[5]);
-  const Vector d22=delta(cvin.pos[7],cvin.pos[6]);
-
-  Vector dd20,dd21,dd22;
-  PLMD::Torsion t2;
-  const double phi2 = t2.compute( d20, d21, d22, dd20, dd21, dd22 );
-
-  // Calculate value
-  const double diff = phi2 - phi1;
-  cvout.values[0] = 0.5*(1.+std::cos(diff));
-  // Derivatives wrt phi1
-  const double dval = 0.5*std::sin(diff);
-  dd10 *= dval;
-  dd11 *= dval;
-  dd12 *= dval;
-  // And add
-  cvout.derivs[0][0]=dd10;
-  cvout.derivs[0][1]=dd11-dd10;
-  cvout.derivs[0][2]=dd12-dd11;
-  cvout.derivs[0][3]=-dd12;
-  // Derivative wrt phi2
-  dd20 *= -dval;
-  dd21 *= -dval;
-  dd22 *= -dval;
-  // And add
-  cvout.derivs[0][4]=dd20;
-  cvout.derivs[0][5]=dd21-dd20;
-  cvout.derivs[0][6]=dd22-dd21;
-  cvout.derivs[0][7]=-dd22;
-  cvout.virial.set( 0,
-                    -(extProduct(d10,dd10)+extProduct(d11,dd11)+extProduct(d12,dd12))
-                    - (extProduct(d20,dd20)+extProduct(d21,dd21)+extProduct(d22,dd22)) );
-}
-
 }
 }
diff --git a/src/colvar/DihedralCorrelation.h b/src/colvar/DihedralCorrelation.h
new file mode 100644
index 0000000000..0c6310ec4f
--- /dev/null
+++ b/src/colvar/DihedralCorrelation.h
@@ -0,0 +1,154 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2013-2020 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "Colvar.h"
+#include "ColvarInput.h"
+#include "tools/Torsion.h"
+
+#include <string>
+#include <cmath>
+
+namespace PLMD {
+namespace colvar {
+
+class DihedralCorrelation : public Colvar {
+private:
+  bool pbc;
+  std::vector<double> value;
+  std::vector<double> derivs;
+public:
+  static void registerKeywords( Keywords& keys );
+  explicit DihedralCorrelation(const ActionOptions&);
+  static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
+  static unsigned getModeAndSetupValues( ActionWithValue* av );
+  void calculate() override;
+  static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout );
+};
+
+void DihedralCorrelation::registerKeywords( Keywords& keys ) {
+  Colvar::registerKeywords( keys );
+  keys.setDisplayName("DIHEDRAL_CORRELATION");
+  keys.add("atoms","ATOMS","the set of 8 atoms that are being used to calculate this quantity");
+  keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
+  keys.setValueDescription("scalar/vector","the DIHEDRAL_CORRELATION for these atoms");
+  keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
+}
+
+DihedralCorrelation::DihedralCorrelation(const ActionOptions&ao):
+  PLUMED_COLVAR_INIT(ao),
+  pbc(true),
+  value(1) {
+  std::vector<AtomNumber> atoms;
+  parseAtomList(-1,atoms,this);
+  if( atoms.size()!=8 ) {
+    error("Number of specified atoms should be 8");
+  }
+
+  bool nopbc=!pbc;
+  parseFlag("NOPBC",nopbc);
+  pbc=!nopbc;
+
+  if(pbc) {
+    log.printf("  using periodic boundary conditions\n");
+  } else {
+    log.printf("  without periodic boundary conditions\n");
+  }
+  addValueWithDerivatives();
+  setNotPeriodic();
+}
+
+void DihedralCorrelation::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
+  aa->parseAtomList("ATOMS",num,t);
+  if( num<0 && t.size()!=8 ) {
+    aa->error("Number of specified atoms should be 8");
+  }
+  if( t.size()==8 ) {
+    aa->log.printf("  correlation between dihedral angle for atoms %d %d %d %d and atoms %d %d %d %d\n",
+                   t[0].serial(),t[1].serial(),t[2].serial(),t[3].serial(),t[4].serial(),t[5].serial(),t[6].serial(),t[7].serial());
+  }
+}
+
+unsigned DihedralCorrelation::getModeAndSetupValues( ActionWithValue* av ) {
+  av->addValueWithDerivatives();
+  av->setNotPeriodic();
+  return 0;
+}
+
+void DihedralCorrelation::calculate() {
+
+  if(pbc) {
+    makeWhole();
+  }
+  ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this);
+  calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout );
+  setValue( value[0] );
+  for(unsigned i=0; i<getPositions().size(); ++i) {
+    setAtomsDerivatives( i, cvout.getAtomDerivatives(0, i) );
+  }
+  setBoxDerivatives( cvout.virial[0] );
+}
+
+void DihedralCorrelation::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
+  const Vector d10=delta(cvin.pos[1],cvin.pos[0]);
+  const Vector d11=delta(cvin.pos[2],cvin.pos[1]);
+  const Vector d12=delta(cvin.pos[3],cvin.pos[2]);
+
+  Vector dd10,dd11,dd12;
+  PLMD::Torsion t1;
+  const double phi1  = t1.compute(d10,d11,d12,dd10,dd11,dd12);
+
+  const Vector d20=delta(cvin.pos[5],cvin.pos[4]);
+  const Vector d21=delta(cvin.pos[6],cvin.pos[5]);
+  const Vector d22=delta(cvin.pos[7],cvin.pos[6]);
+
+  Vector dd20,dd21,dd22;
+  PLMD::Torsion t2;
+  const double phi2 = t2.compute( d20, d21, d22, dd20, dd21, dd22 );
+
+  // Calculate value
+  const double diff = phi2 - phi1;
+  cvout.values[0] = 0.5*(1.+std::cos(diff));
+  // Derivatives wrt phi1
+  const double dval = 0.5*std::sin(diff);
+  dd10 *= dval;
+  dd11 *= dval;
+  dd12 *= dval;
+  // And add
+  cvout.derivs[0][0]=dd10;
+  cvout.derivs[0][1]=dd11-dd10;
+  cvout.derivs[0][2]=dd12-dd11;
+  cvout.derivs[0][3]=-dd12;
+  // Derivative wrt phi2
+  dd20 *= -dval;
+  dd21 *= -dval;
+  dd22 *= -dval;
+  // And add
+  cvout.derivs[0][4]=dd20;
+  cvout.derivs[0][5]=dd21-dd20;
+  cvout.derivs[0][6]=dd22-dd21;
+  cvout.derivs[0][7]=-dd22;
+  cvout.virial.set( 0,
+                    -(extProduct(d10,dd10)+extProduct(d11,dd11)+extProduct(d12,dd12))
+                    - (extProduct(d20,dd20)+extProduct(d21,dd21)+extProduct(d22,dd22)) );
+}
+
+}
+}
diff --git a/src/colvar/Dipole.cpp b/src/colvar/Dipole.cpp
index f7955b7f46..1f367015b1 100644
--- a/src/colvar/Dipole.cpp
+++ b/src/colvar/Dipole.cpp
@@ -19,10 +19,7 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#ifdef __PLUMED_HAS_OPENACC
-#define __PLUMED_USE_OPENACC 1
-#endif //__PLUMED_HAS_OPENACC
-#include "Colvar.h"
+#include "Dipole.h"
 #include "ColvarShortcut.h"
 #include "MultiColvarTemplate.h"
 #include "core/ActionRegister.h"
@@ -91,165 +88,11 @@ PRINT FILE=output STRIDE=5 ARG=d
 */
 //+ENDPLUMEDOC
 
-class Dipole : public Colvar {
-  std::vector<AtomNumber> ga_lista;
-  bool components;
-  bool nopbc;
-  std::vector<double> value;
-  std::vector<double> derivs;
-  Value* valuex=nullptr;
-  Value* valuey=nullptr;
-  Value* valuez=nullptr;
-public:
-  explicit Dipole(const ActionOptions&);
-  static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
-  static unsigned getModeAndSetupValues( ActionWithValue* av );
-  void calculate() override;
-  static void registerKeywords(Keywords& keys);
-  static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout );
-};
-
 typedef ColvarShortcut<Dipole> DipoleShortcut;
 PLUMED_REGISTER_ACTION(DipoleShortcut,"DIPOLE")
 PLUMED_REGISTER_ACTION(Dipole,"DIPOLE_SCALAR")
 typedef MultiColvarTemplate<Dipole> DipoleMulti;
 PLUMED_REGISTER_ACTION(DipoleMulti,"DIPOLE_VECTOR")
 
-void Dipole::registerKeywords(Keywords& keys) {
-  Colvar::registerKeywords(keys);
-  keys.setDisplayName("DIPOLE");
-  keys.add("atoms","GROUP","the group of atoms we are calculating the dipole moment for");
-  keys.addFlag("COMPONENTS",false,"calculate the x, y and z components of the dipole separately and store them as label.x, label.y and label.z");
-  keys.addOutputComponent("x","COMPONENTS","scalar/vector","the x-component of the dipole");
-  keys.addOutputComponent("y","COMPONENTS","scalar/vector","the y-component of the dipole");
-  keys.addOutputComponent("z","COMPONENTS","scalar/vector","the z-component of the dipole");
-  keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
-  keys.setValueDescription("scalar/vector","the DIPOLE for these atoms");
-  keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
-}
-
-Dipole::Dipole(const ActionOptions&ao):
-  PLUMED_COLVAR_INIT(ao),
-  components(false),
-  value(1) {
-  parseAtomList(-1,ga_lista,this);
-  components=(getModeAndSetupValues(this)==1);
-  if( components ) {
-    value.resize(3);
-    valuex=getPntrToComponent("x");
-    valuey=getPntrToComponent("y");
-    valuez=getPntrToComponent("z");
-  } else {
-    derivs.resize(1,ga_lista.size());
-  }
-  parseFlag("NOPBC",nopbc);
-  checkRead();
-
-  if(nopbc) {
-    log.printf("  without periodic boundary conditions\n");
-  } else {
-    log.printf("  using periodic boundary conditions\n");
-  }
-
-  requestAtoms(ga_lista);
-}
-
-void Dipole::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
-  aa->parseAtomList("GROUP",num,t);
-  if( t.size()>0 ) {
-    aa->log.printf("  of %u atoms\n",static_cast<unsigned>(t.size()));
-    for(unsigned int i=0; i<t.size(); ++i) {
-      aa->log.printf("  %d", t[i].serial());
-    }
-    aa->log.printf("  \n");
-  }
-}
-
-unsigned Dipole::getModeAndSetupValues( ActionWithValue* av ) {
-  bool c;
-  av->parseFlag("COMPONENTS",c);
-  if( c ) {
-    av->addComponentWithDerivatives("x");
-    av->componentIsNotPeriodic("x");
-    av->addComponentWithDerivatives("y");
-    av->componentIsNotPeriodic("y");
-    av->addComponentWithDerivatives("z");
-    av->componentIsNotPeriodic("z");
-    return 1;
-  }
-  av->addValueWithDerivatives();
-  av->setNotPeriodic();
-  return 0;
-}
-
-// calculator
-void Dipole::calculate() {
-  if( !chargesWereSet ) {
-    error("charges were not set by MD code");
-  }
-
-  if(!nopbc) {
-    makeWhole();
-  }
-  unsigned N=getNumberOfAtoms();
-
-  if(!components) {
-    ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) );
-    calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout );
-    for(unsigned i=0; i<N; i++) {
-      setAtomsDerivatives(i,cvout.getAtomDerivatives(0,i));
-    }
-    setBoxDerivatives(cvout.virial[0]);
-    setValue(value[0]);
-  } else {
-    ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) );
-    calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), cvout );
-    for(unsigned i=0; i<N; i++) {
-      setAtomsDerivatives(valuex,i,cvout.getAtomDerivatives(0,i));
-      setAtomsDerivatives(valuey,i,cvout.getAtomDerivatives(1,i));
-      setAtomsDerivatives(valuez,i,cvout.getAtomDerivatives(2,i));
-    }
-    setBoxDerivatives(valuex,cvout.virial[0]);
-    setBoxDerivatives(valuey,cvout.virial[1]);
-    setBoxDerivatives(valuez,cvout.virial[2]);
-    valuex->set(value[0]);
-    valuey->set(value[1]);
-    valuez->set(value[2]);
-  }
-}
-
-void Dipole::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
-  unsigned N=cvin.pos.size();
-  double ctot=0.;
-  for(unsigned i=0; i<N; ++i) {
-    ctot += cvin.charges[i];
-  }
-  ctot/=(double)N;
-
-  Vector dipje;
-  for(unsigned i=0; i<N; ++i) {
-    dipje += (cvin.charges[i]-ctot)*Vector(cvin.pos[i][0],cvin.pos[i][1],cvin.pos[i][2]);
-  }
-
-  if( cvin.mode==1 ) {
-    for(unsigned i=0; i<N; i++) {
-      cvout.derivs[0][i]=(cvin.charges[i]-ctot)*Vector(1.0,0.0,0.0);
-      cvout.derivs[1][i]=(cvin.charges[i]-ctot)*Vector(0.0,1.0,0.0);
-      cvout.derivs[2][i]=(cvin.charges[i]-ctot)*Vector(0.0,0.0,1.0);
-    }
-    for(unsigned i=0; i<3; ++i ) {
-      cvout.values[i] = dipje[i];
-    }
-  } else {
-    cvout.values[0] = dipje.modulo();
-    double idip = 1./cvout.values[0];
-    for(unsigned i=0; i<N; i++) {
-      double dfunc=(cvin.charges[i]-ctot)*idip;
-      cvout.derivs[0][i] = dfunc*dipje;
-    }
-  }
-  ColvarInput::setBoxDerivativesNoPbc( cvin, cvout );
-}
-
 }
 }
diff --git a/src/colvar/Dipole.h b/src/colvar/Dipole.h
new file mode 100644
index 0000000000..96bd2193d9
--- /dev/null
+++ b/src/colvar/Dipole.h
@@ -0,0 +1,184 @@
+
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2012-2023 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "Colvar.h"
+#include "ColvarInput.h"
+
+namespace PLMD {
+namespace colvar {
+
+class Dipole : public Colvar {
+  std::vector<AtomNumber> ga_lista;
+  bool components;
+  bool nopbc;
+  std::vector<double> value;
+  std::vector<double> derivs;
+  Value* valuex=nullptr;
+  Value* valuey=nullptr;
+  Value* valuez=nullptr;
+public:
+  explicit Dipole(const ActionOptions&);
+  static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
+  static unsigned getModeAndSetupValues( ActionWithValue* av );
+  void calculate() override;
+  static void registerKeywords(Keywords& keys);
+  static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout );
+};
+
+void Dipole::registerKeywords(Keywords& keys) {
+  Colvar::registerKeywords(keys);
+  keys.setDisplayName("DIPOLE");
+  keys.add("atoms","GROUP","the group of atoms we are calculating the dipole moment for");
+  keys.addFlag("COMPONENTS",false,"calculate the x, y and z components of the dipole separately and store them as label.x, label.y and label.z");
+  keys.addOutputComponent("x","COMPONENTS","scalar/vector","the x-component of the dipole");
+  keys.addOutputComponent("y","COMPONENTS","scalar/vector","the y-component of the dipole");
+  keys.addOutputComponent("z","COMPONENTS","scalar/vector","the z-component of the dipole");
+  keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
+  keys.setValueDescription("scalar/vector","the DIPOLE for these atoms");
+  keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
+}
+
+Dipole::Dipole(const ActionOptions&ao):
+  PLUMED_COLVAR_INIT(ao),
+  components(false),
+  value(1) {
+  parseAtomList(-1,ga_lista,this);
+  components=(getModeAndSetupValues(this)==1);
+  if( components ) {
+    value.resize(3);
+    valuex=getPntrToComponent("x");
+    valuey=getPntrToComponent("y");
+    valuez=getPntrToComponent("z");
+  } else {
+    derivs.resize(1,ga_lista.size());
+  }
+  parseFlag("NOPBC",nopbc);
+  checkRead();
+
+  if(nopbc) {
+    log.printf("  without periodic boundary conditions\n");
+  } else {
+    log.printf("  using periodic boundary conditions\n");
+  }
+
+  requestAtoms(ga_lista);
+}
+
+void Dipole::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
+  aa->parseAtomList("GROUP",num,t);
+  if( t.size()>0 ) {
+    aa->log.printf("  of %u atoms\n",static_cast<unsigned>(t.size()));
+    for(unsigned int i=0; i<t.size(); ++i) {
+      aa->log.printf("  %d", t[i].serial());
+    }
+    aa->log.printf("  \n");
+  }
+}
+
+unsigned Dipole::getModeAndSetupValues( ActionWithValue* av ) {
+  bool c;
+  av->parseFlag("COMPONENTS",c);
+  if( c ) {
+    av->addComponentWithDerivatives("x");
+    av->componentIsNotPeriodic("x");
+    av->addComponentWithDerivatives("y");
+    av->componentIsNotPeriodic("y");
+    av->addComponentWithDerivatives("z");
+    av->componentIsNotPeriodic("z");
+    return 1;
+  }
+  av->addValueWithDerivatives();
+  av->setNotPeriodic();
+  return 0;
+}
+
+// calculator
+void Dipole::calculate() {
+  if( !chargesWereSet ) {
+    error("charges were not set by MD code");
+  }
+
+  if(!nopbc) {
+    makeWhole();
+  }
+  unsigned N=getNumberOfAtoms();
+
+  if(!components) {
+    ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) );
+    calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout );
+    for(unsigned i=0; i<N; i++) {
+      setAtomsDerivatives(i,cvout.getAtomDerivatives(0,i));
+    }
+    setBoxDerivatives(cvout.virial[0]);
+    setValue(value[0]);
+  } else {
+    ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) );
+    calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), cvout );
+    for(unsigned i=0; i<N; i++) {
+      setAtomsDerivatives(valuex,i,cvout.getAtomDerivatives(0,i));
+      setAtomsDerivatives(valuey,i,cvout.getAtomDerivatives(1,i));
+      setAtomsDerivatives(valuez,i,cvout.getAtomDerivatives(2,i));
+    }
+    setBoxDerivatives(valuex,cvout.virial[0]);
+    setBoxDerivatives(valuey,cvout.virial[1]);
+    setBoxDerivatives(valuez,cvout.virial[2]);
+    valuex->set(value[0]);
+    valuey->set(value[1]);
+    valuez->set(value[2]);
+  }
+}
+
+void Dipole::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
+  unsigned N=cvin.pos.size();
+  double ctot=0.;
+  for(unsigned i=0; i<N; ++i) {
+    ctot += cvin.charges[i];
+  }
+  ctot/=(double)N;
+
+  Vector dipje;
+  for(unsigned i=0; i<N; ++i) {
+    dipje += (cvin.charges[i]-ctot)*Vector(cvin.pos[i][0],cvin.pos[i][1],cvin.pos[i][2]);
+  }
+
+  if( cvin.mode==1 ) {
+    for(unsigned i=0; i<N; i++) {
+      cvout.derivs[0][i]=(cvin.charges[i]-ctot)*Vector(1.0,0.0,0.0);
+      cvout.derivs[1][i]=(cvin.charges[i]-ctot)*Vector(0.0,1.0,0.0);
+      cvout.derivs[2][i]=(cvin.charges[i]-ctot)*Vector(0.0,0.0,1.0);
+    }
+    for(unsigned i=0; i<3; ++i ) {
+      cvout.values[i] = dipje[i];
+    }
+  } else {
+    cvout.values[0] = dipje.modulo();
+    double idip = 1./cvout.values[0];
+    for(unsigned i=0; i<N; i++) {
+      double dfunc=(cvin.charges[i]-ctot)*idip;
+      cvout.derivs[0][i] = dfunc*dipje;
+    }
+  }
+  ColvarInput::setBoxDerivativesNoPbc( cvin, cvout );
+}
+
+}
+}
diff --git a/src/colvar/Plane.cpp b/src/colvar/Plane.cpp
index a88f2ea845..1a4b48ecc5 100644
--- a/src/colvar/Plane.cpp
+++ b/src/colvar/Plane.cpp
@@ -19,18 +19,12 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#ifdef __PLUMED_HAS_OPENACC
-#define __PLUMED_USE_OPENACC 1
-#endif //__PLUMED_HAS_OPENACC
-#include "Colvar.h"
+#include "Plane.h"
 #include "ColvarShortcut.h"
 #include "core/ActionRegister.h"
 #include "MultiColvarTemplate.h"
 #include "tools/Pbc.h"
 
-#include <string>
-#include <cmath>
-
 //+PLUMEDOC COLVAR PLANE
 /*
 Calculate the plane perpendicular to two vectors in order to represent the orientation of a planar molecule.
@@ -78,126 +72,13 @@ nucleation of molecular crystals such as [SMAC](SMAC.md).
 namespace PLMD {
 namespace colvar {
 
-class Plane : public Colvar {
-private:
-  bool pbc;
-  std::vector<double> value;
-  std::vector<double> derivs;
-public:
-  static void registerKeywords( Keywords& keys );
-  explicit Plane(const ActionOptions&);
-  static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
-  static unsigned getModeAndSetupValues( ActionWithValue* av );
-// active methods:
-  void calculate() override;
-  static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout );
-};
-
 typedef ColvarShortcut<Plane> PlaneShortcut;
 PLUMED_REGISTER_ACTION(PlaneShortcut,"PLANE")
 PLUMED_REGISTER_ACTION(Plane,"PLANE_SCALAR")
 typedef MultiColvarTemplate<Plane> PlaneMulti;
 PLUMED_REGISTER_ACTION(PlaneMulti,"PLANE_VECTOR")
 
-void Plane::registerKeywords( Keywords& keys ) {
-  Colvar::registerKeywords( keys );
-  keys.setDisplayName("PLANE");
-  keys.add("atoms","ATOMS","the three or four atoms whose plane we are computing");
-  keys.addOutputComponent("x","default","scalar/vector","the x-component of the vector that is normal to the plane containing the atoms");
-  keys.addOutputComponent("y","default","scalar/vector","the y-component of the vector that is normal to the plane containing the atoms");
-  keys.addOutputComponent("z","default","scalar/vector","the z-component of the vector that is normal to the plane containing the atoms");
-  keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
-  keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
-}
-
-void Plane::parseAtomList( const int& num, std::vector<AtomNumber>& atoms, ActionAtomistic* aa ) {
-  aa->parseAtomList("ATOMS",num,atoms);
-  if(atoms.size()==3) {
-    aa->log.printf("  containing atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial());
-    atoms.resize(4);
-    atoms[3]=atoms[2];
-    atoms[2]=atoms[1];
-  } else if(atoms.size()==4) {
-    aa->log.printf("  containing lines %d-%d and %d-%d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial());
-  } else if( num<0 || atoms.size()>0 ) {
-    aa->error("Number of specified atoms should be either 3 or 4");
-  }
-}
-
-unsigned Plane::getModeAndSetupValues( ActionWithValue* av ) {
-  av->addComponentWithDerivatives("x");
-  av->componentIsNotPeriodic("x");
-  av->addComponentWithDerivatives("y");
-  av->componentIsNotPeriodic("y");
-  av->addComponentWithDerivatives("z");
-  av->componentIsNotPeriodic("z");
-  return 0;
-}
 
-Plane::Plane(const ActionOptions&ao):
-  PLUMED_COLVAR_INIT(ao),
-  pbc(true),
-  value(3) {
-  std::vector<AtomNumber> atoms;
-  parseAtomList(-1,atoms,this);
-  bool nopbc=!pbc;
-  parseFlag("NOPBC",nopbc);
-  pbc=!nopbc;
-
-  if(pbc) {
-    log.printf("  using periodic boundary conditions\n");
-  } else {
-    log.printf("  without periodic boundary conditions\n");
-  }
-
-  unsigned mode = getModeAndSetupValues( this );
-  requestAtoms(atoms);
-  checkRead();
 }
-
-void Plane::calculate() {
-
-  if(pbc) {
-    makeWhole();
-  }
-  ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this);
-  calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout );
-  setValue( value[0] );
-  for(unsigned i=0; i<getPositions().size(); ++i) {
-    setAtomsDerivatives( i, cvout.getAtomDerivatives(0,i) );
-  }
-  setBoxDerivatives( cvout.virial[0] );
 }
 
-void Plane::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
-  Vector d1=delta( cvin.pos[1], cvin.pos[0] );
-  Vector d2=delta( cvin.pos[2], cvin.pos[3] );
-  Vector cp = crossProduct( d1, d2 );
-
-  cvout.derivs[0][0] = crossProduct( Vector(-1.0,0,0), d2 );
-  cvout.derivs[0][1] = crossProduct( Vector(+1.0,0,0), d2 );
-  cvout.derivs[0][2] = crossProduct( Vector(-1.0,0,0), d1 );
-  cvout.derivs[0][3] = crossProduct( Vector(+1.0,0,0), d1 );
-  cvout.virial.set( 0, Tensor(d1,crossProduct(Vector(+1.0,0,0), d2)) + Tensor( d2, crossProduct(Vector(-1.0,0,0), d1)) );
-  cvout.values[0] = cp[0];
-
-  cvout.derivs[1][0] = crossProduct( Vector(0,-1.0,0), d2 );
-  cvout.derivs[1][1] = crossProduct( Vector(0,+1.0,0), d2 );
-  cvout.derivs[1][2] = crossProduct( Vector(0,-1.0,0), d1 );
-  cvout.derivs[1][3] = crossProduct( Vector(0,+1.0,0), d1 );
-  cvout.virial.set(1, Tensor(d1,crossProduct(Vector(0,+1.0,0), d2)) + Tensor( d2, crossProduct(Vector(0,-1.0,0), d1)) );
-  cvout.values[1] = cp[1];
-
-  cvout.derivs[2][0] = crossProduct( Vector(0,0,-1.0), d2 );
-  cvout.derivs[2][1] = crossProduct( Vector(0,0,+1.0), d2 );
-  cvout.derivs[2][2] = crossProduct( Vector(0,0,-1.0), d1 );
-  cvout.derivs[2][3] = crossProduct( Vector(0,0,+1.0), d1 );
-  cvout.virial.set(2, Tensor(d1,crossProduct(Vector(0,0,+1.0), d2)) + Tensor( d2, crossProduct(Vector(0,0,-1.0), d1)) );
-  cvout.values[2] = cp[2];
-}
-
-}
-}
-
-
-
diff --git a/src/colvar/Plane.h b/src/colvar/Plane.h
new file mode 100644
index 0000000000..431b3c1e14
--- /dev/null
+++ b/src/colvar/Plane.h
@@ -0,0 +1,146 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2011-2023 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "Colvar.h"
+#include "ColvarInput.h"
+#include "tools/Pbc.h"
+
+#include <string>
+#include <cmath>
+
+namespace PLMD {
+namespace colvar {
+
+class Plane : public Colvar {
+private:
+  bool pbc;
+  std::vector<double> value;
+  std::vector<double> derivs;
+public:
+  static void registerKeywords( Keywords& keys );
+  explicit Plane(const ActionOptions&);
+  static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
+  static unsigned getModeAndSetupValues( ActionWithValue* av );
+// active methods:
+  void calculate() override;
+  static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout );
+};
+
+void Plane::registerKeywords( Keywords& keys ) {
+  Colvar::registerKeywords( keys );
+  keys.setDisplayName("PLANE");
+  keys.add("atoms","ATOMS","the three or four atoms whose plane we are computing");
+  keys.addOutputComponent("x","default","scalar/vector","the x-component of the vector that is normal to the plane containing the atoms");
+  keys.addOutputComponent("y","default","scalar/vector","the y-component of the vector that is normal to the plane containing the atoms");
+  keys.addOutputComponent("z","default","scalar/vector","the z-component of the vector that is normal to the plane containing the atoms");
+  keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
+  keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
+}
+
+void Plane::parseAtomList( const int& num, std::vector<AtomNumber>& atoms, ActionAtomistic* aa ) {
+  aa->parseAtomList("ATOMS",num,atoms);
+  if(atoms.size()==3) {
+    aa->log.printf("  containing atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial());
+    atoms.resize(4);
+    atoms[3]=atoms[2];
+    atoms[2]=atoms[1];
+  } else if(atoms.size()==4) {
+    aa->log.printf("  containing lines %d-%d and %d-%d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial());
+  } else if( num<0 || atoms.size()>0 ) {
+    aa->error("Number of specified atoms should be either 3 or 4");
+  }
+}
+
+unsigned Plane::getModeAndSetupValues( ActionWithValue* av ) {
+  av->addComponentWithDerivatives("x");
+  av->componentIsNotPeriodic("x");
+  av->addComponentWithDerivatives("y");
+  av->componentIsNotPeriodic("y");
+  av->addComponentWithDerivatives("z");
+  av->componentIsNotPeriodic("z");
+  return 0;
+}
+
+Plane::Plane(const ActionOptions&ao):
+  PLUMED_COLVAR_INIT(ao),
+  pbc(true),
+  value(3) {
+  std::vector<AtomNumber> atoms;
+  parseAtomList(-1,atoms,this);
+  bool nopbc=!pbc;
+  parseFlag("NOPBC",nopbc);
+  pbc=!nopbc;
+
+  if(pbc) {
+    log.printf("  using periodic boundary conditions\n");
+  } else {
+    log.printf("  without periodic boundary conditions\n");
+  }
+
+  unsigned mode = getModeAndSetupValues( this );
+  requestAtoms(atoms);
+  checkRead();
+}
+
+void Plane::calculate() {
+
+  if(pbc) {
+    makeWhole();
+  }
+  ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this);
+  calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout );
+  setValue( value[0] );
+  for(unsigned i=0; i<getPositions().size(); ++i) {
+    setAtomsDerivatives( i, cvout.getAtomDerivatives(0,i) );
+  }
+  setBoxDerivatives( cvout.virial[0] );
+}
+
+void Plane::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
+  Vector d1=delta( cvin.pos[1], cvin.pos[0] );
+  Vector d2=delta( cvin.pos[2], cvin.pos[3] );
+  Vector cp = crossProduct( d1, d2 );
+
+  cvout.derivs[0][0] = crossProduct( Vector(-1.0,0,0), d2 );
+  cvout.derivs[0][1] = crossProduct( Vector(+1.0,0,0), d2 );
+  cvout.derivs[0][2] = crossProduct( Vector(-1.0,0,0), d1 );
+  cvout.derivs[0][3] = crossProduct( Vector(+1.0,0,0), d1 );
+  cvout.virial.set( 0, Tensor(d1,crossProduct(Vector(+1.0,0,0), d2)) + Tensor( d2, crossProduct(Vector(-1.0,0,0), d1)) );
+  cvout.values[0] = cp[0];
+
+  cvout.derivs[1][0] = crossProduct( Vector(0,-1.0,0), d2 );
+  cvout.derivs[1][1] = crossProduct( Vector(0,+1.0,0), d2 );
+  cvout.derivs[1][2] = crossProduct( Vector(0,-1.0,0), d1 );
+  cvout.derivs[1][3] = crossProduct( Vector(0,+1.0,0), d1 );
+  cvout.virial.set(1, Tensor(d1,crossProduct(Vector(0,+1.0,0), d2)) + Tensor( d2, crossProduct(Vector(0,-1.0,0), d1)) );
+  cvout.values[1] = cp[1];
+
+  cvout.derivs[2][0] = crossProduct( Vector(0,0,-1.0), d2 );
+  cvout.derivs[2][1] = crossProduct( Vector(0,0,+1.0), d2 );
+  cvout.derivs[2][2] = crossProduct( Vector(0,0,-1.0), d1 );
+  cvout.derivs[2][3] = crossProduct( Vector(0,0,+1.0), d1 );
+  cvout.virial.set(2, Tensor(d1,crossProduct(Vector(0,0,+1.0), d2)) + Tensor( d2, crossProduct(Vector(0,0,-1.0), d1)) );
+  cvout.values[2] = cp[2];
+}
+
+}
+}
+
diff --git a/src/colvar/Position.cpp b/src/colvar/Position.cpp
index f9e0cd0352..99875c6424 100644
--- a/src/colvar/Position.cpp
+++ b/src/colvar/Position.cpp
@@ -19,14 +19,10 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#ifdef __PLUMED_HAS_OPENACC
-#define __PLUMED_USE_OPENACC 1
-#endif //__PLUMED_HAS_OPENACC
-#include "Colvar.h"
+#include "Position.h"
 #include "ColvarShortcut.h"
 #include "MultiColvarTemplate.h"
 #include "core/ActionRegister.h"
-#include "tools/Pbc.h"
 
 namespace PLMD {
 namespace colvar {
@@ -106,164 +102,12 @@ PRINT ARG=p.x,p.y,p.z
 */
 //+ENDPLUMEDOC
 
-class Position : public Colvar {
-  bool scaled_components;
-  bool pbc;
-  std::vector<double> value;
-  std::vector<double> derivs;
-public:
-  static void registerKeywords( Keywords& keys );
-  explicit Position(const ActionOptions&);
-  static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
-  static unsigned getModeAndSetupValues( ActionWithValue* av );
-// active methods:
-  void calculate() override;
-  static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout );
-};
-
 typedef ColvarShortcut<Position> PositionShortcut;
 PLUMED_REGISTER_ACTION(PositionShortcut,"POSITION")
 PLUMED_REGISTER_ACTION(Position,"POSITION_SCALAR")
 typedef MultiColvarTemplate<Position> PositionMulti;
 PLUMED_REGISTER_ACTION(PositionMulti,"POSITION_VECTOR")
 
-void Position::registerKeywords( Keywords& keys ) {
-  Colvar::registerKeywords( keys );
-  keys.setDisplayName("POSITION");
-  keys.add("atoms","ATOM","the atom number");
-  keys.add("atoms","ATOMS","the atom numbers that you would like to use the positions of");
-  keys.addFlag("WHOLEMOLECULES",false,"if this is a vector of positions do you want to make the positions into a whole before");
-  keys.addFlag("SCALED_COMPONENTS",false,"calculate the a, b and c scaled components of the position separately and store them as label.a, label.b and label.c");
-  keys.addOutputComponent("x","default","scalar/vector","the x-component of the atom position");
-  keys.addOutputComponent("y","default","scalar/vector","the y-component of the atom position");
-  keys.addOutputComponent("z","default","scalar/vector","the z-component of the atom position");
-  keys.addOutputComponent("a","SCALED_COMPONENTS","scalar/vector","the normalized projection on the first lattice vector of the atom position");
-  keys.addOutputComponent("b","SCALED_COMPONENTS","scalar/vector","the normalized projection on the second lattice vector of the atom position");
-  keys.addOutputComponent("c","SCALED_COMPONENTS","scalar/vector","the normalized projection on the third lattice vector of the atom position");
-  keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
-  keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
-}
-
-Position::Position(const ActionOptions&ao):
-  PLUMED_COLVAR_INIT(ao),
-  scaled_components(false),
-  pbc(true),
-  value(3) {
-  std::vector<AtomNumber> atoms;
-  parseAtomList(-1,atoms,this);
-  unsigned mode=getModeAndSetupValues(this);
-  if( mode==1 ) {
-    scaled_components=true;
-  }
-
-  bool nopbc=!pbc;
-  parseFlag("NOPBC",nopbc);
-  pbc=!nopbc;
-  checkRead();
-
-  if(pbc) {
-    log.printf("  using periodic boundary conditions\n");
-  } else {
-    log.printf("  without periodic boundary conditions\n");
-  }
-
-  requestAtoms(atoms);
-}
-
-void Position::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
-  aa->parseAtomList("ATOM",num,t);
-  if( t.size()==1 ) {
-    aa->log.printf("  for atom %d\n",t[0].serial());
-  } else if( num<0 || t.size()!=0 ) {
-    aa->error("Number of specified atoms should be 1");
-  }
-}
-
-unsigned Position::getModeAndSetupValues( ActionWithValue* av ) {
-  bool sc;
-  av->parseFlag("SCALED_COMPONENTS",sc);
-  if(sc) {
-    av->addComponentWithDerivatives("a");
-    av->componentIsPeriodic("a","-0.5","+0.5");
-    av->addComponentWithDerivatives("b");
-    av->componentIsPeriodic("b","-0.5","+0.5");
-    av->addComponentWithDerivatives("c");
-    av->componentIsPeriodic("c","-0.5","+0.5");
-    return 1;
-  }
-  av->addComponentWithDerivatives("x");
-  av->componentIsNotPeriodic("x");
-  av->addComponentWithDerivatives("y");
-  av->componentIsNotPeriodic("y");
-  av->addComponentWithDerivatives("z");
-  av->componentIsNotPeriodic("z");
-  av->log<<"  WARNING: components will not have the proper periodicity - see manual\n";
-  return 0;
-}
-
-// calculator
-void Position::calculate() {
-
-  std::vector<Vector> distance(1);
-  if(pbc) {
-    distance[0]=pbcDistance(Vector(0.0,0.0,0.0),getPosition(0));
-  } else {
-    distance[0]=delta(Vector(0.0,0.0,0.0),getPosition(0));
-  }
-
-  ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this);
-  if(scaled_components) {
-    calculateCV( ColvarInput::createColvarInput( 1, distance, this ), cvout );
-    Value* valuea=getPntrToComponent("a");
-    Value* valueb=getPntrToComponent("b");
-    Value* valuec=getPntrToComponent("c");
-    setAtomsDerivatives (valuea,0,cvout.getAtomDerivatives(0,0));
-    valuea->set(value[0]);
-    setAtomsDerivatives (valueb,0,cvout.getAtomDerivatives(1,0));
-    valueb->set(value[1]);
-    setAtomsDerivatives (valuec,0,cvout.getAtomDerivatives(2,0));
-    valuec->set(value[2]);
-  } else {
-    calculateCV( ColvarInput::createColvarInput( 0, distance, this ), cvout );
-    Value* valuex=getPntrToComponent("x");
-    Value* valuey=getPntrToComponent("y");
-    Value* valuez=getPntrToComponent("z");
-
-    setAtomsDerivatives (valuex,0,cvout.getAtomDerivatives(0,0));
-    setBoxDerivatives   (valuex,cvout.virial[0]);
-    valuex->set(value[0]);
-
-    setAtomsDerivatives (valuey,0,cvout.getAtomDerivatives(1,0));
-    setBoxDerivatives   (valuey,cvout.virial[1]);
-    valuey->set(value[1]);
-
-    setAtomsDerivatives (valuez,0,cvout.getAtomDerivatives(2,0));
-    setBoxDerivatives   (valuez,cvout.virial[2]);
-    valuez->set(value[2]);
-  }
-}
-
-void Position::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
-  if( cvin.mode==1 ) {
-    Vector d=cvin.pbc.realToScaled(Vector(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]));
-    cvout.values[0]=Tools::pbc(d[0]);
-    cvout.values[1]=Tools::pbc(d[1]);
-    cvout.values[2]=Tools::pbc(d[2]);
-    cvout.derivs[0][0]=matmul(cvin.pbc.getInvBox(),Vector(+1,0,0));
-    cvout.derivs[1][0]=matmul(cvin.pbc.getInvBox(),Vector(0,+1,0));
-    cvout.derivs[2][0]=matmul(cvin.pbc.getInvBox(),Vector(0,0,+1));
-  } else {
-    for(unsigned i=0; i<3; ++i) {
-      cvout.values[i]=cvin.pos[0][i];
-    }
-    cvout.derivs[0][0]=Vector(+1,0,0);
-    cvout.derivs[1][0]=Vector(0,+1,0);
-    cvout.derivs[2][0]=Vector(0,0,+1);
-    cvout.virial.set(0, Tensor(Vector(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]),Vector(-1,0,0)) );
-    cvout.virial.set(1, Tensor(Vector(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]),Vector(0,-1,0)) );
-    cvout.virial.set(2, Tensor(Vector(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]),Vector(0,0,-1)) );
-  }
-}
 
 }
 }
diff --git a/src/colvar/Position.h b/src/colvar/Position.h
new file mode 100644
index 0000000000..8325886d32
--- /dev/null
+++ b/src/colvar/Position.h
@@ -0,0 +1,186 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2014-2023 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "Colvar.h"
+#include "tools/Pbc.h"
+#include "ColvarInput.h"
+
+namespace PLMD {
+namespace colvar {
+
+class Position : public Colvar {
+  bool scaled_components;
+  bool pbc;
+  std::vector<double> value;
+  std::vector<double> derivs;
+public:
+  static void registerKeywords( Keywords& keys );
+  explicit Position(const ActionOptions&);
+  static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
+  static unsigned getModeAndSetupValues( ActionWithValue* av );
+// active methods:
+  void calculate() override;
+  static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout );
+};
+
+void Position::registerKeywords( Keywords& keys ) {
+  Colvar::registerKeywords( keys );
+  keys.setDisplayName("POSITION");
+  keys.add("atoms","ATOM","the atom number");
+  keys.add("atoms","ATOMS","the atom numbers that you would like to use the positions of");
+  keys.addFlag("WHOLEMOLECULES",false,"if this is a vector of positions do you want to make the positions into a whole before");
+  keys.addFlag("SCALED_COMPONENTS",false,"calculate the a, b and c scaled components of the position separately and store them as label.a, label.b and label.c");
+  keys.addOutputComponent("x","default","scalar/vector","the x-component of the atom position");
+  keys.addOutputComponent("y","default","scalar/vector","the y-component of the atom position");
+  keys.addOutputComponent("z","default","scalar/vector","the z-component of the atom position");
+  keys.addOutputComponent("a","SCALED_COMPONENTS","scalar/vector","the normalized projection on the first lattice vector of the atom position");
+  keys.addOutputComponent("b","SCALED_COMPONENTS","scalar/vector","the normalized projection on the second lattice vector of the atom position");
+  keys.addOutputComponent("c","SCALED_COMPONENTS","scalar/vector","the normalized projection on the third lattice vector of the atom position");
+  keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
+  keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
+}
+
+Position::Position(const ActionOptions&ao):
+  PLUMED_COLVAR_INIT(ao),
+  scaled_components(false),
+  pbc(true),
+  value(3) {
+  std::vector<AtomNumber> atoms;
+  parseAtomList(-1,atoms,this);
+  unsigned mode=getModeAndSetupValues(this);
+  if( mode==1 ) {
+    scaled_components=true;
+  }
+
+  bool nopbc=!pbc;
+  parseFlag("NOPBC",nopbc);
+  pbc=!nopbc;
+  checkRead();
+
+  if(pbc) {
+    log.printf("  using periodic boundary conditions\n");
+  } else {
+    log.printf("  without periodic boundary conditions\n");
+  }
+
+  requestAtoms(atoms);
+}
+
+void Position::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
+  aa->parseAtomList("ATOM",num,t);
+  if( t.size()==1 ) {
+    aa->log.printf("  for atom %d\n",t[0].serial());
+  } else if( num<0 || t.size()!=0 ) {
+    aa->error("Number of specified atoms should be 1");
+  }
+}
+
+unsigned Position::getModeAndSetupValues( ActionWithValue* av ) {
+  bool sc;
+  av->parseFlag("SCALED_COMPONENTS",sc);
+  if(sc) {
+    av->addComponentWithDerivatives("a");
+    av->componentIsPeriodic("a","-0.5","+0.5");
+    av->addComponentWithDerivatives("b");
+    av->componentIsPeriodic("b","-0.5","+0.5");
+    av->addComponentWithDerivatives("c");
+    av->componentIsPeriodic("c","-0.5","+0.5");
+    return 1;
+  }
+  av->addComponentWithDerivatives("x");
+  av->componentIsNotPeriodic("x");
+  av->addComponentWithDerivatives("y");
+  av->componentIsNotPeriodic("y");
+  av->addComponentWithDerivatives("z");
+  av->componentIsNotPeriodic("z");
+  av->log<<"  WARNING: components will not have the proper periodicity - see manual\n";
+  return 0;
+}
+
+// calculator
+void Position::calculate() {
+
+  std::vector<Vector> distance(1);
+  if(pbc) {
+    distance[0]=pbcDistance(Vector(0.0,0.0,0.0),getPosition(0));
+  } else {
+    distance[0]=delta(Vector(0.0,0.0,0.0),getPosition(0));
+  }
+
+  ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this);
+  if(scaled_components) {
+    calculateCV( ColvarInput::createColvarInput( 1, distance, this ), cvout );
+    Value* valuea=getPntrToComponent("a");
+    Value* valueb=getPntrToComponent("b");
+    Value* valuec=getPntrToComponent("c");
+    setAtomsDerivatives (valuea,0,cvout.getAtomDerivatives(0,0));
+    valuea->set(value[0]);
+    setAtomsDerivatives (valueb,0,cvout.getAtomDerivatives(1,0));
+    valueb->set(value[1]);
+    setAtomsDerivatives (valuec,0,cvout.getAtomDerivatives(2,0));
+    valuec->set(value[2]);
+  } else {
+    calculateCV( ColvarInput::createColvarInput( 0, distance, this ), cvout );
+    Value* valuex=getPntrToComponent("x");
+    Value* valuey=getPntrToComponent("y");
+    Value* valuez=getPntrToComponent("z");
+
+    setAtomsDerivatives (valuex,0,cvout.getAtomDerivatives(0,0));
+    setBoxDerivatives   (valuex,cvout.virial[0]);
+    valuex->set(value[0]);
+
+    setAtomsDerivatives (valuey,0,cvout.getAtomDerivatives(1,0));
+    setBoxDerivatives   (valuey,cvout.virial[1]);
+    valuey->set(value[1]);
+
+    setAtomsDerivatives (valuez,0,cvout.getAtomDerivatives(2,0));
+    setBoxDerivatives   (valuez,cvout.virial[2]);
+    valuez->set(value[2]);
+  }
+}
+
+void Position::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
+  if( cvin.mode==1 ) {
+    Vector d=cvin.pbc.realToScaled(Vector(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]));
+    cvout.values[0]=Tools::pbc(d[0]);
+    cvout.values[1]=Tools::pbc(d[1]);
+    cvout.values[2]=Tools::pbc(d[2]);
+    cvout.derivs[0][0]=matmul(cvin.pbc.getInvBox(),Vector(+1,0,0));
+    cvout.derivs[1][0]=matmul(cvin.pbc.getInvBox(),Vector(0,+1,0));
+    cvout.derivs[2][0]=matmul(cvin.pbc.getInvBox(),Vector(0,0,+1));
+  } else {
+    for(unsigned i=0; i<3; ++i) {
+      cvout.values[i]=cvin.pos[0][i];
+    }
+    cvout.derivs[0][0]=Vector(+1,0,0);
+    cvout.derivs[1][0]=Vector(0,+1,0);
+    cvout.derivs[2][0]=Vector(0,0,+1);
+    cvout.virial.set(0, Tensor(Vector(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]),Vector(-1,0,0)) );
+    cvout.virial.set(1, Tensor(Vector(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]),Vector(0,-1,0)) );
+    cvout.virial.set(2, Tensor(Vector(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]),Vector(0,0,-1)) );
+  }
+}
+
+}
+}
+
+
+
diff --git a/src/colvar/Torsion.cpp b/src/colvar/Torsion.cpp
index a55cb3dedb..ca7cbddba5 100644
--- a/src/colvar/Torsion.cpp
+++ b/src/colvar/Torsion.cpp
@@ -19,15 +19,10 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-//this is temporary:
-#ifdef __PLUMED_HAS_OPENACC
-#define __PLUMED_USE_OPENACC 1
-#endif //__PLUMED_HAS_OPENACC
-#include "Colvar.h"
+#include "Torsion.h"
 #include "ColvarShortcut.h"
 #include "MultiColvarTemplate.h"
 #include "core/ActionRegister.h"
-#include "tools/Torsion.h"
 
 namespace PLMD {
 namespace colvar {
@@ -115,185 +110,12 @@ PRINT ARG=t FILE=COLVAR
 */
 //+ENDPLUMEDOC
 
-class Torsion : public Colvar {
-  bool pbc;
-  bool do_cosine;
-
-  std::vector<double> value;
-  std::vector<double> derivs;
-public:
-  explicit Torsion(const ActionOptions&);
-  static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
-  static unsigned getModeAndSetupValues( ActionWithValue* av );
-// active methods:
-  void calculate() override;
-  static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout );
-  static void registerKeywords(Keywords& keys);
-};
-
 typedef ColvarShortcut<Torsion> TorsionShortcut;
 PLUMED_REGISTER_ACTION(TorsionShortcut,"TORSION")
 PLUMED_REGISTER_ACTION(Torsion,"TORSION_SCALAR")
 typedef MultiColvarTemplate<Torsion> TorsionMulti;
 PLUMED_REGISTER_ACTION(TorsionMulti,"TORSION_VECTOR")
 
-void Torsion::registerKeywords(Keywords& keys) {
-  Colvar::registerKeywords( keys );
-  keys.setDisplayName("TORSION");
-  keys.add("atoms-1","ATOMS","the four atoms involved in the torsional angle");
-  keys.add("atoms-2","AXIS","two atoms that define an axis.  You can use this to find the angle in the plane perpendicular to the axis between the vectors specified using the VECTORA and VECTORB keywords.");
-  keys.add("atoms-2","VECTORA","two atoms that define a vector.  You can use this in combination with VECTOR2 and AXIS");
-  keys.add("atoms-2","VECTORB","two atoms that define a vector.  You can use this in combination with VECTOR1 and AXIS");
-  keys.addDeprecatedKeyword("VECTOR1","VECTORA");
-  keys.addDeprecatedKeyword("VECTOR2","VECTORB");
-  keys.addFlag("COSINE",false,"calculate cosine instead of dihedral");
-  keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
-  keys.setValueDescription("scalar/vector","the TORSION involving these atoms");
-  keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
-}
-
-Torsion::Torsion(const ActionOptions&ao):
-  PLUMED_COLVAR_INIT(ao),
-  pbc(true),
-  do_cosine(false),
-  value(1) {
-  std::vector<AtomNumber> atoms;
-  std::vector<AtomNumber> v1;
-  ActionAtomistic::parseAtomList("VECTOR1",v1);
-  if( v1.size()>0 ) {
-    std::vector<AtomNumber> v2;
-    ActionAtomistic::parseAtomList("VECTOR2",v2);
-    std::vector<AtomNumber> axis;
-    ActionAtomistic::parseAtomList("AXIS",axis);
-    if( !(v1.size()==2 && v2.size()==2 && axis.size()==2)) {
-      error("VECTOR1, VECTOR2 and AXIS should specify 2 atoms each");
-    }
-    atoms.resize(6);
-    atoms[0]=v1[1];
-    atoms[1]=v1[0];
-    atoms[2]=axis[0];
-    atoms[3]=axis[1];
-    atoms[4]=v2[0];
-    atoms[5]=v2[1];
-    log.printf("  between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n",
-               v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial());
-  } else {
-    parseAtomList(-1,atoms,this);
-  }
-  unsigned mode=getModeAndSetupValues(this);
-  if( mode==1 ) {
-    do_cosine=true;
-  }
-
-  bool nopbc=!pbc;
-  parseFlag("NOPBC",nopbc);
-  pbc=!nopbc;
-  checkRead();
-
-  if(pbc) {
-    log.printf("  using periodic boundary conditions\n");
-  } else {
-    log.printf("  without periodic boundary conditions\n");
-  }
-  requestAtoms(atoms);
-}
-
-void Torsion::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
-  std::vector<AtomNumber> v1,v2,axis;
-  aa->parseAtomList("ATOMS",num,t);
-  aa->parseAtomList("VECTORA",num,v1);
-  aa->parseAtomList("VECTORB",num,v2);
-  aa->parseAtomList("AXIS",num,axis);
-
-  if(t.size()==4) {
-    if(!(v1.empty() && v2.empty() && axis.empty())) {
-      aa->error("ATOMS keyword is not compatible with VECTORA, VECTORB and AXIS keywords");
-    }
-    aa->log.printf("  between atoms %d %d %d %d\n",t[0].serial(),t[1].serial(),t[2].serial(),t[3].serial());
-    t.resize(6);
-    t[5]=t[3];
-    t[4]=t[2];
-    t[3]=t[2];
-    t[2]=t[1];
-  } else if(t.empty()) {
-    if( num>0 && v1.empty() && v2.empty() && axis.empty() ) {
-      return;
-    }
-    if(!(v1.size()==2 && v2.size()==2 && axis.size()==2)) {
-      aa->error("VECTORA, VECTORB and AXIS should specify 2 atoms each");
-    }
-    aa->log.printf("  between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n",
-                   v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial());
-    t.resize(6);
-    t[0]=v1[1];
-    t[1]=v1[0];
-    t[2]=axis[0];
-    t[3]=axis[1];
-    t[4]=v2[0];
-    t[5]=v2[1];
-  } else if( t.size()!=4 ) {
-    aa->error("ATOMS should specify 4 atoms");
-  }
-}
-
-unsigned Torsion::getModeAndSetupValues( ActionWithValue* av ) {
-  bool do_cos;
-  av->parseFlag("COSINE",do_cos);
-  if(do_cos) {
-    av->log.printf("  calculating cosine instead of torsion\n");
-  }
-
-  av->addValueWithDerivatives();
-  if(!do_cos) {
-    av->setPeriodic("-pi","pi");
-    return 0;
-  }
-  av->setNotPeriodic();
-  return 1;
-}
-
-// calculator
-void Torsion::calculate() {
-  if(pbc) {
-    makeWhole();
-  }
-  ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this);
-  if(do_cosine) {
-    calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), cvout );
-  } else {
-    calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout );
-  }
-  for(unsigned i=0; i<6; ++i) {
-    setAtomsDerivatives(i,cvout.getAtomDerivatives(0,i) );
-  }
-  setValue(value[0]);
-  setBoxDerivatives( cvout.virial[0] );
-}
-
-void Torsion::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
-  Vector d0=delta(cvin.pos[1],cvin.pos[0]);
-  Vector d1=delta(cvin.pos[3],cvin.pos[2]);
-  Vector d2=delta(cvin.pos[5],cvin.pos[4]);
-  Vector dd0,dd1,dd2;
-  PLMD::Torsion t;
-  cvout.values[0] = t.compute(d0,d1,d2,dd0,dd1,dd2);
-  if(cvin.mode==1) {
-    dd0 *= -std::sin(cvout.values[0]);
-    dd1 *= -std::sin(cvout.values[0]);
-    dd2 *= -std::sin(cvout.values[0]);
-    cvout.values[0] = std::cos(cvout.values[0]);
-  }
-  cvout.derivs[0][0] = dd0;
-  cvout.derivs[0][1] = -dd0;
-  cvout.derivs[0][2] = dd1;
-  cvout.derivs[0][3] = -dd1;
-  cvout.derivs[0][4] = dd2;
-  cvout.derivs[0][5] = -dd2;
-  ColvarInput::setBoxDerivativesNoPbc( cvin, cvout );
-}
-
 }
 }
 
-
-
diff --git a/src/colvar/Torsion.h b/src/colvar/Torsion.h
new file mode 100644
index 0000000000..61475dde0f
--- /dev/null
+++ b/src/colvar/Torsion.h
@@ -0,0 +1,204 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2011-2023 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "Colvar.h"
+#include "tools/Torsion.h"
+#include "ColvarInput.h"
+
+namespace PLMD {
+namespace colvar {
+
+class Torsion : public Colvar {
+  bool pbc;
+  bool do_cosine;
+
+  std::vector<double> value;
+  std::vector<double> derivs;
+public:
+  explicit Torsion(const ActionOptions&);
+  static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
+  static unsigned getModeAndSetupValues( ActionWithValue* av );
+// active methods:
+  void calculate() override;
+  static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout );
+  static void registerKeywords(Keywords& keys);
+};
+
+void Torsion::registerKeywords(Keywords& keys) {
+  Colvar::registerKeywords( keys );
+  keys.setDisplayName("TORSION");
+  keys.add("atoms-1","ATOMS","the four atoms involved in the torsional angle");
+  keys.add("atoms-2","AXIS","two atoms that define an axis.  You can use this to find the angle in the plane perpendicular to the axis between the vectors specified using the VECTORA and VECTORB keywords.");
+  keys.add("atoms-2","VECTORA","two atoms that define a vector.  You can use this in combination with VECTOR2 and AXIS");
+  keys.add("atoms-2","VECTORB","two atoms that define a vector.  You can use this in combination with VECTOR1 and AXIS");
+  keys.addDeprecatedKeyword("VECTOR1","VECTORA");
+  keys.addDeprecatedKeyword("VECTOR2","VECTORB");
+  keys.addFlag("COSINE",false,"calculate cosine instead of dihedral");
+  keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
+  keys.setValueDescription("scalar/vector","the TORSION involving these atoms");
+  keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
+}
+
+Torsion::Torsion(const ActionOptions&ao):
+  PLUMED_COLVAR_INIT(ao),
+  pbc(true),
+  do_cosine(false),
+  value(1) {
+  std::vector<AtomNumber> atoms;
+  std::vector<AtomNumber> v1;
+  ActionAtomistic::parseAtomList("VECTOR1",v1);
+  if( v1.size()>0 ) {
+    std::vector<AtomNumber> v2;
+    ActionAtomistic::parseAtomList("VECTOR2",v2);
+    std::vector<AtomNumber> axis;
+    ActionAtomistic::parseAtomList("AXIS",axis);
+    if( !(v1.size()==2 && v2.size()==2 && axis.size()==2)) {
+      error("VECTOR1, VECTOR2 and AXIS should specify 2 atoms each");
+    }
+    atoms.resize(6);
+    atoms[0]=v1[1];
+    atoms[1]=v1[0];
+    atoms[2]=axis[0];
+    atoms[3]=axis[1];
+    atoms[4]=v2[0];
+    atoms[5]=v2[1];
+    log.printf("  between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n",
+               v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial());
+  } else {
+    parseAtomList(-1,atoms,this);
+  }
+  unsigned mode=getModeAndSetupValues(this);
+  if( mode==1 ) {
+    do_cosine=true;
+  }
+
+  bool nopbc=!pbc;
+  parseFlag("NOPBC",nopbc);
+  pbc=!nopbc;
+  checkRead();
+
+  if(pbc) {
+    log.printf("  using periodic boundary conditions\n");
+  } else {
+    log.printf("  without periodic boundary conditions\n");
+  }
+  requestAtoms(atoms);
+}
+
+void Torsion::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
+  std::vector<AtomNumber> v1,v2,axis;
+  aa->parseAtomList("ATOMS",num,t);
+  aa->parseAtomList("VECTORA",num,v1);
+  aa->parseAtomList("VECTORB",num,v2);
+  aa->parseAtomList("AXIS",num,axis);
+
+  if(t.size()==4) {
+    if(!(v1.empty() && v2.empty() && axis.empty())) {
+      aa->error("ATOMS keyword is not compatible with VECTORA, VECTORB and AXIS keywords");
+    }
+    aa->log.printf("  between atoms %d %d %d %d\n",t[0].serial(),t[1].serial(),t[2].serial(),t[3].serial());
+    t.resize(6);
+    t[5]=t[3];
+    t[4]=t[2];
+    t[3]=t[2];
+    t[2]=t[1];
+  } else if(t.empty()) {
+    if( num>0 && v1.empty() && v2.empty() && axis.empty() ) {
+      return;
+    }
+    if(!(v1.size()==2 && v2.size()==2 && axis.size()==2)) {
+      aa->error("VECTORA, VECTORB and AXIS should specify 2 atoms each");
+    }
+    aa->log.printf("  between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n",
+                   v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial());
+    t.resize(6);
+    t[0]=v1[1];
+    t[1]=v1[0];
+    t[2]=axis[0];
+    t[3]=axis[1];
+    t[4]=v2[0];
+    t[5]=v2[1];
+  } else if( t.size()!=4 ) {
+    aa->error("ATOMS should specify 4 atoms");
+  }
+}
+
+unsigned Torsion::getModeAndSetupValues( ActionWithValue* av ) {
+  bool do_cos;
+  av->parseFlag("COSINE",do_cos);
+  if(do_cos) {
+    av->log.printf("  calculating cosine instead of torsion\n");
+  }
+
+  av->addValueWithDerivatives();
+  if(!do_cos) {
+    av->setPeriodic("-pi","pi");
+    return 0;
+  }
+  av->setNotPeriodic();
+  return 1;
+}
+
+// calculator
+void Torsion::calculate() {
+  if(pbc) {
+    makeWhole();
+  }
+  ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this);
+  if(do_cosine) {
+    calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), cvout );
+  } else {
+    calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout );
+  }
+  for(unsigned i=0; i<6; ++i) {
+    setAtomsDerivatives(i,cvout.getAtomDerivatives(0,i) );
+  }
+  setValue(value[0]);
+  setBoxDerivatives( cvout.virial[0] );
+}
+
+void Torsion::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
+  Vector d0=delta(cvin.pos[1],cvin.pos[0]);
+  Vector d1=delta(cvin.pos[3],cvin.pos[2]);
+  Vector d2=delta(cvin.pos[5],cvin.pos[4]);
+  Vector dd0,dd1,dd2;
+  PLMD::Torsion t;
+  cvout.values[0] = t.compute(d0,d1,d2,dd0,dd1,dd2);
+  if(cvin.mode==1) {
+    dd0 *= -std::sin(cvout.values[0]);
+    dd1 *= -std::sin(cvout.values[0]);
+    dd2 *= -std::sin(cvout.values[0]);
+    cvout.values[0] = std::cos(cvout.values[0]);
+  }
+  cvout.derivs[0][0] = dd0;
+  cvout.derivs[0][1] = -dd0;
+  cvout.derivs[0][2] = dd1;
+  cvout.derivs[0][3] = -dd1;
+  cvout.derivs[0][4] = dd2;
+  cvout.derivs[0][5] = -dd2;
+  ColvarInput::setBoxDerivativesNoPbc( cvin, cvout );
+}
+
+}
+}
+
+
+

From 157feafe08adce73dc7bffd43743d1843e096e96 Mon Sep 17 00:00:00 2001
From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com>
Date: Mon, 1 Sep 2025 16:45:44 +0200
Subject: [PATCH 9/9] added the forgotten compiler guards

---
 src/colvar/Angle.h               | 3 +++
 src/colvar/DihedralCorrelation.h | 3 +++
 src/colvar/Dipole.h              | 3 +++
 src/colvar/Plane.h               | 3 +++
 src/colvar/Position.h            | 6 +++---
 src/colvar/Torsion.h             | 6 +++---
 6 files changed, 18 insertions(+), 6 deletions(-)

diff --git a/src/colvar/Angle.h b/src/colvar/Angle.h
index 6d9d0981be..d652a327f7 100644
--- a/src/colvar/Angle.h
+++ b/src/colvar/Angle.h
@@ -19,6 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_colvar_Angle_h
+#define __PLUMED_colvar_Angle_h
 #include "Colvar.h"
 #include "ColvarInput.h"
 #include "tools/Angle.h"
@@ -123,3 +125,4 @@ void Angle::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
 } //namespace colvar
 } //namespace PLMD
 
+#endif //__PLUMED_colvar_Angle_h
diff --git a/src/colvar/DihedralCorrelation.h b/src/colvar/DihedralCorrelation.h
index 0c6310ec4f..1d10571a17 100644
--- a/src/colvar/DihedralCorrelation.h
+++ b/src/colvar/DihedralCorrelation.h
@@ -19,6 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_colvar_DihedralCorrelation_h
+#define __PLUMED_colvar_DihedralCorrelation_h
 #include "Colvar.h"
 #include "ColvarInput.h"
 #include "tools/Torsion.h"
@@ -152,3 +154,4 @@ void DihedralCorrelation::calculateCV( const ColvarInput& cvin, ColvarOutput& cv
 
 }
 }
+#endif  //__PLUMED_colvar_DihedralCorrelation_h
diff --git a/src/colvar/Dipole.h b/src/colvar/Dipole.h
index 96bd2193d9..21f389574e 100644
--- a/src/colvar/Dipole.h
+++ b/src/colvar/Dipole.h
@@ -20,6 +20,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_colvar_Dipole_h
+#define __PLUMED_colvar_Dipole_h
 #include "Colvar.h"
 #include "ColvarInput.h"
 
@@ -182,3 +184,4 @@ void Dipole::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
 
 }
 }
+#endif //__PLUMED_colvar_Dipole_h
diff --git a/src/colvar/Plane.h b/src/colvar/Plane.h
index 431b3c1e14..872ad4d723 100644
--- a/src/colvar/Plane.h
+++ b/src/colvar/Plane.h
@@ -19,6 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_colvar_Plane_h
+#define __PLUMED_colvar_Plane_h
 #include "Colvar.h"
 #include "ColvarInput.h"
 #include "tools/Pbc.h"
@@ -144,3 +146,4 @@ void Plane::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
 }
 }
 
+#endif // __PLUMED_colvar_Plane_h
diff --git a/src/colvar/Position.h b/src/colvar/Position.h
index 8325886d32..ce56cdc7a8 100644
--- a/src/colvar/Position.h
+++ b/src/colvar/Position.h
@@ -19,6 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_colvar_Position_h
+#define __PLUMED_colvar_Position_h
 #include "Colvar.h"
 #include "tools/Pbc.h"
 #include "ColvarInput.h"
@@ -181,6 +183,4 @@ void Position::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
 
 }
 }
-
-
-
+#endif //__PLUMED_colvar_Position_h
diff --git a/src/colvar/Torsion.h b/src/colvar/Torsion.h
index 61475dde0f..00a2b9778d 100644
--- a/src/colvar/Torsion.h
+++ b/src/colvar/Torsion.h
@@ -19,6 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_colvar_Torsion_h
+#define __PLUMED_colvar_Torsion_h
 #include "Colvar.h"
 #include "tools/Torsion.h"
 #include "ColvarInput.h"
@@ -199,6 +201,4 @@ void Torsion::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) {
 
 }
 }
-
-
-
+#endif // __PLUMED_colvar_Torsion.h