From 7db807933c692a2d87ec0783f67b9b375cadb2d0 Mon Sep 17 00:00:00 2001
From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com>
Date: Thu, 19 Sep 2024 14:01:52 +0200
Subject: [PATCH] Applying the new style to adjmat

---
 src/adjmat/ActionWithInputMatrix.cpp |  54 ++++++++-----
 src/adjmat/ActionWithInputMatrix.h   |   8 +-
 src/adjmat/AdjacencyMatrixBase.cpp   | 105 +++++++++++++++++--------
 src/adjmat/AdjacencyMatrixBase.h     |  11 ++-
 src/adjmat/AdjacencyMatrixVessel.cpp |  77 +++++++++++++------
 src/adjmat/AlignedMatrixBase.cpp     |  36 ++++++---
 src/adjmat/ClusterAnalysisBase.cpp   |  24 ++++--
 src/adjmat/ClusterAnalysisBase.h     |   4 +-
 src/adjmat/ClusterDiameter.cpp       |  22 ++++--
 src/adjmat/ClusterDistribution.cpp   |  49 ++++++++----
 src/adjmat/ClusterProperties.cpp     |  48 ++++++++----
 src/adjmat/ClusterSize.cpp           |  23 ++++--
 src/adjmat/ClusterWithSurface.cpp    |  45 +++++++----
 src/adjmat/ClusteringBase.cpp        |  28 +++++--
 src/adjmat/ContactAlignedMatrix.cpp  |  21 +++--
 src/adjmat/ContactMatrix.cpp         |  22 ++++--
 src/adjmat/DFSClustering.cpp         |  35 ++++++---
 src/adjmat/DumpGraph.cpp             |  38 ++++++---
 src/adjmat/HbondMatrix.cpp           |  62 ++++++++++-----
 src/adjmat/MatrixColumnSums.cpp      |  41 +++++++---
 src/adjmat/MatrixRowSums.cpp         |  34 +++++---
 src/adjmat/OutputCluster.cpp         |  96 ++++++++++++++++-------
 src/adjmat/SMACMatrix.cpp            |  42 +++++++---
 src/adjmat/Sprint.cpp                |  34 +++++---
 src/adjmat/TopologyMatrix.cpp        | 111 ++++++++++++++++++---------
 src/astyle.sh                        |   1 -
 src/config/Config.inc.in             |  24 +++---
 27 files changed, 756 insertions(+), 339 deletions(-)

diff --git a/src/adjmat/ActionWithInputMatrix.cpp b/src/adjmat/ActionWithInputMatrix.cpp
index 449240d64a..41dd1d3241 100644
--- a/src/adjmat/ActionWithInputMatrix.cpp
+++ b/src/adjmat/ActionWithInputMatrix.cpp
@@ -43,20 +43,29 @@ ActionWithInputMatrix::ActionWithInputMatrix(const ActionOptions& ao):
   matsums=true;
   if( keywords.exists("MATRIX") ) {
     std::vector<AtomNumber> fake_atoms;
-    if( !parseMultiColvarAtomList("MATRIX",-1,fake_atoms ) ) error("unable to interpret input matrix");
-    if( mybasemulticolvars.size()!=1 ) error("should be exactly one matrix input");
+    if( !parseMultiColvarAtomList("MATRIX",-1,fake_atoms ) )
+      error("unable to interpret input matrix");
+    if( mybasemulticolvars.size()!=1 )
+      error("should be exactly one matrix input");
 
     // Retrieve the adjacency matrix of interest
     for(unsigned i=0; i<mybasemulticolvars[0]->getNumberOfVessels(); ++i) {
       mymatrix = dynamic_cast<AdjacencyMatrixVessel*>( mybasemulticolvars[0]->getPntrToVessel(i) );
-      if( mymatrix ) break ;
+      if( mymatrix )
+        break ;
     }
-    if( !mymatrix ) error( mybasemulticolvars[0]->getLabel() + " does not calculate an adjacency matrix");
-
-    atom_lab.resize(0); unsigned nnodes; // Delete all the atom labels that have been created
-    if( mymatrix->undirectedGraph() ) nnodes = (mymatrix->function)->ablocks[0].size();
-    else nnodes = (mymatrix->function)->ablocks[0].size() + (mymatrix->function)->ablocks[1].size();
-    for(unsigned i=0; i<nnodes; ++i) atom_lab.push_back( std::pair<unsigned,unsigned>( 1, i ) );
+    if( !mymatrix )
+      error( mybasemulticolvars[0]->getLabel() + " does not calculate an adjacency matrix");
+
+    // Delete all the atom labels that have been created
+    atom_lab.resize(0);
+    unsigned nnodes;
+    if( mymatrix->undirectedGraph() )
+      nnodes = (mymatrix->function)->ablocks[0].size();
+    else
+      nnodes = (mymatrix->function)->ablocks[0].size() + (mymatrix->function)->ablocks[1].size();
+    for(unsigned i=0; i<nnodes; ++i)
+      atom_lab.push_back( std::pair<unsigned,unsigned>( 1, i ) );
   }
 }
 
