diff --git a/regtest/crystdistrib/rt-dops-shortcut/config b/regtest/crystdistrib/rt-dops-shortcut/config
index 557dccbe32..d09ec688c6 100644
--- a/regtest/crystdistrib/rt-dops-shortcut/config
+++ b/regtest/crystdistrib/rt-dops-shortcut/config
@@ -1,4 +1,4 @@
 type=driver
-plumed_modules=crystdistrib
+plumed_modules="crystdistrib adjmat"
 arg="--plumed plumed.dat --mf_xtc em.xtc --mc mcfile --timestep=0.001 --trajectory-stride=25"
 extra_files="../../trajectories/em.xtc"
diff --git a/regtest/multicolvar/rt-input/Makefile b/regtest/multicolvar/rt-input/Makefile
new file mode 100644
index 0000000000..3703b27cea
--- /dev/null
+++ b/regtest/multicolvar/rt-input/Makefile
@@ -0,0 +1 @@
+include ../../scripts/test.make
diff --git a/regtest/multicolvar/rt-input/config b/regtest/multicolvar/rt-input/config
new file mode 100644
index 0000000000..df1f95bf3e
--- /dev/null
+++ b/regtest/multicolvar/rt-input/config
@@ -0,0 +1 @@
+type=make
diff --git a/regtest/multicolvar/rt-input/main.cpp b/regtest/multicolvar/rt-input/main.cpp
new file mode 100644
index 0000000000..d0dd3fdc5e
--- /dev/null
+++ b/regtest/multicolvar/rt-input/main.cpp
@@ -0,0 +1,201 @@
+#include "plumed/colvar/MultiColvarTemplate.h"
+#include "plumed/tools/Vector.h"
+
+#include <fstream>
+#include <iostream>
+#include <variant>
+#include <vector>
+
+using namespace PLMD;
+using colvar::multiColvars::Input;
+
+template <typename T> void pv(std::ofstream &stream, const std::vector<T> &v) {
+  // todo: put this is a common header
+  std::string div = "";
+  stream << "(" << v.size() << ") ";
+  for (const auto &t : v) {
+    stream << div << t;
+    div = ", ";
+  }
+}
+
+void print(std::ofstream &stream, const std::vector<double> &masses,
+           const std::vector<double> &charges,
+           const std::vector<Vector> positions) {
+  stream << "Masses: ";
+  pv(stream, masses);
+  stream << "\nCharges: ";
+  pv(stream, charges);
+  stream << "\nPositions: ";
+  pv(stream, positions);
+  stream << "\n";
+}
+
+void inputTest(std::ofstream &stream, Input input) {
+  try {
+    const auto &masses = input.masses();
+    stream << " Masses: ";
+    pv(stream, masses);
+    stream << "\n";
+  } catch (std::bad_variant_access &) {
+    stream << " There are no masses\n";
+  }
+
+  try {
+    const auto &charges = input.charges();
+    stream << " Charges: ";
+    pv(stream, charges);
+    stream << "\n";
+  } catch (std::bad_variant_access &) {
+    stream << " There are no charges\n";
+  }
+  try {
+    const auto &positions = input.positions();
+    stream << " Positions: ";
+    pv(stream, positions);
+    stream << "\n";
+  } catch (std::bad_variant_access &) {
+    stream << " There are no positions\n";
+  }
+  stream << std::endl;
+}
+
+void workOnMass(Input input) {
+  auto &masses = input.var_masses();
+  for (auto &m : masses) {
+    m *= 2;
+  }
+}
+
+void workOnCharges(Input input) {
+  auto &charges = input.var_charges();
+  for (auto &c : charges) {
+    c *= 2;
+  }
+}
+
+void workOnPositions(Input input) {
+  auto &pos = input.var_positions();
+  for (auto &p : pos) {
+    p = p * 2.0;
+  }
+}
+
+class testData {
+  using vd = std::vector<double>;
+  using vv = std::vector<Vector>;
+  vv p = {{0, 0, 0}, {1, 1, 1}};
+  vd c = {-1, 1};
+  vd m = {3, 4};
+
+public:
+  testData() = default;
+  vd &masses() { return m; }
+  vd &charges() { return c; }
+  vv &positions() { return p; }
+  const vd &c_masses() const { return m; }
+  const vd &c_charges() const { return c; }
+  const vv &c_positions() const { return p; }
+};
+
+class test {
+  using vd = std::vector<double>;
+  using vv = std::vector<Vector>;
+
+public:
+  test() = default;
+  void dotestWithReference(std::ofstream &output) {
+    testData td;
+    output << "*dotestWithReference:" << std::endl;
+    output << "*Passing full input:\n";
+    inputTest(output, Input(td.masses(), td.charges(), td.positions()));
+    output << "*Passing only positions:\n";
+    inputTest(output, Input().positions(td.positions()));
+    output << "*Passing only masses:\n";
+    inputTest(output, Input().masses(td.masses()));
+    output << "*Passing only charges:\n";
+    inputTest(output, Input().charges(td.charges()));
+  }
+  void dotestWithConstReference(std::ofstream &output) {
+    testData td;
+    output << "*dotestWithConstReference:" << std::endl;
+    output << "*Passing full input:\n";
+    inputTest(output, Input(td.c_masses(), td.c_charges(), td.c_positions()));
+    output << "*Passing only positions:\n";
+    inputTest(output, Input().positions(td.c_positions()));
+    output << "*Passing only masses:\n";
+    inputTest(output, Input().masses(td.c_masses()));
+    output << "*Passing only charges:\n";
+    inputTest(output, Input().charges(td.c_charges()));
+  }
+
+  void testUsingInputAsBuffer_masses(std::ofstream &output) {
+    testData td;
+    output << "*testUsingInputAsBuffer_masses:" << std::endl;
+    output << "*Passing constant_reference:\n";
+    try {
+      workOnMass(Input(td.c_masses(), td.c_charges(), td.c_positions()));
+    } catch (std::bad_variant_access &) {
+      output << " should throw\n";
+    }
+    output << "*Passing reference:\n";
+    output << " Initial values\n Masses:";
+    pv(output, td.masses());
+    workOnMass(Input(td.masses(), td.charges(), td.positions()));
+    output << "\n";
+    // this extra line is to remember that something has changed!!!
+    output << " Values after call\n Masses:";
+    pv(output, td.masses());
+    output << std::endl;
+  }
+  void testUsingInputAsBuffer_charges(std::ofstream &output) {
+    testData td;
+    output << "*testUsingInputAsBuffer_charges:" << std::endl;
+    output << "*Passing constant_reference:\n";
+    try {
+      workOnCharges(Input(td.c_masses(), td.c_charges(), td.c_positions()));
+    } catch (std::bad_variant_access &) {
+      output << " should throw\n";
+    }
+    output << "*Passing reference:\n";
+    output << " Initial values\n Charges:";
+    pv(output, td.c_charges());
+    workOnCharges(Input(td.masses(), td.charges(), td.positions()));
+    output << "\n";
+    // this extra line is to remember that something has changed!!!
+    output << " Values after call\n Charges:";
+    pv(output, td.c_charges());
+    output << std::endl;
+  }
+  void testUsingInputAsBuffer_positions(std::ofstream &output) {
+    testData td;
+    output << "*testUsingInputAsBuffer_positions:" << std::endl;
+    output << "*Passing constant_reference:\n";
+    try {
+      workOnPositions(Input(td.c_masses(), td.c_charges(), td.c_positions()));
+    } catch (std::bad_variant_access &) {
+      output << " should throw\n";
+    }
+    output << "*Passing reference:\n";
+    output << " Initial values\n Positions:";
+    pv(output, td.positions());
+    workOnPositions(Input(td.masses(), td.charges(), td.positions()));
+    output << "\n";
+    // this extra line is to remember that something has changed!!!
+    output << " Values after call\n Positions:";
+    pv(output, td.positions());
+    output << std::endl;
+  }
+};
+
+int main() {
+  std::ofstream output("output");
+
+  test t;
+  t.dotestWithReference(output);
+  t.dotestWithConstReference(output);
+  t.testUsingInputAsBuffer_masses(output);
+  t.testUsingInputAsBuffer_charges(output);
+  t.testUsingInputAsBuffer_positions(output);
+  return 0;
+}
diff --git a/regtest/multicolvar/rt-input/output.reference b/regtest/multicolvar/rt-input/output.reference
new file mode 100644
index 0000000000..43c6532d27
--- /dev/null
+++ b/regtest/multicolvar/rt-input/output.reference
@@ -0,0 +1,66 @@
+*dotestWithReference:
+*Passing full input:
+ Masses: (2) 3, 4
+ Charges: (2) -1, 1
+ Positions: (2) 0 0 0, 1 1 1
+
+*Passing only positions:
+ There are no masses
+ There are no charges
+ Positions: (2) 0 0 0, 1 1 1
+
+*Passing only masses:
+ Masses: (2) 3, 4
+ There are no charges
+ There are no positions
+
+*Passing only charges:
+ There are no masses
+ Charges: (2) -1, 1
+ There are no positions
+
+*dotestWithConstReference:
+*Passing full input:
+ Masses: (2) 3, 4
+ Charges: (2) -1, 1
+ Positions: (2) 0 0 0, 1 1 1
+
+*Passing only positions:
+ There are no masses
+ There are no charges
+ Positions: (2) 0 0 0, 1 1 1
+
+*Passing only masses:
+ Masses: (2) 3, 4
+ There are no charges
+ There are no positions
+
+*Passing only charges:
+ There are no masses
+ Charges: (2) -1, 1
+ There are no positions
+
+*testUsingInputAsBuffer_masses:
+*Passing constant_reference:
+ should throw
+*Passing reference:
+ Initial values
+ Masses:(2) 3, 4
+ Values after call
+ Masses:(2) 6, 8
+*testUsingInputAsBuffer_charges:
+*Passing constant_reference:
+ should throw
+*Passing reference:
+ Initial values
+ Charges:(2) -1, 1
+ Values after call
+ Charges:(2) -2, 2
+*testUsingInputAsBuffer_positions:
+*Passing constant_reference:
+ should throw
+*Passing reference:
+ Initial values
+ Positions:(2) 0 0 0, 1 1 1
+ Values after call
+ Positions:(2) 0 0 0, 2 2 2
diff --git a/src/colvar/Angle.cpp b/src/colvar/Angle.cpp
index 846e0432f7..7f458330df 100644
--- a/src/colvar/Angle.cpp
+++ b/src/colvar/Angle.cpp
@@ -100,19 +100,16 @@ Calculate multiple angles.
 
 class Angle : public Colvar {
   bool pbc;
-  std::vector<double> value, masses, charges;
+  std::vector<double> value;
   std::vector<std::vector<Vector> > derivs;
   std::vector<Tensor> virial;
 public:
+  MULTICOLVAR_DEFAULT(multiColvars::emptyMode);
   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 unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
-                           const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
+
 };
 
 typedef ColvarShortcut<Angle> AngleShortcut;
