diff --git a/src/colvar/Angle.cpp b/src/colvar/Angle.cpp
index 282d881404..ec109ac1f9 100644
--- a/src/colvar/Angle.cpp
+++ b/src/colvar/Angle.cpp
@@ -110,8 +110,7 @@ class Angle : public Colvar {
   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, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
+  static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs, std::vector<Tensor>& virial );
 };
 
 typedef ColvarShortcut<Angle> AngleShortcut;
@@ -168,14 +167,13 @@ Angle::Angle(const ActionOptions&ao):
 void Angle::calculate() {
 
   if(pbc) makeWhole();
-  calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial, this );
+  calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial );
   setValue( value[0] );
   for(unsigned i=0; i<derivs[0].size(); ++i) setAtomsDerivatives( i, derivs[0][i] );
   setBoxDerivatives( virial[0] );
 }
 
-void Angle::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                         std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
+void Angle::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs, std::vector<Tensor>& virial ) {
   Vector dij,dik;
   dij=delta(cvin.pos[2],cvin.pos[3]);
   dik=delta(cvin.pos[1],cvin.pos[0]);
diff --git a/src/colvar/DihedralCorrelation.cpp b/src/colvar/DihedralCorrelation.cpp
index 58d2783ab8..19b3b09d01 100644
--- a/src/colvar/DihedralCorrelation.cpp
+++ b/src/colvar/DihedralCorrelation.cpp
@@ -73,8 +73,7 @@ class DihedralCorrelation : public Colvar {
   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, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
+  static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs, std::vector<Tensor>& virial );
 };
 
 typedef ColvarShortcut<DihedralCorrelation> DihedralCorrelationShortcut;
@@ -125,14 +124,13 @@ unsigned DihedralCorrelation::getModeAndSetupValues( ActionWithValue* av ) {
 void DihedralCorrelation::calculate() {
 
   if(pbc) makeWhole();
-  calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial, this );
+  calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial );
   setValue( value[0] );
   for(unsigned i=0; i<derivs[0].size(); ++i) setAtomsDerivatives( i, derivs[0][i] );
   setBoxDerivatives( virial[0] );
 }
 
-void DihedralCorrelation::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                                       std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
+void DihedralCorrelation::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs, std::vector<Tensor>& virial ) {
   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]);
diff --git a/src/colvar/Dipole.cpp b/src/colvar/Dipole.cpp
index 02a61ef588..e1ce420351 100644
--- a/src/colvar/Dipole.cpp
+++ b/src/colvar/Dipole.cpp
@@ -94,8 +94,7 @@ class Dipole : public Colvar {
   static unsigned getModeAndSetupValues( ActionWithValue* av );
   void calculate() override;
   static void registerKeywords(Keywords& keys);
-  static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
+  static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs, std::vector<Tensor>& virial );
 };
 
 typedef ColvarShortcut<Dipole> DipoleShortcut;
@@ -122,7 +121,7 @@ Dipole::Dipole(const ActionOptions&ao):
   derivs(1),
   virial(1)
 {
-  parseAtomList(-1,ga_lista,this); 
+  parseAtomList(-1,ga_lista,this);
   components=(getModeAndSetupValues(this)==1);
   if( components ) {
     value.resize(3); derivs.resize(3); virial.resize(3);
@@ -171,12 +170,12 @@ void Dipole::calculate()
   unsigned N=getNumberOfAtoms();
 
   if(!components) {
-    calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial, this );
+    calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial );
     for(unsigned i=0; i<N; i++) setAtomsDerivatives(i,derivs[0][i]);
     setBoxDerivatives(virial[0]);
     setValue(value[0]);
   } else {
-    calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), value, derivs, virial, this );
+    calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), value, derivs, virial );
     for(unsigned i=0; i<N; i++) {
       setAtomsDerivatives(valuex,i,derivs[0][i]);
       setAtomsDerivatives(valuey,i,derivs[1][i]);
@@ -191,8 +190,7 @@ void Dipole::calculate()
   }
 }
 