@@ -77,7 +86,8 @@ AtomNumber ActionWithInputMatrix::getAbsoluteIndexOfCentralAtom(const unsigned&
 }
 
 double ActionWithInputMatrix::retrieveConnectionValue( const unsigned& i, const unsigned& j, std::vector<double>& vals ) const {
-  if( !mymatrix->matrixElementIsActive( i, j ) ) return 0;
+  if( !mymatrix->matrixElementIsActive( i, j ) )
+    return 0;
   unsigned myelem = mymatrix->getStoreIndexFromMatrixIndices( i, j );
 
   // unsigned vi; double df;
@@ -87,18 +97,22 @@ double ActionWithInputMatrix::retrieveConnectionValue( const unsigned& i, const
 
 void ActionWithInputMatrix::getInputData( const unsigned& ind, const bool& normed, const multicolvar::AtomValuePack& myatoms, std::vector<double>& orient0 ) const {
   if( (mymatrix->function)->mybasemulticolvars.size()==0  ) {
-    std::vector<double> tvals( mymatrix->getNumberOfComponents() ); orient0.assign(orient0.size(),0);
+    std::vector<double> tvals( mymatrix->getNumberOfComponents() );
+    orient0.assign(orient0.size(),0);
     for(unsigned i=0; i<mymatrix->getNumberOfColumns(); ++i) {
-      if( mymatrix->undirectedGraph() && ind==i ) continue;
+      if( mymatrix->undirectedGraph() && ind==i )
+        continue;
       orient0[1]+=retrieveConnectionValue( ind, i, tvals );
     }
-    orient0[0]=1.0; return;
+    orient0[0]=1.0;
+    return;
   }
   (mymatrix->function)->getInputData( ind, normed, myatoms, orient0 );
 }
 
 void ActionWithInputMatrix::addConnectionDerivatives( const unsigned& i, const unsigned& j, MultiValue& myvals, MultiValue& myvout ) const {
-  if( !mymatrix->matrixElementIsActive( i, j ) ) return;
+  if( !mymatrix->matrixElementIsActive( i, j ) )
+    return;
   unsigned myelem = mymatrix->getStoreIndexFromMatrixIndices( i, j );
   // Get derivatives and add
   mymatrix->retrieveDerivatives( myelem, false, myvals );
@@ -117,22 +131,26 @@ MultiValue& ActionWithInputMatrix::getInputDerivatives( const unsigned& ind, con
     myder.clearAll();
     MultiValue myvals( (mymatrix->function)->getNumberOfQuantities(), (mymatrix->function)->getNumberOfDerivatives() );
     for(unsigned i=0; i<mymatrix->getNumberOfColumns(); ++i) {
-      if( mymatrix->undirectedGraph() && ind==i ) continue;
+      if( mymatrix->undirectedGraph() && ind==i )
+        continue;
       addConnectionDerivatives( ind, i, myvals, myder );
     }
-    myder.updateDynamicList(); return myder;
+    myder.updateDynamicList();
+    return myder;
   }
   return (mymatrix->function)->getInputDerivatives( ind, normed, myatoms );
 }
 
 unsigned ActionWithInputMatrix::getNumberOfNodeTypes() const {
   unsigned size = (mymatrix->function)->mybasemulticolvars.size();
-  if( size==0 ) return 1;
+  if( size==0 )
+    return 1;
   return size;
 }
 
 unsigned ActionWithInputMatrix::getNumberOfQuantities() const {
-  if( (mymatrix->function)->mybasemulticolvars.size()==0 ) return 2;
+  if( (mymatrix->function)->mybasemulticolvars.size()==0 )
+    return 2;
   return (mymatrix->function)->mybasemulticolvars[0]->getNumberOfQuantities();
 }
 
diff --git a/src/adjmat/ActionWithInputMatrix.h b/src/adjmat/ActionWithInputMatrix.h
index 9d5c44dcdd..e21270468b 100644
--- a/src/adjmat/ActionWithInputMatrix.h
+++ b/src/adjmat/ActionWithInputMatrix.h
@@ -57,12 +57,16 @@ class ActionWithInputMatrix : public multicolvar::MultiColvarBase {
   unsigned getNumberOfDerivatives() override;
 ///  Get the number of rows/cols in the adjacency matrix vessel
   virtual unsigned getNumberOfNodes() const;
-  bool isPeriodic() override { return false; }
+  bool isPeriodic() override {
+    return false;
+  }
   unsigned getNumberOfQuantities() const override;
 ///
   AtomNumber getAbsoluteIndexOfCentralAtom(const unsigned& i) const override;
 /// No loop over tasks for ActionWithInputMatrix
-  double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const override { plumed_error(); }
+  double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const override {
+    plumed_error();
+  }
 ///
   Vector getPositionOfAtomForLinkCells( const unsigned& iatom ) const override;
 };
diff --git a/src/adjmat/AdjacencyMatrixBase.cpp b/src/adjmat/AdjacencyMatrixBase.cpp
index 87c47ca732..5a1acaf62c 100644
--- a/src/adjmat/AdjacencyMatrixBase.cpp
+++ b/src/adjmat/AdjacencyMatrixBase.cpp
@@ -31,7 +31,8 @@ namespace adjmat {
 
 void AdjacencyMatrixBase::registerKeywords( Keywords& keys ) {
   multicolvar::MultiColvarBase::registerKeywords( keys );
-  keys.remove("LOWMEM"); keys.use("HIGHMEM");
+  keys.remove("LOWMEM");
+  keys.use("HIGHMEM");
 }
 
 AdjacencyMatrixBase::AdjacencyMatrixBase(const ActionOptions& ao):
@@ -48,23 +49,28 @@ void AdjacencyMatrixBase::parseConnectionDescriptions( const std::string& key, c
   if( getNumberOfNodeTypes()==1 || (getNumberOfNodeTypes()==2 && nrow_t==1) ) {
     std::vector<std::string> sw;
     if( !multiple ) {
-      sw.resize(1); parse(key,sw[0]);
-      if(sw[0].length()==0) error("could not find " + key + " keyword");
+      sw.resize(1);
+      parse(key,sw[0]);
+      if(sw[0].length()==0)
+        error("could not find " + key + " keyword");
     } else {
       std::string input;
       for(int i=1;; i++) {
-        if( !parseNumbered(key, i, input ) ) break;
+        if( !parseNumbered(key, i, input ) )
+          break;
         sw.push_back( input );
       }
     }
     setupConnector( connect_id, 0, 0, sw );
   } else {
-    if( multiple ) error("keyword " + key + " does not work with multiple input strings");
+    if( multiple )
+      error("keyword " + key + " does not work with multiple input strings");
     unsigned nr, nc;
     if( nrow_t==0 ) {
       nr=nc=getNumberOfNodeTypes();
     } else {
-      nr=nrow_t; nc = getNumberOfNodeTypes() - nr;
+      nr=nrow_t;
+      nc = getNumberOfNodeTypes() - nr;
     }
     for(unsigned i=0; i<nr; ++i) {
       // Retrieve the base number
@@ -78,9 +84,11 @@ void AdjacencyMatrixBase::parseConnectionDescriptions( const std::string& key, c
       }
 
       for(unsigned j=i; j<nc; ++j) {
-        std::vector<std::string> sw(1); parseNumbered(key,ibase+j+1,sw[0]);
+        std::vector<std::string> sw(1);
+        parseNumbered(key,ibase+j+1,sw[0]);
         if(sw[0].length()==0) {
-          std::string num; Tools::convert(ibase+j+1,num);
+          std::string num;
+          Tools::convert(ibase+j+1,num);
           error("could not find " + key + num + " keyword. Need one " + key + " keyword for each distinct base-multicolvar-pair type");
         }
         setupConnector( connect_id, i, j, sw );
@@ -91,18 +99,21 @@ void AdjacencyMatrixBase::parseConnectionDescriptions( const std::string& key, c
 }
 
 unsigned AdjacencyMatrixBase::getSizeOfInputVectors() const {
-  if( mybasemulticolvars.size()==0 ) return 2;
+  if( mybasemulticolvars.size()==0 )
+    return 2;
 
   unsigned nq = mybasemulticolvars[0]->getNumberOfQuantities();
   for(unsigned i=1; i<mybasemulticolvars.size(); ++i) {
-    if( mybasemulticolvars[i]->getNumberOfQuantities()!=nq ) error("mismatch between vectors in base colvars");
+    if( mybasemulticolvars[i]->getNumberOfQuantities()!=nq )
+      error("mismatch between vectors in base colvars");
   }
   return nq;
 }
 
 unsigned AdjacencyMatrixBase::getNumberOfNodeTypes() const {
   unsigned size=mybasemulticolvars.size();
-  if( size==0 ) return 1;
+  if( size==0 )
+    return 1;
   return size;
 }
 
@@ -110,57 +121,85 @@ void AdjacencyMatrixBase::retrieveTypeDimensions( unsigned& nrows, unsigned& nco
   bool allsame=(ablocks[0].size()==ablocks[1].size());
   if( allsame ) {
     for(unsigned i=0; i<ablocks[0].size(); ++i) {
-      if( ablocks[0][i]!=ablocks[1][i] ) allsame=false;
+      if( ablocks[0][i]!=ablocks[1][i] )
+        allsame=false;
     }
   }
 
   if( allsame ) {
-    std::vector<unsigned> types(1); types[0]=atom_lab[ablocks[0][0]].first;
+    //this could directly be intialized like std::vector<unsigned> types={atom_lab[ablocks[0][0]].first}; faster, and clearer
+    std::vector<unsigned> types(1);
+    types[0]=atom_lab[ablocks[0][0]].first;
     for(unsigned i=1; i<ablocks[0].size(); ++i) {
       bool found = false;
       for(unsigned j=0; j<types.size(); ++j) {
-        if( atom_lab[ablocks[0][i]].first==types[j] ) { found=true; break; }
+        if( atom_lab[ablocks[0][i]].first==types[j] ) {
+          found=true;
+          break;
+        }
       }
-      if( !found ) types.push_back( atom_lab[ablocks[0][i]].first );
+      if( !found )
+        types.push_back( atom_lab[ablocks[0][i]].first );
     }
-    ntype=0; nrows=ncols=types.size();
+    ntype=0;
+    nrows=ncols=types.size();
   } else {
-    std::vector<unsigned> types(1); types[0]=atom_lab[ablocks[0][0]].first;
+    std::vector<unsigned> types(1);
+    types[0]=atom_lab[ablocks[0][0]].first;
     for(unsigned i=1; i<ablocks[0].size(); ++i) {
       bool found = false;
       for(unsigned j=0; j<types.size(); ++j) {
-        if( atom_lab[ablocks[0][i]].first==types[j] ) { found=true; break; }
+        if( atom_lab[ablocks[0][i]].first==types[j] ) {
+          found=true;
+          break;
+        }
       }
-      if( !found ) types.push_back( atom_lab[ablocks[0][i]].first );
+      if( !found )
+        types.push_back( atom_lab[ablocks[0][i]].first );
     }
     nrows=ntype=types.size();
     for(unsigned i=0; i<ablocks[1].size(); ++i) {
       bool found = false;
       for(unsigned j=0; j<types.size(); ++j) {
-        if( atom_lab[ablocks[1][i]].first==types[j] ) { found=true; break; }
+        if( atom_lab[ablocks[1][i]].first==types[j] ) {
+          found=true;
+          break;
+        }
       }
-      if( !found ) types.push_back( atom_lab[ablocks[1][i]].first );
+      if( !found )
+        types.push_back( atom_lab[ablocks[1][i]].first );
     }
-    if( types.size()==nrows ) { ntype=0; ncols=1; plumed_assert( types.size()==1 && atom_lab[ablocks[0][0]].first==0 ); }
-    else ncols = types.size() - ntype;
+    if( types.size()==nrows ) {
+      ntype=0;
+      ncols=1;
+      plumed_assert( types.size()==1 && atom_lab[ablocks[0][0]].first==0 );
+    }
+    else
+      ncols = types.size() - ntype;
   }
 }
 
 void AdjacencyMatrixBase::finishMatrixSetup( const bool& symmetric, const std::vector<AtomNumber>& all_atoms ) {
   std::string param;
-  if( symmetric && ablocks[0].size()==ablocks[1].size() ) param="SYMMETRIC";
+  if( symmetric && ablocks[0].size()==ablocks[1].size() )
+    param="SYMMETRIC";
   if( !symmetric ) {
     bool usehbonds=( ablocks[0].size()==ablocks[1].size() );
     if( usehbonds ) {
       for(unsigned i=0; i<ablocks[0].size(); ++i) {
-        if( ablocks[0][i]!=ablocks[1][i] ) { usehbonds = false; break; }
+        if( ablocks[0][i]!=ablocks[1][i] ) {
+          usehbonds = false;
+          break;
+        }
       }
-      if( usehbonds ) param="HBONDS";
+      if( usehbonds )
+        param="HBONDS";
     }
   }
 
   vesselbase::VesselOptions da("","",0,param,this);
-  Keywords keys; AdjacencyMatrixVessel::registerKeywords( keys );
+  Keywords keys;
+  AdjacencyMatrixVessel::registerKeywords( keys );
   vesselbase::VesselOptions da2(da,keys);
   auto ves=Tools::make_unique<AdjacencyMatrixVessel>(da2);
   addVessel( std::move( ves ) );
@@ -168,12 +207,14 @@ void AdjacencyMatrixBase::finishMatrixSetup( const bool& symmetric, const std::v
 }
 
 void AdjacencyMatrixBase::readMaxTwoSpeciesMatrix( const std::string& key0, const std::string& key1, const std::string& key2, const bool& symmetric ) {
-  std::vector<AtomNumber> all_atoms; readTwoGroups( key0, key1, key2, all_atoms );
+  std::vector<AtomNumber> all_atoms;
+  readTwoGroups( key0, key1, key2, all_atoms );
   finishMatrixSetup( symmetric, all_atoms );
 }
 
 void AdjacencyMatrixBase::readMaxThreeSpeciesMatrix( const std::string& key0, const std::string& key1, const std::string& key2, const std::string& keym, const bool& symmetric ) {
-  std::vector<AtomNumber> all_atoms; readGroupKeywords( key0, key1, key2, keym, true, symmetric, all_atoms );
+  std::vector<AtomNumber> all_atoms;
+  readGroupKeywords( key0, key1, key2, keym, true, symmetric, all_atoms );
   finishMatrixSetup( symmetric, all_atoms );
 }
 
@@ -185,10 +226,12 @@ void AdjacencyMatrixBase::readMaxThreeSpeciesMatrix( const std::string& key0, co
 // }
 
 void AdjacencyMatrixBase::recalculateMatrixElement( const unsigned& myelem, MultiValue& myvals ) {
-  std::vector<unsigned> myatoms; decodeIndexToAtoms( getTaskCode( myelem ), myatoms );
+  std::vector<unsigned> myatoms;
+  decodeIndexToAtoms( getTaskCode( myelem ), myatoms );
   unsigned i=myatoms[0], j=myatoms[1];
   for(unsigned k=bookeeping(i,j).first; k<bookeeping(i,j).second; ++k) {
-    if( !taskIsCurrentlyActive(k) ) continue;
+    if( !taskIsCurrentlyActive(k) )
+      continue;
     performTask( k, getTaskCode(k), myvals );  // This may not accumulate as we would like  GAT
   }
 }
diff --git a/src/adjmat/AdjacencyMatrixBase.h b/src/adjmat/AdjacencyMatrixBase.h
index 6efff8a50f..293ebb6900 100644
--- a/src/adjmat/AdjacencyMatrixBase.h
+++ b/src/adjmat/AdjacencyMatrixBase.h
@@ -73,8 +73,12 @@ class AdjacencyMatrixBase : public multicolvar::MultiColvarBase {
 /// Create the connection object
   virtual void setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) = 0;
 /// None of these things are allowed
-  bool isPeriodic() { return false; }
-  Vector getCentralAtom() { plumed_merror("cannot find central atoms for adjacency matrix actions"); }
+  bool isPeriodic() {
+    return false;
+  }
+  Vector getCentralAtom() {
+    plumed_merror("cannot find central atoms for adjacency matrix actions");
+  }
 /// Get the atom number
   AtomNumber getAbsoluteIndexOfCentralAtom( const unsigned& i ) const ;
 };
@@ -86,7 +90,8 @@ AdjacencyMatrixVessel* AdjacencyMatrixBase::getAdjacencyVessel() {
 
 inline
 unsigned AdjacencyMatrixBase::getBaseColvarNumber( const unsigned& inum ) const {
-  if( atom_lab[inum].first>0 ) return atom_lab[inum].first-1;
+  if( atom_lab[inum].first>0 )
+    return atom_lab[inum].first-1;
   return 0;
 }
 
diff --git a/src/adjmat/AdjacencyMatrixVessel.cpp b/src/adjmat/AdjacencyMatrixVessel.cpp
index b4136f0f53..fd3b581f85 100644
--- a/src/adjmat/AdjacencyMatrixVessel.cpp
+++ b/src/adjmat/AdjacencyMatrixVessel.cpp
@@ -37,10 +37,14 @@ AdjacencyMatrixVessel::AdjacencyMatrixVessel( const vesselbase::VesselOptions& d
 {
   function=dynamic_cast<AdjacencyMatrixBase*>( getAction() );
   plumed_assert( function );
-  parseFlag("SYMMETRIC",symmetric); parseFlag("HBONDS",hbonds);
-  if( symmetric && hbonds ) error("matrix should be either symmetric or hbonds");
-  if( symmetric && function->ablocks[0].size()!=function->ablocks[1].size() ) error("matrix is supposed to be symmetric but nrows!=ncols");
-  if( hbonds &&  function->ablocks[0].size()!=function->ablocks[1].size() ) error("matrix is supposed to be hbonds but nrows!=ncols");
+  parseFlag("SYMMETRIC",symmetric);
+  parseFlag("HBONDS",hbonds);
+  if( symmetric && hbonds )
+    error("matrix should be either symmetric or hbonds");
+  if( symmetric && function->ablocks[0].size()!=function->ablocks[1].size() )
+    error("matrix is supposed to be symmetric but nrows!=ncols");
+  if( hbonds &&  function->ablocks[0].size()!=function->ablocks[1].size() )
+    error("matrix is supposed to be hbonds but nrows!=ncols");
 }
 
 bool AdjacencyMatrixVessel::isSymmetric() const {
@@ -64,13 +68,16 @@ bool AdjacencyMatrixVessel::matrixElementIsActive( const unsigned& ielem, const
 }
 
 unsigned AdjacencyMatrixVessel::getStoreIndexFromMatrixIndices( const unsigned& ielem, const unsigned& jelem ) const {
-  if( !symmetric && !hbonds ) return (function->ablocks[1].size())*ielem + jelem;
+  if( !symmetric && !hbonds )
+    return (function->ablocks[1].size())*ielem + jelem;
   if( !symmetric ) {
     plumed_dbg_assert( ielem!=jelem );
-    if( jelem<ielem ) return (function->ablocks[1].size()-1)*ielem + jelem;
+    if( jelem<ielem )
+      return (function->ablocks[1].size()-1)*ielem + jelem;
     return (function->ablocks[1].size()-1)*ielem + jelem - 1;
   }
-  if( ielem>jelem ) return 0.5*ielem*(ielem-1)+jelem;
+  if( ielem>jelem )
+    return 0.5*ielem*(ielem-1)+jelem;
   return 0.5*jelem*(jelem-1) + ielem;
 }
 
@@ -79,22 +86,32 @@ AdjacencyMatrixBase* AdjacencyMatrixVessel::getMatrixAction() {
 }
 
 void AdjacencyMatrixVessel::getMatrixIndices( const unsigned& code, unsigned& i, unsigned& j ) const {
-  std::vector<unsigned> myatoms; function->decodeIndexToAtoms( function->getTaskCode(code), myatoms );
-  i=myatoms[0]; j=myatoms[1];
-  if( !undirectedGraph() ) j -= function->ablocks[0].size(); // Have to remove number of columns as returns number in ablocks[1]
+  std::vector<unsigned> myatoms;
+  function->decodeIndexToAtoms( function->getTaskCode(code), myatoms );
+  i=myatoms[0];
+  j=myatoms[1];
+  if( !undirectedGraph() ) {
+    // Have to remove number of columns as returns number in ablocks[1]
+    j -= function->ablocks[0].size();
+  }
 }
 
 void AdjacencyMatrixVessel::retrieveMatrix( DynamicList<unsigned>& myactive_elements, Matrix<double>& mymatrix ) {
-  myactive_elements.deactivateAll(); std::vector<double> vals( getNumberOfComponents() );
+  myactive_elements.deactivateAll();
+  std::vector<double> vals( getNumberOfComponents() );
   for(unsigned i=0; i<getNumberOfStoredValues(); ++i) {
     retrieveSequentialValue( i, false, vals );
-    if( vals[0]<epsilon ) continue ;
+    if( vals[0]<epsilon )
+      continue ;
 
     myactive_elements.activate(i);
-    unsigned j, k; getMatrixIndices( function->getPositionInFullTaskList(i), k, j );
+    unsigned j, k;
+    getMatrixIndices( function->getPositionInFullTaskList(i), k, j );
 
-    if( symmetric ) mymatrix(k,j)=mymatrix(j,k)=vals[0]*vals[1];
-    else mymatrix(k,j)=vals[0]*vals[1];
+    if( symmetric )
+      mymatrix(k,j)=mymatrix(j,k)=vals[0]*vals[1];
+    else
+      mymatrix(k,j)=vals[0]*vals[1];
   }
   myactive_elements.updateActiveMembers();
 }
@@ -102,34 +119,43 @@ void AdjacencyMatrixVessel::retrieveMatrix( DynamicList<unsigned>& myactive_elem
 void AdjacencyMatrixVessel::retrieveAdjacencyLists( std::vector<unsigned>& nneigh, Matrix<unsigned>& adj_list ) {
   plumed_dbg_assert( undirectedGraph() );
   // Currently everything has zero neighbors
-  for(unsigned i=0; i<nneigh.size(); ++i) nneigh[i]=0;
+  for(unsigned i=0; i<nneigh.size(); ++i)
+    nneigh[i]=0;
 
   // And set up the adjacency list
   std::vector<double> myvals( getNumberOfComponents() );
   for(unsigned i=0; i<getNumberOfStoredValues(); ++i) {
     // Check if atoms are connected
     retrieveSequentialValue( i, false, myvals );
-    if( myvals[0]<epsilon || myvals[1]<epsilon ) continue ;
+    if( myvals[0]<epsilon || myvals[1]<epsilon )
+      continue ;
 
-    unsigned j, k; getMatrixIndices( function->getPositionInFullTaskList(i), k, j );
+    unsigned j, k;
+    getMatrixIndices( function->getPositionInFullTaskList(i), k, j );
 
-    if( nneigh[j]>=adj_list.ncols() || nneigh[k]>=adj_list.ncols() ) error("adjacency lists are not large enough, increase maxconnections");
+    if( nneigh[j]>=adj_list.ncols() || nneigh[k]>=adj_list.ncols() )
+      error("adjacency lists are not large enough, increase maxconnections");
     // Store if atoms are connected
     // unsigned j, k; getMatrixIndices( i, k, j );
-    adj_list(k,nneigh[k])=j; nneigh[k]++;
-    adj_list(j,nneigh[j])=k; nneigh[j]++;
+    adj_list(k,nneigh[k])=j;
+    nneigh[k]++;
+    adj_list(j,nneigh[j])=k;
+    nneigh[j]++;
   }
 }
 
 void AdjacencyMatrixVessel::retrieveEdgeList( unsigned& nedge, std::vector<std::pair<unsigned,unsigned> >& edge_list ) {
-  plumed_dbg_assert( undirectedGraph() ); nedge=0;
+  plumed_dbg_assert( undirectedGraph() );
+  nedge=0;
   std::vector<double> myvals( getNumberOfComponents() );
-  if( getNumberOfStoredValues()>edge_list.size() ) error("adjacency lists are not large enough, increase maxconnections");
+  if( getNumberOfStoredValues()>edge_list.size() )
+    error("adjacency lists are not large enough, increase maxconnections");
 
   for(unsigned i=0; i<getNumberOfStoredValues(); ++i) {
     // Check if atoms are connected
     retrieveSequentialValue( i, false, myvals );
-    if( myvals[0]<epsilon || myvals[1]<epsilon ) continue ;
+    if( myvals[0]<epsilon || myvals[1]<epsilon )
+      continue ;
 
     getMatrixIndices( function->getPositionInFullTaskList(i), edge_list[nedge].first, edge_list[nedge].second );
     nedge++;
@@ -137,7 +163,8 @@ void AdjacencyMatrixVessel::retrieveEdgeList( unsigned& nedge, std::vector<std::
 }
 
 bool AdjacencyMatrixVessel::nodesAreConnected( const unsigned& iatom, const unsigned& jatom ) const {
-  if( !matrixElementIsActive( iatom, jatom ) ) return false;
+  if( !matrixElementIsActive( iatom, jatom ) )
+    return false;
   unsigned ind=getStoreIndexFromMatrixIndices( iatom, jatom );
 
   std::vector<double> myvals( getNumberOfComponents() );
diff --git a/src/adjmat/AlignedMatrixBase.cpp b/src/adjmat/AlignedMatrixBase.cpp
index 7d865e964d..ba8596f441 100644
--- a/src/adjmat/AlignedMatrixBase.cpp
+++ b/src/adjmat/AlignedMatrixBase.cpp
@@ -51,9 +51,12 @@ AlignedMatrixBase::AlignedMatrixBase( const ActionOptions& ao ):
 {
   // Read in the atomic positions
   readMaxTwoSpeciesMatrix( "ATOMS", "ATOMSA", "ATOMSB", true );
-  unsigned nrows, ncols; retrieveTypeDimensions( nrows, ncols, ncol_t );
-  if( mybasemulticolvars.size()==0 ) error("cannot use atom indices as input to this variable / input was not specified");
-  if( getSizeOfInputVectors()<3 ) error("base multicolvars do not calculate an orientation");
+  unsigned nrows, ncols;
+  retrieveTypeDimensions( nrows, ncols, ncol_t );
+  if( mybasemulticolvars.size()==0 )
+    error("cannot use atom indices as input to this variable / input was not specified");
+  if( getSizeOfInputVectors()<3 )
+    error("base multicolvars do not calculate an orientation");
   // Read in the switching function
   switchingFunction.resize( nrows, ncols );
   parseConnectionDescriptions("SWITCH",false,ncol_t);
@@ -63,7 +66,8 @@ AlignedMatrixBase::AlignedMatrixBase( const ActionOptions& ao ):
   for(unsigned i=0; i<getNumberOfNodeTypes(); ++i) {
     for(unsigned j=0; j<getNumberOfNodeTypes(); ++j) {
       double tsf=switchingFunction(i,j).get_dmax();
-      if( tsf>sfmax ) sfmax=tsf;
+      if( tsf>sfmax )
+        sfmax=tsf;
     }
   }
   // And set the link cell cutoff
@@ -73,9 +77,13 @@ AlignedMatrixBase::AlignedMatrixBase( const ActionOptions& ao ):
 void AlignedMatrixBase::setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) {
   plumed_assert( id<2 );
   if( id==0 ) {
-    plumed_assert( desc.size()==1 ); std::string errors; switchingFunction(j,i).set(desc[0],errors);
-    if( errors.length()!=0 ) error("problem reading switching function in SWITCH keywrd description " + errors);
-    if( j!=i) switchingFunction(i,j).set(desc[0],errors);
+    plumed_assert( desc.size()==1 );
+    std::string errors;
+    switchingFunction(j,i).set(desc[0],errors);
+    if( errors.length()!=0 )
+      error("problem reading switching function in SWITCH keywrd description " + errors);
+    if( j!=i)
+      switchingFunction(i,j).set(desc[0],errors);
     log.printf("  %u th and %u th multicolvar groups must be within %s\n",i+1,j+1,(switchingFunction(i,j).description()).c_str() );
   } else if( id==1 ) {
     readOrientationConnector( i, j, desc );
@@ -84,15 +92,18 @@ void AlignedMatrixBase::setupConnector( const unsigned& id, const unsigned& i, c
 
 double AlignedMatrixBase::calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
   Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
-  if( distance.modulo2()<switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) )-ncol_t ).get_dmax2() ) return 1.0;
+  if( distance.modulo2()<switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) )-ncol_t ).get_dmax2() )
+    return 1.0;
   return 0.0;
 }
 
 double AlignedMatrixBase::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
-  unsigned ncomp=getSizeOfInputVectors(); Vector ddistance;
+  unsigned ncomp=getSizeOfInputVectors();
+  Vector ddistance;
   std::vector<double> orient0(ncomp), orient1(ncomp), dorient0(ncomp), dorient1(ncomp);
   Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
-  getInputData( 0, true, myatoms, orient0 ); getInputData( 1, true, myatoms, orient1 );
+  getInputData( 0, true, myatoms, orient0 );
+  getInputData( 1, true, myatoms, orient1 );
   double f_dot = computeVectorFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) )-ncol_t,
                                         distance, orient0, orient1, ddistance, dorient0, dorient1 );
 
@@ -105,7 +116,10 @@ double AlignedMatrixBase::compute( const unsigned& tindex, multicolvar::AtomValu
     myatoms.addBoxDerivatives( 1, (-dfunc)*f_dot*Tensor(distance,distance) - sw*extProduct( distance, ddistance ) );
 
     // Add derivatives of orientation
-    for(unsigned k=2; k<orient0.size(); ++k) { dorient0[k]*=sw; dorient1[k]*=sw; }
+    for(unsigned k=2; k<orient0.size(); ++k) {
+      dorient0[k]*=sw;
+      dorient1[k]*=sw;
+    }
     mergeInputDerivatives( 1, 2, orient0.size(), 0, dorient0, getInputDerivatives( 0, true, myatoms ), myatoms );
     mergeInputDerivatives( 1, 2, orient1.size(), 1, dorient1, getInputDerivatives( 1, true, myatoms ), myatoms );
   }
diff --git a/src/adjmat/ClusterAnalysisBase.cpp b/src/adjmat/ClusterAnalysisBase.cpp
index f8f3ee959d..f6ef4d4283 100644
--- a/src/adjmat/ClusterAnalysisBase.cpp
+++ b/src/adjmat/ClusterAnalysisBase.cpp
@@ -37,23 +37,31 @@ ClusterAnalysisBase::ClusterAnalysisBase(const ActionOptions& ao):
   myclusters(NULL)
 {
   // This makes these colvars behave appropriately with dump and analysis
-  matsums=usespecies=true; std::vector<AtomNumber> fake_atoms;
+  matsums=usespecies=true;
+  std::vector<AtomNumber> fake_atoms;
   // Find what action we are taking the clusters from
-  if( !parseMultiColvarAtomList("CLUSTERS",-1,fake_atoms ) ) error("unable to interpret input CLUSTERS" );
-  if( mybasemulticolvars.size()!=1 ) error("should be exactly one multicolvar input");
-  atom_lab.resize(0); myclusters = dynamic_cast<ClusteringBase*>( mybasemulticolvars[0] );
-  if( !myclusters ) error("input label is not that of a DFS object");
+  if( !parseMultiColvarAtomList("CLUSTERS",-1,fake_atoms ) )
+    error("unable to interpret input CLUSTERS" );
+  if( mybasemulticolvars.size()!=1 )
+    error("should be exactly one multicolvar input");
+  atom_lab.resize(0);
+  myclusters = dynamic_cast<ClusteringBase*>( mybasemulticolvars[0] );
+  if( !myclusters )
+    error("input label is not that of a DFS object");
   // Setup the atom pack
   myfatoms.setNumberOfAtoms( myclusters->getNumberOfNodes() );
   myfvals.getIndices().resize( myclusters->getNumberOfNodes() );
-  for(unsigned i=0; i<myclusters->getNumberOfNodes(); ++i) myfatoms.setAtomIndex( i, i );
+  for(unsigned i=0; i<myclusters->getNumberOfNodes(); ++i)
+    myfatoms.setAtomIndex( i, i );
 }
 
 void ClusterAnalysisBase::turnOnDerivatives() {
   // Check for dubious vessels
   for(unsigned i=0; i<getNumberOfVessels(); ++i) {
-    if( getPntrToVessel(i)->getName()=="MEAN" ) error("MEAN of cluster is not differentiable");
-    if( getPntrToVessel(i)->getName()=="VMEAN" ) error("VMEAN of cluster is not differentiable");
+    if( getPntrToVessel(i)->getName()=="MEAN" )
+      error("MEAN of cluster is not differentiable");
+    if( getPntrToVessel(i)->getName()=="VMEAN" )
+      error("VMEAN of cluster is not differentiable");
   }
   MultiColvarBase::turnOnDerivatives();
 }
diff --git a/src/adjmat/ClusterAnalysisBase.h b/src/adjmat/ClusterAnalysisBase.h
index 97593afb5d..b43cb465cd 100644
--- a/src/adjmat/ClusterAnalysisBase.h
+++ b/src/adjmat/ClusterAnalysisBase.h
@@ -50,7 +50,9 @@ class ClusterAnalysisBase : public multicolvar::MultiColvarBase {
   void turnOnDerivatives() override;
   void setupActiveTaskSet( std::vector<unsigned>& active_tasks, const std::string& input_label ) {}
   Vector getPositionOfAtomForLinkCells( const unsigned& ) const override;
-  double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const override { plumed_error(); }
+  double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const override {
+    plumed_error();
+  }
 };
 
 }
diff --git a/src/adjmat/ClusterDiameter.cpp b/src/adjmat/ClusterDiameter.cpp
index aa67dce394..cedd2b4373 100644
--- a/src/adjmat/ClusterDiameter.cpp
+++ b/src/adjmat/ClusterDiameter.cpp
@@ -91,15 +91,20 @@ ClusterDiameter::ClusterDiameter(const ActionOptions&ao):
   // Find out which cluster we want
   parse("CLUSTER",clustr);
 
-  if( clustr<1 ) error("cannot look for a cluster larger than the largest cluster");
-  if( clustr>getNumberOfNodes() ) error("cluster selected is invalid - too few atoms in system");
+  if( clustr<1 )
+    error("cannot look for a cluster larger than the largest cluster");
+  if( clustr>getNumberOfNodes() )
+    error("cluster selected is invalid - too few atoms in system");
 
   // Create the task list
   for(unsigned  i=0; i<getNumberOfNodes(); ++i) {
-    for(unsigned j=0; j<getNumberOfNodes(); ++j) addTaskToList( i*getNumberOfNodes() + j );
+    for(unsigned j=0; j<getNumberOfNodes(); ++j)
+      addTaskToList( i*getNumberOfNodes() + j );
   }
   // Now create a highest vessel
-  addVessel("HIGHEST", "", -1); std::vector<AtomNumber> fake_atoms; setupMultiColvarBase( fake_atoms );
+  addVessel("HIGHEST", "", -1);
+  std::vector<AtomNumber> fake_atoms;
+  setupMultiColvarBase( fake_atoms );
 }
 
 void ClusterDiameter::turnOnDerivatives() {
@@ -108,11 +113,13 @@ void ClusterDiameter::turnOnDerivatives() {
 
 void ClusterDiameter::calculate() {
   // Retrieve the atoms in the largest cluster
-  std::vector<unsigned> myatoms; retrieveAtomsInCluster( clustr, myatoms );
+  std::vector<unsigned> myatoms;
+  retrieveAtomsInCluster( clustr, myatoms );
   // Activate the relevant tasks
   deactivateAllTasks();
   for(unsigned i=1; i<myatoms.size(); ++i) {
-    for(unsigned j=0; j<i; ++j) taskFlags[ myatoms[i]*getNumberOfNodes() + myatoms[j] ] = 1;
+    for(unsigned j=0; j<i; ++j)
+      taskFlags[ myatoms[i]*getNumberOfNodes() + myatoms[j] ] = 1;
   }
   lockContributors();
   // Now do the calculation
@@ -123,7 +130,8 @@ void ClusterDiameter::performTask( const unsigned& task_index, const unsigned& c
   unsigned iatom=std::floor(current/getNumberOfNodes()), jatom = current - iatom*getNumberOfNodes();
   Vector distance=getSeparation( getPosition(iatom), getPosition(jatom) );
   double dd = distance.modulo();
-  myvals.setValue( 0, 1.0 ); myvals.setValue( 1, dd );
+  myvals.setValue( 0, 1.0 );
+  myvals.setValue( 1, dd );
 }
 
 }
diff --git a/src/adjmat/ClusterDistribution.cpp b/src/adjmat/ClusterDistribution.cpp
index 5733eb8405..a6e0f44620 100644
--- a/src/adjmat/ClusterDistribution.cpp
+++ b/src/adjmat/ClusterDistribution.cpp
@@ -87,8 +87,13 @@ void ClusterDistribution::registerKeywords( Keywords& keys ) {
   keys.add("compulsory","TRANSFORM","none","the switching function to use to convert the crystallinity parameter to a number between zero and one");
   keys.addFlag("INVERSE_TRANSFORM",false,"when TRANSFORM appears alone the input symmetry functions, \\f$x\\f$ are transformed used \\f$1-s(x)\\f$ "
                "where \\f$s(x)\\f$ is a switching function.  When this option is used you instead transform using \\f$s(x)\\f$ only.");
-  keys.use("MORE_THAN"); keys.use("LESS_THAN"); keys.use("BETWEEN");
-  keys.use("HISTOGRAM"); keys.use("ALT_MIN"); keys.use("MIN"); keys.use("MAX");
+  keys.use("MORE_THAN");
+  keys.use("LESS_THAN");
+  keys.use("BETWEEN");
+  keys.use("HISTOGRAM");
+  keys.use("ALT_MIN");
+  keys.use("MIN");
+  keys.use("MAX");
 }
 
 ClusterDistribution::ClusterDistribution(const ActionOptions&ao):
@@ -97,35 +102,44 @@ ClusterDistribution::ClusterDistribution(const ActionOptions&ao):
   nderivatives(0)
 {
   use_switch=false;
-  std::string input, errors; parse("TRANSFORM",input);
+  std::string input, errors;
+  parse("TRANSFORM",input);
   if( input!="none" ) {
-    use_switch=true; sf.set( input, errors );
-    if( errors.length()!=0 ) error("problem reading SWITCH keyword : " + errors );
+    use_switch=true;
+    sf.set( input, errors );
+    if( errors.length()!=0 )
+      error("problem reading SWITCH keyword : " + errors );
   }
   parseFlag("INVERSE_TRANSFORM",inverse);
-  if( inverse && !use_switch ) error("INVERSE_TRANSFORM option was specified but no TRANSOFRM switching function was given");
+  if( inverse && !use_switch )
+    error("INVERSE_TRANSFORM option was specified but no TRANSOFRM switching function was given");
 
   // Create all tasks by copying those from underlying DFS object (which is actually MultiColvar)
-  for(unsigned i=0; i<getNumberOfNodes(); ++i) addTaskToList(i);
+  for(unsigned i=0; i<getNumberOfNodes(); ++i)
+    addTaskToList(i);
 
   // And now finish the setup of everything in the base
-  std::vector<AtomNumber> fake_atoms; setupMultiColvarBase( fake_atoms );
+  std::vector<AtomNumber> fake_atoms;
+  setupMultiColvarBase( fake_atoms );
 }
 
 void ClusterDistribution::calculate() {
   // Activate the relevant tasks
   nderivatives = getNumberOfDerivatives();
   deactivateAllTasks();
-  for(unsigned i=0; i<getNumberOfClusters(); ++i) taskFlags[i]=1;
+  for(unsigned i=0; i<getNumberOfClusters(); ++i)
+    taskFlags[i]=1;
   lockContributors();
   // Now do the calculation
   runAllTasks();
 }
 
 void ClusterDistribution::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
-  std::vector<unsigned> myatoms; retrieveAtomsInCluster( current+1, myatoms );
+  std::vector<unsigned> myatoms;
+  retrieveAtomsInCluster( current+1, myatoms );
   // This deals with filters
-  if( myatoms.size()==1 && !nodeIsActive(myatoms[0]) ) return ;
+  if( myatoms.size()==1 && !nodeIsActive(myatoms[0]) )
+    return ;
 
   std::vector<double> vals( getNumberOfQuantities() );
   MultiValue tvals( getNumberOfQuantities(), nderivatives );
@@ -137,12 +151,16 @@ void ClusterDistribution::performTask( const unsigned& task_index, const unsigne
     getPropertiesOfNode( i, vals );
     if( use_switch && !inverse ) {
       vv = 1.0 - sf.calculate( vals[1], df );
-      tval += vals[0]*vv; df=-df*vals[1];
+      tval += vals[0]*vv;
+      df=-df*vals[1];
     } else if( use_switch ) {
       vv = sf.calculate( vals[1], df );
-      tval += vals[0]*vv; df=df*vals[1];
+      tval += vals[0]*vv;
+      df=df*vals[1];
     } else {
-      tval += vals[0]*vals[1]; df=1.; vv=vals[1];
+      tval += vals[0]*vals[1];
+      df=1.;
+      vv=vals[1];
     }
     if( !doNotCalculateDerivatives() ) {
       getNodePropertyDerivatives( i, tvals );
@@ -153,7 +171,8 @@ void ClusterDistribution::performTask( const unsigned& task_index, const unsigne
       tvals.clearAll();
     }
   }
-  myvals.setValue( 0, 1.0 ); myvals.setValue( 1, tval );
+  myvals.setValue( 0, 1.0 );
+  myvals.setValue( 1, tval );
 }
 
 }
diff --git a/src/adjmat/ClusterProperties.cpp b/src/adjmat/ClusterProperties.cpp
index 5739fd443b..860d20da21 100644
--- a/src/adjmat/ClusterProperties.cpp
+++ b/src/adjmat/ClusterProperties.cpp
@@ -76,11 +76,22 @@ PLUMED_REGISTER_ACTION(ClusterProperties,"CLUSTER_PROPERTIES")
 void ClusterProperties::registerKeywords( Keywords& keys ) {
   ClusterAnalysisBase::registerKeywords( keys );
   keys.add("compulsory","CLUSTER","1","which cluster would you like to look at 1 is the largest cluster, 2 is the second largest, 3 is the the third largest and so on.");
-  keys.use("MEAN"); keys.use("MORE_THAN"); keys.use("LESS_THAN");
-  if( keys.reserved("VMEAN") ) keys.use("VMEAN");
-  if( keys.reserved("VSUM") ) keys.use("VSUM");
-  keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS"); keys.use("ALT_MIN");
-  keys.use("MIN"); keys.use("MAX"); keys.use("SUM"); keys.use("LOWEST"); keys.use("HIGHEST");
+  keys.use("MEAN");
+  keys.use("MORE_THAN");
+  keys.use("LESS_THAN");
+  if( keys.reserved("VMEAN") )
+    keys.use("VMEAN");
+  if( keys.reserved("VSUM") )
+    keys.use("VSUM");
+  keys.use("BETWEEN");
+  keys.use("HISTOGRAM");
+  keys.use("MOMENTS");
+  keys.use("ALT_MIN");
+  keys.use("MIN");
+  keys.use("MAX");
+  keys.use("SUM");
+  keys.use("LOWEST");
+  keys.use("HIGHEST");
 }
 
 ClusterProperties::ClusterProperties(const ActionOptions&ao):
@@ -90,31 +101,40 @@ ClusterProperties::ClusterProperties(const ActionOptions&ao):
   // Find out which cluster we want
   parse("CLUSTER",clustr);
 
-  if( clustr<1 ) error("cannot look for a cluster larger than the largest cluster");
-  if( clustr>getNumberOfNodes() ) error("cluster selected is invalid - too few atoms in system");
+  if( clustr<1 )
+    error("cannot look for a cluster larger than the largest cluster");
+  if( clustr>getNumberOfNodes() )
+    error("cluster selected is invalid - too few atoms in system");
 
   // Create all tasks by copying those from underlying DFS object (which is actually MultiColvar)
-  for(unsigned i=0; i<getNumberOfNodes(); ++i) addTaskToList(i);
+  for(unsigned i=0; i<getNumberOfNodes(); ++i)
+    addTaskToList(i);
 
   // And now finish the setup of everything in the base
-  std::vector<AtomNumber> fake_atoms; setupMultiColvarBase( fake_atoms );
+  std::vector<AtomNumber> fake_atoms;
+  setupMultiColvarBase( fake_atoms );
 }
 
 void ClusterProperties::calculate() {
   // Retrieve the atoms in the largest cluster
-  std::vector<unsigned> myatoms; retrieveAtomsInCluster( clustr, myatoms );
+  std::vector<unsigned> myatoms;
+  retrieveAtomsInCluster( clustr, myatoms );
   // Activate the relevant tasks
   deactivateAllTasks();
-  for(unsigned i=0; i<myatoms.size(); ++i) taskFlags[myatoms[i]]=1;
+  for(unsigned i=0; i<myatoms.size(); ++i)
+    taskFlags[myatoms[i]]=1;
   lockContributors();
   // Now do the calculation
   runAllTasks();
 }
 
 void ClusterProperties::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
-  std::vector<double> vals( myvals.getNumberOfValues() ); getPropertiesOfNode( current, vals );
-  if( !doNotCalculateDerivatives() ) getNodePropertyDerivatives( current, myvals );
-  for(unsigned k=0; k<vals.size(); ++k) myvals.setValue( k, vals[k] );
+  std::vector<double> vals( myvals.getNumberOfValues() );
+  getPropertiesOfNode( current, vals );
+  if( !doNotCalculateDerivatives() )
+    getNodePropertyDerivatives( current, myvals );
+  for(unsigned k=0; k<vals.size(); ++k)
+    myvals.setValue( k, vals[k] );
 }
 
 }
diff --git a/src/adjmat/ClusterSize.cpp b/src/adjmat/ClusterSize.cpp
index 1f69058325..df577d2a7f 100644
--- a/src/adjmat/ClusterSize.cpp
+++ b/src/adjmat/ClusterSize.cpp
@@ -71,7 +71,9 @@ class ClusterSize : public ClusterAnalysisBase {
 ///
   void calculate() override;
 ///
-  void performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const override { plumed_error(); }
+  void performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const override {
+    plumed_error();
+  }
 ///
   void turnOnDerivatives() override;
 };
@@ -90,14 +92,19 @@ ClusterSize::ClusterSize(const ActionOptions&ao):
   // Find out which cluster we want
   parse("CLUSTER",clustr);
 
-  if( clustr<1 ) error("cannot look for a cluster larger than the largest cluster");
-  if( clustr>getNumberOfNodes() ) error("cluster selected is invalid - too few atoms in system");
+  if( clustr<1 )
+    error("cannot look for a cluster larger than the largest cluster");
+  if( clustr>getNumberOfNodes() )
+    error("cluster selected is invalid - too few atoms in system");
 
   // Create all tasks by copying those from underlying DFS object (which is actually MultiColvar)
-  for(unsigned i=0; i<getNumberOfNodes(); ++i) addTaskToList(i);
+  for(unsigned i=0; i<getNumberOfNodes(); ++i)
+    addTaskToList(i);
   // And now finish the setup of everything in the base
-  std::vector<AtomNumber> fake_atoms; setupMultiColvarBase( fake_atoms );
-  addValue(); setNotPeriodic();
+  std::vector<AtomNumber> fake_atoms;
+  setupMultiColvarBase( fake_atoms );
+  addValue();
+  setNotPeriodic();
 }
 
 void ClusterSize::turnOnDerivatives() {
@@ -106,7 +113,9 @@ void ClusterSize::turnOnDerivatives() {
 
 void ClusterSize::calculate() {
   // Retrieve the atoms in the largest cluster
-  std::vector<unsigned> myatoms; retrieveAtomsInCluster( clustr, myatoms ); setValue( myatoms.size() );
+  std::vector<unsigned> myatoms;
+  retrieveAtomsInCluster( clustr, myatoms );
+  setValue( myatoms.size() );
 }
 
 }
diff --git a/src/adjmat/ClusterWithSurface.cpp b/src/adjmat/ClusterWithSurface.cpp
index 780dd7ab95..5e1e6b663e 100644
--- a/src/adjmat/ClusterWithSurface.cpp
+++ b/src/adjmat/ClusterWithSurface.cpp
@@ -110,16 +110,22 @@ ClusterWithSurface::ClusterWithSurface(const ActionOptions&ao):
   ClusteringBase(ao)
 {
   std::vector<AtomNumber> fake_atoms;
-  if( !parseMultiColvarAtomList("CLUSTERS",-1,fake_atoms ) ) error("unable to find CLUSTERS input");
-  if( mybasemulticolvars.size()!=1 ) error("should be exactly one multicolvar input");
+  if( !parseMultiColvarAtomList("CLUSTERS",-1,fake_atoms ) )
+    error("unable to find CLUSTERS input");
+  if( mybasemulticolvars.size()!=1 )
+    error("should be exactly one multicolvar input");
 
   // Retrieve the adjacency matrix of interest
-  atom_lab.resize(0); myclusters = dynamic_cast<ClusteringBase*>( mybasemulticolvars[0] );
-  if( !myclusters ) error( mybasemulticolvars[0]->getLabel() + " does not calculate clusters");
+  atom_lab.resize(0);
+  myclusters = dynamic_cast<ClusteringBase*>( mybasemulticolvars[0] );
+  if( !myclusters )
+    error( mybasemulticolvars[0]->getLabel() + " does not calculate clusters");
 
   // Setup switching function for surface atoms
-  double rcut_surf; parse("RCUT_SURF",rcut_surf);
-  if( rcut_surf>0 ) log.printf("  counting surface atoms that are within %f of the cluster atoms \n",rcut_surf);
+  double rcut_surf;
+  parse("RCUT_SURF",rcut_surf);
+  if( rcut_surf>0 )
+    log.printf("  counting surface atoms that are within %f of the cluster atoms \n",rcut_surf);
   rcut_surf2=rcut_surf*rcut_surf;
 
   // And now finish the setup of everything in the base
@@ -152,35 +158,46 @@ unsigned ClusterWithSurface::getNumberOfQuantities() const {
 
 double  ClusterWithSurface::getCutoffForConnection() const {
   double tcut = myclusters->getCutoffForConnection();
-  if( tcut>std::sqrt(rcut_surf2) ) return tcut;
+  if( tcut>std::sqrt(rcut_surf2) )
+    return tcut;
   return std::sqrt(rcut_surf2);
 }
 
 void ClusterWithSurface::retrieveAtomsInCluster( const unsigned& clust, std::vector<unsigned>& myatoms ) const {
-  std::vector<unsigned> tmpat; myclusters->retrieveAtomsInCluster( clust, tmpat );
+  std::vector<unsigned> tmpat;
+  myclusters->retrieveAtomsInCluster( clust, tmpat );
 
   // Prevent double counting
   std::vector<bool> incluster( getNumberOfNodes(), false );
-  for(unsigned i=0; i<tmpat.size(); ++i) incluster[tmpat[i]]=true;
+  for(unsigned i=0; i<tmpat.size(); ++i)
+    incluster[tmpat[i]]=true;
 
   // Find the atoms in the the clusters
   std::vector<bool> surface_atom( getNumberOfNodes(), false );
   for(unsigned i=0; i<tmpat.size(); ++i) {
     for(unsigned j=0; j<getNumberOfNodes(); ++j) {
-      if( incluster[j] ) continue;
+      if( incluster[j] )
+        continue;
       double dist2=getSeparation( getPosition(tmpat[i]), getPosition(j) ).modulo2();
-      if( dist2<rcut_surf2 ) { surface_atom[j]=true; }
+      if( dist2<rcut_surf2 ) {
+        surface_atom[j]=true;
+      }
     }
   }
   unsigned nsurf_at=0;
   for(unsigned j=0; j<getNumberOfNodes(); ++j) {
-    if( surface_atom[j] ) nsurf_at++;
+    if( surface_atom[j] )
+      nsurf_at++;
   }
   myatoms.resize( nsurf_at + tmpat.size() );
-  for(unsigned i=0; i<tmpat.size(); ++i) myatoms[i]=tmpat[i];
+  for(unsigned i=0; i<tmpat.size(); ++i)
+    myatoms[i]=tmpat[i];
   unsigned nn=tmpat.size();
   for(unsigned j=0; j<getNumberOfNodes(); ++j) {
-    if( surface_atom[j] ) { myatoms[nn]=j; nn++; }
+    if( surface_atom[j] ) {
+      myatoms[nn]=j;
+      nn++;
+    }
   }
   plumed_assert( nn==myatoms.size() );
 }
diff --git a/src/adjmat/ClusteringBase.cpp b/src/adjmat/ClusteringBase.cpp
index f4ccaba15a..afc327c1c2 100644
--- a/src/adjmat/ClusteringBase.cpp
+++ b/src/adjmat/ClusteringBase.cpp
@@ -36,19 +36,24 @@ ClusteringBase::ClusteringBase(const ActionOptions&ao):
   number_of_cluster(-1)
 {
   if( getAdjacencyVessel() ) {
-    cluster_sizes.resize(getNumberOfNodes()); which_cluster.resize(getNumberOfNodes());
-    if( getNumberOfNodeTypes()!=1 ) error("should only be running clustering with one base multicolvar in function");
-    if( !getAdjacencyVessel()->undirectedGraph() ) error("input contact matrix is incompatible with clustering");
+    cluster_sizes.resize(getNumberOfNodes());
+    which_cluster.resize(getNumberOfNodes());
+    if( getNumberOfNodeTypes()!=1 )
+      error("should only be running clustering with one base multicolvar in function");
+    if( !getAdjacencyVessel()->undirectedGraph() )
+      error("input contact matrix is incompatible with clustering");
   }
   if( keywords.exists("MATRIX") ) {
-    std::vector<AtomNumber> fake_atoms; setupMultiColvarBase( fake_atoms );
+    std::vector<AtomNumber> fake_atoms;
+    setupMultiColvarBase( fake_atoms );
   }
 }
 
 void ClusteringBase::turnOnDerivatives() {
   // Check base multicolvar isn't density probably other things shouldn't be allowed here as well
   if( (getAdjacencyVessel()->getMatrixAction())->getNumberOfBaseMultiColvars()>0 ) {
-    if( getBaseMultiColvar(0)->isDensity() ) error("DFS clustering cannot be differentiated if base multicolvar is DENSITY");
+    if( getBaseMultiColvar(0)->isDensity() )
+      error("DFS clustering cannot be differentiated if base multicolvar is DENSITY");
   }
 
   // Ensure that derivatives are turned on in base classes
@@ -57,7 +62,10 @@ void ClusteringBase::turnOnDerivatives() {
 
 void ClusteringBase::calculate() {
   // All the clusters have zero size initially
-  for(unsigned i=0; i<cluster_sizes.size(); ++i) { cluster_sizes[i].first=0; cluster_sizes[i].second=i; }
+  for(unsigned i=0; i<cluster_sizes.size(); ++i) {
+    cluster_sizes[i].first=0;
+    cluster_sizes[i].second=i;
+  }
   // Do the clustering bit
   performClustering();
   // Order the clusters in the system by size (this returns ascending order )
@@ -65,9 +73,13 @@ void ClusteringBase::calculate() {
 }
 
 void ClusteringBase::retrieveAtomsInCluster( const unsigned& clust, std::vector<unsigned>& myatoms ) const {
-  unsigned n=0; myatoms.resize( cluster_sizes[cluster_sizes.size() - clust].first );
+  unsigned n=0;
+  myatoms.resize( cluster_sizes[cluster_sizes.size() - clust].first );
   for(unsigned i=0; i<getNumberOfNodes(); ++i) {
-    if( which_cluster[i]==cluster_sizes[cluster_sizes.size() - clust].second ) { myatoms[n]=i; n++; }
+    if( which_cluster[i]==cluster_sizes[cluster_sizes.size() - clust].second ) {
+      myatoms[n]=i;
+      n++;
+    }
   }
 }
 
diff --git a/src/adjmat/ContactAlignedMatrix.cpp b/src/adjmat/ContactAlignedMatrix.cpp
index 4a5669b955..5f8afdded2 100644
--- a/src/adjmat/ContactAlignedMatrix.cpp
+++ b/src/adjmat/ContactAlignedMatrix.cpp
@@ -95,24 +95,33 @@ ContactAlignedMatrix::ContactAlignedMatrix( const ActionOptions& ao ):
   Action(ao),
   AlignedMatrixBase(ao)
 {
-  unsigned nrows, ncols, ig; retrieveTypeDimensions( nrows, ncols, ig );
+  unsigned nrows, ncols, ig;
+  retrieveTypeDimensions( nrows, ncols, ig );
   sf.resize( nrows, ncols );
   parseConnectionDescriptions("ORIENTATION_SWITCH",false,0);
 }
 
 void ContactAlignedMatrix::readOrientationConnector( const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) {
-  plumed_assert( desc.size()==1 ); std::string errors; sf(j,i).set(desc[0],errors);
-  if( j!=i ) sf(i,j).set(desc[0],errors);
+  plumed_assert( desc.size()==1 );
+  std::string errors;
+  sf(j,i).set(desc[0],errors);
+  if( j!=i )
+    sf(i,j).set(desc[0],errors);
   log.printf("  vectors in %u th and %u th groups must have a dot product that is greater than %s \n",i+1,j+1,(sf(i,j).description()).c_str() );
 }
 
 double ContactAlignedMatrix::computeVectorFunction( const unsigned& iv, const unsigned& jv,
     const Vector& conn, const std::vector<double>& vec1, const std::vector<double>& vec2,
     Vector& dconn, std::vector<double>& dvec1, std::vector<double>& dvec2 ) const {
-  double dot_df, dot=0; dconn.zero();
-  for(unsigned k=2; k<vec1.size(); ++k) dot+=vec1[k]*vec2[k];
+  double dot_df, dot=0;
+  dconn.zero();
+  for(unsigned k=2; k<vec1.size(); ++k)
+    dot+=vec1[k]*vec2[k];
   double f_dot = sf(iv,jv).calculate( dot, dot_df );
-  for(unsigned k=2; k<vec1.size(); ++k) { dvec1[k]=dot_df*vec2[k]; dvec2[k]=dot_df*vec1[k]; }
+  for(unsigned k=2; k<vec1.size(); ++k) {
+    dvec1[k]=dot_df*vec2[k];
+    dvec2[k]=dot_df*vec1[k];
+  }
   return f_dot;
 }
 
diff --git a/src/adjmat/ContactMatrix.cpp b/src/adjmat/ContactMatrix.cpp
index 6e4b68496f..6fddb5ba9a 100644
--- a/src/adjmat/ContactMatrix.cpp
+++ b/src/adjmat/ContactMatrix.cpp
@@ -93,7 +93,8 @@ void ContactMatrix::registerKeywords( Keywords& keys ) {
            "The following provides information on the \\ref switchingfunction that are available. ");
 // I added these keywords so I can test the results I get for column and row sums against output from COORDINATIONNUMBERS
 /// These  should never be used in production as I think they will be much slower than COORDINATIONNUMBERS
-  keys.add("hidden","ATOMSA",""); keys.add("hidden","ATOMSB","");
+  keys.add("hidden","ATOMSA","");
+  keys.add("hidden","ATOMSB","");
 }
 
 ContactMatrix::ContactMatrix( const ActionOptions& ao ):
@@ -102,7 +103,8 @@ ContactMatrix::ContactMatrix( const ActionOptions& ao ):
 {
   // Read in the atoms and setup the matrix
   readMaxTwoSpeciesMatrix( "ATOMS", "ATOMSA", "ATOMSB", true );
-  unsigned nrows, ncols; retrieveTypeDimensions( nrows, ncols, ncol_t );
+  unsigned nrows, ncols;
+  retrieveTypeDimensions( nrows, ncols, ncol_t );
   switchingFunction.resize( nrows, ncols );
   // Read in the switching functions
   parseConnectionDescriptions("SWITCH",false,ncol_t);
@@ -112,7 +114,8 @@ ContactMatrix::ContactMatrix( const ActionOptions& ao ):
   for(unsigned i=0; i<switchingFunction.nrows(); ++i) {
     for(unsigned j=0; j<switchingFunction.ncols(); ++j) {
       double tsf=switchingFunction(i,j).get_dmax();
-      if( tsf>sfmax ) sfmax=tsf;
+      if( tsf>sfmax )
+        sfmax=tsf;
     }
   }
   // And set the link cell cutoff
@@ -120,15 +123,20 @@ ContactMatrix::ContactMatrix( const ActionOptions& ao ):
 }
 
 void ContactMatrix::setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) {
-  plumed_assert( id==0 && desc.size()==1 ); std::string errors; switchingFunction(j,i).set(desc[0],errors);
-  if( errors.length()!=0 ) error("problem reading switching function description " + errors);
-  if( j!=i) switchingFunction(i,j).set(desc[0],errors);
+  plumed_assert( id==0 && desc.size()==1 );
+  std::string errors;
+  switchingFunction(j,i).set(desc[0],errors);
+  if( errors.length()!=0 )
+    error("problem reading switching function description " + errors);
+  if( j!=i)
+    switchingFunction(i,j).set(desc[0],errors);
   log.printf("  %u th and %u th multicolvar groups must be within %s\n",i+1,j+1,(switchingFunction(i,j).description()).c_str() );
 }
 
 double ContactMatrix::calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
   Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
-  if( distance.modulo2()<switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) - ncol_t ).get_dmax2() ) return 1.0;
+  if( distance.modulo2()<switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) - ncol_t ).get_dmax2() )
+    return 1.0;
   return 0.0;
 }
 
diff --git a/src/adjmat/DFSClustering.cpp b/src/adjmat/DFSClustering.cpp
index aea152c059..2cca84221c 100644
--- a/src/adjmat/DFSClustering.cpp
+++ b/src/adjmat/DFSClustering.cpp
@@ -105,21 +105,28 @@ DFSClustering::DFSClustering(const ActionOptions&ao):
   Action(ao),
   ClusteringBase(ao)
 {
-  unsigned maxconnections; parse("MAXCONNECT",maxconnections);
+  unsigned maxconnections;
+  parse("MAXCONNECT",maxconnections);
 #ifdef __PLUMED_HAS_BOOST_GRAPH
-  if( maxconnections>0 ) edge_list.resize( getNumberOfNodes()*maxconnections );
-  else edge_list.resize(0.5*getNumberOfNodes()*(getNumberOfNodes()-1));
+  if( maxconnections>0 )
+    edge_list.resize( getNumberOfNodes()*maxconnections );
+  else
+    edge_list.resize(0.5*getNumberOfNodes()*(getNumberOfNodes()-1));
 #else
-  nneigh.resize( getNumberOfNodes() ); color.resize(getNumberOfNodes());
-  if( maxconnections>0 ) adj_list.resize(getNumberOfNodes(),maxconnections);
-  else adj_list.resize(getNumberOfNodes(),getNumberOfNodes());
+  nneigh.resize( getNumberOfNodes() );
+  color.resize(getNumberOfNodes());
+  if( maxconnections>0 )
+    adj_list.resize(getNumberOfNodes(),maxconnections);
+  else
+    adj_list.resize(getNumberOfNodes(),getNumberOfNodes());
 #endif
 }
 
 void DFSClustering::performClustering() {
 #ifdef __PLUMED_HAS_BOOST_GRAPH
   // Get the list of edges
-  unsigned nedges=0; getAdjacencyVessel()->retrieveEdgeList( nedges, edge_list );
+  unsigned nedges=0;
+  getAdjacencyVessel()->retrieveEdgeList( nedges, edge_list );
 
   // Build the graph using boost
   boost::adjacency_list<boost::vecS,boost::vecS,boost::undirectedS> sg(&edge_list[0],&edge_list[nedges],getNumberOfNodes());
@@ -128,15 +135,20 @@ void DFSClustering::performClustering() {
   number_of_cluster=boost::connected_components(sg,&which_cluster[0]) - 1;
 
   // And work out the size of each cluster
-  for(unsigned i=0; i<which_cluster.size(); ++i) cluster_sizes[which_cluster[i]].first++;
+  for(unsigned i=0; i<which_cluster.size(); ++i)
+    cluster_sizes[which_cluster[i]].first++;
 #else
   // Get the adjacency matrix
   getAdjacencyVessel()->retrieveAdjacencyLists( nneigh, adj_list );
 
   // Perform clustering
-  number_of_cluster=-1; color.assign(color.size(),0);
+  number_of_cluster=-1;
+  color.assign(color.size(),0);
   for(unsigned i=0; i<getNumberOfNodes(); ++i) {
-    if( color[i]==0 ) { number_of_cluster++; color[i]=explore(i); }
+    if( color[i]==0 ) {
+      number_of_cluster++;
+      color[i]=explore(i);
+    }
   }
 #endif
 }
@@ -147,7 +159,8 @@ int DFSClustering::explore( const unsigned& index ) {
   color[index]=1;
   for(unsigned i=0; i<nneigh[index]; ++i) {
     unsigned j=adj_list(index,i);
-    if( color[j]==0 ) color[j]=explore(j);
+    if( color[j]==0 )
+      color[j]=explore(j);
   }
 
   // Count the size of the cluster
diff --git a/src/adjmat/DumpGraph.cpp b/src/adjmat/DumpGraph.cpp
index e7755ee8fd..6321d70425 100644
--- a/src/adjmat/DumpGraph.cpp
+++ b/src/adjmat/DumpGraph.cpp
@@ -63,7 +63,8 @@ class DumpGraph : public ActionPilot {
 PLUMED_REGISTER_ACTION(DumpGraph,"DUMPGRAPH")
 
 void DumpGraph::registerKeywords( Keywords& keys ) {
-  Action::registerKeywords( keys ); ActionPilot::registerKeywords( keys );
+  Action::registerKeywords( keys );
+  ActionPilot::registerKeywords( keys );
   keys.add("compulsory","MATRIX","the action that calculates the adjacency matrix vessel we would like to analyze");
   keys.add("compulsory","STRIDE","1","the frequency with which you would like to output the graph");
   keys.add("compulsory","FILE","the name of the file on which to output the data");
@@ -79,33 +80,46 @@ DumpGraph::DumpGraph( const ActionOptions& ao):
   ActionPilot(ao),
   mymatrix(NULL)
 {
-  parse("MAXCONNECT",maxconnections); std::string mstring; parse("MATRIX",mstring);
+  parse("MAXCONNECT",maxconnections);
+  std::string mstring;
+  parse("MATRIX",mstring);
   AdjacencyMatrixBase* mm = plumed.getActionSet().selectWithLabel<AdjacencyMatrixBase*>( mstring );
-  if( !mm ) error("found no action in set with label " + mstring + " that calculates matrix");
+  if( !mm )
+    error("found no action in set with label " + mstring + " that calculates matrix");
   log.printf("  printing graph for matrix calculated by action %s\n", mm->getLabel().c_str() );
 
   // Retrieve the adjacency matrix of interest
   for(unsigned i=0; i<mm->getNumberOfVessels(); ++i) {
     mymatrix = dynamic_cast<AdjacencyMatrixVessel*>( mm->getPntrToVessel(i) );
-    if( mymatrix ) break ;
+    if( mymatrix )
+      break ;
   }
-  if( !mymatrix ) error( mm->getLabel() + " does not calculate an adjacency matrix");
-  if( !mymatrix->isSymmetric() ) error("input contact matrix must be symmetric");
-  if( maxconnections==0 ) maxconnections=mymatrix->getNumberOfRows();
+  if( !mymatrix )
+    error( mm->getLabel() + " does not calculate an adjacency matrix");
+  if( !mymatrix->isSymmetric() )
+    error("input contact matrix must be symmetric");
+  if( maxconnections==0 )
+    maxconnections=mymatrix->getNumberOfRows();
   parse("FILE",filename);
   log.printf("  printing graph to file named %s \n",filename.c_str() );
   checkRead();
 }
 
 void DumpGraph::update() {
-  OFile ofile; ofile.link(*this); ofile.setBackupString("graph");
-  ofile.open( filename ); ofile.printf("graph G { \n");
+  OFile ofile;
+  ofile.link(*this);
+  ofile.setBackupString("graph");
+  ofile.open( filename );
+  ofile.printf("graph G { \n");
   // Print all nodes
-  for(unsigned i=0; i<mymatrix->getNumberOfRows(); ++i) ofile.printf("%u [label=\"%u\"];\n",i,i);
+  for(unsigned i=0; i<mymatrix->getNumberOfRows(); ++i)
+    ofile.printf("%u [label=\"%u\"];\n",i,i);
   // Now retrieve connectivitives
-  unsigned nedge; std::vector<std::pair<unsigned,unsigned> > edge_list( mymatrix->getNumberOfRows()*maxconnections );
+  unsigned nedge;
+  std::vector<std::pair<unsigned,unsigned> > edge_list( mymatrix->getNumberOfRows()*maxconnections );
   mymatrix->retrieveEdgeList( nedge, edge_list );
-  for(unsigned i=0; i<nedge; ++i) ofile.printf("%u -- %u \n", edge_list[i].first, edge_list[i].second );
+  for(unsigned i=0; i<nedge; ++i)
+    ofile.printf("%u -- %u \n", edge_list[i].first, edge_list[i].second );
   ofile.printf("} \n");
 }
 
diff --git a/src/adjmat/HbondMatrix.cpp b/src/adjmat/HbondMatrix.cpp
index 887ec94ff3..89b0208f4d 100644
--- a/src/adjmat/HbondMatrix.cpp
+++ b/src/adjmat/HbondMatrix.cpp
@@ -140,8 +140,11 @@ HBondMatrix::HBondMatrix( const ActionOptions& ao ):
   AdjacencyMatrixBase(ao)
 {
   readMaxThreeSpeciesMatrix( "ATOMS", "DONORS", "ACCEPTORS", "HYDROGENS", false );
-  unsigned nrows, ncols; retrieveTypeDimensions( nrows, ncols, ndonor_types );
-  distanceOOSwitch.resize( nrows, ncols ); distanceOHSwitch.resize( nrows, ncols ); angleSwitch.resize( nrows, ncols );
+  unsigned nrows, ncols;
+  retrieveTypeDimensions( nrows, ncols, ndonor_types );
+  distanceOOSwitch.resize( nrows, ncols );
+  distanceOHSwitch.resize( nrows, ncols );
+  angleSwitch.resize( nrows, ncols );
   parseConnectionDescriptions("SWITCH",false,ndonor_types);
   parseConnectionDescriptions("HSWITCH",false,ndonor_types);
   parseConnectionDescriptions("ASWITCH",false,ndonor_types);
@@ -151,7 +154,8 @@ HBondMatrix::HBondMatrix( const ActionOptions& ao ):
   for(unsigned i=0; i<getNumberOfNodeTypes(); ++i) {
     for(unsigned j=0; j<getNumberOfNodeTypes(); ++j) {
       double tsf=distanceOOSwitch(i,j).get_dmax();
-      if( tsf>sfmax ) sfmax=tsf;
+      if( tsf>sfmax )
+        sfmax=tsf;
     }
   }
   // Set the link cell cutoff
@@ -161,45 +165,60 @@ HBondMatrix::HBondMatrix( const ActionOptions& ao ):
 void HBondMatrix::setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) {
   plumed_assert( id<3 && desc.size()==1 );
   if( id==0 ) {
-    std::string errors; distanceOOSwitch(j,i).set(desc[0],errors);
-    if( errors.length()!=0 ) error("problem reading switching function description " + errors);
-    if( j!=i) distanceOOSwitch(i,j).set(desc[0],errors);
+    std::string errors;
+    distanceOOSwitch(j,i).set(desc[0],errors);
+    if( errors.length()!=0 )
+      error("problem reading switching function description " + errors);
+    if( j!=i)
+      distanceOOSwitch(i,j).set(desc[0],errors);
     log.printf("  atoms of type %u and %u must be within %s\n",i+1,j+1,(distanceOOSwitch(i,j).description()).c_str() );
   } else if( id==1 ) {
-    std::string errors; distanceOHSwitch(j,i).set(desc[0],errors);
-    if( errors.length()!=0 ) error("problem reading switching function description " + errors);
-    if( j!=i) distanceOHSwitch(i,j).set(desc[0],errors);
+    std::string errors;
+    distanceOHSwitch(j,i).set(desc[0],errors);
+    if( errors.length()!=0 )
+      error("problem reading switching function description " + errors);
+    if( j!=i)
+      distanceOHSwitch(i,j).set(desc[0],errors);
     log.printf("  for atoms of type %u and %u the OH distance must be less than %s \n",i+1,j+1,(distanceOHSwitch(i,j).description()).c_str() );
   } else if( id==2 ) {
-    std::string errors; angleSwitch(j,i).set(desc[0],errors);
-    if( errors.length()!=0 ) error("problem reading switching function description " + errors);
-    if( j!=i) angleSwitch(i,j).set(desc[0],errors);
+    std::string errors;
+    angleSwitch(j,i).set(desc[0],errors);
+    if( errors.length()!=0 )
+      error("problem reading switching function description " + errors);
+    if( j!=i)
+      angleSwitch(i,j).set(desc[0],errors);
     log.printf("  for atoms of type %u and %u the OOH angle must be less than %s \n",i+1,j+1,(angleSwitch(i,j).description()).c_str() );
   }
 }
 
 double HBondMatrix::calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
   // Ensure we skip diagonal elements of square matrix
-  if( myatoms.getIndex(0)==myatoms.getIndex(1) ) return 0.0;
+  if( myatoms.getIndex(0)==myatoms.getIndex(1) )
+    return 0.0;
 
   Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
-  if( distance.modulo2()<distanceOOSwitch( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).get_dmax2() ) return 1.0;
+  if( distance.modulo2()<distanceOOSwitch( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).get_dmax2() )
+    return 1.0;
   return 0.0;
 }
 
 double HBondMatrix::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
-  Vector ood = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) ); double ood_l = ood.modulo(); // acceptor - donor
+  Vector ood = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
+  double ood_l = ood.modulo(); // acceptor - donor
   double ood_df, ood_sw=distanceOOSwitch( getBaseColvarNumber( myatoms.getIndex(0) ),
                                           getBaseColvarNumber( myatoms.getIndex(1) ) ).calculate( ood_l, ood_df );
 
   // Get the base colvar numbers
   unsigned ano, dno = getBaseColvarNumber( myatoms.getIndex(0) );
-  if( ndonor_types==0 ) ano = getBaseColvarNumber( myatoms.getIndex(1) );
-  else ano = getBaseColvarNumber( myatoms.getIndex(1) ) - ndonor_types;
+  if( ndonor_types==0 )
+    ano = getBaseColvarNumber( myatoms.getIndex(1) );
+  else
+    ano = getBaseColvarNumber( myatoms.getIndex(1) ) - ndonor_types;
 
   double value=0;
   if( myatoms.getNumberOfAtoms()>3 ) {
-    for(unsigned i=2; i<myatoms.getNumberOfAtoms(); ++i) value+=calculateForThree( i, ano, dno, ood, ood_df, ood_sw,  myatoms );
+    for(unsigned i=2; i<myatoms.getNumberOfAtoms(); ++i)
+      value+=calculateForThree( i, ano, dno, ood, ood_df, ood_sw,  myatoms );
   } else {
     plumed_dbg_assert( myatoms.getNumberOfAtoms()==3 );
     value=calculateForThree( 2, ano, dno, ood, ood_df, ood_sw, myatoms );
@@ -209,11 +228,14 @@ double HBondMatrix::compute( const unsigned& tindex, multicolvar::AtomValuePack&
 
 double HBondMatrix::calculateForThree( const unsigned& iat, const unsigned& ano, const unsigned& dno, const Vector& ood,
                                        const double& ood_df, const double& ood_sw, multicolvar::AtomValuePack& myatoms ) const {
-  Vector ohd=getSeparation( myatoms.getPosition(0), myatoms.getPosition(iat) ); double ohd_l=ohd.modulo();
+  Vector ohd=getSeparation( myatoms.getPosition(0), myatoms.getPosition(iat) );
+  double ohd_l=ohd.modulo();
   double ohd_df, ohd_sw=distanceOHSwitch( getBaseColvarNumber( myatoms.getIndex(0) ),
                                           getBaseColvarNumber( myatoms.getIndex(1) ) ).calculate( ohd_l, ohd_df );
 
-  Angle a; Vector ood_adf, ohd_adf; double angle=a.compute( ood, ohd, ood_adf, ohd_adf );
+  Angle a;
+  Vector ood_adf, ohd_adf;
+  double angle=a.compute( ood, ohd, ood_adf, ohd_adf );
   double angle_df, angle_sw=angleSwitch( getBaseColvarNumber( myatoms.getIndex(0) ),
                                          getBaseColvarNumber( myatoms.getIndex(1) ) ).calculate( angle, angle_df );
 
diff --git a/src/adjmat/MatrixColumnSums.cpp b/src/adjmat/MatrixColumnSums.cpp
index a3775c2f41..6f3e3beeb6 100644
--- a/src/adjmat/MatrixColumnSums.cpp
+++ b/src/adjmat/MatrixColumnSums.cpp
@@ -76,34 +76,50 @@ PLUMED_REGISTER_ACTION(MatrixColumnSums,"COLUMNSUMS")
 
 void MatrixColumnSums::registerKeywords( Keywords& keys ) {
   ActionWithInputMatrix::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");
+  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");
 }
 
 MatrixColumnSums::MatrixColumnSums(const ActionOptions& ao):
   Action(ao),
   ActionWithInputMatrix(ao)
 {
-  if( (mymatrix->getMatrixAction())->mybasemulticolvars.size()>0 ) error("matrix row sums should only be calculated when inputs are atoms");
+  if( (mymatrix->getMatrixAction())->mybasemulticolvars.size()>0 )
+    error("matrix row sums should only be calculated when inputs are atoms");
   // Setup the tasks
   unsigned ncols = mymatrix->getNumberOfColumns();
-  ablocks.resize(1); ablocks[0].resize( ncols );
-  for(unsigned i=0; i<ncols; ++i) addTaskToList( i );
+  ablocks.resize(1);
+  ablocks[0].resize( ncols );
+  for(unsigned i=0; i<ncols; ++i)
+    addTaskToList( i );
   // Set the positions - this is only used when getting positions for central atoms
   if( mymatrix->undirectedGraph() ) {
-    for(unsigned i=0; i<ncols; ++i) ablocks[0][i]=i;
+    for(unsigned i=0; i<ncols; ++i)
+      ablocks[0][i]=i;
   } else {
-    for(unsigned i=0; i<ncols; ++i) ablocks[0][i]=mymatrix->getNumberOfRows() + i;
+    for(unsigned i=0; i<ncols; ++i)
+      ablocks[0][i]=mymatrix->getNumberOfRows() + i;
   }
-  std::vector<AtomNumber> fake_atoms; setupMultiColvarBase( fake_atoms );
+  std::vector<AtomNumber> fake_atoms;
+  setupMultiColvarBase( fake_atoms );
 }
 
 double MatrixColumnSums::compute( const unsigned& tinded, multicolvar::AtomValuePack& myatoms ) const {
-  double sum=0.0; std::vector<double> tvals( mymatrix->getNumberOfComponents() );
+  double sum=0.0;
+  std::vector<double> tvals( mymatrix->getNumberOfComponents() );
   unsigned nrows = mymatrix->getNumberOfRows();
   for(unsigned i=0; i<nrows; ++i) {
-    if( mymatrix->undirectedGraph() && tinded==i ) continue;
+    if( mymatrix->undirectedGraph() && tinded==i )
+      continue;
     sum+=retrieveConnectionValue( i, tinded, tvals );
   }
 
@@ -111,7 +127,8 @@ double MatrixColumnSums::compute( const unsigned& tinded, multicolvar::AtomValue
     MultiValue myvals( mymatrix->getNumberOfComponents(), myatoms.getNumberOfDerivatives() );
     MultiValue& myvout=myatoms.getUnderlyingMultiValue();
     for(unsigned i=0; i<nrows; ++i) {
-      if( mymatrix->isSymmetric() && tinded==i ) continue ;
+      if( mymatrix->isSymmetric() && tinded==i )
+        continue ;
       addConnectionDerivatives( i, tinded, myvals, myvout );
     }
   }
diff --git a/src/adjmat/MatrixRowSums.cpp b/src/adjmat/MatrixRowSums.cpp
index 4905f9b0aa..e96943b3b4 100644
--- a/src/adjmat/MatrixRowSums.cpp
+++ b/src/adjmat/MatrixRowSums.cpp
@@ -76,29 +76,45 @@ PLUMED_REGISTER_ACTION(MatrixRowSums,"ROWSUMS")
 
 void MatrixRowSums::registerKeywords( Keywords& keys ) {
   ActionWithInputMatrix::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");
+  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");
 }
 
 MatrixRowSums::MatrixRowSums(const ActionOptions& ao):
   Action(ao),
   ActionWithInputMatrix(ao)
 {
-  if( (mymatrix->getMatrixAction())->mybasemulticolvars.size()>0 ) warning("matrix row may be problematic when inputs are not atoms");
+  if( (mymatrix->getMatrixAction())->mybasemulticolvars.size()>0 )
+    warning("matrix row may be problematic when inputs are not atoms");
   // Setup the tasks
   unsigned nrows = mymatrix->getNumberOfRows();
-  ablocks.resize(1); ablocks[0].resize( nrows );
-  for(unsigned i=0; i<nrows; ++i) { ablocks[0][i]=i; addTaskToList( i ); }
-  std::vector<AtomNumber> fake_atoms; setupMultiColvarBase( fake_atoms );
+  ablocks.resize(1);
+  ablocks[0].resize( nrows );
+  for(unsigned i=0; i<nrows; ++i) {
+    ablocks[0][i]=i;
+    addTaskToList( i );
+  }
+  std::vector<AtomNumber> fake_atoms;
+  setupMultiColvarBase( fake_atoms );
 }
 
 double MatrixRowSums::compute( const unsigned& tinded, multicolvar::AtomValuePack& myatoms ) const {
   std::vector<double> tvals( mymatrix->getNumberOfComponents() );
-  getInputData( tinded, false, myatoms, tvals ); double fval=tvals[1];
+  getInputData( tinded, false, myatoms, tvals );
+  double fval=tvals[1];
 
   if( !doNotCalculateDerivatives() ) {
-    tvals.assign( tvals.size(), 0 ); tvals[1]=1.0;
+    tvals.assign( tvals.size(), 0 );
+    tvals[1]=1.0;
     mergeInputDerivatives( 1, 1, 2, tinded, tvals, getInputDerivatives( tinded, false, myatoms ), myatoms );
   }
   return fval;
diff --git a/src/adjmat/OutputCluster.cpp b/src/adjmat/OutputCluster.cpp
index 7fe7a3cd27..cb29341051 100644
--- a/src/adjmat/OutputCluster.cpp
+++ b/src/adjmat/OutputCluster.cpp
@@ -106,30 +106,44 @@ OutputCluster::OutputCluster(const ActionOptions& ao):
   myclusters(NULL)
 {
   // Setup output file
-  ofile.link(*this); std::string file; parse("FILE",file);
-  if( file.length()==0 ) error("output file name was not specified");
+  ofile.link(*this);
+  std::string file;
+  parse("FILE",file);
+  if( file.length()==0 )
+    error("output file name was not specified");
   // Search for xyz extension
   output_xyz=false;
   if( file.find(".")!=std::string::npos ) {
     std::size_t dot=file.find_first_of('.');
-    if( file.substr(dot+1)=="xyz" ) output_xyz=true;
+    if( file.substr(dot+1)=="xyz" )
+      output_xyz=true;
   }
 
-  ofile.open(file); log.printf("  on file %s \n",file.c_str());
-  parseFlag("MAKE_WHOLE",makewhole); parse("MAXDEPTH",maxdepth); parse("MAXGOES",maxgoes);
-  if( makewhole && !output_xyz) error("MAKE_WHOLE flag is not compatible with output of non-xyz files");
+  ofile.open(file);
+  log.printf("  on file %s \n",file.c_str());
+  parseFlag("MAKE_WHOLE",makewhole);
+  parse("MAXDEPTH",maxdepth);
+  parse("MAXGOES",maxgoes);
+  if( makewhole && !output_xyz)
+    error("MAKE_WHOLE flag is not compatible with output of non-xyz files");
 
   // Find what action we are taking the clusters from
-  std::vector<std::string> matname(1); parse("CLUSTERS",matname[0]);
+  std::vector<std::string> matname(1);
+  parse("CLUSTERS",matname[0]);
   myclusters = plumed.getActionSet().selectWithLabel<ClusteringBase*>( matname[0] );
-  if( !myclusters ) error( matname[0] + " does not calculate perform a clustering of the atomic positions");
+  if( !myclusters )
+    error( matname[0] + " does not calculate perform a clustering of the atomic positions");
   // N.B. the +0.3 is a fudge factor.  Reconstrucing PBC doesnt work without this GAT
-  addDependency( myclusters ); double rcut=myclusters->getCutoffForConnection() + 0.3; rcut2=rcut*rcut;
+  addDependency( myclusters );
+  double rcut=myclusters->getCutoffForConnection() + 0.3;
+  rcut2=rcut*rcut;
 
   // Read in the cluster we are calculating
   parse("CLUSTER",clustr);
-  if( clustr<1 ) error("cannot look for a cluster larger than the largest cluster");
-  if( clustr>myclusters->getNumberOfNodes() ) error("cluster selected is invalid - too few atoms in system");
+  if( clustr<1 )
+    error("cannot look for a cluster larger than the largest cluster");
+  if( clustr>myclusters->getNumberOfNodes() )
+    error("cluster selected is invalid - too few atoms in system");
   log.printf("  outputting atoms in %u th largest cluster found by %s \n",clustr,matname[0].c_str() );
 }
 
@@ -141,35 +155,56 @@ void OutputCluster::update() {
     if( makewhole ) {
       // Retrieve the atom positions
       atomsin.resize( myatoms.size() );
-      for(unsigned i=0; i<myatoms.size(); ++i) atomsin[i]=myclusters->getPositionOfAtomForLinkCells( myatoms[i] );
+      for(unsigned i=0; i<myatoms.size(); ++i)
+        atomsin[i]=myclusters->getPositionOfAtomForLinkCells( myatoms[i] );
       // Build a connectivity matrix neglecting the pbc
-      nneigh.resize( myatoms.size(), 0 ); adj_list.resize( myatoms.size(), myatoms.size() );
+      nneigh.resize( myatoms.size(), 0 );
+      adj_list.resize( myatoms.size(), myatoms.size() );
       for(unsigned i=1; i<myatoms.size(); ++i) {
         for(unsigned j=0; j<i; ++j) {
-          if( delta( atomsin[i], atomsin[j] ).modulo2()<=rcut2 ) { adj_list(i,nneigh[i])=j; adj_list(j,nneigh[j])=i; nneigh[i]++; nneigh[j]++; }
+          if( delta( atomsin[i], atomsin[j] ).modulo2()<=rcut2 ) {
+            adj_list(i,nneigh[i])=j;
+            adj_list(j,nneigh[j])=i;
+            nneigh[i]++;
+            nneigh[j]++;
+          }
         }
       }
       // Use DFS to find the largest cluster not broken by PBC
-      number_of_cluster=-1; visited.resize( myatoms.size(), false );
-      cluster_sizes.resize( myatoms.size() ); which_cluster.resize( myatoms.size() );
-      for(unsigned i=0; i<cluster_sizes.size(); ++i) { cluster_sizes[i].first=0; cluster_sizes[i].second=i; }
+      number_of_cluster=-1;
+      visited.resize( myatoms.size(), false );
+      cluster_sizes.resize( myatoms.size() );
+      which_cluster.resize( myatoms.size() );
+      for(unsigned i=0; i<cluster_sizes.size(); ++i) {
+        cluster_sizes[i].first=0;
+        cluster_sizes[i].second=i;
+      }
 
       for(unsigned i=0; i<myatoms.size(); ++i) {
-        if( !visited[i] ) { number_of_cluster++; visited[i]=explore_dfs(i); }
+        if( !visited[i] ) {
+          number_of_cluster++;
+          visited[i]=explore_dfs(i);
+        }
       }
       std::sort( cluster_sizes.begin(), cluster_sizes.end() );
 
       // Now set visited so that only those atoms in largest cluster will be start points for PBCing
       visited.assign( visited.size(), false );
       for(unsigned i=0; i<myatoms.size(); ++i) {
-        if( which_cluster[i]==cluster_sizes[cluster_sizes.size()-1].second ) visited[i]=true;
+        if( which_cluster[i]==cluster_sizes[cluster_sizes.size()-1].second )
+          visited[i]=true;
       }
 
       // Now retrieve the original connectivity matrix (including pbc)
       nneigh.assign( nneigh.size(), 0 );
       for(unsigned i=1; i<myatoms.size(); ++i) {
         for(unsigned j=0; j<i; ++j) {
-          if( myclusters->areConnected( myatoms[i], myatoms[j] ) ) { adj_list(i,nneigh[i])=j; adj_list(j,nneigh[j])=i; nneigh[i]++; nneigh[j]++; }
+          if( myclusters->areConnected( myatoms[i], myatoms[j] ) ) {
+            adj_list(i,nneigh[i])=j;
+            adj_list(j,nneigh[j])=i;
+            nneigh[i]++;
+            nneigh[j]++;
+          }
         }
       }
 
@@ -177,18 +212,21 @@ void OutputCluster::update() {
       for(unsigned jj=0; jj<maxgoes; ++jj) {
 
         for(unsigned j=0; j<myatoms.size(); ++j) {
-          if( !visited[j] ) continue;
+          if( !visited[j] )
+            continue;
 
           for(unsigned k=0; k<nneigh[j]; ++k) {
             if( delta( atomsin[j],atomsin[adj_list(j,k)] ).modulo2()>rcut2 ) {
               visited[j]=true;
-              for(unsigned depth=0; depth<=maxdepth; ++depth) explore( j, depth );
+              for(unsigned depth=0; depth<=maxdepth; ++depth)
+                explore( j, depth );
             }
           }
         }
       }
       // And print final positions
-      for(unsigned i=0; i<myatoms.size(); ++i) ofile.printf( "X %f %f %f \n", atomsin[i][0], atomsin[i][1], atomsin[i][2] );
+      for(unsigned i=0; i<myatoms.size(); ++i)
+        ofile.printf( "X %f %f %f \n", atomsin[i][0], atomsin[i][1], atomsin[i][2] );
     } else {
       for(unsigned i=0; i<myatoms.size(); ++i) {
         Vector pos=myclusters->getPositionOfAtomForLinkCells( myatoms[i] );
@@ -198,16 +236,19 @@ void OutputCluster::update() {
   } else {
     ofile.printf("CLUSTERING RESULTS AT TIME %f : NUMBER OF ATOMS IN %u TH LARGEST CLUSTER EQUALS %u \n",getTime(),clustr,static_cast<unsigned>(myatoms.size()) );
     ofile.printf("INDICES OF ATOMS : ");
-    for(unsigned i=0; i<myatoms.size(); ++i) ofile.printf("%d ",(myclusters->getAbsoluteIndexOfCentralAtom(myatoms[i])).index());
+    for(unsigned i=0; i<myatoms.size(); ++i)
+      ofile.printf("%d ",(myclusters->getAbsoluteIndexOfCentralAtom(myatoms[i])).index());
     ofile.printf("\n");
   }
 }
 
 void OutputCluster::explore( const unsigned& index, const unsigned& depth ) {
-  if( depth==0 ) return ;
+  if( depth==0 )
+    return ;
 
   for(unsigned i=0; i<nneigh[index]; ++i) {
-    unsigned j=adj_list(index,i); visited[j]=true;
+    unsigned j=adj_list(index,i);
+    visited[j]=true;
     Vector svec=myclusters->pbcDistance( atomsin[index], atomsin[j] );
     atomsin[j] = atomsin[index] + svec;
     explore( j, depth-1 );
@@ -218,7 +259,8 @@ bool OutputCluster::explore_dfs( const unsigned& index ) {
   visited[index]=true;
   for(unsigned i=0; i<nneigh[index]; ++i) {
     unsigned j=adj_list(index,i);
-    if( !visited[j] ) visited[j]=explore_dfs(j);
+    if( !visited[j] )
+      visited[j]=explore_dfs(j);
   }
 
   // Count the size of the cluster
diff --git a/src/adjmat/SMACMatrix.cpp b/src/adjmat/SMACMatrix.cpp
index bd4dede889..41ad1a22c3 100644
--- a/src/adjmat/SMACMatrix.cpp
+++ b/src/adjmat/SMACMatrix.cpp
@@ -103,47 +103,65 @@ SMACMatrix::SMACMatrix( const ActionOptions& ao ):
   Action(ao),
   AlignedMatrixBase(ao)
 {
-  unsigned nrows, ncols, ig; retrieveTypeDimensions( nrows, ncols, ig );
-  kernels.resize( nrows, ncols ); parseConnectionDescriptions("KERNEL",true,0);
+  unsigned nrows, ncols, ig;
+  retrieveTypeDimensions( nrows, ncols, ig );
+  kernels.resize( nrows, ncols );
+  parseConnectionDescriptions("KERNEL",true,0);
 }
 
 void SMACMatrix::readOrientationConnector( const unsigned& iv, const unsigned& jv, const std::vector<std::string>& desc ) {
   for(int i=0; i<desc.size(); i++) {
     KernelFunctions mykernel( desc[i] );
     kernels(iv,jv).push_back( mykernel );
-    if( jv!=iv ) kernels(jv,iv).push_back( mykernel );
+    if( jv!=iv )
+      kernels(jv,iv).push_back( mykernel );
   }
-  if( kernels(iv,jv).size()==0 ) error("no kernels defined");
+  if( kernels(iv,jv).size()==0 )
+    error("no kernels defined");
 }
 
 double SMACMatrix::computeVectorFunction( const unsigned& iv, const unsigned& jv,
     const Vector& conn, const std::vector<double>& vec1, const std::vector<double>& vec2,
     Vector& dconn, std::vector<double>& dvec1, std::vector<double>& dvec2 ) const {
 
-  unsigned nvectors = ( vec1.size() - 2 ) / 3; plumed_assert( (vec1.size()-2)%3==0 );
-  std::vector<Vector> dv1(nvectors), dv2(nvectors), tdconn(nvectors); Torsion t; std::vector<Vector> v1(nvectors), v2(nvectors);
+  unsigned nvectors = ( vec1.size() - 2 ) / 3;
+  plumed_assert( (vec1.size()-2)%3==0 );
+  std::vector<Vector> dv1(nvectors), dv2(nvectors), tdconn(nvectors);
+  Torsion t;
+  std::vector<Vector> v1(nvectors), v2(nvectors);
   std::vector<std::unique_ptr<Value>> pos;
-  for(unsigned i=0; i<nvectors; ++i) { pos.emplace_back( Tools::make_unique<Value>() ); pos[i]->setDomain( "-pi", "pi" ); }
+  for(unsigned i=0; i<nvectors; ++i) {
+    pos.emplace_back( Tools::make_unique<Value>() );
+    pos[i]->setDomain( "-pi", "pi" );
+  }
 
   for(unsigned j=0; j<nvectors; ++j) {
     for(unsigned k=0; k<3; ++k) {
-      v1[j][k]=vec1[2+3*j+k]; v2[j][k]=vec2[2+3*j+k];
+      v1[j][k]=vec1[2+3*j+k];
+      v2[j][k]=vec2[2+3*j+k];
     }
     double angle = t.compute( v1[j], conn, v2[j], dv1[j], tdconn[j], dv2[j] );
     pos[j]->set( angle );
   }
 
-  double ans=0; std::vector<double> deriv( nvectors ), df( nvectors, 0 );
+  double ans=0;
+  std::vector<double> deriv( nvectors ), df( nvectors, 0 );
 
   auto pos_ptr=Tools::unique2raw(pos);
 
   for(unsigned i=0; i<kernels(iv,jv).size(); ++i) {
     ans += kernels(iv,jv)[i].evaluate( pos_ptr, deriv );
-    for(unsigned j=0; j<nvectors; ++j) df[j] += deriv[j];
+    for(unsigned j=0; j<nvectors; ++j)
+      df[j] += deriv[j];
   }
-  dconn.zero(); for(unsigned j=0; j<nvectors; ++j) dconn += df[j]*tdconn[j];
+  dconn.zero();
+  for(unsigned j=0; j<nvectors; ++j)
+    dconn += df[j]*tdconn[j];
   for(unsigned j=0; j<nvectors; ++j) {
-    for(unsigned k=0; k<3; ++k) { dvec1[2+3*j+k]=df[j]*dv1[j][k]; dvec2[2+3*j+k]=df[j]*dv2[j][k]; }
+    for(unsigned k=0; k<3; ++k) {
+      dvec1[2+3*j+k]=df[j]*dv1[j][k];
+      dvec2[2+3*j+k]=df[j]*dv2[j][k];
+    }
   }
   return ans;
 }
diff --git a/src/adjmat/Sprint.cpp b/src/adjmat/Sprint.cpp
index d0c135883c..224b50227b 100644
--- a/src/adjmat/Sprint.cpp
+++ b/src/adjmat/Sprint.cpp
@@ -131,13 +131,16 @@ Sprint::Sprint(const ActionOptions&ao):
   //    // if( !getBaseMultiColvar(i)->hasDifferentiableOrientation() ) error("cannot use multicolvar of type " + getBaseMultiColvar(i)->getName() );
   // }
 
-  if( !getAdjacencyVessel()->isSymmetric() ) error("input contact matrix is not symmetric");
-  std::vector<AtomNumber> fake_atoms; setupMultiColvarBase( fake_atoms );
+  if( !getAdjacencyVessel()->isSymmetric() )
+    error("input contact matrix is not symmetric");
+  std::vector<AtomNumber> fake_atoms;
+  setupMultiColvarBase( fake_atoms );
 
   // Create all the values
   sqrtn = std::sqrt( static_cast<double>( getNumberOfNodes() ) );
   for(unsigned i=0; i<getNumberOfNodes(); ++i) {
-    std::string num; Tools::convert(i,num);
+    std::string num;
+    Tools::convert(i,num);
     addComponentWithDerivatives("coord-"+num);
     componentIsNotPeriodic("coord-"+num);
     getPntrToComponent(i)->resizeDerivatives( getNumberOfDerivatives() );
@@ -145,7 +148,8 @@ Sprint::Sprint(const ActionOptions&ao):
 
   // Setup the dynamic list to hold all the tasks
   unsigned ntriangle = 0.5*getNumberOfNodes()*(getNumberOfNodes()-1);
-  for(unsigned i=0; i<ntriangle; ++i) active_elements.addIndexToList( i );
+  for(unsigned i=0; i<ntriangle; ++i)
+    active_elements.addIndexToList( i );
 }
 
 void Sprint::calculate() {
@@ -179,16 +183,24 @@ void Sprint::calculate() {
 
   // Parallelism
   unsigned rank, stride;
-  if( serialCalculation() ) { stride=1; rank=0; }
-  else { rank=comm.Get_rank(); stride=comm.Get_size(); }
+  if( serialCalculation() ) {
+    stride=1;
+    rank=0;
+  }
+  else {
+    rank=comm.Get_rank();
+    stride=comm.Get_size();
+  }
 
   // Derivatives
   MultiValue myvals( 2, getNumberOfDerivatives() );
   Matrix<double> mymat_ders( getNumberOfComponents(), getNumberOfDerivatives() );
   // std::vector<unsigned> catoms(2);
-  unsigned nval = getNumberOfNodes(); mymat_ders=0;
+  unsigned nval = getNumberOfNodes();
+  mymat_ders=0;
   for(unsigned i=rank; i<active_elements.getNumberActive(); i+=stride) {
-    unsigned j, k; getAdjacencyVessel()->getMatrixIndices( active_elements[i], j, k );
+    unsigned j, k;
+    getAdjacencyVessel()->getMatrixIndices( active_elements[i], j, k );
     double tmp1 = 2 * eigenvecs(nval-1,j)*eigenvecs(nval-1,k);
     for(int icomp=0; icomp<getNumberOfComponents(); ++icomp) {
       double tmp2 = 0.;
@@ -203,11 +215,13 @@ void Sprint::calculate() {
       }
     }
   }
-  if( !serialCalculation() ) comm.Sum( mymat_ders );
+  if( !serialCalculation() )
+    comm.Sum( mymat_ders );
 
   for(int j=0; j<getNumberOfComponents(); ++j) {
     Value* val=getPntrToComponent(j);
-    for(unsigned i=0; i<getNumberOfDerivatives(); ++i) val->addDerivative( i, mymat_ders(j,i) );
+    for(unsigned i=0; i<getNumberOfDerivatives(); ++i)
+      val->addDerivative( i, mymat_ders(j,i) );
   }
 }
 
diff --git a/src/adjmat/TopologyMatrix.cpp b/src/adjmat/TopologyMatrix.cpp
index 07abbd501a..043989ae03 100644
--- a/src/adjmat/TopologyMatrix.cpp
+++ b/src/adjmat/TopologyMatrix.cpp
@@ -102,18 +102,26 @@ TopologyMatrix::TopologyMatrix( const ActionOptions& ao ):
   maxbins(0)
 {
   readMaxThreeSpeciesMatrix("NODES", "FAKE", "FAKE", "ATOMS", true );
-  unsigned nrows, ncols, ndonor_types; retrieveTypeDimensions( nrows, ncols, ndonor_types );
-  switchingFunction.resize( nrows, ncols ); parseConnectionDescriptions("SWITCH",false,ndonor_types);
-  cylinder_sw.resize( nrows, ncols ); parseConnectionDescriptions("RADIUS",false,ndonor_types);
-  low_sf.resize( nrows, ncols ); parseConnectionDescriptions("CYLINDER_SWITCH",false,ndonor_types);
-  binw_mat.resize( nrows, ncols ); cell_volume.resize( nrows, ncols );
+  unsigned nrows, ncols, ndonor_types;
+  retrieveTypeDimensions( nrows, ncols, ndonor_types );
+  switchingFunction.resize( nrows, ncols );
+  parseConnectionDescriptions("SWITCH",false,ndonor_types);
+  cylinder_sw.resize( nrows, ncols );
+  parseConnectionDescriptions("RADIUS",false,ndonor_types);
+  low_sf.resize( nrows, ncols );
+  parseConnectionDescriptions("CYLINDER_SWITCH",false,ndonor_types);
+  binw_mat.resize( nrows, ncols );
+  cell_volume.resize( nrows, ncols );
   parseConnectionDescriptions("BIN_SIZE",false,ndonor_types);
   // Read in stuff for grid
-  parse("SIGMA",sigma); parse("KERNEL",kerneltype);
+  parse("SIGMA",sigma);
+  parse("KERNEL",kerneltype);
   // Read in threshold for density cutoff
-  std::string errors, thresh_sw_str; parse("DENSITY_THRESHOLD",thresh_sw_str);
+  std::string errors, thresh_sw_str;
+  parse("DENSITY_THRESHOLD",thresh_sw_str);
   threshold_switch.set(thresh_sw_str, errors );
-  if( errors.length()>0 ) error("errors in DENSITY_THRESHOLD switching function : " + errors );
+  if( errors.length()>0 )
+    error("errors in DENSITY_THRESHOLD switching function : " + errors );
   log.printf("  threshold on density of atoms in cylinder equals %s\n",threshold_switch.description().c_str() );
 
   for(unsigned i=0; i<getNumberOfNodeTypes(); ++i) {
@@ -130,16 +138,21 @@ TopologyMatrix::TopologyMatrix( const ActionOptions& ao ):
   for(unsigned i=0; i<getNumberOfNodeTypes(); ++i) {
     for(unsigned j=0; j<getNumberOfNodeTypes(); ++j) {
       double tsf=switchingFunction(i,j).get_dmax();
-      if( tsf>sfmax ) sfmax=tsf;
+      if( tsf>sfmax )
+        sfmax=tsf;
       double rsf=cylinder_sw(i,j).get_dmax();
-      if( rsf>rfmax ) rfmax=rsf;
+      if( rsf>rfmax )
+        rfmax=rsf;
       double lsf=low_sf(i,j).get_dmax();
-      if( lsf>lsfmax ) lsfmax=lsf;
+      if( lsf>lsfmax )
+        lsfmax=lsf;
     }
   }
   // Get the width of the bead
-  HistogramBead bead; bead.isNotPeriodic();
-  bead.setKernelType( kerneltype ); bead.set( 0.0, 1.0, sigma );
+  HistogramBead bead;
+  bead.isNotPeriodic();
+  bead.setKernelType( kerneltype );
+  bead.set( 0.0, 1.0, sigma );
   beadrad = bead.getCutoff();
 
   // Set the link cell cutoff
@@ -149,7 +162,8 @@ TopologyMatrix::TopologyMatrix( const ActionOptions& ao ):
   double maxsize=0;
   for(unsigned i=0; i<getNumberOfNodeTypes(); ++i) {
     for(unsigned j=0; j<getNumberOfNodeTypes(); ++j) {
-      if( binw_mat(i,j)>maxsize ) maxsize=binw_mat(i,j);
+      if( binw_mat(i,j)>maxsize )
+        maxsize=binw_mat(i,j);
     }
   }
   // Set the maximum number of bins that we will need to compute
@@ -166,55 +180,75 @@ unsigned TopologyMatrix::getNumberOfQuantities() const {
 void TopologyMatrix::setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) {
   plumed_assert( id<4 );
   if( id==0 ) {
-    std::string errors; switchingFunction(j,i).set(desc[0],errors);
-    if( errors.length()!=0 ) error("problem reading switching function description " + errors);
-    if( j!=i) switchingFunction(i,j).set(desc[0],errors);
+    std::string errors;
+    switchingFunction(j,i).set(desc[0],errors);
+    if( errors.length()!=0 )
+      error("problem reading switching function description " + errors);
+    if( j!=i)
+      switchingFunction(i,j).set(desc[0],errors);
     log.printf("  %u th and %u th multicolvar groups must be within %s\n",i+1,j+1,(switchingFunction(i,j).description()).c_str() );
   } else if( id==1 ) {
-    std::string errors; cylinder_sw(j,i).set(desc[0],errors);
-    if( errors.length()!=0 ) error("problem reading switching function description " + errors);
-    if( j!=i) cylinder_sw(i,j).set(desc[0],errors);
+    std::string errors;
+    cylinder_sw(j,i).set(desc[0],errors);
+    if( errors.length()!=0 )
+      error("problem reading switching function description " + errors);
+    if( j!=i)
+      cylinder_sw(i,j).set(desc[0],errors);
     log.printf("  there must be not atoms within the cylinder connections atoms of multicolvar groups %u th and %u th.  This cylinder has radius %s \n",i+1,j+1,(cylinder_sw(i,j).description()).c_str() );
   } else if( id==2 ) {
-    std::string errors; low_sf(j,i).set(desc[0],errors);
-    if( errors.length()!=0 ) error("problem reading switching function description " + errors);
-    if( j!=i ) low_sf(i,j).set(desc[0],errors);
+    std::string errors;
+    low_sf(j,i).set(desc[0],errors);
+    if( errors.length()!=0 )
+      error("problem reading switching function description " + errors);
+    if( j!=i )
+      low_sf(i,j).set(desc[0],errors);
     log.printf("  %u th and %u th multicolvar groups must be further apart than %s\n",i+1,j+1,(low_sf(j,i).description()).c_str() );
   } else if( id==3 ) {
     Tools::convert( desc[0], binw_mat(j,i) );
-    if( i!=j ) binw_mat(i,j)=binw_mat(j,i);
+    if( i!=j )
+      binw_mat(i,j)=binw_mat(j,i);
     log.printf("  cylinder for %u th and %u th multicolvar groups is split into bins of length %f \n",i,j,binw_mat(i,j) );
   }
 }
 
 double TopologyMatrix::calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
   Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
-  if( distance.modulo2()<switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).get_dmax2() ) return 1.0;
+  if( distance.modulo2()<switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).get_dmax2() )
+    return 1.0;
   return 0.0;
 }
 
 double TopologyMatrix::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
-  HistogramBead bead; bead.isNotPeriodic(); bead.setKernelType( kerneltype );
+  HistogramBead bead;
+  bead.isNotPeriodic();
+  bead.setKernelType( kerneltype );
 
   // Initialise to zero density on all bins
-  for(unsigned bin=0; bin<maxbins; ++bin) myatoms.setValue(bin+1,0);
+  for(unsigned bin=0; bin<maxbins; ++bin)
+    myatoms.setValue(bin+1,0);
   // Calculate whether or not atoms 1 and 2 are within cutoff (can use delta here as pbc are done in atom setup)
-  Vector d1 = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) ); double d1_len = d1.modulo();
+  Vector d1 = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
+  double d1_len = d1.modulo();
   d1 = d1 / d1_len;  // Convert vector into director
   AtomNumber a1 = myatoms.getAbsoluteIndex( 0 );
   AtomNumber a2 = myatoms.getAbsoluteIndex( 1 );
   for(unsigned i=2; i<myatoms.getNumberOfAtoms(); ++i) {
     AtomNumber a3 = myatoms.getAbsoluteIndex( i );
-    if( a3!=a1 && a3!=a2 ) calculateForThreeAtoms( i, d1, d1_len, bead, myatoms );
+    if( a3!=a1 && a3!=a2 )
+      calculateForThreeAtoms( i, d1, d1_len, bead, myatoms );
   }
   // std::vector<double> binvals( 1+maxbins ); for(unsigned i=1;i<maxbins;++i) binvals[i]=myatoms.getValue(i);
   // unsigned ii; double fdf;
   //std::cout<<"HELLO DENSITY "<<myatoms.getIndex(0)<<" "<<myatoms.getIndex(1)<<" "<<transformStoredValues( binvals, ii, fdf )<<std::endl;
 
   // Now find the element for which the density is maximal
-  unsigned vout=2; double max=myatoms.getValue( 2 );
+  unsigned vout=2;
+  double max=myatoms.getValue( 2 );
   for(unsigned i=3; i<myatoms.getUnderlyingMultiValue().getNumberOfValues()-1; ++i) {
-    if( myatoms.getValue(i)>max ) { max=myatoms.getValue(i); vout=i; }
+    if( myatoms.getValue(i)>max ) {
+      max=myatoms.getValue(i);
+      vout=i;
+    }
   }
   // Calculate value and derivative of switching function between atoms 1 and 2
   double dfuncl, sw = switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ),
@@ -248,12 +282,14 @@ void TopologyMatrix::calculateForThreeAtoms( const unsigned& iat, const Vector&
   // This tells us if we are outside the end of the cylinder
   double excess = proj - d1_len;
   // Return if we are outside of the cylinder as calculated based on excess
-  if( excess>low_sf( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).get_dmax() ) return;
+  if( excess>low_sf( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).get_dmax() )
+    return;
   // Find the length of the cylinder
   double binw = binw_mat( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) );
   double lcylinder = (std::floor( d1_len / binw ) + 1)*binw;
   // Return if the projection is outside the length of interest
-  if( proj<-bead.getCutoff() || proj>(lcylinder+bead.getCutoff()) ) return;
+  if( proj<-bead.getCutoff() || proj>(lcylinder+bead.getCutoff()) )
+    return;
 
   // Calculate the excess swiching function
   double edf, eval = low_sf( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).calculate( excess, edf );
@@ -297,11 +333,14 @@ void TopologyMatrix::calculateForThreeAtoms( const unsigned& iat, const Vector&
 
     Vector pos1 = myatoms.getPosition(0) + d1_len*d1;
     Vector pos2 = myatoms.getPosition(0) + d2;
-    Vector g1derivf,g2derivf,lderivf; Tensor vir;
+    Vector g1derivf,g2derivf,lderivf;
+    Tensor vir;
     for(unsigned bin=0; bin<maxbins; ++bin) {
       bead.set( bin*binw, (bin+1)*binw, sigma );
-      if( proj<(bin*binw-bead.getCutoff()) || proj>binw*(bin+1)+bead.getCutoff() ) continue;
-      double der, contr=bead.calculateWithCutoff( proj, der ) / cellv; der /= cellv;
+      if( proj<(bin*binw-bead.getCutoff()) || proj>binw*(bin+1)+bead.getCutoff() )
+        continue;
+      double der, contr=bead.calculateWithCutoff( proj, der ) / cellv;
+      der /= cellv;
       myatoms.addValue( 2+bin, contr*val*eval );
 
       if( !doNotCalculateDerivatives() ) {
diff --git a/src/astyle.sh b/src/astyle.sh
index 7ff80b99b1..463066547e 100755
--- a/src/astyle.sh
+++ b/src/astyle.sh
@@ -8,7 +8,6 @@ echo "$DIRS"
 # shellcheck disable=SC2125
 test -z "$DIRS" && DIRS=*
 use_legacy=(
-  adjmat
   analysis
   annfunc
   asmjit
diff --git a/src/config/Config.inc.in b/src/config/Config.inc.in
index 0a6f99e9cf..9c55128e68 100644
--- a/src/config/Config.inc.in
+++ b/src/config/Config.inc.in
@@ -29,20 +29,20 @@ namespace PLMD {
 namespace config {
 
 namespace {
-  /// local tool that, given a string, returns a new string which is:
-  /// - enclosed in single quotes (')
-  /// - with all single quotes escaped
-  std::string escapeForSingleQuote(const std::string& input) {
-    std::string escaped;
-    for (char c : input) {
-        if (c == '\'') {
-            escaped += "'\\''";
-        } else {
-            escaped += c;
-        }
+/// local tool that, given a string, returns a new string which is:
+/// - enclosed in single quotes (')
+/// - with all single quotes escaped
+std::string escapeForSingleQuote(const std::string& input) {
+  std::string escaped;
+  for (char c : input) {
+    if (c == '\'') {
+      escaped += "'\\''";
+    } else {
+      escaped += c;
     }
-    return "'" + escaped + "'";
   }
+  return "'" + escaped + "'";
+}
 }
 
 // This is a fix to allow conda to correctly replace paths in binary files.