From e61f5f16daf22806460b1813aa851327daef63b9 Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Fri, 17 Mar 2017 19:36:58 +0000
Subject: [PATCH] Removed some small and unhelpful classes

---
 src/adjmat/ActionWithInputMatrix.cpp          |   4 +-
 src/adjmat/ActionWithInputMatrix.h            |   3 +-
 src/crystallization/BondOrientation.cpp       |  11 +-
 .../InterMolecularTorsions.cpp                |  11 +-
 src/crystallization/MoleculeOrientation.cpp   |  12 +-
 src/crystallization/MoleculePlane.cpp         |   6 +-
 src/crystallization/OrientationSphere.cpp     |  12 +-
 src/crystallization/OrientationSphere.h       |   5 +-
 src/crystallization/VectorMultiColvar.cpp     |   8 +-
 src/crystallization/VectorMultiColvar.h       |   9 +-
 src/multicolvar/AlphaBeta.cpp                 |  22 ++--
 src/multicolvar/Angles.cpp                    |  21 +++-
 src/multicolvar/Bridge.cpp                    |  11 +-
 .../BridgedMultiColvarFunction.cpp            |   4 +-
 src/multicolvar/BridgedMultiColvarFunction.h  |   2 +-
 src/multicolvar/CenterOfMultiColvar.cpp       |   2 +-
 src/multicolvar/Density.cpp                   |  10 +-
 src/multicolvar/DihedralCorrelation.cpp       |  22 ++--
 src/multicolvar/Distances.cpp                 |  21 +++-
 src/multicolvar/InPlaneDistances.cpp          |  10 +-
 src/multicolvar/LocalAverage.cpp              |  27 ++--
 src/multicolvar/MultiColvar.cpp               |  89 --------------
 src/multicolvar/MultiColvar.h                 |  56 ---------
 src/multicolvar/MultiColvarBase.cpp           | 115 ++++++++++++++----
 src/multicolvar/MultiColvarBase.h             |  12 +-
 src/multicolvar/MultiColvarCombine.cpp        |  15 ++-
 src/multicolvar/MultiColvarFunction.cpp       |  76 ------------
 src/multicolvar/MultiColvarFunction.h         |  58 ---------
 src/multicolvar/NumberOfLinks.cpp             |  11 +-
 src/multicolvar/Torsions.cpp                  |  21 +++-
 src/multicolvar/VolumeGradientBase.cpp        |   6 +-
 src/multicolvar/XAngle.cpp                    |  21 +++-
 src/multicolvar/XDistances.cpp                |  21 +++-
 src/multicolvar/XYDistances.cpp               |  24 ++--
 src/multicolvar/XYTorsion.cpp                 |  21 +++-
 35 files changed, 330 insertions(+), 449 deletions(-)
 delete mode 100644 src/multicolvar/MultiColvar.cpp
 delete mode 100644 src/multicolvar/MultiColvar.h
 delete mode 100644 src/multicolvar/MultiColvarFunction.cpp
 delete mode 100644 src/multicolvar/MultiColvarFunction.h