-void Dipole::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                          std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
+void Dipole::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs, std::vector<Tensor>& virial ) {
   unsigned N=cvin.pos.size(); double ctot=0.;
   for(unsigned i=0; i<N; ++i) ctot += cvin.charges[i];
   ctot/=(double)N;
diff --git a/src/colvar/Distance.cpp b/src/colvar/Distance.cpp
index 3e41f3bba0..bc9298a04d 100644
--- a/src/colvar/Distance.cpp
+++ b/src/colvar/Distance.cpp
@@ -138,8 +138,7 @@ class Distance : public Colvar {
   static unsigned getModeAndSetupValues( ActionWithValue* av );
 // active methods:
   void calculate() override;
-  static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
+  static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs, std::vector<Tensor>& virial );
 };
 
 typedef ColvarShortcut<Distance> DistanceShortcut;
@@ -226,7 +225,7 @@ void Distance::calculate() {
   if(pbc) makeWhole();
 
   if( components ) {
-    calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), value, derivs, virial, this );
+    calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), value, derivs, virial );
     Value* valuex=getPntrToComponent("x");
     Value* valuey=getPntrToComponent("y");
     Value* valuez=getPntrToComponent("z");
@@ -243,7 +242,7 @@ void Distance::calculate() {
     setBoxDerivatives(valuez,virial[2]);
     valuez->set(value[2]);
   } else if( scaled_components ) {
-    calculateCV( ColvarInput::createColvarInput( 2, getPositions(), this ), value, derivs, virial, this );
+    calculateCV( ColvarInput::createColvarInput( 2, getPositions(), this ), value, derivs, virial );
 
     Value* valuea=getPntrToComponent("a");
     Value* valueb=getPntrToComponent("b");
@@ -255,15 +254,14 @@ void Distance::calculate() {
     for(unsigned i=0; i<2; ++i) setAtomsDerivatives(valuec,i,derivs[2][i] );
     valuec->set(value[2]);
   } else  {
-    calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial, this );
+    calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial );
     for(unsigned i=0; i<2; ++i) setAtomsDerivatives(i,derivs[0][i] );
     setBoxDerivatives(virial[0]);
     setValue           (value[0]);
   }
 }
 
-void Distance::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                            std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
+void Distance::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs, std::vector<Tensor>& virial ) {
   Vector distance=delta(cvin.pos[0],cvin.pos[1]);
   const double value=distance.modulo();
   const double invvalue=1.0/value;
@@ -282,15 +280,15 @@ void Distance::calculateCV( const ColvarInput& cvin, std::vector<double>& vals,
     vals[2] = distance[2];
     setBoxDerivativesNoPbc( cvin.pos, derivs, virial );
   } else if(cvin.mode==2) {
-    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));
+    Vector d=cvin.pbc.realToScaled(distance);
+    derivs[0][0] = matmul(cvin.pbc.getInvBox(),Vector(-1,0,0));
+    derivs[0][1] = matmul(cvin.pbc.getInvBox(),Vector(+1,0,0));
     vals[0] = Tools::pbc(d[0]);
-    derivs[1][0] = matmul(aa->getPbc().getInvBox(),Vector(0,-1,0));
-    derivs[1][1] = matmul(aa->getPbc().getInvBox(),Vector(0,+1,0));
+    derivs[1][0] = matmul(cvin.pbc.getInvBox(),Vector(0,-1,0));
+    derivs[1][1] = matmul(cvin.pbc.getInvBox(),Vector(0,+1,0));
     vals[1] = Tools::pbc(d[1]);