@@ -128,7 +125,7 @@ void Angle::registerKeywords( Keywords& keys ) {
   keys.setValueDescription("the ANGLE involving these atoms");
 }
 
-void Angle::parseAtomList( const int& num, std::vector<AtomNumber>& atoms, ActionAtomistic* aa ) {
+void Angle::parseAtomList( int const 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());
@@ -140,8 +137,10 @@ void Angle::parseAtomList( const int& num, std::vector<AtomNumber>& atoms, Actio
   } 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::Modetype Angle::getModeAndSetupValues( ActionWithValue* av ) {
+  av->addValueWithDerivatives();
+  av->setNotPeriodic();
+  return {};
 }
 
 Angle::Angle(const ActionOptions&ao):
@@ -169,19 +168,26 @@ Angle::Angle(const ActionOptions&ao):
 void Angle::calculate() {
 
   if(pbc) makeWhole();
-  calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
+  calculateCV( {}, multiColvars::Input().positions(getPositions()), multiColvars::Ouput(value, derivs, virial), this );
   setValue( value[0] );
-  for(unsigned i=0; i<derivs[0].size(); ++i) setAtomsDerivatives( i, derivs[0][i] );
+  for(unsigned i=0; i<derivs[0].size(); ++i)
+    setAtomsDerivatives( i, derivs[0][i] );
   setBoxDerivatives( virial[0] );
 }
 
-void Angle::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
-                         const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                         std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
+void Angle::calculateCV( Modetype /*mode*/,
+                         multiColvars::Input const in,
+                         multiColvars::Ouput out,
+                         const ActionAtomistic* aa ) {
+  auto & vals=out.vals();
+  auto & derivs=out.derivs();
+  auto & virial=out.virial();
+  const auto & pos=in.positions();
   Vector dij,dik;
   dij=delta(pos[2],pos[3]);
   dik=delta(pos[1],pos[0]);
-  Vector ddij,ddik; PLMD::Angle a;
+  Vector ddij,ddik;
+  PLMD::Angle a;
   vals[0]=a.compute(dij,dik,ddij,ddik);
   derivs[0][0]=ddik; derivs[0][1]=-ddik;
   derivs[0][2]=-ddij; derivs[0][3]=ddij;
diff --git a/src/colvar/DihedralCorrelation.cpp b/src/colvar/DihedralCorrelation.cpp
index 19e08feb1b..41ac7388bd 100644
--- a/src/colvar/DihedralCorrelation.cpp
+++ b/src/colvar/DihedralCorrelation.cpp
@@ -64,18 +64,14 @@ Measure the correlation between a multiple pairs of dihedral angles
 class DihedralCorrelation : public Colvar {
 private:
   bool pbc;
-  std::vector<double> value, masses, charges;
+  std::vector<double> value;
   std::vector<std::vector<Vector> > derivs;
   std::vector<Tensor> virial;
 public:
+  MULTICOLVAR_DEFAULT(multiColvars::emptyMode);
   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 unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
-                           const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
 };
 
 typedef ColvarShortcut<DihedralCorrelation> DihedralCorrelationShortcut;
@@ -110,7 +106,7 @@ DihedralCorrelation::DihedralCorrelation(const ActionOptions&ao):
   else    log.printf("  without periodic boundary conditions\n");
 }
 
-void DihedralCorrelation::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
+void DihedralCorrelation::parseAtomList( int const 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 ) {
@@ -119,22 +115,28 @@ void DihedralCorrelation::parseAtomList( const int& num, std::vector<AtomNumber>
   }
 }
 
-unsigned DihedralCorrelation::getModeAndSetupValues( ActionWithValue* av ) {
-  av->addValueWithDerivatives(); av->setNotPeriodic(); return 0;
+DihedralCorrelation::Modetype DihedralCorrelation::getModeAndSetupValues( ActionWithValue* av ) {
+  av->addValueWithDerivatives();
+  av->setNotPeriodic();
+  return {};
 }
 
 void DihedralCorrelation::calculate() {
 
   if(pbc) makeWhole();
-  calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
+  calculateCV( {}, multiColvars::Input().positions(getPositions()), multiColvars::Ouput(value, derivs, virial), this );
   setValue( value[0] );
   for(unsigned i=0; i<derivs[0].size(); ++i) setAtomsDerivatives( i, derivs[0][i] );
   setBoxDerivatives( virial[0] );
 }
 
-void DihedralCorrelation::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
-                                       const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                                       std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
+void DihedralCorrelation::calculateCV(Modetype  /*mode*/,
+                                      multiColvars::Input const in,
+                                      multiColvars::Ouput out, const ActionAtomistic* aa ) {
+  auto & vals=out.vals();
+  auto & derivs=out.derivs();
+  auto & virial=out.virial();
+  const auto& pos = in.positions();
   const Vector d10=delta(pos[1],pos[0]);
   const Vector d11=delta(pos[2],pos[1]);
   const Vector d12=delta(pos[3],pos[2]);
diff --git a/src/colvar/Dipole.cpp b/src/colvar/Dipole.cpp
index e0a9212455..2ef27f7736 100644
--- a/src/colvar/Dipole.cpp
+++ b/src/colvar/Dipole.cpp
@@ -79,24 +79,22 @@ Calculate a vector of dipole moments for a set of groups of atoms.
 //+ENDPLUMEDOC
 
 class Dipole : public Colvar {
+public:
+  explicit Dipole(const ActionOptions&);
+  void calculate() override;
+  static void registerKeywords(Keywords& keys);
+  MULTICOLVAR_DEFAULT(multiColvars::components);
+private:
   std::vector<AtomNumber> ga_lista;
-  bool components;
-  bool nopbc;
-  std::vector<double> value, masses, charges;
+  std::vector<double> value;
+  std::vector<double> charges;
   std::vector<std::vector<Vector> > derivs;
   std::vector<Tensor> virial;
   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 unsigned& mode, const std::vector<double>& masses, std::vector<double>& charges,
-                           const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
+  Modetype components;
+  bool nopbc=false;
 };
 
 typedef ColvarShortcut<Dipole> DipoleShortcut;
@@ -118,20 +116,24 @@ void Dipole::registerKeywords(Keywords& keys) {
 
 Dipole::Dipole(const ActionOptions&ao):
   PLUMED_COLVAR_INIT(ao),
-  components(false),
+
   value(1),
   derivs(1),
-  virial(1)
+  virial(1),
+  components(getModeAndSetupValues(this))
 {
-  parseAtomList(-1,ga_lista,this); charges.resize(ga_lista.size());
-  components=(getModeAndSetupValues(this)==1);
-  if( components ) {
+  parseAtomList(-1,ga_lista,this);
+  charges.resize(ga_lista.size());
+  // components=getModeAndSetupValues(this);
+  if( components == Modetype::withCompontents ) {
     value.resize(3); derivs.resize(3); virial.resize(3);
     valuex=getPntrToComponent("x");
     valuey=getPntrToComponent("y");
     valuez=getPntrToComponent("z");
   }
-  for(unsigned i=0; i<derivs.size(); ++i) derivs[i].resize( ga_lista.size() );
+  for(unsigned i=0; i<derivs.size(); ++i) {
+    derivs[i].resize( ga_lista.size() );
+  }
   parseFlag("NOPBC",nopbc);
   checkRead();
 
@@ -141,7 +143,7 @@ Dipole::Dipole(const ActionOptions&ao):
   requestAtoms(ga_lista);
 }
 
-void Dipole::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
+void Dipole::parseAtomList( 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()));
@@ -152,15 +154,17 @@ void Dipole::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAt
   }
 }
 
-unsigned Dipole::getModeAndSetupValues( ActionWithValue* av ) {
-  bool c; av->parseFlag("COMPONENTS",c);
+Dipole::Modetype 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;
+    return Modetype::withCompontents;
   }
-  av->addValueWithDerivatives(); av->setNotPeriodic(); return 0;
+  av->addValueWithDerivatives(); av->setNotPeriodic();
+  return Modetype::noComponents;
 }
 
 // calculator
@@ -169,14 +173,14 @@ void Dipole::calculate()
   if(!nopbc) makeWhole();
   unsigned N=getNumberOfAtoms();
   for(unsigned i=0; i<N; ++i) charges[i]=getCharge(i);
-
-  if(!components) {
-    calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
+  auto inputs = multiColvars::Input().positions(getPositions())
+                .charges(charges);
+  calculateCV( components, inputs, multiColvars::Ouput(value, derivs, virial), this );
+  if(components == Modetype::noComponents) {
     for(unsigned i=0; i<N; i++) setAtomsDerivatives(i,derivs[0][i]);
     setBoxDerivatives(virial[0]);
     setValue(value[0]);
   } else {
-    calculateCV( 1, masses, charges, getPositions(), value, derivs, virial, this );
     for(unsigned i=0; i<N; i++) {
       setAtomsDerivatives(valuex,i,derivs[0][i]);
       setAtomsDerivatives(valuey,i,derivs[1][i]);
@@ -191,25 +195,40 @@ void Dipole::calculate()
   }
 }
 
-void Dipole::calculateCV( const unsigned& mode, const std::vector<double>& masses, std::vector<double>& charges,
-                          const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                          std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
-  unsigned N=pos.size(); double ctot=0.;
-  for(unsigned i=0; i<N; ++i) ctot += charges[i];
+void Dipole::calculateCV( Modetype mode,
+                          multiColvars::Input const in,
+                          multiColvars::Ouput out,
+                          const ActionAtomistic* aa ) {
+  auto & vals=out.vals();
+  auto & derivs=out.derivs();
+  auto & virial=out.virial();
+  const auto & pos=in.positions();
+  auto & charges=in.var_charges();
+  unsigned N=pos.size();
+  double ctot=0.;
+  for(unsigned i=0; i<N; ++i) {
+    ctot += charges[i];
+  }
+  //check if
+  //std::accumulate(charges.begin(),charges.end(),ctot);
+  // is faster
   ctot/=(double)N;
 
   Vector dipje;
   for(unsigned i=0; i<N; ++i) {
-    charges[i]-=ctot; dipje += charges[i]*pos[i];
+    charges[i]-=ctot;
+    dipje += charges[i]*pos[i];
   }
 
-  if( mode==1 ) {
+  if( mode==Modetype::withCompontents ) {
     for(unsigned i=0; i<N; i++) {
-      derivs[0][i]=charges[i]*Vector(1.0,0.0,0.0);
-      derivs[1][i]=charges[i]*Vector(0.0,1.0,0.0);
-      derivs[2][i]=charges[i]*Vector(0.0,0.0,1.0);
+      derivs[0][i]=Vector(charges[i],0.0,0.0);
+      derivs[1][i]=Vector(0.0,charges[i],0.0);
+      derivs[2][i]=Vector(0.0,0.0,charges[i]);
+    }
+    for(unsigned i=0; i<3; ++i ) {
+      vals[i] = dipje[i];
     }
-    for(unsigned i=0; i<3; ++i ) vals[i] = dipje[i];
   } else {
     vals[0] = dipje.modulo();
     double idip = 1./vals[0];
diff --git a/src/colvar/Distance.cpp b/src/colvar/Distance.cpp
index b2d75a27fe..0dd403753c 100644
--- a/src/colvar/Distance.cpp
+++ b/src/colvar/Distance.cpp
@@ -124,23 +124,18 @@ Calculate a vector containing the distances between various pairs of atoms
 //+ENDPLUMEDOC
 
 class Distance : public Colvar {
-  bool components;
-  bool scaled_components;
-  bool pbc;
-
-  std::vector<double> value, masses, charges;
-  std::vector<std::vector<Vector> > derivs;
-  std::vector<Tensor> virial;
 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 unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
-                           const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
+  MULTICOLVAR_DEFAULT(multiColvars::scaledComponents);
+private:
+  std::vector<double> value{{0.0}};
+  std::vector<std::vector<Vector> > derivs{std::vector<Vector>{Vector{},Vector{}}};
+  std::vector<Tensor> virial{Tensor()};
+  Modetype mode_;
+  bool pbc=true;
 };
 
 typedef ColvarShortcut<Distance> DistanceShortcut;
@@ -165,13 +160,7 @@ void Distance::registerKeywords( Keywords& keys ) {
 }
 
 Distance::Distance(const ActionOptions&ao):
-  PLUMED_COLVAR_INIT(ao),
-  components(false),
-  scaled_components(false),
-  pbc(true),
-  value(1),
-  derivs(1),
-  virial(1)
+  PLUMED_COLVAR_INIT(ao)
 {
   derivs[0].resize(2);
   std::vector<AtomNumber> atoms;
@@ -186,91 +175,120 @@ Distance::Distance(const ActionOptions&ao):
   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); derivs.resize(3); virial.resize(3);
-    for(unsigned i=0; i<3; ++i) derivs[i].resize(2);
+  mode_ = getModeAndSetupValues( this );
+  if(mode_==Modetype::withCompontents || mode_==Modetype::scaledComponents) {
+    value.resize(3);
+    derivs.resize(3);
+    virial.resize(3);
+    for(unsigned i=0; i<3; ++i)
+      derivs[i].resize(2);
   }
   requestAtoms(atoms);
 }
 
-void Distance::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
+void Distance::parseAtomList( int const 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");
-
+Distance::Modetype 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;
+    return Modetype::withCompontents;
   } 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;
+    return Modetype::scaledComponents;
   }
   av->addValueWithDerivatives(); av->setNotPeriodic();
-  return 0;
+  return Modetype::noComponents;
 }
 
 // calculator
 void Distance::calculate() {
 
   if(pbc) makeWhole();
+  calculateCV( mode_, multiColvars::Input().positions(getPositions()), multiColvars::Ouput(value, derivs, virial), this );
 
-  if( components ) {
-    calculateCV( 1, masses, charges, getPositions(), value, derivs, virial, this );
+  switch (mode_) {
+  case Modetype::withCompontents: {
     Value* valuex=getPntrToComponent("x");
     Value* valuey=getPntrToComponent("y");
     Value* valuez=getPntrToComponent("z");
 
-    for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuex,i,derivs[0][i] );
+    for(unsigned i=0; i<2; ++i) {
+      setAtomsDerivatives(valuex,i,derivs[0][i] );
+    }
     setBoxDerivatives(valuex,virial[0]);
     valuex->set(value[0]);
 
-    for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuey,i,derivs[1][i] );
+    for(unsigned i=0; i<2; ++i) {
+      setAtomsDerivatives(valuey,i,derivs[1][i] );
+    }
     setBoxDerivatives(valuey,virial[1]);
     valuey->set(value[1]);
 
-    for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuez,i,derivs[2][i] );
+    for(unsigned i=0; i<2; ++i) {
+      setAtomsDerivatives(valuez,i,derivs[2][i] );
+    }
     setBoxDerivatives(valuez,virial[2]);
     valuez->set(value[2]);
-  } else if( scaled_components ) {
-    calculateCV( 2, masses, charges, getPositions(), value, derivs, virial, this );
+  }
+  break;
 
+  case Modetype::scaledComponents: {
     Value* valuea=getPntrToComponent("a");
     Value* valueb=getPntrToComponent("b");
     Value* valuec=getPntrToComponent("c");
-    for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuea,i,derivs[0][i] );
+    for(unsigned i=0; i<2; ++i) {
+      setAtomsDerivatives(valuea,i,derivs[0][i] );
+    }
     valuea->set(value[0]);
-    for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valueb,i,derivs[1][i] );
+    for(unsigned i=0; i<2; ++i) {
+      setAtomsDerivatives(valueb,i,derivs[1][i] );
+    }
     valueb->set(value[1]);
-    for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuec,i,derivs[2][i] );
+    for(unsigned i=0; i<2; ++i) {
+      setAtomsDerivatives(valuec,i,derivs[2][i] );
+    }
     valuec->set(value[2]);
-  } else  {
-    calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
-    for(unsigned i=0; i<2; ++i) setAtomsDerivatives(i,derivs[0][i] );
+  }
+  break;
+
+  case Modetype::noComponents: {
+    for(unsigned i=0; i<2; ++i) {
+      setAtomsDerivatives(i,derivs[0][i] );
+    }
     setBoxDerivatives(virial[0]);
     setValue           (value[0]);
   }