diff --git a/src/adjmat/ActionWithInputMatrix.cpp b/src/adjmat/ActionWithInputMatrix.cpp
index 553c2ae25..b0f6551d0 100644
--- a/src/adjmat/ActionWithInputMatrix.cpp
+++ b/src/adjmat/ActionWithInputMatrix.cpp
@@ -30,14 +30,14 @@ namespace PLMD {
 namespace adjmat {
 
 void ActionWithInputMatrix::registerKeywords( Keywords& keys ){
-  MultiColvarFunction::registerKeywords( keys ); keys.remove("DATA");
+  MultiColvarBase::registerKeywords( keys ); 
   keys.add("compulsory","MATRIX","the action that calcualtes the adjacency matrix vessel we would like to analyse"); 
 }
 
 
 ActionWithInputMatrix::ActionWithInputMatrix(const ActionOptions& ao):
 Action(ao),
-MultiColvarFunction(ao),
+MultiColvarBase(ao),
 mymatrix(NULL)
 {
   matsums=true; 
diff --git a/src/adjmat/ActionWithInputMatrix.h b/src/adjmat/ActionWithInputMatrix.h
index bbc400ec1..e08124bd2 100644
--- a/src/adjmat/ActionWithInputMatrix.h
+++ b/src/adjmat/ActionWithInputMatrix.h
@@ -25,14 +25,13 @@
 #include "core/ActionWithValue.h"
 #include "core/ActionAtomistic.h"
 #include "multicolvar/MultiColvarBase.h"
-#include "multicolvar/MultiColvarFunction.h"
 
 namespace PLMD {
 namespace adjmat {
 
 class AdjacencyMatrixVessel;
 
-class ActionWithInputMatrix : public multicolvar::MultiColvarFunction {
+class ActionWithInputMatrix : public multicolvar::MultiColvarBase {
 protected:
 /// The vessel that holds the adjacency matrix
   AdjacencyMatrixVessel* mymatrix;
diff --git a/src/crystallization/BondOrientation.cpp b/src/crystallization/BondOrientation.cpp
index 2076bbfd2..9369c5117 100644
--- a/src/crystallization/BondOrientation.cpp
+++ b/src/crystallization/BondOrientation.cpp
@@ -21,6 +21,7 @@
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
 #include "tools/SwitchingFunction.h"
 #include "core/ActionRegister.h"
+#include "multicolvar/AtomValuePack.h"
 #include "VectorMultiColvar.h"
 
 //+PLUMEDOC MCOLVAR BOND_DIRECTIONS
@@ -50,7 +51,12 @@ PLUMED_REGISTER_ACTION(BondOrientation,"BOND_DIRECTIONS")
 
 void BondOrientation::registerKeywords( Keywords& keys ){
   VectorMultiColvar::registerKeywords( keys );
-  keys.use("ATOMS");
+  keys.add("numbered","ATOMS","the atoms involved in each of the vectors you wish to calculate. "
+                               "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one vector will be "
+                               "calculated for each ATOM keyword you specify (all ATOM keywords should "
+                               "specify the indices of two atoms).  The eventual number of quantities calculated by this "
+                               "action will depend on what functions of the distribution you choose to calculate.");
+  keys.reset_style("ATOMS","atoms");
   keys.add("atoms-1","GROUP","Calculate the distance between each distinct pair of atoms in the group");
   keys.add("atoms-2","GROUPA","Calculate the distances between all the atoms in GROUPA and all "
                               "the atoms in GROUPB. This must be used in conjuction with GROUPB.");
@@ -74,7 +80,8 @@ VectorMultiColvar(ao)
   weightHasDerivatives=true;
   std::vector<AtomNumber> all_atoms; 
   readTwoGroups( "GROUP", "GROUPA", "GROUPB", all_atoms );
-  int natoms=2; readAtoms( natoms, all_atoms ); 
+  if( atom_lab.size()==0 ) readAtomsLikeKeyword( "ATOMS", 2, all_atoms ); 
+  setupMultiColvarBase( all_atoms );
   // Read in the switching function
   std::string sw, errors; parse("SWITCH",sw);
   if(sw.length()>0){
diff --git a/src/crystallization/InterMolecularTorsions.cpp b/src/crystallization/InterMolecularTorsions.cpp
index e2432e133..14df0dd23 100644
--- a/src/crystallization/InterMolecularTorsions.cpp
+++ b/src/crystallization/InterMolecularTorsions.cpp
@@ -19,7 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "multicolvar/MultiColvarFunction.h"
+#include "multicolvar/MultiColvarBase.h"
+#include "multicolvar/AtomValuePack.h"
 #include "core/ActionRegister.h"
 #include "tools/SwitchingFunction.h"
 #include "tools/Torsion.h" 
@@ -50,7 +51,7 @@ the torsional angle between an orientation vector for molecule \f$i\f$ and molec
 namespace PLMD {
 namespace crystallization {
 
-class InterMolecularTorsions : public multicolvar::MultiColvarFunction {
+class InterMolecularTorsions : public multicolvar::MultiColvarBase {
 private:
 /// The switching function that tells us if atoms are close enough together
   SwitchingFunction switchingFunction;
@@ -69,7 +70,7 @@ public:
 PLUMED_REGISTER_ACTION(InterMolecularTorsions,"INTERMOLECULARTORSIONS")
 
 void InterMolecularTorsions::registerKeywords( Keywords& keys ){
-  MultiColvarFunction::registerKeywords( keys );
+  MultiColvarBase::registerKeywords( keys );
   keys.add("atoms","MOLS","The molecules you would like to calculate the torsional angles between. This should be the label/s of \\ref MOLECULES or \\ref PLANES actions");
   keys.add("atoms-1","MOLSA","In this version of the input the torsional angles between all pairs of atoms including one atom from MOLA one atom from MOLB will be computed. " 
                              "This should be the label/s of \\ref MOLECULES or \\ref PLANES actions");
@@ -83,13 +84,13 @@ void InterMolecularTorsions::registerKeywords( Keywords& keys ){
                                "The following provides information on the \\ref switchingfunction that are available. "
                                "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
   // Use actionWithDistributionKeywords
-  keys.remove("LOWMEM"); keys.remove("DATA"); 
+  keys.remove("LOWMEM"); 
   keys.addFlag("LOWMEM",false,"lower the memory requirements");
 }
 
 InterMolecularTorsions::InterMolecularTorsions(const ActionOptions& ao):
 Action(ao),
-MultiColvarFunction(ao)
+MultiColvarBase(ao)
 {
   for(unsigned i=0;i<getNumberOfBaseMultiColvars();++i){
       if( getBaseMultiColvar(i)->getNumberOfQuantities()!=5 ) error("input multicolvar does not calculate molecular orientations");
diff --git a/src/crystallization/MoleculeOrientation.cpp b/src/crystallization/MoleculeOrientation.cpp
index ab0a04c6d..5102013eb 100644
--- a/src/crystallization/MoleculeOrientation.cpp
+++ b/src/crystallization/MoleculeOrientation.cpp
@@ -75,16 +75,16 @@ MoleculeOrientation::MoleculeOrientation( const ActionOptions& ao ):
 Action(ao),
 VectorMultiColvar(ao)
 {
-  int natoms=-1; std::vector<AtomNumber> all_atoms;
-  readAtomsLikeKeyword("MOL",natoms,all_atoms); 
-  nvectors = std::floor( natoms / 2 );
-  if( natoms%2!=0 && 2*nvectors+1!=natoms ) error("number of atoms in molecule specification is wrong.  Should be two or three.");
+  std::vector<AtomNumber> all_atoms;
+  readAtomsLikeKeyword("MOL",-1,all_atoms); 
+  nvectors = std::floor( ablocks.size() / 2 );
+  if( ablocks.size()%2!=0 && 2*nvectors+1!=ablocks.size() ) error("number of atoms in molecule specification is wrong.  Should be two or three.");
 
   if( all_atoms.size()==0 ) error("No atoms were specified");
   setVectorDimensionality( 3*nvectors ); setupMultiColvarBase( all_atoms );
 
-  if( natoms==2*nvectors+1  ){
-    std::vector<bool> catom_ind(natoms, false); catom_ind[natoms-1]=true;
+  if( ablocks.size()==2*nvectors+1  ){
+    std::vector<bool> catom_ind(ablocks.size(), false); catom_ind[ablocks.size()-1]=true;
     setAtomsForCentralAtom( catom_ind );
   } 
 }
diff --git a/src/crystallization/MoleculePlane.cpp b/src/crystallization/MoleculePlane.cpp
index 618b4edb4..24fcdf540 100644
--- a/src/crystallization/MoleculePlane.cpp
+++ b/src/crystallization/MoleculePlane.cpp
@@ -61,9 +61,9 @@ MoleculePlane::MoleculePlane( const ActionOptions& ao ):
 Action(ao),
 VectorMultiColvar(ao)
 {
-  int natoms=-1; std::vector<AtomNumber> all_atoms;
-  readAtomsLikeKeyword("MOL",natoms,all_atoms); 
-  if( natoms!=3 && natoms!=4 ) error("number of atoms in molecule specification is wrong.  Should be three or four.");
+  std::vector<AtomNumber> all_atoms;
+  readAtomsLikeKeyword("MOL",-1,all_atoms); 
+  if( ablocks.size()!=3 && ablocks.size()!=4 ) error("number of atoms in molecule specification is wrong.  Should be three or four.");
 
   if( all_atoms.size()==0 ) error("No atoms were specified");
   setVectorDimensionality( 3 ); setupMultiColvarBase( all_atoms );
diff --git a/src/crystallization/OrientationSphere.cpp b/src/crystallization/OrientationSphere.cpp
index cd568037a..94cb9dac2 100644
--- a/src/crystallization/OrientationSphere.cpp
+++ b/src/crystallization/OrientationSphere.cpp
@@ -28,7 +28,7 @@ namespace PLMD{
 namespace crystallization {
 
 void OrientationSphere::registerKeywords( Keywords& keys ){
-  multicolvar::MultiColvarFunction::registerKeywords( keys );
+  multicolvar::MultiColvarBase::registerKeywords( keys );
   keys.add("compulsory","NN","6","The n parameter of the switching function ");
   keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
   keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
@@ -40,12 +40,12 @@ void OrientationSphere::registerKeywords( Keywords& keys ){
   keys.use("SPECIES"); keys.use("SPECIESA"); keys.use("SPECIESB"); 
   keys.use("MEAN"); keys.use("MORE_THAN"); keys.use("LESS_THAN"); 
   keys.use("MIN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
-  keys.use("LOWEST"); keys.use("HIGHEST"); keys.remove("DATA");
+  keys.use("LOWEST"); keys.use("HIGHEST"); 
 }
 
 OrientationSphere::OrientationSphere(const ActionOptions&ao):
 Action(ao),
-MultiColvarFunction(ao)
+MultiColvarBase(ao)
 {
   if( getNumberOfBaseMultiColvars()>1 ) warning("not sure if orientation sphere works with more than one base multicolvar - check numerical derivatives");
   // Read in the switching function
@@ -73,9 +73,7 @@ double OrientationSphere::compute( const unsigned& tindex, multicolvar::AtomValu
    std::vector<double> this_der( ncomponents ), catom_der( ncomponents ); 
 
    getInputData( 0, true, myatoms, catom_orient ); 
-   multicolvar::CatomPack atom0; 
    MultiValue& myder0=getInputDerivatives( 0, true, myatoms ); 
-   if( !doNotCalculateDerivatives() ) atom0=getCentralAtomPackFromInput( myatoms.getIndex(0) );
 
    for(unsigned i=1;i<myatoms.getNumberOfAtoms();++i){
       Vector& distance=myatoms.getPosition(i);  
@@ -94,12 +92,12 @@ double OrientationSphere::compute( const unsigned& tindex, multicolvar::AtomValu
              MultiValue& myder1=getInputDerivatives( i, true, myatoms );
              mergeInputDerivatives( 1, 2, this_orient.size(), 0, catom_der, myder0, myatoms );
              mergeInputDerivatives( 1, 2, catom_der.size(), i, this_der, myder1, myatoms );
-             myatoms.addComDerivatives( 1, f_dot*(-dfunc)*distance - sw*ddistance, atom0 );
+             addAtomDerivatives( 1, 0, f_dot*(-dfunc)*distance - sw*ddistance, myatoms );
              addAtomDerivatives( 1, i, f_dot*(dfunc)*distance + sw*ddistance, myatoms );
              myatoms.addBoxDerivatives( 1, (-dfunc)*f_dot*Tensor(distance,distance) - sw*extProduct(distance,ddistance) );
              myder1.clearAll();
               
-             myatoms.addComDerivatives( -1, (-dfunc)*distance, atom0 );
+             addAtomDerivatives( -1, 0, (-dfunc)*distance, myatoms );
              addAtomDerivatives( -1, i, (dfunc)*distance, myatoms );
              myatoms.addTemporyBoxDerivatives( (-dfunc)*Tensor(distance,distance) );
 
diff --git a/src/crystallization/OrientationSphere.h b/src/crystallization/OrientationSphere.h
index fd480417c..1b8ce7fae 100644
--- a/src/crystallization/OrientationSphere.h
+++ b/src/crystallization/OrientationSphere.h
@@ -22,13 +22,14 @@
 #ifndef __PLUMED_crystallization_OrientationSphere_h
 #define __PLUMED_crystallization_OrientationSphere_h
 
-#include "multicolvar/MultiColvarFunction.h"
+#include "multicolvar/MultiColvarBase.h"
+#include "multicolvar/AtomValuePack.h"
 #include "tools/SwitchingFunction.h"
 
 namespace PLMD {
 namespace crystallization {
 
-class OrientationSphere : public multicolvar::MultiColvarFunction {
+class OrientationSphere : public multicolvar::MultiColvarBase {
 private:
   double rcut2;
   SwitchingFunction switchingFunction;
diff --git a/src/crystallization/VectorMultiColvar.cpp b/src/crystallization/VectorMultiColvar.cpp
index a7b6596ea..8ca7f03cc 100644
--- a/src/crystallization/VectorMultiColvar.cpp
+++ b/src/crystallization/VectorMultiColvar.cpp
@@ -20,18 +20,18 @@
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
 #include "VectorMultiColvar.h"
-#include "multicolvar/MultiColvarFunction.h"
 #include "multicolvar/BridgedMultiColvarFunction.h"
 
 namespace PLMD {
 namespace crystallization {
 
 void VectorMultiColvar::registerKeywords( Keywords& keys ){
-  MultiColvar::registerKeywords( keys );
+  MultiColvarBase::registerKeywords( keys );
 }
 
 VectorMultiColvar::VectorMultiColvar(const ActionOptions& ao):
-PLUMED_MULTICOLVAR_INIT(ao),
+Action(ao),
+MultiColvarBase(ao),
 store_director(true)
 {
   setLowMemOption(true);
@@ -45,7 +45,7 @@ void VectorMultiColvar::doNotCalculateDirector(){
   store_director=false;    // Need a sanity check in here  so that you don't use the same instance of Q4 to calcualte vectors and directors 
 }
 
-double VectorMultiColvar::doCalculation( const unsigned& taskIndex, multicolvar::AtomValuePack& myatoms ) const {
+double VectorMultiColvar::compute( const unsigned& taskIndex, multicolvar::AtomValuePack& myatoms ) const {
   // Now calculate the vector
   calculateVector( myatoms );
   // Sort out the active derivatives
diff --git a/src/crystallization/VectorMultiColvar.h b/src/crystallization/VectorMultiColvar.h
index f0ad3d552..98a391b4e 100644
--- a/src/crystallization/VectorMultiColvar.h
+++ b/src/crystallization/VectorMultiColvar.h
@@ -23,12 +23,13 @@
 #define __PLUMED_crystallization_VectorMultiColvar_h
 
 #include "tools/Matrix.h"
-#include "multicolvar/MultiColvar.h"
+#include "multicolvar/MultiColvarBase.h"
+#include "multicolvar/AtomValuePack.h"
 
 namespace PLMD {
 namespace crystallization {
 
-class VectorMultiColvar : public multicolvar::MultiColvar {
+class VectorMultiColvar : public multicolvar::MultiColvarBase {
 friend class OrientationSphere;
 friend class VolumeGradientBase;
 private:
@@ -50,9 +51,9 @@ public:
 /// The norm of a vector is not periodic
   virtual bool isPeriodic(){ return false; }
 /// Calculate the multicolvar
-  double doCalculation( const unsigned& taskIndex, multicolvar::AtomValuePack& myatoms ) const ;
+//  double doCalculation( const unsigned& taskIndex, multicolvar::AtomValuePack& myatoms ) const ;
 /// This shouldn't do anything
-  double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const { plumed_error(); }
+  double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const ;
 /// Calculate the vector
   virtual void calculateVector( multicolvar::AtomValuePack& myatoms ) const=0;
 /// Get the number of components in the vector
diff --git a/src/multicolvar/AlphaBeta.cpp b/src/multicolvar/AlphaBeta.cpp
index 71a957319..f68fe750e 100644
--- a/src/multicolvar/AlphaBeta.cpp
+++ b/src/multicolvar/AlphaBeta.cpp
@@ -19,7 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvar.h"
+#include "MultiColvarBase.h"
+#include "AtomValuePack.h"
 #include "tools/Torsion.h"
 #include "core/ActionRegister.h"
 
@@ -92,7 +93,7 @@ Similarly \@psi-4 tells plumed that you want to calculate the \f$\psi\f$ angle o
 */
 //+ENDPLUMEDOC
 
-class AlphaBeta : public MultiColvar {
+class AlphaBeta : public MultiColvarBase {
 private:
   std::vector<double> target;
 public:
@@ -105,19 +106,26 @@ public:
 PLUMED_REGISTER_ACTION(AlphaBeta,"ALPHABETA")
 
 void AlphaBeta::registerKeywords( Keywords& keys ){
-  MultiColvar::registerKeywords( keys );
-  keys.use("ATOMS");
+  MultiColvarBase::registerKeywords( keys );
+  keys.add("numbered","ATOMS","the atoms involved in each of the alpha-beta variables you wish to calculate. "
+                               "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one alpha-beta values will be "
+                               "calculated for each ATOM keyword you specify (all ATOM keywords should "
+                               "specify the indices of four atoms).  The eventual number of quantities calculated by this "
+                               "action will depend on what functions of the distribution you choose to calculate.");
+  keys.reset_style("ATOMS","atoms");
   keys.add("numbered","REFERENCE","the reference values for each of the torsional angles.  If you use a single REFERENCE value the " 
                                   "same reference value is used for all torsions");
   keys.reset_style("REFERENCE","compulsory");
 }
 
 AlphaBeta::AlphaBeta(const ActionOptions&ao):
-PLUMED_MULTICOLVAR_INIT(ao)
+Action(ao),
+MultiColvarBase(ao)
 {
   // Read in the atoms
-  int natoms=4; std::vector<AtomNumber> all_atoms;
-  readAtoms( natoms, all_atoms );
+  std::vector<AtomNumber> all_atoms;
+  readAtomsLikeKeyword( "ATOMS", 4, all_atoms );
+  setupMultiColvarBase( all_atoms );
   // Resize target
   target.resize( getFullNumberOfTasks() );
   // Setup central atom indices
diff --git a/src/multicolvar/Angles.cpp b/src/multicolvar/Angles.cpp
index f630c6352..bf4a81907 100644
--- a/src/multicolvar/Angles.cpp
+++ b/src/multicolvar/Angles.cpp
@@ -19,7 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvar.h"
+#include "MultiColvarBase.h"
+#include "AtomValuePack.h"
 #include "tools/Angle.h"
 #include "tools/SwitchingFunction.h"
 #include "core/ActionRegister.h"
@@ -83,7 +84,7 @@ PRINT ARG=a1.* FILE=colvar
 */
 //+ENDPLUMEDOC
 
-class Angles : public MultiColvar {
+class Angles : public MultiColvarBase {
 private:
   bool use_sf;
   double rcut2_1, rcut2_2;
@@ -102,10 +103,16 @@ public:
 PLUMED_REGISTER_ACTION(Angles,"ANGLES")
 
 void Angles::registerKeywords( Keywords& keys ){
-  MultiColvar::registerKeywords( keys );
-  keys.use("ATOMS"); keys.use("MEAN"); keys.use("LESS_THAN"); 
+  MultiColvarBase::registerKeywords( keys );
+  keys.use("MEAN"); keys.use("LESS_THAN"); 
   keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MORE_THAN");
   // Could also add Region here in theory
+  keys.add("numbered","ATOMS","the atoms involved in each of the angles you wish to calculate. "
+                               "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one angle will be "
+                               "calculated for each ATOM keyword you specify (all ATOM keywords should "
+                               "provide the indices of three atoms).  The eventual number of quantities calculated by this "
+                               "action will depend on what functions of the distribution you choose to calculate.");
+  keys.reset_style("ATOMS","atoms");
   keys.add("atoms-1","GROUP","Calculate angles for each distinct set of three atoms in the group");
   keys.add("atoms-2","GROUPA","A group of central atoms about which angles should be calculated");
   keys.add("atoms-2","GROUPB","When used in conjuction with GROUPA this keyword instructs plumed "
@@ -125,7 +132,8 @@ void Angles::registerKeywords( Keywords& keys ){
 }
 
 Angles::Angles(const ActionOptions&ao):
-PLUMED_MULTICOLVAR_INIT(ao),
+Action(ao),
+MultiColvarBase(ao),
 use_sf(false)
 {
   std::string sfinput,errors; parse("SWITCH",sfinput);
@@ -155,7 +163,8 @@ use_sf(false)
   // Read in the atoms
   std::vector<AtomNumber> all_atoms;
   readGroupKeywords( "GROUP", "GROUPA", "GROUPB", "GROUPC", false, true, all_atoms );
-  int natoms=3; readAtoms( natoms, all_atoms );
+  if( atom_lab.size()==0 ) readAtomsLikeKeyword( "ATOMS", 3, all_atoms );
+  setupMultiColvarBase( all_atoms );
   // Set cutoff for link cells 
   if( use_sf ){ 
     setLinkCellCutoff( sf1.get_dmax() ); 
diff --git a/src/multicolvar/Bridge.cpp b/src/multicolvar/Bridge.cpp
index c53c38594..4d091f063 100644
--- a/src/multicolvar/Bridge.cpp
+++ b/src/multicolvar/Bridge.cpp
@@ -19,7 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvar.h"
+#include "MultiColvarBase.h"
+#include "AtomValuePack.h"
 #include "tools/SwitchingFunction.h"
 #include "core/ActionRegister.h"
 
@@ -58,7 +59,7 @@ PRINT ARG=a1.mean FILE=colvar
 */
 //+ENDPLUMEDOC
 
-class Bridge : public MultiColvar {
+class Bridge : public MultiColvarBase {
 private:
   Vector dij, dik;
   SwitchingFunction sf1;
@@ -74,8 +75,7 @@ public:
 PLUMED_REGISTER_ACTION(Bridge,"BRIDGE")
 
 void Bridge::registerKeywords( Keywords& keys ){
-  MultiColvar::registerKeywords( keys );
-  keys.use("ATOMS"); 
+  MultiColvarBase::registerKeywords( keys );
   keys.add("atoms-2","BRIDGING_ATOMS","The list of atoms that can form the bridge between the two interesting parts "
                                       "of the structure.");
   keys.add("atoms-2","GROUPA","The list of atoms that are in the first interesting part of the structure"); 
@@ -88,7 +88,8 @@ void Bridge::registerKeywords( Keywords& keys ){
 }
 
 Bridge::Bridge(const ActionOptions&ao):
-PLUMED_MULTICOLVAR_INIT(ao)
+Action(ao),
+MultiColvarBase(ao)
 {
   // Read in the atoms
   std::vector<AtomNumber> all_atoms;
diff --git a/src/multicolvar/BridgedMultiColvarFunction.cpp b/src/multicolvar/BridgedMultiColvarFunction.cpp
index 374ffec99..170eee30d 100644
--- a/src/multicolvar/BridgedMultiColvarFunction.cpp
+++ b/src/multicolvar/BridgedMultiColvarFunction.cpp
@@ -115,8 +115,8 @@ void BridgedMultiColvarFunction::deactivate_task( const unsigned& taskno ){
   plumed_merror("This should never be called");
 }
 
-CatomPack BridgedMultiColvarFunction::getCentralAtomPack( const unsigned& basn, const unsigned& curr ){
-  return mycolv->getCentralAtomPack( basn, curr );
+void BridgedMultiColvarFunction::getCentralAtomPack( const unsigned& basn, const unsigned& curr, CatomPack& mypack ){
+  return mycolv->getCentralAtomPack( basn, curr, mypack );
 }
 
 }
diff --git a/src/multicolvar/BridgedMultiColvarFunction.h b/src/multicolvar/BridgedMultiColvarFunction.h
index 75751605a..74a02b764 100644
--- a/src/multicolvar/BridgedMultiColvarFunction.h
+++ b/src/multicolvar/BridgedMultiColvarFunction.h
@@ -87,7 +87,7 @@ public:
   void getIndexList( const unsigned& ntotal, const unsigned& jstore, const unsigned& maxder, std::vector<unsigned>& indices );
   void applyBridgeForces( const std::vector<double>& bb );
   Vector getCentralAtomPos( const unsigned& curr );
-  CatomPack getCentralAtomPack( const unsigned& basn, const unsigned& curr );
+  void getCentralAtomPack( const unsigned& basn, const unsigned& curr, CatomPack& mypack );
 };
 
 inline
diff --git a/src/multicolvar/CenterOfMultiColvar.cpp b/src/multicolvar/CenterOfMultiColvar.cpp
index aa3878a08..5dc543385 100644
--- a/src/multicolvar/CenterOfMultiColvar.cpp
+++ b/src/multicolvar/CenterOfMultiColvar.cpp
@@ -170,7 +170,7 @@ void CenterOfMultiColvar::calculate(){
           }
       }
       // Get the central atom pack 
-      CatomPack mypack( mycolv->getCentralAtomPack( 0, mycolv->getPositionInFullTaskList(i) ) );
+      CatomPack mypack; mycolv->getCentralAtomPack( 0, mycolv->getPositionInFullTaskList(i), mypack );
       for(unsigned j=0;j<mypack.getNumberOfAtomsWithDerivatives();++j){
           unsigned jder=3*mypack.getIndex(j);
           // Derivatives of sine
diff --git a/src/multicolvar/Density.cpp b/src/multicolvar/Density.cpp
index af192ec6e..4041c651e 100644
--- a/src/multicolvar/Density.cpp
+++ b/src/multicolvar/Density.cpp
@@ -19,7 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvar.h"
+#include "MultiColvarBase.h"
+#include "AtomValuePack.h"
 #include "core/ActionRegister.h"
 
 #include <string>
@@ -49,7 +50,7 @@ PRINT ARG=d1.* FILE=colvar1 FMT=%8.4f
 //+ENDPLUMEDOC
 
 
-class Density : public MultiColvar {
+class Density : public MultiColvarBase {
 public:
   static void registerKeywords( Keywords& keys );
   explicit Density(const ActionOptions&);
@@ -69,12 +70,13 @@ public:
 PLUMED_REGISTER_ACTION(Density,"DENSITY")
 
 void Density::registerKeywords( Keywords& keys ){
-  MultiColvar::registerKeywords( keys );
+  MultiColvarBase::registerKeywords( keys );
   keys.use("SPECIES"); 
 }
 
 Density::Density(const ActionOptions&ao):
-PLUMED_MULTICOLVAR_INIT(ao)
+Action(ao),
+MultiColvarBase(ao)
 {
   std::vector<AtomNumber> all_atoms; parseMultiColvarAtomList("SPECIES", -1, all_atoms);
   ablocks.resize(1); ablocks[0].resize( atom_lab.size() ); 
diff --git a/src/multicolvar/DihedralCorrelation.cpp b/src/multicolvar/DihedralCorrelation.cpp
index 955836635..d8781c8f3 100644
--- a/src/multicolvar/DihedralCorrelation.cpp
+++ b/src/multicolvar/DihedralCorrelation.cpp
@@ -19,7 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvar.h"
+#include "MultiColvarBase.h"
+#include "AtomValuePack.h"
 #include "tools/Torsion.h"
 #include "core/ActionRegister.h"
 
@@ -80,7 +81,7 @@ Similarly \@psi-4 tells plumed that you want to calculate the \f$\psi\f$ angle o
 */
 //+ENDPLUMEDOC
 
-class DihedralCorrelation : public MultiColvar {
+class DihedralCorrelation : public MultiColvarBase {
 private:
 public:
   static void registerKeywords( Keywords& keys );
@@ -92,16 +93,23 @@ public:
 PLUMED_REGISTER_ACTION(DihedralCorrelation,"DIHCOR")
 
 void DihedralCorrelation::registerKeywords( Keywords& keys ){
-  MultiColvar::registerKeywords( keys );
-  keys.use("ATOMS");
+  MultiColvarBase::registerKeywords( keys );
+  keys.add("numbered","ATOMS","the atoms involved in each of the dihedral correlation values you wish to calculate. "
+                               "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one dihedral correlation will be "
+                               "calculated for each ATOM keyword you specify (all ATOM keywords should "
+                               "specify the indices of 8 atoms).  The eventual number of quantities calculated by this "
+                               "action will depend on what functions of the distribution you choose to calculate.");
+  keys.reset_style("ATOMS","atoms");
 }
 
 DihedralCorrelation::DihedralCorrelation(const ActionOptions&ao):
-PLUMED_MULTICOLVAR_INIT(ao)
+Action(ao),
+MultiColvarBase(ao)
 {
   // Read in the atoms
-  int natoms=8; std::vector<AtomNumber> all_atoms;
-  readAtoms( natoms, all_atoms );
+  std::vector<AtomNumber> all_atoms;
+  readAtomsLikeKeyword( "ATOMS", 8, all_atoms );
+  setupMultiColvarBase( all_atoms );
   // Stuff for central atoms
   std::vector<bool> catom_ind(8, false); 
   catom_ind[1]=catom_ind[2]=catom_ind[5]=catom_ind[6]=true;
diff --git a/src/multicolvar/Distances.cpp b/src/multicolvar/Distances.cpp
index 3f76a8771..ac0c9a3bd 100644
--- a/src/multicolvar/Distances.cpp
+++ b/src/multicolvar/Distances.cpp
@@ -19,7 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvar.h"
+#include "MultiColvarBase.h"
+#include "AtomValuePack.h"
 #include "core/ActionRegister.h"
 #include "vesselbase/LessThan.h"
 #include "vesselbase/Between.h"
@@ -122,7 +123,7 @@ calculated at every step.
 //+ENDPLUMEDOC
 
 
-class Distances : public MultiColvar {
+class Distances : public MultiColvarBase {
 private:
 public:
   static void registerKeywords( Keywords& keys );
@@ -136,10 +137,16 @@ public:
 PLUMED_REGISTER_ACTION(Distances,"DISTANCES")
 
 void Distances::registerKeywords( Keywords& keys ){
-  MultiColvar::registerKeywords( keys );
-  keys.use("ATOMS"); keys.use("ALT_MIN"); keys.use("LOWEST"); keys.use("HIGHEST"); 
+  MultiColvarBase::registerKeywords( keys );
+  keys.use("ALT_MIN"); keys.use("LOWEST"); keys.use("HIGHEST"); 
   keys.use("MEAN"); keys.use("MIN"); keys.use("MAX"); keys.use("LESS_THAN"); // keys.use("DHENERGY");
   keys.use("MORE_THAN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
+  keys.add("numbered","ATOMS","the atoms involved in each of the distances you wish to calculate. "
+                               "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one distance will be "
+                               "calculated for each ATOM keyword you specify (all ATOM keywords should "
+                               "specify the indices of two atoms).  The eventual number of quantities calculated by this "
+                               "action will depend on what functions of the distribution you choose to calculate.");
+  keys.reset_style("ATOMS","atoms");
   keys.add("atoms-1","GROUP","Calculate the distance between each distinct pair of atoms in the group");
   keys.add("atoms-2","GROUPA","Calculate the distances between all the atoms in GROUPA and all "
                               "the atoms in GROUPB. This must be used in conjuction with GROUPB.");
@@ -148,12 +155,14 @@ void Distances::registerKeywords( Keywords& keys ){
 }
 
 Distances::Distances(const ActionOptions&ao):
-PLUMED_MULTICOLVAR_INIT(ao)
+Action(ao),
+MultiColvarBase(ao)
 {
   // Read in the atoms
   std::vector<AtomNumber> all_atoms; 
   readTwoGroups( "GROUP", "GROUPA", "GROUPB", all_atoms );
-  int natoms=2; readAtoms( natoms, all_atoms );
+  if( atom_lab.size()==0 ) readAtomsLikeKeyword( "ATOMS", 2, all_atoms );
+  setupMultiColvarBase( all_atoms );
   // And check everything has been read in correctly
   checkRead();
 
diff --git a/src/multicolvar/InPlaneDistances.cpp b/src/multicolvar/InPlaneDistances.cpp
index fc9f07a63..9c3941bfb 100644
--- a/src/multicolvar/InPlaneDistances.cpp
+++ b/src/multicolvar/InPlaneDistances.cpp
@@ -19,7 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvar.h"
+#include "MultiColvarBase.h"
+#include "AtomValuePack.h"
 #include "tools/SwitchingFunction.h"
 #include "core/ActionRegister.h"
 #include "vesselbase/LessThan.h"
@@ -43,7 +44,7 @@ Calculate distances in the plane perpendicular to an axis
 */
 //+ENDPLUMEDOC
 
-class InPlaneDistances : public MultiColvar {
+class InPlaneDistances : public MultiColvarBase {
 public:
   static void registerKeywords( Keywords& keys );
   explicit InPlaneDistances(const ActionOptions&);
@@ -55,7 +56,7 @@ public:
 PLUMED_REGISTER_ACTION(InPlaneDistances,"INPLANEDISTANCES")
 
 void InPlaneDistances::registerKeywords( Keywords& keys ){
-  MultiColvar::registerKeywords( keys );
+  MultiColvarBase::registerKeywords( keys );
   keys.use("ALT_MIN"); keys.use("LOWEST"); keys.use("HIGHEST");
   keys.use("MEAN"); keys.use("MIN"); keys.use("MAX"); keys.use("LESS_THAN"); 
   keys.use("MORE_THAN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
@@ -65,7 +66,8 @@ void InPlaneDistances::registerKeywords( Keywords& keys ){
 }
 
 InPlaneDistances::InPlaneDistances(const ActionOptions&ao):
-PLUMED_MULTICOLVAR_INIT(ao)
+Action(ao),
+MultiColvarBase(ao)
 {
   // Read in the atoms
   std::vector<AtomNumber> all_atoms;
diff --git a/src/multicolvar/LocalAverage.cpp b/src/multicolvar/LocalAverage.cpp
index 4eeb6d491..ac7a40f20 100644
--- a/src/multicolvar/LocalAverage.cpp
+++ b/src/multicolvar/LocalAverage.cpp
@@ -19,7 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvarFunction.h"
+#include "MultiColvarBase.h"
+#include "AtomValuePack.h"
 #include "core/ActionRegister.h"
 #include "tools/SwitchingFunction.h"
 
@@ -80,7 +81,7 @@ PRINT ARG=la.* FILE=colvar
 namespace PLMD {
 namespace multicolvar {
 
-class LocalAverage : public MultiColvarFunction {
+class LocalAverage : public MultiColvarBase {
 private:
 /// Cutoff
   double rcut2;
@@ -100,7 +101,7 @@ public:
 PLUMED_REGISTER_ACTION(LocalAverage,"LOCAL_AVERAGE")
 
 void LocalAverage::registerKeywords( Keywords& keys ){
-  MultiColvarFunction::registerKeywords( keys );
+  MultiColvarBase::registerKeywords( keys );
   keys.add("compulsory","NN","6","The n parameter of the switching function ");
   keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
   keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
@@ -111,7 +112,7 @@ void LocalAverage::registerKeywords( Keywords& keys ){
   // Use actionWithDistributionKeywords
   keys.use("SPECIES"); keys.use("SPECIESA"); keys.use("SPECIESB");
   keys.remove("LOWMEM"); keys.use("MEAN"); keys.use("MORE_THAN"); keys.use("LESS_THAN");
-  keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS"); keys.remove("DATA");
+  keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS"); 
   keys.addFlag("LOWMEM",false,"lower the memory requirements");
   if( keys.reserved("VMEAN") ) keys.use("VMEAN");
   if( keys.reserved("VSUM") ) keys.use("VSUM");
@@ -119,7 +120,7 @@ void LocalAverage::registerKeywords( Keywords& keys ){
 
 LocalAverage::LocalAverage(const ActionOptions& ao):
 Action(ao),
-MultiColvarFunction(ao)
+MultiColvarBase(ao)
 {
   if( getNumberOfBaseMultiColvars()>1 ) error("local average with more than one base colvar makes no sense");
   // Read in the switching function
@@ -144,7 +145,7 @@ unsigned LocalAverage::getNumberOfQuantities() const {
 }
 
 double LocalAverage::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
-  double d2, sw, dfunc; CatomPack atom0, atom1; MultiValue& myvals = myatoms.getUnderlyingMultiValue();
+  double d2, sw, dfunc; MultiValue& myvals = myatoms.getUnderlyingMultiValue();
   std::vector<double> values( getBaseMultiColvar(0)->getNumberOfQuantities() ); 
 
   getInputData( 0, false, myatoms, values );
@@ -156,7 +157,6 @@ double LocalAverage::compute( const unsigned& tindex, AtomValuePack& myatoms ) c
   }
 
   if( !doNotCalculateDerivatives() ){
-      atom0=getCentralAtomPackFromInput( myatoms.getIndex(0) );
       MultiValue& myder=getInputDerivatives( 0, false, myatoms );
       if( values.size()>2 ){
           for(unsigned j=0;j<myder.getNumberActive();++j){
@@ -198,7 +198,6 @@ double LocalAverage::compute( const unsigned& tindex, AtomValuePack& myatoms ) c
          if( !doNotCalculateDerivatives() ){
              Tensor vir(distance,distance);
              MultiValue& myder=getInputDerivatives( i, false, myatoms );
-             atom1=getCentralAtomPackFromInput( myatoms.getIndex(i) );
              if( values.size()>2 ){
                  for(unsigned j=0;j<myder.getNumberActive();++j){
                      unsigned jder=myder.getActiveIndex(j);
@@ -208,8 +207,8 @@ double LocalAverage::compute( const unsigned& tindex, AtomValuePack& myatoms ) c
                      } 
                  }
                  for(unsigned k=2;k<values.size();++k){
-                     myatoms.addComDerivatives( k, (-dfunc)*values[0]*values[k]*distance, atom0 );
-                     myatoms.addComDerivatives( k, (+dfunc)*values[0]*values[k]*distance, atom1 );
+                     addAtomDerivatives( k, 0, (-dfunc)*values[0]*values[k]*distance, myatoms ); 
+                     addAtomDerivatives( k, i, (+dfunc)*values[0]*values[k]*distance, myatoms );
                      myatoms.addBoxDerivatives( k, (-dfunc)*values[0]*values[k]*vir );
                  }
              } else {
@@ -218,13 +217,13 @@ double LocalAverage::compute( const unsigned& tindex, AtomValuePack& myatoms ) c
                      myatoms.addDerivative( 1, jder, sw*values[0]*myder.getDerivative(1,jder) ); 
                      myatoms.addDerivative( 1, jder, sw*values[1]*myder.getDerivative(0,jder) );
                  }
-                 myatoms.addComDerivatives( 1, (-dfunc)*values[0]*values[1]*distance, atom0 );
-                 myatoms.addComDerivatives( 1, (+dfunc)*values[0]*values[1]*distance, atom1 );
+                 addAtomDerivatives( 1, 0, (-dfunc)*values[0]*values[1]*distance, myatoms ); 
+                 addAtomDerivatives( 1, i, (+dfunc)*values[0]*values[1]*distance, myatoms );
                  myatoms.addBoxDerivatives( 1, (-dfunc)*values[0]*values[1]*vir );
              }
              // And the bit we use to average the vector
-             myatoms.addComDerivatives( -1, (-dfunc)*values[0]*distance, atom0 );
-             myatoms.addComDerivatives( -1, (+dfunc)*values[0]*distance, atom1 );
+             addAtomDerivatives( -1, 0, (-dfunc)*values[0]*distance, myatoms );
+             addAtomDerivatives( -1, i, (+dfunc)*values[0]*distance, myatoms ); 
              for(unsigned j=0;j<myder.getNumberActive();++j){
                  unsigned jder=myder.getActiveIndex(j); myvals.addTemporyDerivative( jder, sw*myder.getDerivative(0, jder) );
              }
diff --git a/src/multicolvar/MultiColvar.cpp b/src/multicolvar/MultiColvar.cpp
deleted file mode 100644
index bc8db06eb..000000000
--- a/src/multicolvar/MultiColvar.cpp
+++ /dev/null
@@ -1,89 +0,0 @@
-/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-   Copyright (c) 2012-2016 The plumed team
-   (see the PEOPLE file at the root of the distribution for a list of names)
-
-   See http://www.plumed.org for more information.
-
-   This file is part of plumed, version 2.
-
-   plumed is free software: you can redistribute it and/or modify
-   it under the terms of the GNU Lesser General Public License as published by
-   the Free Software Foundation, either version 3 of the License, or
-   (at your option) any later version.
-
-   plumed is distributed in the hope that it will be useful,
-   but WITHOUT ANY WARRANTY; without even the implied warranty of
-   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-   GNU Lesser General Public License for more details.
-
-   You should have received a copy of the GNU Lesser General Public License
-   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
-+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvar.h"
-#include "core/PlumedMain.h"
-#include "core/ActionSet.h"
-#include "core/SetupMolInfo.h"
-#include "vesselbase/Vessel.h"
-#include "tools/Pbc.h"
-#include <vector>
-#include <string>
-
-using namespace std;
-namespace PLMD{
-namespace multicolvar{
-
-void MultiColvar::registerKeywords( Keywords& keys ){
-  MultiColvarBase::registerKeywords( keys );
-  keys.reserve("numbered","ATOMS","the atoms involved in each of the collective variables you wish to calculate. "
-                               "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one CV will be "
-                               "calculated for each ATOM keyword you specify (all ATOM keywords should "
-                               "define the same number of atoms).  The eventual number of quantities calculated by this "
-                               "action will depend on what functions of the distribution you choose to calculate."); 
-  keys.reset_style("ATOMS","atoms");
-} 
-
-MultiColvar::MultiColvar(const ActionOptions&ao):
-Action(ao),
-MultiColvarBase(ao)
-{
-}
-
-void MultiColvar::readAtoms( int& natoms, std::vector<AtomNumber> all_atoms ){
-  if( atom_lab.size()==0 && keywords.exists("ATOMS")  ) readAtomsLikeKeyword( "ATOMS", natoms, all_atoms );
-  // Setup the multicolvar base
-  setupMultiColvarBase( all_atoms );
-}
-
-void MultiColvar::readAtomsLikeKeyword( const std::string & key, int& natoms, std::vector<AtomNumber>& all_atoms ){ 
-  plumed_assert( !usespecies );
-  if( all_atoms.size()>0 ) return; 
-
-  std::vector<AtomNumber> t; 
-  for(int i=1;;++i ){
-     parseAtomList(key, i, t );
-     if( t.empty() ) break;
-
-     log.printf("  Colvar %d is calculated from atoms : ", i);
-     for(unsigned j=0;j<t.size();++j) log.printf("%d ",t[j].serial() );
-     log.printf("\n"); 
-
-     if( i==1 && natoms<0 ){ natoms=t.size(); ablocks.resize(natoms); }
-     else if( i==1 ) ablocks.resize(natoms);
-     if( t.size()!=natoms ){
-        std::string ss; Tools::convert(i,ss); 
-        error(key + ss + " keyword has the wrong number of atoms"); 
-     }
-     for(unsigned j=0;j<natoms;++j){ 
-        ablocks[j].push_back( natoms*(i-1)+j ); all_atoms.push_back( t[j] ); 
-        atom_lab.push_back( std::pair<unsigned,unsigned>( 0, natoms*(i-1)+j ) );
-     }
-     t.resize(0); 
-  }
-  if( all_atoms.size()>0 ){
-     nblock=0; 
-     for(unsigned i=0;i<ablocks[0].size();++i) addTaskToList( i );
-  }
-}
-     
-}
-}
diff --git a/src/multicolvar/MultiColvar.h b/src/multicolvar/MultiColvar.h
deleted file mode 100644
index 4544d285e..000000000
--- a/src/multicolvar/MultiColvar.h
+++ /dev/null
@@ -1,56 +0,0 @@
-/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-   Copyright (c) 2012-2016 The plumed team
-   (see the PEOPLE file at the root of the distribution for a list of names)
-
-   See http://www.plumed.org for more information.
-
-   This file is part of plumed, version 2.
-
-   plumed is free software: you can redistribute it and/or modify
-   it under the terms of the GNU Lesser General Public License as published by
-   the Free Software Foundation, either version 3 of the License, or
-   (at your option) any later version.
-
-   plumed is distributed in the hope that it will be useful,
-   but WITHOUT ANY WARRANTY; without even the implied warranty of
-   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-   GNU Lesser General Public License for more details.
-
-   You should have received a copy of the GNU Lesser General Public License
-   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
-+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#ifndef __PLUMED_multicolvar_MultiColvar_h
-#define __PLUMED_multicolvar_MultiColvar_h
-
-#include "MultiColvarBase.h"
-#include "AtomValuePack.h"
-#include "tools/SwitchingFunction.h"
-#include <vector>
-
-#define PLUMED_MULTICOLVAR_INIT(ao) Action(ao),MultiColvar(ao)
-
-namespace PLMD {
-namespace multicolvar {
-
-/**
-\ingroup INHERIT
-This is the abstract base class to use for creating distributions of colvars and functions
-thereof, whtin it there is \ref AddingAMultiColvar "information" as to how to go implementing these types of actions.
-*/
-
-class MultiColvar : public MultiColvarBase {
-protected:
-/// Read in all the keywords that can be used to define atoms
-  void readAtoms( int& natoms, std::vector<AtomNumber> all_atoms );
-/// Read in ATOMS keyword
-  void readAtomsLikeKeyword( const std::string & key, int& natoms, std::vector<AtomNumber>& all_atoms );
-public:
-  explicit MultiColvar(const ActionOptions&);
-  ~MultiColvar(){}
-  static void registerKeywords( Keywords& keys );
-};
-
-}
-}
-
-#endif
diff --git a/src/multicolvar/MultiColvarBase.cpp b/src/multicolvar/MultiColvarBase.cpp
index 939227a29..a955c841e 100644
--- a/src/multicolvar/MultiColvarBase.cpp
+++ b/src/multicolvar/MultiColvarBase.cpp
@@ -28,8 +28,6 @@
 #include "core/ActionSet.h"
 #include "tools/Pbc.h"
 #include "AtomValuePack.h"
-#include "CatomPack.h"
-#include "CatomPack.h"
 #include <vector>
 #include <string>
 
@@ -96,6 +94,37 @@ nblock(0)
   if( keywords.exists("SPECIESA") ){ matsums=usespecies=true; }
 }
 
+void MultiColvarBase::readAtomsLikeKeyword( const std::string & key, const int& natoms, std::vector<AtomNumber>& all_atoms ){
+  plumed_assert( !usespecies );
+  if( all_atoms.size()>0 ) return;
+
+  std::vector<AtomNumber> t;
+  for(int i=1;;++i ){
+     parseAtomList(key, i, t );
+     if( t.empty() ) break;
+
+     log.printf("  Colvar %d is calculated from atoms : ", i);
+     for(unsigned j=0;j<t.size();++j) log.printf("%d ",t[j].serial() );
+     log.printf("\n");
+
+     if( i==1 && natoms<0 ){ ablocks.resize(t.size()); }
+     else if( i==1 ) ablocks.resize(natoms);
+     if( t.size()!=ablocks.size() ){
+        std::string ss; Tools::convert(i,ss);
+        error(key + 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] );
+        atom_lab.push_back( std::pair<unsigned,unsigned>( 0, ablocks.size()*(i-1)+j ) );
+     }
+     t.resize(0);
+  }
+  if( all_atoms.size()>0 ){
+     nblock=0;
+     for(unsigned i=0;i<ablocks[0].size();++i) addTaskToList( i );
+  }
+}
+
 bool MultiColvarBase::parseMultiColvarAtomList(const std::string& key, const int& num, std::vector<AtomNumber>& t){
   std::vector<std::string> mlabs; 
   if( num<0 ) parseVector(key,mlabs);
@@ -310,6 +339,34 @@ void MultiColvarBase::readThreeGroups( const std::string& key1, const std::strin
   }
 }
 
+void MultiColvarBase::buildSets(){
+  std::vector<AtomNumber> fake_atoms;
+  if( !parseMultiColvarAtomList("DATA",-1,fake_atoms) ) error("missing DATA keyword");
+  if( fake_atoms.size()>0 ) error("no atoms should appear in the specification for this object.  Input should be other multicolvars");
+
+  nblock = mybasemulticolvars[0]->getFullNumberOfTasks();
+  for(unsigned i=0;i<mybasemulticolvars.size();++i){
+     if( mybasemulticolvars[i]->getFullNumberOfTasks()!=nblock ){
+          error("mismatch between numbers of tasks in various base multicolvars");
+     }
+  }
+  ablocks.resize( mybasemulticolvars.size() ); usespecies=false;
+  for(unsigned i=0;i<mybasemulticolvars.size();++i){
+      ablocks[i].resize( nblock );
+      for(unsigned j=0;j<nblock;++j) ablocks[i][j]=i*nblock+j;
+  }
+  for(unsigned i=0;i<nblock;++i){
+      if( mybasemulticolvars.size()<4 ){
+          unsigned cvcode=0, tmpc=1;
+          for(unsigned j=0;j<ablocks.size();++j){ cvcode += i*tmpc; tmpc *= nblock; }
+          addTaskToList( cvcode );
+      } else {
+          addTaskToList( i );
+      }
+  }
+  mybasedata[0]->resizeTemporyMultiValues( mybasemulticolvars.size() ); setupMultiColvarBase( fake_atoms );
+}
+
 void MultiColvarBase::addTaskToList( const unsigned& taskCode ){
   plumed_assert( getNumberOfVessels()==0 );
   ActionWithVessel::addTaskToList( taskCode );
@@ -388,7 +445,11 @@ void MultiColvarBase::setupMultiColvarBase( const std::vector<AtomNumber>& atoms
          }  
      }
   }
-  if( mybasemulticolvars.size()>0 ){ for(unsigned i=0;i<mybasedata.size();++i) mybasedata[i]->resizeTemporyMultiValues(2); }
+  if( mybasemulticolvars.size()>0 ){ 
+      for(unsigned i=0;i<mybasedata.size();++i){
+          mybasedata[i]->resizeTemporyMultiValues(2); mybasemulticolvars[i]->my_tmp_capacks.resize(2);
+      }
+  }
 
   // Copy lists of atoms involved from base multicolvars 
   std::vector<AtomNumber> tmp_atoms;
@@ -786,10 +847,12 @@ void MultiColvarBase::addComDerivatives( const int& ival, const unsigned& iatom,
   unsigned katom = myatoms.getIndex( iatom ); plumed_dbg_assert( atom_lab[katom].first>0 );
   // Find base colvar
   unsigned mmc = atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
+  if( usespecies && iatom==0 ){ myatoms.addComDerivatives( ival, der, mybasemulticolvars[mmc]->my_tmp_capacks[0] ); return; }
+
   // Get start of indices for this atom
   unsigned basen=0; for(unsigned i=0;i<mmc;++i) basen+=(mybasemulticolvars[i]->getNumberOfDerivatives() - 9) / 3;
-  multicolvar::CatomPack atom0=mybasemulticolvars[mmc]->getCentralAtomPack( basen, atom_lab[katom].second );
-  myatoms.addComDerivatives( ival, der, atom0 );
+  mybasemulticolvars[mmc]->getCentralAtomPack( basen, atom_lab[katom].second, mybasemulticolvars[mmc]->my_tmp_capacks[1] );
+  myatoms.addComDerivatives( ival, der, mybasemulticolvars[mmc]->my_tmp_capacks[1] );
 }
 
 void MultiColvarBase::getInputData( const unsigned& ind, const bool& normed, 
@@ -919,26 +982,26 @@ void MultiColvarBase::performTask( const unsigned& task_index, const unsigned& c
          myder.clearAll();
       }
   }
+  // Retrieve derivative stuff for central atom
+  if( !doNotCalculateDerivatives() ){
+      if( usespecies && mybasemulticolvars.size()>0 && atom_lab[myatoms.getIndex(0)].first>0 ){
+          unsigned mmc = atom_lab[0].first - 1;
+          MultiValue& myder=mybasedata[mmc]->getTemporyMultiValue(0);
+          if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ||
+              myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ){
+                  myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
+          }
+          mybasedata[mmc]->retrieveDerivatives( atom_lab[myatoms.getIndex(0)].second, false, myder );
+          unsigned basen=0; for(unsigned i=0;i<mmc;++i) basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
+          mybasemulticolvars[mmc]->getCentralAtomPack( basen, atom_lab[myatoms.getIndex(0)].second,  mybasemulticolvars[mmc]->my_tmp_capacks[0] );
+      }
+  }
   // Compute everything
-  double vv=doCalculation( task_index, myatoms ); 
+  double vv=compute( task_index, myatoms ); updateActiveAtoms( myatoms );
   myatoms.setValue( 1, vv );
   return;
 }
 
-double MultiColvarBase::doCalculation( const unsigned& taskIndex, AtomValuePack& myatoms ) const {
-  if( usespecies && mybasemulticolvars.size()>0 && atom_lab[myatoms.getIndex(0)].first>0 ){  
-      unsigned mmc = atom_lab[0].first - 1; 
-      MultiValue& myder=mybasedata[mmc]->getTemporyMultiValue(0);
-      if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ||
-          myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ){
-              myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
-      }
-      mybasedata[mmc]->retrieveDerivatives( atom_lab[myatoms.getIndex(0)].second, false, myder );  
-  }
-  double val=compute( taskIndex, myatoms ); updateActiveAtoms( myatoms );
-  return val;
-}
-
 void MultiColvarBase::updateActiveAtoms( AtomValuePack& myatoms ) const {
   if( mybasemulticolvars.size()==0 ) myatoms.updateUsingIndices();
   else myatoms.updateDynamicList();
@@ -966,16 +1029,16 @@ Vector MultiColvarBase::getCentralAtomPos( const unsigned& taskIndex ){
   }
 }
 
-CatomPack MultiColvarBase::getCentralAtomPack( const unsigned& basn, const unsigned& taskIndex ){
+void MultiColvarBase::getCentralAtomPack( const unsigned& basn, const unsigned& taskIndex, CatomPack& mypack ){
   unsigned curr=getTaskCode( taskIndex );
 
-  CatomPack mypack;
   if(usespecies){
-     mypack.resize(1);
+     if( mypack.getNumberOfAtomsWithDerivatives()!=1 ) mypack.resize(1);
      mypack.setIndex( 0, basn + curr );
      mypack.setDerivative( 0, Tensor::identity() );
   } else if( nblock>0 ){
-     mypack.resize(ncentral); unsigned k=0;
+     if( mypack.getNumberOfAtomsWithDerivatives()!=ncentral ) mypack.resize(ncentral); 
+     unsigned k=0;
      std::vector<unsigned> atoms( ablocks.size() ); decodeIndexToAtoms( curr, atoms );
      for(unsigned i=0;i<ablocks.size();++i){
          if( use_for_central_atom[i] ){
@@ -985,7 +1048,8 @@ CatomPack MultiColvarBase::getCentralAtomPack( const unsigned& basn, const unsig
          }
      }
   } else {
-     mypack.resize(ncentral); unsigned k=0;
+     if( mypack.getNumberOfAtomsWithDerivatives()!=ncentral ) mypack.resize(ncentral); 
+     unsigned k=0;
      for(unsigned i=0;i<ablocks.size();++i){
          if( use_for_central_atom[i] ){
              mypack.setIndex( k, basn + ablocks[i][curr] );
@@ -994,7 +1058,6 @@ CatomPack MultiColvarBase::getCentralAtomPack( const unsigned& basn, const unsig
          }
      }
   }
-  return mypack;
 } 
 
 Vector MultiColvarBase::getSeparation( const Vector& vec1, const Vector& vec2 ) const {
diff --git a/src/multicolvar/MultiColvarBase.h b/src/multicolvar/MultiColvarBase.h
index ca5e93fb1..5dc5b82e5 100644
--- a/src/multicolvar/MultiColvarBase.h
+++ b/src/multicolvar/MultiColvarBase.h
@@ -28,13 +28,13 @@
 #include "tools/LinkCells.h"
 #include "vesselbase/StoreDataVessel.h"
 #include "vesselbase/ActionWithVessel.h"
+#include "CatomPack.h"
 #include <vector>
 
 namespace PLMD {
 namespace multicolvar {
 
 class AtomValuePack;
-class CatomPack;
 class BridgedMultiColvarFunction;
 class ActionVolume;
 
@@ -66,6 +66,8 @@ private:
   unsigned ncentral;
 /// Bool vector telling us which atoms are required to calculate central atom position
   std::vector<bool> use_for_central_atom;
+/// Vector of tempory holders for central atom values
+  std::vector<CatomPack> my_tmp_capacks;
 /// 1/number of atoms involved in central atoms
   double numberForCentralAtom;
 /// Ensures that setup is only performed once per loop
@@ -134,6 +136,8 @@ protected:
   void decodeIndexToAtoms( const unsigned& taskCode, std::vector<unsigned>& atoms ) const ;
 /// Read in some atoms
   bool parseMultiColvarAtomList(const std::string& key, const int& num, std::vector<AtomNumber>& t);
+/// Read in ATOMS keyword
+  void readAtomsLikeKeyword( const std::string & key, const int& natoms, std::vector<AtomNumber>& all_atoms );
 /// Read in two groups of atoms and setup multicolvars to calculate
   void readTwoGroups( const std::string& key0, const std::string& key1, const std::string& key2, std::vector<AtomNumber>& all_atoms );
 /// Read in three groups of atoms
@@ -142,6 +146,8 @@ protected:
 /// Read in three groups of atoms and construct CVs involving at least one
   void readThreeGroups( const std::string& key1, const std::string& key2, const std::string& key3,
                         const bool& allow2, const bool& no_third_dim_accum, std::vector<AtomNumber>& all_atoms );
+/// Build sets by taking one multicolvar from each base
+  void buildSets();
 public:
   explicit MultiColvarBase(const ActionOptions&);
   ~MultiColvarBase(){}
@@ -170,8 +176,6 @@ public:
   virtual Vector getPositionOfAtomForLinkCells( const unsigned& iatom ) const ;
 /// Returns the position where we should assume the center is for link cell calculations
   virtual Vector getLinkCellPosition( const std::vector<unsigned>& atoms ) const ;
-/// And a virtual function which actually computes the colvar
-  virtual double doCalculation( const unsigned& tindex, AtomValuePack& myatoms ) const ;  
 /// Get the absolute index of the central atom
   virtual AtomNumber getAbsoluteIndexOfCentralAtom( const unsigned& i ) const ;
 /// This is replaced once we have a function to calculate the cv
@@ -183,7 +187,7 @@ public:
 /// Checks if an task is being performed at the present time
   virtual bool isCurrentlyActive( const unsigned& code );
 ///
-  virtual CatomPack getCentralAtomPack( const unsigned& basn, const unsigned& curr );
+  virtual void getCentralAtomPack( const unsigned& basn, const unsigned& curr, CatomPack& mypack);
 /// Get the index where the central atom is stored
   virtual Vector getCentralAtomPos( const unsigned& curr );
 /// You can use this to screen contributions that are very small so we can avoid expensive (and pointless) calculations
diff --git a/src/multicolvar/MultiColvarCombine.cpp b/src/multicolvar/MultiColvarCombine.cpp
index 1b0a1c3ea..8e5058851 100644
--- a/src/multicolvar/MultiColvarCombine.cpp
+++ b/src/multicolvar/MultiColvarCombine.cpp
@@ -19,7 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvarFunction.h"
+#include "MultiColvarBase.h"
+#include "AtomValuePack.h"
 #include "core/ActionRegister.h"
 
 #include <string>
@@ -37,7 +38,7 @@ Calculate linear combinations of multiple multicolvars
 namespace PLMD {
 namespace multicolvar {
 
-class MultiColvarCombine : public MultiColvarFunction {
+class MultiColvarCombine : public MultiColvarBase {
 private:
   std::vector<double> coeff;
 public:
@@ -52,16 +53,18 @@ public:
 PLUMED_REGISTER_ACTION(MultiColvarCombine,"MCOLV_COMBINE")
 
 void MultiColvarCombine::registerKeywords( Keywords& keys ){
-  MultiColvarFunction::registerKeywords( keys );
+  MultiColvarBase::registerKeywords( keys );
+  keys.add("compulsory","DATA","the multicolvars you are calculating linear combinations for");
   keys.add("compulsory","COEFFICIENTS","1.0","the coeficients to use for the various multicolvars");
-  keys.use("MEAN"); keys.use("MORE_THAN"); keys.use("SUM"); keys.use("LESS_THAN"); keys.use("HISTOGRAM"); keys.use("HISTOGRAM");
+  keys.use("MEAN"); keys.use("MORE_THAN"); keys.use("SUM"); keys.use("LESS_THAN"); keys.use("HISTOGRAM"); keys.use("HISTOGRAM"); 
   keys.use("MIN"); keys.use("MAX"); keys.use("LOWEST"); keys.use("HIGHEST"); keys.use("ALT_MIN"); keys.use("BETWEEN"); keys.use("MOMENTS");
 }
 
 MultiColvarCombine::MultiColvarCombine(const ActionOptions& ao):
 Action(ao),
-MultiColvarFunction(ao)
+MultiColvarBase(ao)
 {
+  buildSets();
   for(unsigned i=0;i<getNumberOfBaseMultiColvars();++i){
       if( mybasemulticolvars[i]->weightWithDerivatives() ) error("cannot combine multicolvars with weights");
   }
@@ -69,7 +72,7 @@ MultiColvarFunction(ao)
   parseVector("COEFFICIENTS",coeff);
   log.printf("  coefficients of multicolvars %f", coeff[0] );
   for(unsigned i=1;i<coeff.size();++i) log.printf(", %f", coeff[i] );
-  log.printf("\n"); buildSets();
+  log.printf("\n"); 
 }
 
 double MultiColvarCombine::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
diff --git a/src/multicolvar/MultiColvarFunction.cpp b/src/multicolvar/MultiColvarFunction.cpp
deleted file mode 100644
index b392fb559..000000000
--- a/src/multicolvar/MultiColvarFunction.cpp
+++ /dev/null
@@ -1,76 +0,0 @@
-/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-   Copyright (c) 2013-2016 The plumed team
-   (see the PEOPLE file at the root of the distribution for a list of names)
-
-   See http://www.plumed.org for more information.
-
-   This file is part of plumed, version 2.
-
-   plumed is free software: you can redistribute it and/or modify
-   it under the terms of the GNU Lesser General Public License as published by
-   the Free Software Foundation, either version 3 of the License, or
-   (at your option) any later version.
-
-   plumed is distributed in the hope that it will be useful,
-   but WITHOUT ANY WARRANTY; without even the implied warranty of
-   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-   GNU Lesser General Public License for more details.
-
-   You should have received a copy of the GNU Lesser General Public License
-   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
-+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvarFunction.h"
-#include "core/PlumedMain.h"
-#include "core/ActionSet.h"
-#include "tools/Pbc.h"
-#include "MultiColvarBase.h"
-
-namespace PLMD {
-namespace multicolvar { 
-
-void MultiColvarFunction::registerKeywords( Keywords& keys ){
-  MultiColvarBase::registerKeywords( keys );
-  keys.add("compulsory","DATA","the labels of the action that calculates the multicolvars we are interested in");
-  keys.remove("NUMERICAL_DERIVATIVES");
-}
-
-MultiColvarFunction::MultiColvarFunction(const ActionOptions& ao):
-Action(ao),
-MultiColvarBase(ao)
-{
-  // Read in the arguments
-  if( keywords.exists("DATA") ){ 
-      std::vector<AtomNumber> fake_atoms; 
-      if( !parseMultiColvarAtomList("DATA",-1,fake_atoms) ) error("missing DATA keyword");
-      if( fake_atoms.size()>0 ) error("no atoms should appear in the specification for this object.  Input should be other multicolvars");
-  }
-}
-
-void MultiColvarFunction::buildSets(){
-  nblock = mybasemulticolvars[0]->getFullNumberOfTasks();
-  for(unsigned i=0;i<mybasemulticolvars.size();++i){
-     if( mybasemulticolvars[i]->getFullNumberOfTasks()!=nblock ){
-          error("mismatch between numbers of tasks in various base multicolvars");
-     }
-  }
-  ablocks.resize( mybasemulticolvars.size() ); usespecies=false; 
-  for(unsigned i=0;i<mybasemulticolvars.size();++i){
-      ablocks[i].resize( nblock ); 
-      for(unsigned j=0;j<nblock;++j) ablocks[i][j]=i*nblock+j;  
-  }
-  for(unsigned i=0;i<nblock;++i){
-      if( mybasemulticolvars.size()<4 ){
-          unsigned cvcode=0, tmpc=1;
-          for(unsigned j=0;j<ablocks.size();++j){ cvcode += i*tmpc; tmpc *= nblock; }
-          addTaskToList( cvcode );
-      } else {
-          addTaskToList( i );
-      }
-  }
-  mybasedata[0]->resizeTemporyMultiValues( mybasemulticolvars.size() );
-  std::vector<AtomNumber> fake_atoms; setupMultiColvarBase( fake_atoms );
-}
-
-}
-}
-
diff --git a/src/multicolvar/MultiColvarFunction.h b/src/multicolvar/MultiColvarFunction.h
deleted file mode 100644
index 60877a639..000000000
--- a/src/multicolvar/MultiColvarFunction.h
+++ /dev/null
@@ -1,58 +0,0 @@
-/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-   Copyright (c) 2013-2016 The plumed team
-   (see the PEOPLE file at the root of the distribution for a list of names)
-
-   See http://www.plumed.org for more information.
-
-   This file is part of plumed, version 2.
-
-   plumed is free software: you can redistribute it and/or modify
-   it under the terms of the GNU Lesser General Public License as published by
-   the Free Software Foundation, either version 3 of the License, or
-   (at your option) any later version.
-
-   plumed is distributed in the hope that it will be useful,
-   but WITHOUT ANY WARRANTY; without even the implied warranty of
-   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-   GNU Lesser General Public License for more details.
-
-   You should have received a copy of the GNU Lesser General Public License
-   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
-+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#ifndef __PLUMED_multicolvar_MultiColvarFunction_h
-#define __PLUMED_multicolvar_MultiColvarFunction_h
-
-#include "tools/Matrix.h"
-#include "MultiColvarBase.h"
-#include "AtomValuePack.h"
-#include "CatomPack.h"
-#include "vesselbase/StoreDataVessel.h"
-
-namespace PLMD {
-namespace multicolvar {
-
-class MultiColvarFunction : public MultiColvarBase {
-private:
-/// A tempory vector that is used for retrieving vectors
-  std::vector<double> tvals;
-protected:
-/// Get the derivatives for the central atom with index ind
-  CatomPack getCentralAtomPackFromInput( const unsigned& ind ) const ;
-/// Build sets by taking one multicolvar from each base
-  void buildSets();
-public:
-  explicit MultiColvarFunction(const ActionOptions&);
-  static void registerKeywords( Keywords& keys );
-};
-
-inline
-CatomPack MultiColvarFunction::getCentralAtomPackFromInput( const unsigned& ind ) const {
-  plumed_dbg_assert( atom_lab[ind].first>0 ); unsigned mmc=atom_lab[ind].first-1;
-  unsigned basen=0;
-  for(unsigned i=0;i<mmc;++i) basen+=mybasemulticolvars[i]->getNumberOfAtoms();
-  return mybasemulticolvars[mmc]->getCentralAtomPack( basen, atom_lab[ind].second );
-}
-
-}
-}
-#endif
diff --git a/src/multicolvar/NumberOfLinks.cpp b/src/multicolvar/NumberOfLinks.cpp
index 060481d08..b5c351344 100644
--- a/src/multicolvar/NumberOfLinks.cpp
+++ b/src/multicolvar/NumberOfLinks.cpp
@@ -19,7 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvarFunction.h"
+#include "MultiColvarBase.h"
+#include "AtomValuePack.h"
 #include "core/ActionRegister.h"
 #include "tools/SwitchingFunction.h"
 
@@ -65,7 +66,7 @@ PRINT ARG=dd FILE=colvar
 namespace PLMD {
 namespace multicolvar {
 
-class NumberOfLinks : public MultiColvarFunction {
+class NumberOfLinks : public MultiColvarBase {
 private:
 /// The values of the quantities in the dot products
   std::vector<double> orient0, orient1; 
@@ -85,7 +86,7 @@ public:
 PLUMED_REGISTER_ACTION(NumberOfLinks,"NLINKS")
 
 void NumberOfLinks::registerKeywords( Keywords& keys ){
-  MultiColvarFunction::registerKeywords( keys );
+  MultiColvarBase::registerKeywords( keys );
   keys.add("atoms","GROUP","");
   keys.add("atoms-1","GROUPA","");
   keys.add("atoms-1","GROUPB","");
@@ -97,13 +98,13 @@ void NumberOfLinks::registerKeywords( Keywords& keys ){
                                "The following provides information on the \\ref switchingfunction that are available. "
                                "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
   // Use actionWithDistributionKeywords
-  keys.remove("LOWMEM"); keys.remove("DATA");
+  keys.remove("LOWMEM"); 
   keys.addFlag("LOWMEM",false,"lower the memory requirements");
 }
 
 NumberOfLinks::NumberOfLinks(const ActionOptions& ao):
 Action(ao),
-MultiColvarFunction(ao)
+MultiColvarBase(ao)
 {
   // The weight of this does have derivatives
   weightHasDerivatives=true;
diff --git a/src/multicolvar/Torsions.cpp b/src/multicolvar/Torsions.cpp
index eefd0e79b..57af60a91 100644
--- a/src/multicolvar/Torsions.cpp
+++ b/src/multicolvar/Torsions.cpp
@@ -19,7 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvar.h"
+#include "MultiColvarBase.h"
+#include "AtomValuePack.h"
 #include "tools/Torsion.h"
 #include "core/ActionRegister.h"
 
@@ -71,7 +72,7 @@ Similarly \@psi-4 tells plumed that you want to calculate the \f$\psi\f$ angle o
 */
 //+ENDPLUMEDOC
 
-class Torsions : public MultiColvar {
+class Torsions : public MultiColvarBase {
 public:
   static void registerKeywords( Keywords& keys );
   explicit Torsions(const ActionOptions&);
@@ -83,16 +84,24 @@ public:
 PLUMED_REGISTER_ACTION(Torsions,"TORSIONS")
 
 void Torsions::registerKeywords( Keywords& keys ){
-  MultiColvar::registerKeywords( keys );
-  keys.use("ATOMS"); keys.use("BETWEEN"); keys.use("HISTOGRAM");
+  MultiColvarBase::registerKeywords( keys );
+  keys.reserve("numbered","ATOMS","the atoms involved in each of the torsion angles you wish to calculate. "
+                             "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one torsion will be "
+                             "calculated for each ATOM keyword you specify (all ATOM keywords should "
+                             "provide the indices of four atoms).  The eventual number of quantities calculated by this "
+                             "action will depend on what functions of the distribution you choose to calculate.");
+  keys.reset_style("ATOMS","atoms");
+  keys.use("BETWEEN"); keys.use("HISTOGRAM");
 }
 
 Torsions::Torsions(const ActionOptions&ao):
-PLUMED_MULTICOLVAR_INIT(ao)
+Action(ao),
+MultiColvarBase(ao)
 {
   // Read in the atoms
   int natoms=4; std::vector<AtomNumber> all_atoms; 
-  readAtoms( natoms, all_atoms );
+  readAtomsLikeKeyword( "ATOMS", natoms, all_atoms );
+  setupMultiColvarBase( all_atoms );
   std::vector<bool> catom_ind(4, false); 
   catom_ind[1]=catom_ind[2]=true;
   setAtomsForCentralAtom( catom_ind );
diff --git a/src/multicolvar/VolumeGradientBase.cpp b/src/multicolvar/VolumeGradientBase.cpp
index e30592754..ff2bf534f 100644
--- a/src/multicolvar/VolumeGradientBase.cpp
+++ b/src/multicolvar/VolumeGradientBase.cpp
@@ -65,7 +65,7 @@ void VolumeGradientBase::setNumberInVolume( const unsigned& ivol, const unsigned
   if( !mcolv->weightHasDerivatives ){
       outvals.setValue(ivol, weight ); 
       if( derivativesAreRequired() ){
-         CatomPack catom( mcolv->getCentralAtomPack( 0, curr ) );
+         CatomPack catom; mcolv->getCentralAtomPack( 0, curr, catom );
          for(unsigned i=0;i<catom.getNumberOfAtomsWithDerivatives();++i){
              unsigned jatom=3*catom.getIndex(i);
              outvals.addDerivative( ivol, jatom+0, catom.getDerivative(i,0,wdf) );
@@ -86,7 +86,7 @@ void VolumeGradientBase::setNumberInVolume( const unsigned& ivol, const unsigned
       double ww=outvals.get(0); outvals.setValue(ivol,ww*weight);
       if( derivativesAreRequired() ){
          plumed_merror("This needs testing");
-         CatomPack catom( mcolv->getCentralAtomPack( 0, curr ) );
+         CatomPack catom; mcolv->getCentralAtomPack( 0, curr, catom );
          for(unsigned i=0;i<catom.getNumberOfAtomsWithDerivatives();++i){
              unsigned jatom=3*catom.getIndex(i);
              outvals.addDerivative( ivol, jatom+0, weight*outvals.getDerivative(ivol,jatom+0) + ww*catom.getDerivative(i,0,wdf) );
@@ -106,7 +106,7 @@ void VolumeGradientBase::setNumberInVolume( const unsigned& ivol, const unsigned
       double ww=outvals.get(0); outvals.setValue(ivol,ww*weight);
       if( derivativesAreRequired() ){
          plumed_merror("This needs testing");
-         CatomPack catom( mcolv->getCentralAtomPack( 0, curr ) ); 
+         CatomPack catom; mcolv->getCentralAtomPack( 0, curr, catom ); 
          for(unsigned i=0;i<catom.getNumberOfAtomsWithDerivatives();++i){
              unsigned jatom=3*catom.getIndex(i);
              outvals.addDerivative( ivol, jatom+0, ww*catom.getDerivative(i,0,wdf) );
diff --git a/src/multicolvar/XAngle.cpp b/src/multicolvar/XAngle.cpp
index f1bc4759b..865f96ab7 100644
--- a/src/multicolvar/XAngle.cpp
+++ b/src/multicolvar/XAngle.cpp
@@ -19,7 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvar.h"
+#include "MultiColvarBase.h"
+#include "AtomValuePack.h"
 #include "core/ActionRegister.h"
 #include "tools/Angle.h"
 #include "tools/SwitchingFunction.h"
@@ -82,7 +83,7 @@ PRINT ARG=d1.min
 
 
 
-class XAngles : public MultiColvar {
+class XAngles : public MultiColvarBase {
 private:
   bool use_sf;
   unsigned myc; 
@@ -102,11 +103,17 @@ PLUMED_REGISTER_ACTION(XAngles,"YANGLES")
 PLUMED_REGISTER_ACTION(XAngles,"ZANGLES")
 
 void XAngles::registerKeywords( Keywords& keys ){
-  MultiColvar::registerKeywords( keys );
-  keys.use("ATOMS");  keys.use("MAX"); keys.use("ALT_MIN"); 
+  MultiColvarBase::registerKeywords( keys );
+  keys.use("MAX"); keys.use("ALT_MIN"); 
   keys.use("MEAN"); keys.use("MIN"); keys.use("LESS_THAN");
   keys.use("LOWEST"); keys.use("HIGHEST"); 
   keys.use("MORE_THAN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
+  keys.add("numbered","ATOMS","the atoms involved in each of the angles you wish to calculate. "
+                              "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one angle will be "
+                              "calculated for each ATOM keyword you specify (all ATOM keywords should "
+                              "specify the indices of two atoms).  The eventual number of quantities calculated by this "
+                              "action will depend on what functions of the distribution you choose to calculate.");
+  keys.reset_style("ATOMS","atoms");
   keys.add("atoms-1","GROUP","Calculate the distance between each distinct pair of atoms in the group");
   keys.add("atoms-2","GROUPA","Calculate the distances between all the atoms in GROUPA and all "
                               "the atoms in GROUPB. This must be used in conjuction with GROUPB.");
@@ -117,7 +124,8 @@ void XAngles::registerKeywords( Keywords& keys ){
 }
 
 XAngles::XAngles(const ActionOptions&ao):
-PLUMED_MULTICOLVAR_INIT(ao),
+Action(ao),
+MultiColvarBase(ao),
 use_sf(false)
 {
   if( getName().find("X")!=std::string::npos) myc=0;
@@ -138,7 +146,8 @@ use_sf(false)
   // Read in the atoms
   std::vector<AtomNumber> all_atoms;
   readTwoGroups( "GROUP", "GROUPA", "GROUPB", all_atoms );
-  int natoms=2; readAtoms( natoms, all_atoms );
+  if( atom_lab.size()==0 ) readAtomsLikeKeyword( "ATOMS", 2, all_atoms );
+  setupMultiColvarBase( all_atoms );
   // And check everything has been read in correctly
   checkRead();
 }
diff --git a/src/multicolvar/XDistances.cpp b/src/multicolvar/XDistances.cpp
index 4843ea5de..3aedb3f17 100644
--- a/src/multicolvar/XDistances.cpp
+++ b/src/multicolvar/XDistances.cpp
@@ -19,7 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvar.h"
+#include "MultiColvarBase.h"
+#include "AtomValuePack.h"
 #include "core/ActionRegister.h"
 
 #include <string>
@@ -108,7 +109,7 @@ rather than the x component.
 //+ENDPLUMEDOC
 
 
-class XDistances : public MultiColvar {
+class XDistances : public MultiColvarBase {
 private:
   unsigned myc;
 public:
@@ -125,11 +126,17 @@ PLUMED_REGISTER_ACTION(XDistances,"YDISTANCES")
 PLUMED_REGISTER_ACTION(XDistances,"ZDISTANCES")
 
 void XDistances::registerKeywords( Keywords& keys ){
-  MultiColvar::registerKeywords( keys );
-  keys.use("ATOMS");  keys.use("MAX"); keys.use("ALT_MIN"); 
+  MultiColvarBase::registerKeywords( keys );
+  keys.use("MAX"); keys.use("ALT_MIN"); 
   keys.use("MEAN"); keys.use("MIN"); keys.use("LESS_THAN");
   keys.use("LOWEST"); keys.use("HIGHEST"); 
   keys.use("MORE_THAN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
+  keys.add("numbered","ATOMS","the atoms involved in each of the distances you wish to calculate. "
+                               "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one distance will be "
+                               "calculated for each ATOM keyword you specify (all ATOM keywords should "
+                               "specify the indices of two atoms).  The eventual number of quantities calculated by this "
+                               "action will depend on what functions of the distribution you choose to calculate.");
+  keys.reset_style("ATOMS","atoms");
   keys.add("atoms-1","GROUP","Calculate the distance between each distinct pair of atoms in the group");
   keys.add("atoms-2","GROUPA","Calculate the distances between all the atoms in GROUPA and all "
                               "the atoms in GROUPB. This must be used in conjuction with GROUPB.");
@@ -138,7 +145,8 @@ void XDistances::registerKeywords( Keywords& keys ){
 }
 
 XDistances::XDistances(const ActionOptions&ao):
-PLUMED_MULTICOLVAR_INIT(ao)
+Action(ao),
+MultiColvarBase(ao)
 {
   if( getName().find("X")!=std::string::npos) myc=0;
   else if( getName().find("Y")!=std::string::npos) myc=1;
@@ -148,7 +156,8 @@ PLUMED_MULTICOLVAR_INIT(ao)
   // Read in the atoms
   std::vector<AtomNumber> all_atoms;
   readTwoGroups( "GROUP", "GROUPA", "GROUPB", all_atoms );
-  int natoms=2; readAtoms( natoms, all_atoms );
+  if( atom_lab.size()==0 ) readAtomsLikeKeyword( "ATOMS", 2, all_atoms );
+  setupMultiColvarBase( all_atoms );
   // And check everything has been read in correctly
   checkRead();
 }
diff --git a/src/multicolvar/XYDistances.cpp b/src/multicolvar/XYDistances.cpp
index a800b57ff..e06c49f00 100644
--- a/src/multicolvar/XYDistances.cpp
+++ b/src/multicolvar/XYDistances.cpp
@@ -19,7 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvar.h"
+#include "MultiColvarBase.h"
+#include "AtomValuePack.h"
 #include "core/ActionRegister.h"
 
 #include <string>
@@ -82,7 +83,7 @@ rather than the xy component.
 //+ENDPLUMEDOC
 
 
-class XYDistances : public MultiColvar {
+class XYDistances : public MultiColvarBase {
 private:
   unsigned myc1, myc2;
 public:
@@ -99,11 +100,16 @@ PLUMED_REGISTER_ACTION(XYDistances,"XZDISTANCES")
 PLUMED_REGISTER_ACTION(XYDistances,"YZDISTANCES")
 
 void XYDistances::registerKeywords( Keywords& keys ){
-  MultiColvar::registerKeywords( keys );
-  keys.use("ATOMS"); keys.use("MAX"); keys.use("ALT_MIN"); 
-  keys.use("MEAN"); keys.use("MIN"); keys.use("LESS_THAN");
-  keys.use("LOWEST"); keys.use("HIGHEST"); 
+  MultiColvarBase::registerKeywords( keys );
+  keys.use("MAX"); keys.use("ALT_MIN"); 
+  keys.use("MEAN"); keys.use("MIN"); keys.use("LESS_THAN"); keys.use("LOWEST"); keys.use("HIGHEST"); 
   keys.use("MORE_THAN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
+  keys.add("numbered","ATOMS","the atoms involved in each of the distances you wish to calculate. "
+                               "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one distance will be "
+                               "calculated for each ATOM keyword you specify (all ATOM keywords should "
+                               "specify the incides of two atoms).  The eventual number of quantities calculated by this "
+                               "action will depend on what functions of the distribution you choose to calculate.");
+  keys.reset_style("ATOMS","atoms");
   keys.add("atoms-1","GROUP","Calculate the distance between each distinct pair of atoms in the group");
   keys.add("atoms-2","GROUPA","Calculate the distances between all the atoms in GROUPA and all "
                               "the atoms in GROUPB. This must be used in conjuction with GROUPB.");
@@ -112,7 +118,8 @@ void XYDistances::registerKeywords( Keywords& keys ){
 }
 
 XYDistances::XYDistances(const ActionOptions&ao):
-PLUMED_MULTICOLVAR_INIT(ao)
+Action(ao),
+MultiColvarBase(ao)
 {
   if( getName().find("XY")!=std::string::npos){
       myc1=0; myc2=1;
@@ -125,7 +132,8 @@ PLUMED_MULTICOLVAR_INIT(ao)
   // Read in the atoms
   std::vector<AtomNumber> all_atoms;
   readTwoGroups( "GROUP", "GROUPA", "GROUPB", all_atoms );
-  int natoms=2; readAtoms( natoms, all_atoms );
+  if( atom_lab.size()==0 ) readAtomsLikeKeyword( "ATOMS", 2, all_atoms );
+  setupMultiColvarBase( all_atoms );
   // And check everything has been read in correctly
   checkRead();
 }
diff --git a/src/multicolvar/XYTorsion.cpp b/src/multicolvar/XYTorsion.cpp
index 3a73dc6e9..e92559a45 100644
--- a/src/multicolvar/XYTorsion.cpp
+++ b/src/multicolvar/XYTorsion.cpp
@@ -19,7 +19,8 @@
    You should have received a copy of the GNU Lesser General Public License
    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "MultiColvar.h"
+#include "MultiColvarBase.h"
+#include "AtomValuePack.h"
 #include "core/ActionRegister.h"
 #include "tools/Torsion.h"
 #include "tools/SwitchingFunction.h"
@@ -131,7 +132,7 @@ PRINT ARG=d1.min
 
 
 
-class XYTorsion : public MultiColvar {
+class XYTorsion : public MultiColvarBase {
 private:
   bool use_sf;
   unsigned myc1, myc2; 
@@ -155,11 +156,17 @@ PLUMED_REGISTER_ACTION(XYTorsion,"ZXTORSIONS")
 PLUMED_REGISTER_ACTION(XYTorsion,"ZYTORSIONS")
 
 void XYTorsion::registerKeywords( Keywords& keys ){
-  MultiColvar::registerKeywords( keys );
-  keys.use("ATOMS");  keys.use("MAX"); keys.use("ALT_MIN"); 
+  MultiColvarBase::registerKeywords( keys );
+  keys.use("MAX"); keys.use("ALT_MIN"); 
   keys.use("MEAN"); keys.use("MIN"); 
   keys.use("LOWEST"); keys.use("HIGHEST"); 
   keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
+  keys.add("numbered","ATOMS","the atoms involved in each of the torsions you wish to calculate. "
+                               "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one torsion will be "
+                               "calculated for each ATOM keyword you specify (all ATOM keywords should "
+                               "specify the incides of two atoms).  The eventual number of quantities calculated by this "
+                               "action will depend on what functions of the distribution you choose to calculate.");
+  keys.reset_style("ATOMS","atoms");
   keys.add("atoms-1","GROUP","Calculate the distance between each distinct pair of atoms in the group");
   keys.add("atoms-2","GROUPA","Calculate the distances between all the atoms in GROUPA and all "
                               "the atoms in GROUPB. This must be used in conjuction with GROUPB.");
@@ -170,7 +177,8 @@ void XYTorsion::registerKeywords( Keywords& keys ){
 }
 
 XYTorsion::XYTorsion(const ActionOptions&ao):
-PLUMED_MULTICOLVAR_INIT(ao),
+Action(ao),
+MultiColvarBase(ao),
 use_sf(false)
 {
   if( getName().find("XY")!=std::string::npos){ myc1=0; myc2=1; }
@@ -194,7 +202,8 @@ use_sf(false)
   // Read in the atoms
   std::vector<AtomNumber> all_atoms;
   readTwoGroups( "GROUP", "GROUPA", "GROUPB", all_atoms );
-  int natoms=2; readAtoms( natoms, all_atoms );
+  if( atom_lab.size()==0 ) readAtomsLikeKeyword( "ATOMS", 2, all_atoms );
+  setupMultiColvarBase( all_atoms );
   // And check everything has been read in correctly
   checkRead();
 }
-- 
GitLab