-    derivs[2][0] = matmul(aa->getPbc().getInvBox(),Vector(0,0,-1));
-    derivs[2][1] = matmul(aa->getPbc().getInvBox(),Vector(0,0,+1));
+    derivs[2][0] = matmul(cvin.pbc.getInvBox(),Vector(0,0,-1));
+    derivs[2][1] = matmul(cvin.pbc.getInvBox(),Vector(0,0,+1));
     vals[2] = Tools::pbc(d[2]);
   } else {
     derivs[0][0] = -invvalue*distance;
diff --git a/src/colvar/MultiColvarTemplate.cpp b/src/colvar/MultiColvarTemplate.cpp
index 4602c337c4..30d6fcd625 100644
--- a/src/colvar/MultiColvarTemplate.cpp
+++ b/src/colvar/MultiColvarTemplate.cpp
@@ -25,16 +25,17 @@
 namespace PLMD {
 namespace colvar {
 
-ColvarInput::ColvarInput( const unsigned& m, const std::vector<Vector>& p, const std::vector<double>& w, const std::vector<double>& q ) :
-mode(m),
-pos(p),
-mass(w),
-charges(q)
+ColvarInput::ColvarInput( const unsigned& m, const std::vector<Vector>& p, const std::vector<double>& w, const std::vector<double>& q, const Pbc& box ) :
+  mode(m),
+  pbc(box),
+  pos(p),
+  mass(w),
+  charges(q)
 {
 }
 
 ColvarInput ColvarInput::createColvarInput( const unsigned& m, const std::vector<Vector>& p, const Colvar* colv ) {
-  return ColvarInput( m, p, colv->getMasses(), colv->getCharges(true) );
+  return ColvarInput( m, p, colv->getMasses(), colv->getCharges(true), colv->getPbc() );
 }
 
 }
diff --git a/src/colvar/MultiColvarTemplate.h b/src/colvar/MultiColvarTemplate.h
index ceea73e8d3..9087975543 100644
--- a/src/colvar/MultiColvarTemplate.h
+++ b/src/colvar/MultiColvarTemplate.h
@@ -34,10 +34,11 @@ namespace colvar {
 class ColvarInput {
 public:
   unsigned mode;
+  const Pbc& pbc;
   const std::vector<Vector>& pos;
   const std::vector<double>& mass;
   const std::vector<double>& charges;
-  ColvarInput( const unsigned& m, const std::vector<Vector>& p, const std::vector<double>& w, const std::vector<double>& q );
+  ColvarInput( const unsigned& m, const std::vector<Vector>& p, const std::vector<double>& w, const std::vector<double>& q, const Pbc& box );
   static ColvarInput createColvarInput( const unsigned& m, const std::vector<Vector>& p, const Colvar* colv );
 };
 
@@ -177,7 +178,7 @@ void MultiColvarTemplate<T>::performTask( const unsigned& task_index, MultiValue
   std::vector<Tensor> & virial( myvals.getFirstAtomVirialVector() );
   std::vector<std::vector<Vector> > & derivs( myvals.getFirstAtomDerivativeVector() );
   // Calculate the CVs using the method in the Colvar
-  T::calculateCV( ColvarInput(mode, fpositions, mass, charge), values, derivs, virial, this );
+  T::calculateCV( ColvarInput(mode, fpositions, mass, charge, getPbc() ), values, derivs, virial );
   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 8cd6164697..0a9968cd62 100644
--- a/src/colvar/Plane.cpp
+++ b/src/colvar/Plane.cpp
@@ -72,8 +72,7 @@ class Plane : public Colvar {
   static unsigned getModeAndSetupValues( ActionWithValue* av );
 // active methods:
   void calculate() override;
-  static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
+  static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs, std::vector<Tensor>& virial );
 };
 
 typedef ColvarShortcut<Plane> PlaneShortcut;
@@ -134,14 +133,13 @@ Plane::Plane(const ActionOptions&ao):
 void Plane::calculate() {
 
   if(pbc) makeWhole();
-  calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial, this );
+  calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial );
   setValue( value[0] );
   for(unsigned i=0; i<derivs[0].size(); ++i) setAtomsDerivatives( i, derivs[0][i] );
   setBoxDerivatives( virial[0] );
 }
 
-void Plane::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                         std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
+void Plane::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs, std::vector<Tensor>& virial ) {
   Vector d1=delta( cvin.pos[1], cvin.pos[0] );
   Vector d2=delta( cvin.pos[2], cvin.pos[3] );
   Vector cp = crossProduct( d1, d2 );
diff --git a/src/colvar/Position.cpp b/src/colvar/Position.cpp
index 02bfdcbb6f..4288429c58 100644
--- a/src/colvar/Position.cpp
+++ b/src/colvar/Position.cpp
@@ -103,8 +103,7 @@ class Position : public Colvar {
   static unsigned getModeAndSetupValues( ActionWithValue* av );
 // active methods:
   void calculate() override;
-  static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
+  static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs, std::vector<Tensor>& virial );
 };
 
 typedef ColvarShortcut<Position> PositionShortcut;