+  }
 }
 
-void Distance::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
-                            const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                            std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
+void Distance::calculateCV( Modetype mode,
+                            multiColvars::Input const in,
+                            multiColvars::Ouput out, const ActionAtomistic* aa ) {
+  auto & vals=out.vals();
+  auto & derivs=out.derivs();
+  auto & virial=out.virial();
+  const auto & pos = in.positions();
   Vector distance=delta(pos[0],pos[1]);
   const double value=distance.modulo();
   const double invvalue=1.0/value;
-
-  if(mode==1) {
+  //the compiler will alert me if I forgot a case :)
+  switch (mode) {
+  case Modetype::withCompontents: {
     derivs[0][0] = Vector(-1,0,0);
     derivs[0][1] = Vector(+1,0,0);
     vals[0] = distance[0];
@@ -283,7 +301,10 @@ void Distance::calculateCV( const unsigned& mode, const std::vector<double>& mas
     derivs[2][1] = Vector(0,0,+1);
     vals[2] = distance[2];
     setBoxDerivativesNoPbc( pos, derivs, virial );
-  } else if(mode==2) {
+  }
+  break;
+
+  case Modetype::scaledComponents: {
     Vector d=aa->getPbc().realToScaled(distance);
     derivs[0][0] = matmul(aa->getPbc().getInvBox(),Vector(-1,0,0));
     derivs[0][1] = matmul(aa->getPbc().getInvBox(),Vector(+1,0,0));
@@ -294,16 +315,18 @@ void Distance::calculateCV( const unsigned& mode, const std::vector<double>& mas
     derivs[2][0] = matmul(aa->getPbc().getInvBox(),Vector(0,0,-1));
     derivs[2][1] = matmul(aa->getPbc().getInvBox(),Vector(0,0,+1));
     vals[2] = Tools::pbc(d[2]);
-  } else {
+  }
+  break;
+
+  case Modetype::noComponents: {
     derivs[0][0] = -invvalue*distance;
     derivs[0][1] = invvalue*distance;
     setBoxDerivativesNoPbc( pos, derivs, virial );
     vals[0] = value;
   }
-}
+  }
 
 }
-}
-
-
+} // namespace colvar
 