@@ -184,7 +183,7 @@ void Position::calculate() {
   }
 
   if(scaled_components) {
-    calculateCV( ColvarInput::createColvarInput( 1, distance, this ), value, derivs, virial, this );
+    calculateCV( ColvarInput::createColvarInput( 1, distance, this ), value, derivs, virial );
     Value* valuea=getPntrToComponent("a");
     Value* valueb=getPntrToComponent("b");
     Value* valuec=getPntrToComponent("c");
@@ -195,7 +194,7 @@ void Position::calculate() {
     setAtomsDerivatives (valuec,0,derivs[2][0]);
     valuec->set(value[2]);
   } else {
-    calculateCV( ColvarInput::createColvarInput( 0, distance, this ), value, derivs, virial, this );
+    calculateCV( ColvarInput::createColvarInput( 0, distance, this ), value, derivs, virial );
     Value* valuex=getPntrToComponent("x");
     Value* valuey=getPntrToComponent("y");
     Value* valuez=getPntrToComponent("z");
@@ -214,14 +213,13 @@ void Position::calculate() {
   }
 }
 
-void Position::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                            std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
+void Position::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs, std::vector<Tensor>& virial ) {
   if( cvin.mode==1 ) {
-    Vector d=aa->getPbc().realToScaled(cvin.pos[0]);
+    Vector d=cvin.pbc.realToScaled(cvin.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));
+    derivs[0][0]=matmul(cvin.pbc.getInvBox(),Vector(+1,0,0));
+    derivs[1][0]=matmul(cvin.pbc.getInvBox(),Vector(0,+1,0));
+    derivs[2][0]=matmul(cvin.pbc.getInvBox(),Vector(0,0,+1));
   } else {
     for(unsigned i=0; i<3; ++i) vals[i]=cvin.pos[0][i];
     derivs[0][0]=Vector(+1,0,0); derivs[1][0]=Vector(0,+1,0); derivs[2][0]=Vector(0,0,+1);
diff --git a/src/colvar/SelectMassCharge.cpp b/src/colvar/SelectMassCharge.cpp
index 50c8079e0a..1e81128d9e 100644
--- a/src/colvar/SelectMassCharge.cpp
+++ b/src/colvar/SelectMassCharge.cpp
@@ -80,9 +80,9 @@ Get the mass of one or multiple atoms
 namespace PLMD {
 namespace colvar {
 
-class SelectMassCharge : 
-public ActionAtomistic,
-public ActionWithValue {
+class SelectMassCharge :
+  public ActionAtomistic,
+  public ActionWithValue {
 public:
   static void registerKeywords( Keywords& keys );
   explicit SelectMassCharge(const ActionOptions&);
@@ -112,13 +112,13 @@ SelectMassCharge::SelectMassCharge(const ActionOptions&ao):
   std::vector<AtomNumber> atoms; parseAtomList("ATOMS",atoms);
   log.printf("  getting %s of atoms : ", getName().c_str() );
   for(unsigned i=0; i<atoms.size(); ++i) {
-      std::pair<std::size_t,std::size_t> p = getValueIndices( atoms[i] ); 
-      if( getName()=="MASS" && !masv[p.first]->isConstant() ) error("cannot deal with non-constant " + getName() + " values");  
-      if( getName()=="CHARGE" && !chargev[p.first]->isConstant() ) error("cannot deal with non-constant " + getName() + " values");
-      log.printf("%d ", atoms[i].serial() );
+    std::pair<std::size_t,std::size_t> p = getValueIndices( atoms[i] );
+    if( getName()=="MASS" && !masv[p.first]->isConstant() ) error("cannot deal with non-constant " + getName() + " values");
+    if( getName()=="CHARGE" && !chargev[p.first]->isConstant() ) error("cannot deal with non-constant " + getName() + " values");
+    log.printf("%d ", atoms[i].serial() );
   }
   log.printf("\n"); requestAtoms(atoms);
-  std::vector<unsigned> shape(1); 
+  std::vector<unsigned> shape(1);
   if(atoms.size()==1) shape.resize(0);
   else shape[0] = atoms.size();
   addValue( shape ); setNotPeriodic();
@@ -127,13 +127,13 @@ SelectMassCharge::SelectMassCharge(const ActionOptions&ao):
 
 // calculator
 void SelectMassCharge::calculate() {
-   Value* myval = getPntrToComponent(0);
-   if( getName()=="CHARGES" ) {
-       if( !chargesWereSet ) error("cannot determine charges are charges were not set");
-       for(unsigned i=0; i<getNumberOfAtoms(); ++i) myval->set( i, getCharge(i) );
-   } else {
-       for(unsigned i=0; i<getNumberOfAtoms(); ++i) myval->set( i, getMass(i) );
-   }
+  Value* myval = getPntrToComponent(0);
+  if( getName()=="CHARGES" ) {
+    if( !chargesWereSet ) error("cannot determine charges are charges were not set");
+    for(unsigned i=0; i<getNumberOfAtoms(); ++i) myval->set( i, getCharge(i) );
+  } else {
+    for(unsigned i=0; i<getNumberOfAtoms(); ++i) myval->set( i, getMass(i) );
+  }
 }
 
 }
diff --git a/src/colvar/Torsion.cpp b/src/colvar/Torsion.cpp
index 3345925626..2ec5124929 100644
--- a/src/colvar/Torsion.cpp
+++ b/src/colvar/Torsion.cpp
@@ -119,8 +119,7 @@ class Torsion : public Colvar {
   static unsigned getModeAndSetupValues( ActionWithValue* av );
 // active methods:
   void calculate() override;
-  static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
+  static void calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs, std::vector<Tensor>& virial );
   static void registerKeywords(Keywords& keys);
 };
 
@@ -224,14 +223,13 @@ unsigned Torsion::getModeAndSetupValues( ActionWithValue* av ) {
 // calculator
 void Torsion::calculate() {
   if(pbc) makeWhole();
-  if(do_cosine) calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), value, derivs, virial, this );
-  else calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial, this );
+  if(do_cosine) calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), value, derivs, virial );
+  else calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial );
   for(unsigned i=0; i<6; ++i) setAtomsDerivatives(i,derivs[0][i] );
   setValue(value[0]); setBoxDerivatives( virial[0] );
 }
 
-void Torsion::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
+void Torsion::calculateCV( const ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs, std::vector<Tensor>& virial ) {
   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]);
diff --git a/src/crystdistrib/Quaternion.cpp b/src/crystdistrib/Quaternion.cpp
index 826f5c7423..290cb015c9 100644
--- a/src/crystdistrib/Quaternion.cpp
+++ b/src/crystdistrib/Quaternion.cpp
@@ -107,8 +107,7 @@ class Quaternion : public Colvar {
   static unsigned getModeAndSetupValues( ActionWithValue* av );
 // active methods:
   void calculate() override;
-  static void calculateCV( const colvar::ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                           std::vector<Tensor>& virial, const ActionAtomistic* aa );
+  static void calculateCV( const colvar::ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs, std::vector<Tensor>& virial );
 };
 
 typedef colvar::ColvarShortcut<Quaternion> QuaternionShortcut;
@@ -165,7 +164,7 @@ unsigned Quaternion::getModeAndSetupValues( ActionWithValue* av ) {
 void Quaternion::calculate() {
   if(pbc) makeWhole();
 
-  calculateCV( colvar::ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial, this );
+  calculateCV( colvar::ColvarInput::createColvarInput( 0, getPositions(), this ), value, derivs, virial );
   for(unsigned j=0; j<4; ++j) {
     Value* valuej=getPntrToComponent(j);
     for(unsigned i=0; i<3; ++i) setAtomsDerivatives(valuej,i,derivs[j][i] );
@@ -175,8 +174,7 @@ void Quaternion::calculate() {
 }
 
 // calculator
-void Quaternion::calculateCV( const colvar::ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
-                              std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
+void Quaternion::calculateCV( const colvar::ColvarInput& cvin, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs, std::vector<Tensor>& virial ) {
   //declarations
   Vector vec1_comp = delta( cvin.pos[0], cvin.pos[1] ); //components between atom 1 and 2
   Vector vec2_comp = delta( cvin.pos[0], cvin.pos[2] ); //components between atom 1 and 3