+} // namespace PLMD
diff --git a/src/colvar/MultiColvarTemplate.h b/src/colvar/MultiColvarTemplate.h
index bd11210de5..fb1dbcbad6 100644
--- a/src/colvar/MultiColvarTemplate.h
+++ b/src/colvar/MultiColvarTemplate.h
@@ -22,16 +22,215 @@
 #ifndef __PLUMED_colvar_MultiColvarTemplate_h
 #define __PLUMED_colvar_MultiColvarTemplate_h
 
+#include <variant>
+#include <type_traits>
+
 #include "core/ActionWithVector.h"
 
 namespace PLMD {
 namespace colvar {
 
-template <class T>
+namespace multiColvars {
+struct emptyMode {};
+enum class components {
+  withCompontents,
+  noComponents
+};
+enum class plainOrScaled {
+  scaled,
+  plain
+};
+enum class scaledComponents {
+  withCompontents,
+  scaledComponents,
+  noComponents
+};
+
+
+//TODO: test this
+class Ouput final {
+  std::vector<double>* vals_=nullptr;
+  std::vector<std::vector<Vector> >* derivs_=nullptr;
+  std::vector<Tensor>* virial_=nullptr;
+  Ouput()=default;
+public:
+//const because the data in the class is not changing: we are modifying data pointed by the const pointers
+  std::vector<double>&vals() const {return *vals_;}
+  std::vector<std::vector<Vector> >&derivs() const {return *derivs_;}
+  std::vector<Tensor>&virial() const {return *virial_;}
+  Ouput(std::vector<double>& vals,
+        std::vector<std::vector<Vector> >& derivs,
+        std::vector<Tensor>& virial)
+    : vals_(&vals),
+      derivs_(&derivs),
+      virial_(&virial) {}
+  Ouput(Ouput const& other)
+    : vals_(other.vals_),
+      derivs_(other.derivs_),
+      virial_(other.virial_) {};
+  Ouput(Ouput && other) noexcept
+    : Ouput() {
+    swap(*this, other);
+  };
+  // the input is initialized as copy or move ctor
+  Ouput& operator=(Ouput other) =delete;
+  //  {
+  //   //Self assignment protection is free[cit]
+  //   swap(*this, other);
+  //   return *this;
+  // };
+  friend void swap(Ouput& a, Ouput& b) noexcept {
+    std::swap(a.vals_, b.vals_);
+    std::swap(a.derivs_, b.derivs_);
+    std::swap(a.virial_, b.virial_);
+  }
+  ~Ouput() {}//this does not own anything
+};
+
+class Input final {
+  using vd=std::vector<double>*;
+  using cvd = const std::vector<double>*;
+  using vv=std::vector<Vector>*;
+  using cvv = const std::vector<Vector>*;
+  std::variant<std::monostate,std::vector<double>*,const std::vector<double>*> masses_;
+  std::variant<std::monostate,std::vector<double>*,const std::vector<double>*> charges_;
+  std::variant<std::monostate,std::vector<Vector>*,const std::vector<Vector>*> positions_;
+public:
+//const because the data in the class is not changing: we are modifying data pointed by the const pointers
+//get_if return nullptr on error
+//get throws on error
+//references to variable data, work only if the data passed is a reference to variable data
+  std::vector<double>&var_masses() const {return *std::get<std::vector<double>*>(masses_);}
+  std::vector<double>&var_charges() const {return *std::get<std::vector<double>*>(charges_);}
+  std::vector<Vector>&var_positions() const {return *std::get<std::vector<Vector>*>(positions_);}
+  //references to constant data, must ALWAYS work
+  const std::vector<double>&masses() const {
+    if(auto *p=std::get_if<cvd>(&masses_)) {
+      return **p;
+    }
+    //get_if takes a ptr, get takes a ref...
+    //this throws if there is no data
+    return *std::get<vd>(masses_);
+  }
+  const std::vector<double>&charges() const {
+    if(auto *p=std::get_if<cvd>(&charges_)) {
+      return **p;
+    }
+    //get_if takes a ptr, get takes a ref...
+    //this throws if there is no data
+    return *std::get<vd>(charges_);
+  }
+  const std::vector<Vector>&positions() const {
+    if(auto *p=std::get_if<cvv>(&positions_)) {
+      return **p;
+    }
+    //get_if takes a ptr, get takes a ref...
+    //this throws if there is no data
+    return *std::get<vv>(positions_);
+  }
+
+  Input()=default;
+  template<typename masses_T, typename charges_T, typename positions_T>
+  Input(masses_T& masses,
+        charges_T& charges,
+        positions_T& positions)
+    : masses_(&masses),
+      charges_(&charges),
+      positions_(&positions) {
+    static_assert(std::is_same_v<masses_T,    std::vector<double>> || std::is_same_v<masses_T,    const std::vector<double>>);
+    static_assert(std::is_same_v<charges_T,   std::vector<double>> || std::is_same_v<charges_T,   const std::vector<double>>);
+    static_assert(std::is_same_v<positions_T, std::vector<Vector>> || std::is_same_v<positions_T, const std::vector<Vector>>);
+  }
+  Input(Input const& other)
+    : masses_(other.masses_),
+      charges_(other.charges_),
+      positions_(other.positions_) {};
+  Input(Input && other) noexcept
+    : Input() {
+    swap(*this, other);
+  };
+  //with the builder pattern I ma not need the ubercomplex ctor  and the just things that I set up...
+  Input& charges(std::vector<double>& charges) {
+    charges_=&charges;
+    return *this;
+  }
+  Input& charges(const std::vector<double>& charges) {
+    charges_=&charges;
+    return *this;
+  }
+  Input& masses(std::vector<double>& masses) {
+    masses_=&masses;
+    return *this;
+  }
+  Input& masses(const std::vector<double>& masses) {
+    masses_=&masses;
+    return *this;
+  }
+  Input& positions(std::vector<Vector>& positions) {
+    positions_=&positions;
+    return *this;
+  }
+  Input& positions(const std::vector<Vector>& positions) {
+    positions_=&positions;
+    return *this;
+  }
+
+  // the input is initialized as copy or move ctor
+  Input& operator=(Input other) =delete;
+  friend void swap(Input& a, Input& b) noexcept {
+    std::swap(a.masses_, b.masses_);
+    std::swap(a.charges_, b.charges_);
+    std::swap(a.positions_, b.positions_);
+  }
+  ~Input() {}//this does not own anything
+};
+
+} // namespace multiColvars
+
+/*MANUAL DRAFT
+To setup a CV compatible with multicolvartemplate you need to add
+ - A type that will be used to pass to calculateCV the calculation settings
+   - this type will be called `Modetype`
+   - using  Modetype=type;
+   - some default types are already defined in MultiColvarTemplate.h
+   - for example `using  Modetype=PLME::colvar::multiColvars::emptyMode;` when no setting inmact the calculation
+   - others predefined types are:
+     - `enum class components {withCompontents, noComponents};`
+     - `enum class plainOrScaled {scaled, plain};`
+     - `enum class scaledComponents {withCompontents, scaledComponents, noComponents};`
+ - some static function that will be called by the MultiColvarTemplate:
+   - these functions will need to match the foloowing signatures (withour the constness of the non const& arguments):
+   - `static void parseAtomList( int num, std::vector<AtomNumber>& t, ActionAtomistic* aa );`
+   - `Modetype getModeAndSetupValues( ActionWithValue* av )`
+   - `static void calculateCV( Modetype mode, PLMD::colvar::multiColvars::Input in, PLMD::colvar::multiColvars::Ouput out, const ActionAtomistic* aa )
+     - the input variable contain the references to position, masses, charges
+     - the output variable contain the references to value, derivs, virial,and acts as a return argument
+     - this function will be called by MultiColvarTemplate to calculate the CV value on the inputs
+     - to avoid code repetitition you should change ::calculate() to call this function
+     - By default all the inputs are const ref, but the constedness can be changed, since the MulticolvarTemplate will pass the plain references
+ - By default you can decorate the public part of your CV wiht `MULTICOLVAR_DEFAULT(modetype here)` to "conver" the declaration of the CV in MultiColvarTemplate-compatible one, you will need to manually declare and implement the methods
+   - If you need to change some constness of the `calculateCV` signature use `MULTICOLVAR_SETTINGS_BASE(modetype here)` and then declare `calculateCV` with the desired constness in the signature (see Dipole.cpp)
+*/
+
+
+#define MULTICOLVAR_SETTINGS_BASE(type) using  Modetype=type; \
+  static void parseAtomList( int num, std::vector<AtomNumber>& t, ActionAtomistic* aa ); \
+  static Modetype getModeAndSetupValues( ActionWithValue* av )
+#define MULTICOLVAR_SETTINGS_CALCULATE_CONST() static void calculateCV( Modetype mode, \
+                          PLMD::colvar::multiColvars::Input in, \
+                          PLMD::colvar::multiColvars::Ouput out, \
+                          const ActionAtomistic* aa )
+
+#define MULTICOLVAR_DEFAULT(type) MULTICOLVAR_SETTINGS_BASE(type); \
+          MULTICOLVAR_SETTINGS_CALCULATE_CONST()
+//^no ';' here so that the macro will "consume" it and not generate a double semicolon error
+
+template <class CV>
 class MultiColvarTemplate : public ActionWithVector {
 private:
+  using Modetype =typename CV::Modetype;
 /// An index that decides what we are calculating
-  unsigned mode;
+  Modetype mode;
 /// Are we using pbc to calculate the CVs
   bool usepbc;
 /// Do we reassemble the molecule
@@ -48,9 +247,9 @@ class MultiColvarTemplate : public ActionWithVector {
   void calculate() override;
 };
 
-template <class T>
-void MultiColvarTemplate<T>::registerKeywords(Keywords& keys ) {
-  T::registerKeywords( keys );
+template <class CV>
+void MultiColvarTemplate<CV>::registerKeywords(Keywords& keys ) {
+  CV::registerKeywords( keys );
   keys.add("optional","MASK","the label for a sparse matrix that should be used to determine which elements of the matrix should be computed");
   unsigned nkeys = keys.size();
   for(unsigned i=0; i<nkeys; ++i) {
@@ -59,80 +258,96 @@ void MultiColvarTemplate<T>::registerKeywords(Keywords& keys ) {
   if( keys.outputComponentExists(".#!value") ) keys.setValueDescription("the " + keys.getDisplayName() + " for each set of specified atoms");
 }
 
-template <class T>
-MultiColvarTemplate<T>::MultiColvarTemplate(const ActionOptions&ao):
+template <class CV>
+MultiColvarTemplate<CV>::MultiColvarTemplate(const ActionOptions&ao):
   Action(ao),
   ActionWithVector(ao),
-  mode(0),
+  mode(),
   usepbc(true),
   wholemolecules(false)
 {
   std::vector<AtomNumber> all_atoms;
   if( getName()=="POSITION_VECTOR" || getName()=="MASS_VECTOR" || getName()=="CHARGE_VECTOR" ) parseAtomList( "ATOMS", all_atoms );
   if( all_atoms.size()>0 ) {
-    ablocks.resize(1); ablocks[0].resize( all_atoms.size() );
-    for(unsigned i=0; i<all_atoms.size(); ++i) ablocks[0][i] = i;
+    ablocks.resize(1);
+    ablocks[0].resize( all_atoms.size() );
+    for(unsigned i=0; i<all_atoms.size(); ++i) {
+      ablocks[0][i] = i;
+    }
   } else {
     std::vector<AtomNumber> t;
     for(int i=1;; ++i ) {
-      T::parseAtomList( i, t, this );
+      CV::parseAtomList( i, t, this );
       if( t.empty() ) break;
 
-      if( i==1 ) { ablocks.resize(t.size()); }
+      if( i==1 ) {
+        ablocks.resize(t.size());
+      }
       if( t.size()!=ablocks.size() ) {
-        std::string ss; Tools::convert(i,ss);
+        std::string ss;
+        Tools::convert(i,ss);
         error("ATOMS" + ss + " keyword has the wrong number of atoms");
       }
       for(unsigned j=0; j<ablocks.size(); ++j) {
-        ablocks[j].push_back( ablocks.size()*(i-1)+j ); all_atoms.push_back( t[j] );
+        ablocks[j].push_back( ablocks.size()*(i-1)+j );
+        all_atoms.push_back( t[j] );
       }
       t.resize(0);
     }
   }
-  if( all_atoms.size()==0 ) error("No atoms have been specified");
+  if( all_atoms.size()==0 ) {
+    error("No atoms have been specified");
+  }
   requestAtoms(all_atoms);
   if( keywords.exists("NOPBC") ) {
-    bool nopbc=!usepbc; parseFlag("NOPBC",nopbc);
+    bool nopbc=!usepbc;
+    parseFlag("NOPBC",nopbc);
     usepbc=!nopbc;
   }
   if( keywords.exists("WHOLEMOLECULES") ) {
     parseFlag("WHOLEMOLECULES",wholemolecules);
-    if( wholemolecules ) usepbc=false;
+    if( wholemolecules ) {
+      usepbc=false;
+    }
   }
   if( usepbc ) log.printf("  using periodic boundary conditions\n");
   else    log.printf("  without periodic boundary conditions\n");
 
   // Setup the values
-  mode = T::getModeAndSetupValues( this );
+  mode = CV::getModeAndSetupValues( this );
 }
 
-template <class T>
-unsigned MultiColvarTemplate<T>::getNumberOfDerivatives() {
+template <class CV>
+unsigned MultiColvarTemplate<CV>::getNumberOfDerivatives() {
   return 3*getNumberOfAtoms()+9;
 }
 
-template <class T>
-void MultiColvarTemplate<T>::calculate() {
+template <class CV>
+void MultiColvarTemplate<CV>::calculate() {
   if( wholemolecules ) makeWhole();
   runAllTasks();
 }
 
-template <class T>
-void MultiColvarTemplate<T>::addValueWithDerivatives( const std::vector<unsigned>& shape ) {
+template <class CV>
+void MultiColvarTemplate<CV>::addValueWithDerivatives( const std::vector<unsigned>& shape ) {
   std::vector<unsigned> s(1); s[0]=ablocks[0].size(); addValue( s );
 }
 
-template <class T>
-void MultiColvarTemplate<T>::addComponentWithDerivatives( const std::string& name, const std::vector<unsigned>& shape ) {
+template <class CV>
+void MultiColvarTemplate<CV>::addComponentWithDerivatives( const std::string& name, const std::vector<unsigned>& shape ) {
   std::vector<unsigned> s(1); s[0]=ablocks[0].size(); addComponent( name, s );
 }
 
-template <class T>
-void MultiColvarTemplate<T>::performTask( const unsigned& task_index, MultiValue& myvals ) const {
+template <class CV>
+void MultiColvarTemplate<CV>::performTask( const unsigned& task_index, MultiValue& myvals ) const {
   // Retrieve the positions
   std::vector<Vector> & fpositions( myvals.getFirstAtomVector() );
-  if( fpositions.size()!=ablocks.size() ) fpositions.resize( ablocks.size() );
-  for(unsigned i=0; i<ablocks.size(); ++i) fpositions[i] = getPosition( ablocks[i][task_index] );
+  if( fpositions.size()!=ablocks.size() ) {
+    fpositions.resize( ablocks.size() );
+  }
+  for(unsigned i=0; i<ablocks.size(); ++i) {
+    fpositions[i] = getPosition( ablocks[i][task_index] );
+  }
   // If we are using pbc make whole
   if( usepbc ) {
     if( fpositions.size()==1 ) {
@@ -143,23 +358,36 @@ void MultiColvarTemplate<T>::performTask( const unsigned& task_index, MultiValue
         second=first+pbcDistance(first,second);
       }
     }
-  } else if( fpositions.size()==1 ) fpositions[0]=delta(Vector(0.0,0.0,0.0),getPosition( ablocks[0][task_index] ) );
+  } else if( fpositions.size()==1 ) {
+    fpositions[0]=delta(Vector(0.0,0.0,0.0),getPosition( ablocks[0][task_index] ) );
+  }
   // Retrieve the masses and charges
   myvals.resizeTemporyVector(2);
   std::vector<double> & mass( myvals.getTemporyVector(0) );
   std::vector<double> & charge( myvals.getTemporyVector(1) );
-  if( mass.size()!=ablocks.size() ) { mass.resize(ablocks.size()); charge.resize(ablocks.size()); }
-  for(unsigned i=0; i<ablocks.size(); ++i) { mass[i]=getMass( ablocks[i][task_index] ); charge[i]=getCharge( ablocks[i][task_index] ); }
+  if( mass.size()!=ablocks.size() ) {
+    mass.resize(ablocks.size());
+    charge.resize(ablocks.size());
+  }
+  for(unsigned i=0; i<ablocks.size(); ++i) {
+    mass[i]=getMass( ablocks[i][task_index] );
+    charge[i]=getCharge( ablocks[i][task_index] );
+  }
   // Make some space to store various things
   std::vector<double> values( getNumberOfComponents() );
   std::vector<Tensor> & virial( myvals.getFirstAtomVirialVector() );
   std::vector<std::vector<Vector> > & derivs( myvals.getFirstAtomDerivativeVector() );
-  if( derivs.size()!=values.size() ) { derivs.resize( values.size() ); virial.resize( values.size() ); }
+  if( derivs.size()!=values.size() ) {
+    derivs.resize( values.size() );
+    virial.resize( values.size() );
+  }
   for(unsigned i=0; i<derivs.size(); ++i) {
-    if( derivs[i].size()<ablocks.size() ) derivs[i].resize( ablocks.size() );
+    if( derivs[i].size()<ablocks.size() ) {
+      derivs[i].resize( ablocks.size() );
+    }
   }
   // Calculate the CVs using the method in the Colvar
-  T::calculateCV( mode, mass, charge, fpositions, values, derivs, virial, this );
+  CV::calculateCV( mode, multiColvars::Input(mass, charge, fpositions), multiColvars::Ouput{values, derivs, virial}, this );
   for(unsigned i=0; i<values.size(); ++i) myvals.setValue( i, values[i] );
   // Finish if there are no derivatives
   if( doNotCalculateDerivatives() ) return;
diff --git a/src/colvar/Plane.cpp b/src/colvar/Plane.cpp
index 9d57a1832d..c94a9293b0 100644
--- a/src/colvar/Plane.cpp
+++ b/src/colvar/Plane.cpp
@@ -62,19 +62,15 @@ namespace colvar {
 class Plane : public Colvar {
 private:
   bool pbc;
-  std::vector<double> value, masses, charges;
+  std::vector<double> value;
   std::vector<std::vector<Vector> > derivs;
   std::vector<Tensor> virial;
 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 unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
-                           const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
+  MULTICOLVAR_DEFAULT(multiColvars::emptyMode);
 };
 
 typedef ColvarShortcut<Plane> PlaneShortcut;
@@ -92,7 +88,7 @@ void Plane::registerKeywords( Keywords& keys ) {
   keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
 }
 
-void Plane::parseAtomList( const int& num, std::vector<AtomNumber>& atoms, ActionAtomistic* aa ) {
+void Plane::parseAtomList( int const 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());
@@ -104,11 +100,11 @@ void Plane::parseAtomList( const int& num, std::vector<AtomNumber>& atoms, Actio
   } else if( num<0 || atoms.size()>0 ) aa->error("Number of specified atoms should be either 3 or 4");
 }
 
-unsigned Plane::getModeAndSetupValues( ActionWithValue* av ) {
+Plane::Modetype Plane::getModeAndSetupValues( ActionWithValue* av ) {
   av->addComponentWithDerivatives("x"); av->componentIsNotPeriodic("x");
   av->addComponentWithDerivatives("y"); av->componentIsNotPeriodic("y");
   av->addComponentWithDerivatives("z"); av->componentIsNotPeriodic("z");
-  return 0;
+  return {};
 }
 
 Plane::Plane(const ActionOptions&ao):
@@ -127,7 +123,7 @@ Plane::Plane(const ActionOptions&ao):
   if(pbc) log.printf("  using periodic boundary conditions\n");
   else    log.printf("  without periodic boundary conditions\n");
 
-  unsigned mode = getModeAndSetupValues( this );
+  /*Modetype mode = */getModeAndSetupValues( this );
   requestAtoms(atoms);
   checkRead();
 }
@@ -135,15 +131,19 @@ Plane::Plane(const ActionOptions&ao):
 void Plane::calculate() {
 
   if(pbc) makeWhole();
-  calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
+  calculateCV( {}, multiColvars::Input().positions(getPositions()), multiColvars::Ouput(value, derivs, virial), this );
   setValue( value[0] );
   for(unsigned i=0; i<derivs[0].size(); ++i) setAtomsDerivatives( i, derivs[0][i] );
   setBoxDerivatives( virial[0] );
 }
 
-void Plane::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
-                         const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                         std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
+void Plane::calculateCV( Modetype /*mode*/,
+                         multiColvars::Input const in,
+                         multiColvars::Ouput out, const ActionAtomistic* aa ) {
+  auto & vals=out.vals();
+  auto & derivs=out.derivs();
+  auto & virial=out.virial();
+  const auto & pos = in.positions();
   Vector d1=delta( pos[1], pos[0] );
   Vector d2=delta( pos[2], pos[3] );
   Vector cp = crossProduct( d1, d2 );
diff --git a/src/colvar/Position.cpp b/src/colvar/Position.cpp
index 20bdc8f234..0abf9c016e 100644
--- a/src/colvar/Position.cpp
+++ b/src/colvar/Position.cpp
@@ -91,21 +91,18 @@ Create a vector that holds the components of the position of a set of atoms.
 //+ENDPLUMEDOC
 
 class Position : public Colvar {
-  bool scaled_components;
-  bool pbc;
-  std::vector<double> value, masses, charges;
-  std::vector<std::vector<Vector> > derivs;
-  std::vector<Tensor> virial;
 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 unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
-                           const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
+  MULTICOLVAR_DEFAULT(multiColvars::plainOrScaled);
+private:
+  Modetype components;
+  bool pbc;
+  std::vector<double> value;
+  std::vector<std::vector<Vector> > derivs;
+  std::vector<Tensor> virial;
 };
 
 typedef ColvarShortcut<Position> PositionShortcut;
@@ -131,7 +128,6 @@ void Position::registerKeywords( Keywords& keys ) {
 
 Position::Position(const ActionOptions&ao):
   PLUMED_COLVAR_INIT(ao),
-  scaled_components(false),
   pbc(true),
   value(3),
   derivs(3),
@@ -139,8 +135,7 @@ Position::Position(const ActionOptions&ao):
 {
   for(unsigned i=0; i<3; ++i) derivs[i].resize(1);
   std::vector<AtomNumber> atoms; parseAtomList(-1,atoms,this);
-  unsigned mode=getModeAndSetupValues(this);
-  if( mode==1 ) scaled_components=true;
+  components=getModeAndSetupValues(this);
 
   bool nopbc=!pbc;
   parseFlag("NOPBC",nopbc);
@@ -153,25 +148,26 @@ Position::Position(const ActionOptions&ao):
   requestAtoms(atoms);
 }
 
-void Position::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
+void Position::parseAtomList(int const 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);
+Position::Modetype 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;
+    return Modetype::scaled;
   }
   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;
+  return Modetype::plain;
 }
 
 // calculator
@@ -183,9 +179,9 @@ void Position::calculate() {
   } else {
     distance[0]=delta(Vector(0.0,0.0,0.0),getPosition(0));
   }
-
-  if(scaled_components) {
-    calculateCV( 1, masses, charges, distance, value, derivs, virial, this );
+  calculateCV( components, multiColvars::Input().positions(distance), multiColvars::Ouput(value, derivs, virial), this );
+  switch (components) {
+  case Modetype::scaled: {
     Value* valuea=getPntrToComponent("a");
     Value* valueb=getPntrToComponent("b");
     Value* valuec=getPntrToComponent("c");
@@ -195,8 +191,10 @@ void Position::calculate() {
     valueb->set(value[1]);
     setAtomsDerivatives (valuec,0,derivs[2][0]);
     valuec->set(value[2]);
-  } else {
-    calculateCV( 0, masses, charges, distance, value, derivs, virial, this );
+  }
+  break;
+
+  case Modetype::plain: {
     Value* valuex=getPntrToComponent("x");
     Value* valuey=getPntrToComponent("y");
     Value* valuez=getPntrToComponent("z");
@@ -213,26 +211,37 @@ void Position::calculate() {
     setBoxDerivatives   (valuez,virial[2]);
     valuez->set(value[2]);
   }
+  }
 }
 
-void Position::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
-                            const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                            std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
-  if( mode==1 ) {
+void Position::calculateCV( Modetype mode,
+                            multiColvars::Input const in,
+                            multiColvars::Ouput out, const ActionAtomistic* aa ) {
+  auto & vals=out.vals();
+  auto & derivs=out.derivs();
+  auto & virial=out.virial();
+  const auto &pos = in.positions();
+  ////////////////////////////////////////////////////////////////////////////
+  //why here there is not the distance from 000 treatment??????
+  ///////////////////////////////////////////////////////////////////////////
+  switch (mode)
+  {
+  case Modetype::scaled: {
     Vector d=aa->getPbc().realToScaled(pos[0]);
     vals[0]=Tools::pbc(d[0]); vals[1]=Tools::pbc(d[1]); vals[2]=Tools::pbc(d[2]);
     derivs[0][0]=matmul(aa->getPbc().getInvBox(),Vector(+1,0,0));
     derivs[1][0]=matmul(aa->getPbc().getInvBox(),Vector(0,+1,0));
     derivs[2][0]=matmul(aa->getPbc().getInvBox(),Vector(0,0,+1));
-  } else {
+  }
+  break;
+
+  case Modetype::plain: {
     for(unsigned i=0; i<3; ++i) vals[i]=pos[0][i];
     derivs[0][0]=Vector(+1,0,0); derivs[1][0]=Vector(0,+1,0); derivs[2][0]=Vector(0,0,+1);
     virial[0]=Tensor(pos[0],Vector(-1,0,0)); virial[1]=Tensor(pos[0],Vector(0,-1,0)); virial[2]=Tensor(pos[0],Vector(0,0,-1));
   }
+  }
 }
 
-}
-}
-
-
-
+} // namespace colvar
+} // namespacs PLMD
diff --git a/src/colvar/SelectMassCharge.cpp b/src/colvar/SelectMassCharge.cpp
index 4d5a4855ea..69a80d4fd9 100644
--- a/src/colvar/SelectMassCharge.cpp
+++ b/src/colvar/SelectMassCharge.cpp
@@ -82,29 +82,40 @@ Get the mass of one or multiple atoms
 namespace PLMD {
 namespace colvar {
 
+namespace
+{
+enum class MC {Mass,Charge};
+} // namespace unnamed
+
+template <MC mq>
 class SelectMassCharge : public Colvar {
 public:
   static void registerKeywords( Keywords& keys );
   explicit SelectMassCharge(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 unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
-                           const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
+  MULTICOLVAR_DEFAULT(multiColvars::emptyMode);
+private:
+  AtomNumber theAtom;
 };
 
-typedef ColvarShortcut<SelectMassCharge> MQShortcut;
-PLUMED_REGISTER_ACTION(MQShortcut,"MASS")
-PLUMED_REGISTER_ACTION(MQShortcut,"CHARGE")
-PLUMED_REGISTER_ACTION(SelectMassCharge,"MASS_SCALAR")
-PLUMED_REGISTER_ACTION(SelectMassCharge,"CHARGE_SCALAR")
-typedef MultiColvarTemplate<SelectMassCharge> MQMulti;
-PLUMED_REGISTER_ACTION(MQMulti,"MASS_VECTOR")
-PLUMED_REGISTER_ACTION(MQMulti,"CHARGE_VECTOR")
+typedef SelectMassCharge<MC::Mass> SelectMass;
+typedef SelectMassCharge<MC::Charge> SelectCharge;
+PLUMED_REGISTER_ACTION(SelectMass,"MASS_SCALAR")
+PLUMED_REGISTER_ACTION(SelectCharge,"CHARGE_SCALAR")
+
+typedef ColvarShortcut<SelectMass> MassShortcut;
+typedef ColvarShortcut<SelectCharge> ChargeShortcut;
+PLUMED_REGISTER_ACTION(MassShortcut,"MASS")
+PLUMED_REGISTER_ACTION(ChargeShortcut,"CHARGE")
+
+typedef MultiColvarTemplate<SelectMass> MassMulti;
+typedef MultiColvarTemplate<SelectCharge> ChargeMulti;
+PLUMED_REGISTER_ACTION(MassMulti,"MASS_VECTOR")
+PLUMED_REGISTER_ACTION(ChargeMulti,"CHARGE_VECTOR")
 
-void SelectMassCharge::registerKeywords( Keywords& keys ) {
+template <MC mq>
+void SelectMassCharge<mq>::registerKeywords( Keywords& keys ) {
   Colvar::registerKeywords( keys );
   keys.add("atoms","ATOM","the atom number");
   keys.add("atoms","ATOMS","the atom numbers that you would like to store the masses and charges of");
@@ -114,49 +125,89 @@ void SelectMassCharge::registerKeywords( Keywords& keys ) {
   keys.setDisplayName( acname.substr(0,und) ); keys.setValueDescription("the " + keys.getDisplayName() + " of the atom");
 }
 
-SelectMassCharge::SelectMassCharge(const ActionOptions&ao):
+template <MC mq>
+SelectMassCharge<mq>::SelectMassCharge(const ActionOptions&ao):
   PLUMED_COLVAR_INIT(ao)
 {
-  std::vector<AtomNumber> atoms; parseAtomList(-1,atoms,this);
-  unsigned mode=getModeAndSetupValues(this);
+  std::vector<AtomNumber> atoms;
+  parseAtomList(-1,atoms,this);
+  theAtom=atoms[0];
+  /*Modetype mode=*/getModeAndSetupValues(this);
   requestAtoms(atoms);
 }
 
-void SelectMassCharge::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
+template <MC mq>
+void SelectMassCharge<mq>::parseAtomList(  int const 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");
+  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 SelectMassCharge::getModeAndSetupValues( ActionWithValue* av ) {
-  av->addValueWithDerivatives(); av->setNotPeriodic(); bool constant=true;
-  ActionAtomistic* aa=dynamic_cast<ActionAtomistic*>( av ); plumed_assert( aa );
+template <MC mq>
+typename SelectMassCharge<mq>::Modetype SelectMassCharge<mq>::getModeAndSetupValues( ActionWithValue* av ) {
+  av->addValueWithDerivatives();
+  av->setNotPeriodic();
+  bool constant=true;
+  ActionAtomistic* aa=dynamic_cast<ActionAtomistic*>( av );
+  plumed_assert( aa );
   for(unsigned i=0; i<aa->getNumberOfAtoms(); ++i) {
     std::pair<std::size_t,std::size_t> p = aa->getValueIndices( aa->getAbsoluteIndex(i) );
-    if( av->getName().find("MASS")!=std::string::npos && !aa->masv[p.first]->isConstant() ) constant=false;
-    if( av->getName().find("CHARGE")!=std::string::npos && !aa->chargev[p.first]->isConstant() ) constant=false;
+    if constexpr( mq == MC::Mass ) {
+      constant = aa->isMassConstant(p.first);
+    } else {
+      constant = aa->isChargeConstant(p.first);
+    }
+  }
+  if( !constant ) {
+    av->error("cannot deal with non-constant " + av->getName() + " values");
   }
-  if( !constant ) av->error("cannot deal with non-constant " + av->getName() + " values");
   (av->copyOutput(0))->setConstant();
-  return 0;
+  return {};
 }
 
 // calculator
-void SelectMassCharge::calculate() {
-  std::vector<double> masses(1), charges(1), vals(1);
-  std::vector<Vector> pos; std::vector<std::vector<Vector> > derivs; std::vector<Tensor> virial;
-  calculateCV( 0, masses, charges, pos, vals, derivs, virial, this ); setValue( vals[0] );
-}
 
-void SelectMassCharge::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
-                                    const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                                    std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
-  if( aa->getName().find("MASSES")!=std::string::npos ) vals[0]=masses[0];
-  else if( aa->chargesWereSet ) vals[0]=charges[0];
-}
+template <MC mq>
+void SelectMassCharge<mq>::calculate() {
 
-}
-}
+  std::vector<std::vector<Vector> > derivsdummy;
+  std::vector<Tensor> virialdummy;
+  std::vector<double> vals(1);
 
+  if constexpr( mq == MC::Mass ) {
+    std::vector<double> mass{getMass(theAtom.index())};
+    calculateCV( {}, multiColvars::Input().masses(mass), multiColvars::Ouput(vals, derivsdummy, virialdummy), this );
+  } else {
+    std::vector<double> charge{getCharge(theAtom.index())};
+    calculateCV( {}, multiColvars::Input().charges(charge), multiColvars::Ouput(vals, derivsdummy, virialdummy), this );
+  }
 
+  setValue( vals[0] );
+  // does the code above give the same result as doing:
+  // if constexpr( mq == MC::Mass ) {
+  //   setValue( getMass(theAtom.index) );
+  // } else {
+  //   setValue( getCharge(theAtom.index) );
+  // }
+  // ???
+  //calculateCV copies the first element from masses or charges into vals[0]
+}
+
+template <MC mq>
+void SelectMassCharge<mq>::calculateCV( Modetype /*mode*/,
+                                        multiColvars::Input const in,
+                                        multiColvars::Ouput out, const ActionAtomistic* aa ) {
+  auto & vals=out.vals();
+  if constexpr(mq == MC::Mass)  {
+    vals[0]= in.masses()[0];
+    // } else if( aa->chargesWereSet ) { this is done in the getModeAndSetupValues by isChargeConstant
+  } else {
+    vals[0]=in.charges()[0];
+  }
+}
 
+} // namespace colvar
+} // namespace PLMD
diff --git a/src/colvar/Torsion.cpp b/src/colvar/Torsion.cpp
index 710b1359eb..83205ba674 100644
--- a/src/colvar/Torsion.cpp
+++ b/src/colvar/Torsion.cpp
@@ -107,22 +107,21 @@ Calculate multiple torsional angles.
 //+ENDPLUMEDOC
 
 class Torsion : public Colvar {
-  bool pbc;
-  bool do_cosine;
-
-  std::vector<double> value, masses, charges;
-  std::vector<std::vector<Vector> > derivs;
-  std::vector<Tensor> virial;
 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 unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
-                           const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
   static void registerKeywords(Keywords& keys);
+  enum class torsionModes {
+    torsion,cosine
+  };
+  MULTICOLVAR_DEFAULT(torsionModes);
+private:
+  std::vector<double> value;
+  std::vector<std::vector<Vector> > derivs;
+  std::vector<Tensor> virial;
+  bool pbc;
+  Modetype mode;
 };
 
 typedef ColvarShortcut<Torsion> TorsionShortcut;
@@ -146,11 +145,10 @@ void Torsion::registerKeywords(Keywords& keys) {
 
 Torsion::Torsion(const ActionOptions&ao):
   PLUMED_COLVAR_INIT(ao),
-  pbc(true),
-  do_cosine(false),
   value(1),
   derivs(1),
-  virial(1)
+  virial(1),
+  pbc(true)
 {
   derivs[0].resize(6); std::vector<AtomNumber> atoms;
   std::vector<AtomNumber> v1; ActionAtomistic::parseAtomList("VECTOR1",v1);
@@ -168,8 +166,7 @@ Torsion::Torsion(const ActionOptions&ao):
     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;
+  mode=getModeAndSetupValues(this);
 
   bool nopbc=!pbc;
   parseFlag("NOPBC",nopbc);
@@ -181,7 +178,7 @@ Torsion::Torsion(const ActionOptions&ao):
   requestAtoms(atoms);
 }
 
-void Torsion::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
+void Torsion::parseAtomList( int const num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
   std::vector<AtomNumber> v1,v2,axis;
   aa->parseAtomList("ATOMS",num,t);
   aa->parseAtomList("VECTORA",num,v1);
@@ -213,34 +210,46 @@ void Torsion::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionA
   } 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");
-
+Torsion::Modetype Torsion::getModeAndSetupValues( ActionWithValue* av ) {
+  bool do_cos;
+  av->parseFlag("COSINE",do_cos);
   av->addValueWithDerivatives();
-  if(!do_cos) { av->setPeriodic("-pi","pi"); return 0; }
-  av->setNotPeriodic(); return 1;
+  if(do_cos) {
+    av->log.printf("  calculating cosine instead of torsion\n");
+    av->setNotPeriodic();
+    return Modetype::cosine;
+  }
+  av->setPeriodic("-pi","pi");
+  return Modetype::torsion;
 }
 
 // calculator
 void Torsion::calculate() {
-  if(pbc) makeWhole();
-  if(do_cosine) calculateCV( 1, masses, charges, getPositions(), value, derivs, virial, this );
-  else calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
-  for(unsigned i=0; i<6; ++i) setAtomsDerivatives(i,derivs[0][i] );
-  setValue(value[0]); setBoxDerivatives( virial[0] );
+  if(pbc) {
+    makeWhole();
+  }
+  calculateCV( mode, multiColvars::Input().positions(getPositions()), multiColvars::Ouput(value, derivs, virial), this );
+  for(unsigned i=0; i<6; ++i) {
+    setAtomsDerivatives(i,derivs[0][i] );
+  }
+  setValue(value[0]);
+  setBoxDerivatives( virial[0] );
 }
 
-void Torsion::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
-                           const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
+void Torsion::calculateCV( Modetype mode,
+                           multiColvars::Input const in,
+                           multiColvars::Ouput out, const ActionAtomistic* aa ) {
+  auto & vals=out.vals();
+  auto & derivs=out.derivs();
+  auto & virial=out.virial();
+  const auto & pos = in.positions();
   Vector d0=delta(pos[1],pos[0]);
   Vector d1=delta(pos[3],pos[2]);
   Vector d2=delta(pos[5],pos[4]);
   Vector dd0,dd1,dd2;
   PLMD::Torsion t;
   vals[0] = t.compute(d0,d1,d2,dd0,dd1,dd2);
-  if(mode==1) {
+  if(mode==Modetype::cosine) {
     dd0 *= -std::sin(vals[0]);
     dd1 *= -std::sin(vals[0]);
     dd2 *= -std::sin(vals[0]);
diff --git a/src/core/ActionAtomistic.h b/src/core/ActionAtomistic.h
index da8c1266b7..f04f08452e 100644
--- a/src/core/ActionAtomistic.h
+++ b/src/core/ActionAtomistic.h
@@ -35,10 +35,6 @@ namespace PLMD {
 class Pbc;
 class PDB;
 
-namespace colvar {
-class SelectMassCharge;
-}
-
 /// \ingroup MULTIINHERIT
 /// Action used to create objects that access the positions of the atoms from the MD code
 class ActionAtomistic :
@@ -46,7 +42,6 @@ class ActionAtomistic :
 {
   friend class Group;
   friend class DomainDecomposition;
-  friend class colvar::SelectMassCharge;
   friend class ActionWithVirtualAtom;
 
   std::vector<AtomNumber> indexes;         // the set of needed atoms
@@ -82,9 +77,9 @@ class ActionAtomistic :
 protected:
   bool                  chargesWereSet;
   void setExtraCV(const std::string &name);
+public:
 /// Used to interpret whether this index is a virtual atom or a real atom
   std::pair<std::size_t, std::size_t> getValueIndices( const AtomNumber& i ) const ;
-public:
 /// Request an array of atoms.
 /// This method is used to ask for a list of atoms. Atoms
 /// should be asked for by number. If this routine is called
@@ -116,8 +111,12 @@ class ActionAtomistic :
   const double & getEnergy()const;
 /// Get mass of i-th atom
   double getMass(int i)const;
+/// Check if the mass of the i-th atom is constant
+  bool isMassConstant(int i) const;
 /// Get charge of i-th atom
   double getCharge(int i)const;
+/// Check if the charge of the i-th atom is constant
+  bool isChargeConstant(int i) const;
 /// Get the force acting on a particular atom
   Vector getForce( const std::pair<std::size_t, std::size_t>& a ) const ;
 /// Add force to an atom
@@ -205,12 +204,23 @@ double ActionAtomistic::getMass(int i)const {
   return masses[i];
 }
 
+inline
+bool ActionAtomistic::isMassConstant(int i)const {
+  return masv[i]->isConstant();
+}
+
 inline
 double ActionAtomistic::getCharge(int i) const {
   if( !chargesWereSet ) error("charges were not passed to plumed");
   return charges[i];
 }
 
+inline
+bool ActionAtomistic::isChargeConstant(int i) const {
+  if( !chargesWereSet ) error("charges were not passed to plumed");
+  return chargev[i]->isConstant();
+}
+
 inline
 const std::vector<AtomNumber> & ActionAtomistic::getAbsoluteIndexes()const {
   return indexes;
diff --git a/src/core/ActionWithVector.cpp b/src/core/ActionWithVector.cpp
index 8eed7e7771..2bc7f7517a 100644
--- a/src/core/ActionWithVector.cpp
+++ b/src/core/ActionWithVector.cpp
@@ -277,12 +277,16 @@ void ActionWithVector::getNumberOfTasks( unsigned& ntasks ) {
 
 void ActionWithVector::runTask( const unsigned& current, MultiValue& myvals ) const {
   const ActionWithMatrix* am = dynamic_cast<const ActionWithMatrix*>(this);
-  myvals.setTaskIndex(current); myvals.vector_call=true; performTask( current, myvals );
+  myvals.setTaskIndex(current);
+  myvals.vector_call=true;
+  performTask( current, myvals );
   for(unsigned i=0; i<getNumberOfComponents(); ++i) {
     const Value* myval = getConstPntrToComponent(i);
-    if( am || myval->hasDerivatives() ) continue;
+    if( am || myval->hasDerivatives() )
+      continue;
     Value* myv = const_cast<Value*>( myval );
-    if( getName()=="RMSD_VECTOR" && myv->getRank()==2 ) continue;
+    if( getName()=="RMSD_VECTOR" && myv->getRank()==2 )
+      continue;
     myv->set( current, myvals.get( i ) );
   }
 }
diff --git a/src/core/Value.cpp b/src/core/Value.cpp
index 18fe337d15..68b71ac860 100644
--- a/src/core/Value.cpp
+++ b/src/core/Value.cpp
@@ -361,7 +361,7 @@ void Value::print( OFile& ofile ) const {
 
 void Value::printForce( OFile& ofile ) const {
   if( shape.size()==0 || getNumberOfValues()==1 ) {
-    ofile.printField( name, getForce(0) ); 
+    ofile.printField( name, getForce(0) );
   } else {
     std::vector<unsigned> indices( shape.size() );
     for(unsigned i=0; i<getNumberOfValues(); ++i) {
@@ -371,7 +371,7 @@ void Value::printForce( OFile& ofile ) const {
       ofile.printField( fname,  getForce(i) );
     }
   }
-} 
+}
 
 unsigned Value::getGoodNumThreads( const unsigned& j, const unsigned& k ) const {
   return OpenMP::getGoodNumThreads( &data[j], (k-j) );
diff --git a/src/crystdistrib/Quaternion.cpp b/src/crystdistrib/Quaternion.cpp
index 700b04d016..02975cecba 100644
--- a/src/crystdistrib/Quaternion.cpp
+++ b/src/crystdistrib/Quaternion.cpp
@@ -97,19 +97,15 @@ See \ref QUATERNION for more details
 class Quaternion : public Colvar {
 private:
   bool pbc;
-  std::vector<double> value, masses, charges;
+  std::vector<double> value;
   std::vector<std::vector<Vector> > derivs;
   std::vector<Tensor> virial;
 public:
   static void registerKeywords( Keywords& keys );
   explicit Quaternion(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 unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
-                           const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
+  MULTICOLVAR_DEFAULT(::PLMD::colvar::multiColvars::emptyMode);
 };
 
 typedef colvar::ColvarShortcut<Quaternion> QuaternionShortcut;
@@ -145,28 +141,28 @@ Quaternion::Quaternion(const ActionOptions&ao):
   parseFlag("NOPBC",nopbc);
   pbc=!nopbc;
 
-  unsigned mode = getModeAndSetupValues( this );
+  /*Modetype mode =*/ getModeAndSetupValues( this );
   requestAtoms(atoms);
 }
 
-void Quaternion::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
+void Quaternion::parseAtomList(int num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
   aa->parseAtomList("ATOMS",num,t);
   if( t.size()==3 ) aa->log.printf("  involving atoms %d %d %d\n",t[0].serial(),t[1].serial(),t[0].serial());
 }
 
-unsigned Quaternion::getModeAndSetupValues( ActionWithValue* av ) {
+Quaternion::Modetype Quaternion::getModeAndSetupValues( ActionWithValue* av ) {
   // This sets up values that we can pass around in PLUMED
   av->addComponentWithDerivatives("w"); av->componentIsNotPeriodic("w");
   av->addComponentWithDerivatives("i"); av->componentIsNotPeriodic("i");
   av->addComponentWithDerivatives("j"); av->componentIsNotPeriodic("j");
   av->addComponentWithDerivatives("k"); av->componentIsNotPeriodic("k");
-  return 0;
+  return {};
 }
 
 void Quaternion::calculate() {
   if(pbc) makeWhole();
 
-  calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
+  calculateCV( {}, PLMD::colvar::multiColvars::Input().positions(getPositions()), PLMD::colvar::multiColvars::Ouput(value, derivs, virial), this );
   for(unsigned j=0; j<4; ++j) {
     Value* valuej=getPntrToComponent(j);
     for(unsigned i=0; i<3; ++i) setAtomsDerivatives(valuej,i,derivs[j][i] );
@@ -176,9 +172,13 @@ void Quaternion::calculate() {
 }
 
 // calculator
-void Quaternion::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
-                              const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                              std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
+void Quaternion::calculateCV( Modetype /*mode*/,
+                              PLMD::colvar::multiColvars::Input const in,
+                              PLMD::colvar::multiColvars::Ouput out, const ActionAtomistic* aa ) {
+  auto & vals=out.vals();
+  auto & derivs=out.derivs();
+  auto & virial=out.virial();
+  const auto & pos = in.positions();
   //declarations
   Vector vec1_comp = delta( pos[0], pos[1] ); //components between atom 1 and 2
   Vector vec2_comp = delta( pos[0], pos[2] ); //components between atom 1 and 3