Removed code that was used for chained version of derivatives that is…

… no longer used for forces
plumed · Sep 20, 2024 · 46f53ba · 46f53ba · PlumedBot · Sep 20, 2024
1 parent 03e8500
commit 46f53ba
Show file tree

Hide file tree

Showing 63 changed files with 325 additions and 1,109 deletions.
diff --git a/src/adjmat/AdjacencyMatrixBase.cpp b/src/adjmat/AdjacencyMatrixBase.cpp
@@ -206,25 +206,6 @@ void AdjacencyMatrixBase::updateNeighbourList() {
   if( read_one_group && maxcol==getConstPntrToComponent(0)->getShape()[1] ) maxcol = maxcol-1;
 }
 
-void AdjacencyMatrixBase::getAdditionalTasksRequired( ActionWithVector* action, std::vector<unsigned>& atasks ) {
-  if( action==this ) return;
-  // Update the neighbour list
-  updateNeighbourList();
-
-  unsigned nactive = atasks.size();
-  std::vector<unsigned> indlist( 1 + ablocks.size() + threeblocks.size() );
-  for(unsigned i=0; i<nactive; ++i) {
-    unsigned num = retrieveNeighbours( atasks[i], indlist );
-    for(unsigned j=0; j<num; ++j) {
-      bool found=false;
-      for(unsigned k=0; k<atasks.size(); ++k ) {
-        if( indlist[j]==atasks[k] ) { found=true; break; }
-      }
-      if( !found ) atasks.push_back( indlist[j] );
-    }
-  }
-}
-
 unsigned AdjacencyMatrixBase::retrieveNeighbours( const unsigned& current, std::vector<unsigned> & indices ) const {
   unsigned natoms=nlist[current]; indices[0]=current;
   unsigned lstart = getConstPntrToComponent(0)->getShape()[0] + current*(1+natoms_per_list); plumed_dbg_assert( nlist[lstart]==current );
@@ -278,10 +259,9 @@ void AdjacencyMatrixBase::performTask( const std::string& controller, const unsi
     // Update dynamic list indices for virial
     unsigned base = 3*getNumberOfAtoms(); for(unsigned j=0; j<9; ++j) myvals.updateIndex( 0, base+j );
     // And the indices for the derivatives of the row of the matrix
-    unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
-    std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
+    unsigned nmat_ind = myvals.getNumberOfMatrixRowDerivatives(); std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices() );
     matrix_indices[nmat_ind+0]=3*index2+0; matrix_indices[nmat_ind+1]=3*index2+1; matrix_indices[nmat_ind+2]=3*index2+2;
-    myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind+3 );
+    myvals.setNumberOfMatrixRowDerivatives( nmat_ind+3 );
   }
 
   // Calculate the components if we need them
@@ -319,38 +299,30 @@ void AdjacencyMatrixBase::performTask( const std::string& controller, const unsi
       myvals.addDerivative( 3, base+1, 0 ); myvals.addDerivative( 3, base+4, 0 ); myvals.addDerivative( 3, base+7, 0 );
       myvals.addDerivative( 3, base+2, -atom[0] ); myvals.addDerivative( 3, base+5, -atom[1] ); myvals.addDerivative( 3, base+8, -atom[2] );
       for(unsigned k=0; k<9; ++k) { myvals.updateIndex( 1, base+k ); myvals.updateIndex( 2, base+k ); myvals.updateIndex( 3, base+k ); }
-      for(unsigned k=1; k<4; ++k) {
-        unsigned nmat = getConstPntrToComponent(k)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
-        std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
-        matrix_indices[nmat_ind+0]=3*index2+0; matrix_indices[nmat_ind+1]=3*index2+1; matrix_indices[nmat_ind+2]=3*index2+2;
-        myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind+3 );
-      }
     }
   }
 }
 
 void AdjacencyMatrixBase::runEndOfRowJobs( const unsigned& ind, const std::vector<unsigned> & indices, MultiValue& myvals ) const {
   if( doNotCalculateDerivatives() ) return;
 
-  for(int k=0; k<getNumberOfComponents(); ++k) {
-    unsigned nmat = getConstPntrToComponent(k)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
-    std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
-    plumed_assert( nmat_ind<matrix_indices.size() );
-    matrix_indices[nmat_ind+0]=3*ind+0;
-    matrix_indices[nmat_ind+1]=3*ind+1;
-    matrix_indices[nmat_ind+2]=3*ind+2;
+  unsigned nmat_ind = myvals.getNumberOfMatrixRowDerivatives();
+  std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices() );
+  plumed_assert( nmat_ind<matrix_indices.size() );
+  matrix_indices[nmat_ind+0]=3*ind+0;
+  matrix_indices[nmat_ind+1]=3*ind+1;
+  matrix_indices[nmat_ind+2]=3*ind+2;
+  nmat_ind+=3;
+  for(unsigned i=myvals.getSplitIndex(); i<myvals.getNumberOfIndices(); ++i) {
+    matrix_indices[nmat_ind+0]=3*indices[i]+0;
+    matrix_indices[nmat_ind+1]=3*indices[i]+1;
+    matrix_indices[nmat_ind+2]=3*indices[i]+2;
     nmat_ind+=3;
-    for(unsigned i=myvals.getSplitIndex(); i<myvals.getNumberOfIndices(); ++i) {
-      matrix_indices[nmat_ind+0]=3*indices[i]+0;
-      matrix_indices[nmat_ind+1]=3*indices[i]+1;
-      matrix_indices[nmat_ind+2]=3*indices[i]+2;
-      nmat_ind+=3;
-    }
-    unsigned virbase = 3*getNumberOfAtoms();
-    for(unsigned i=0; i<9; ++i) matrix_indices[nmat_ind+i]=virbase+i;
-    nmat_ind+=9; plumed_dbg_massert( nmat_ind<=3*getNumberOfAtoms() + 9, "found too many derivatives in " + getLabel() );
-    myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind );
   }
+  unsigned virbase = 3*getNumberOfAtoms();
+  for(unsigned i=0; i<9; ++i) matrix_indices[nmat_ind+i]=virbase+i;
+  nmat_ind+=9; plumed_dbg_massert( nmat_ind<=3*getNumberOfAtoms() + 9, "found too many derivatives in " + getLabel() );
+  myvals.setNumberOfMatrixRowDerivatives( nmat_ind );
 }
 
 }

diff --git a/src/adjmat/AdjacencyMatrixBase.h b/src/adjmat/AdjacencyMatrixBase.h
@@ -54,7 +54,6 @@ class AdjacencyMatrixBase : public ActionWithMatrix {
   unsigned getNumberOfDerivatives() override ;
   unsigned getNumberOfColumns() const override;
   void prepare() override;
-  void getAdditionalTasksRequired( ActionWithVector* action, std::vector<unsigned>& atasks ) override ;
   void setupForTask( const unsigned& current, std::vector<unsigned> & indices, MultiValue& myvals ) const override;
   // void setupCurrentTaskList() override;
   void updateNeighbourList() override ;

diff --git a/src/adjmat/Neighbors.cpp b/src/adjmat/Neighbors.cpp
@@ -42,7 +42,6 @@ class Neighbors : public ActionWithMatrix {
   explicit Neighbors(const ActionOptions&);
   unsigned getNumberOfDerivatives() override;
   unsigned getNumberOfColumns() const override { return number; }
-  bool canBeAfterInChain( ActionWithVector* av ) { return av->getLabel()!=(getPntrToArgument(0)->getPntrToAction())->getLabel(); }
   void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const ;
   void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override;
   void runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const override {}
@@ -64,7 +63,6 @@ Neighbors::Neighbors(const ActionOptions&ao):
 {
   if( getNumberOfArguments()!=1 ) error("found wrong number of arguments in input");
   if( getPntrToArgument(0)->getRank()!=2 ) error("input argument should be a matrix");
-  getPntrToArgument(0)->buildDataStore();
 
   unsigned nlow; parse("NLOWEST",nlow);
   unsigned nhigh; parse("NHIGHEST",nhigh);

diff --git a/src/adjmat/TorsionsMatrix.cpp b/src/adjmat/TorsionsMatrix.cpp
@@ -36,9 +36,6 @@ namespace PLMD {
 namespace adjmat {
 
 class TorsionsMatrix : public ActionWithMatrix {
-private:
-  unsigned nderivatives;
-  bool stored_matrix1, stored_matrix2;
 public:
   static void registerKeywords( Keywords& keys );
   explicit TorsionsMatrix(const ActionOptions&);
@@ -85,10 +82,7 @@ TorsionsMatrix::TorsionsMatrix(const ActionOptions&ao):
   log.printf("\n"); requestAtoms( atoms_a, false );
 
   std::vector<unsigned> shape(2); shape[0]=getPntrToArgument(0)->getShape()[0]; shape[1]=getPntrToArgument(1)->getShape()[1];
-  addValue( shape ); setPeriodic("-pi","pi"); nderivatives = buildArgumentStore(0) + 3*getNumberOfAtoms() + 9;
-  std::string headstr=getFirstActionInChain()->getLabel();
-  stored_matrix1 = getPntrToArgument(0)->ignoreStoredValue( headstr );
-  stored_matrix2 = getPntrToArgument(1)->ignoreStoredValue( headstr );
+  addValue( shape ); setPeriodic("-pi","pi");
 
   if( nm>0 ) {
     unsigned iarg = getNumberOfArguments()-1;
@@ -98,7 +92,7 @@ TorsionsMatrix::TorsionsMatrix(const ActionOptions&ao):
 }
 
 unsigned TorsionsMatrix::getNumberOfDerivatives() {
-  return nderivatives;
+  return 3*getNumberOfAtoms() + 9 + getPntrToArgument(0)->getNumberOfStoredValues() + getPntrToArgument(1)->getNumberOfStoredValues();
 }
 
 unsigned TorsionsMatrix::getNumberOfColumns() const {
@@ -131,8 +125,8 @@ void TorsionsMatrix::performTask( const std::string& controller, const unsigned&
 
   // Get the two vectors
   for(unsigned i=0; i<3; ++i) {
-    v1[i] = getElementOfMatrixArgument( 0, index1, i, myvals );
-    v2[i] = getElementOfMatrixArgument( 1, i, ind2, myvals );
+    v1[i] = getPntrToArgument(0)->get( index1*getPntrToArgument(0)->getNumberOfColumns() + i, false );
+    v2[i] = getPntrToArgument(1)->get( i*getPntrToArgument(1)->getNumberOfColumns() + ind2, false );
   }
   // Evaluate angle
   Torsion t; double angle = t.compute( v1, conn, v2, dv1, dconn, dv2 );
@@ -141,9 +135,12 @@ void TorsionsMatrix::performTask( const std::string& controller, const unsigned&
   if( doNotCalculateDerivatives() ) return;
 
   // Add the derivatives on the matrices
+  unsigned base1 = index1*getPntrToArgument(0)->getNumberOfColumns();
+  unsigned ncols = getPntrToArgument(1)->getNumberOfColumns();
+  unsigned base2 = getPntrToArgument(0)->getNumberOfStoredValues() + ind2;
   for(unsigned i=0; i<3; ++i) {
-    addDerivativeOnMatrixArgument( stored_matrix1, 0, 0, index1, i, dv1[i], myvals );
-    addDerivativeOnMatrixArgument( stored_matrix2, 0, 1, i, ind2, dv2[i], myvals );
+    myvals.addDerivative( 0, base1 + i, dv1[i] ); myvals.updateIndex( 0, base1 + i );
+    myvals.addDerivative( 0, base2 + i*ncols, dv2[i] ); myvals.updateIndex( 0, base2 + i*ncols );
   }
   // And derivatives on positions
   unsigned narg_derivatives = getPntrToArgument(0)->getNumberOfValues() + getPntrToArgument(1)->getNumberOfValues();
@@ -160,20 +157,20 @@ void TorsionsMatrix::runEndOfRowJobs( const unsigned& ival, const std::vector<un
   if( doNotCalculateDerivatives() ) return ;
 
   unsigned mat1s = 3*ival, ss = getPntrToArgument(1)->getShape()[1];
-  unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
+  unsigned nmat_ind = myvals.getNumberOfMatrixRowDerivatives();
   unsigned narg_derivatives = getPntrToArgument(0)->getNumberOfValues() + getPntrToArgument(1)->getNumberOfValues();
-  std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) ); unsigned ntwo_atoms = myvals.getSplitIndex();
+  std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices() ); unsigned ntwo_atoms = myvals.getSplitIndex();
   for(unsigned j=0; j<3; ++j) {
     matrix_indices[nmat_ind] = mat1s + j; nmat_ind++;
     matrix_indices[nmat_ind] = narg_derivatives + mat1s + j; nmat_ind++;
     for(unsigned i=1; i<ntwo_atoms; ++i) {
       unsigned ind2 = indices[i]; if( ind2>=getPntrToArgument(0)->getShape()[0] ) ind2 = indices[i] - getPntrToArgument(0)->getShape()[0];
-      matrix_indices[nmat_ind] = arg_deriv_starts[1] + j*ss + ind2; nmat_ind++;
+      matrix_indices[nmat_ind] = getPntrToArgument(0)->getNumberOfStoredValues() + j*ss + ind2; nmat_ind++;
       matrix_indices[nmat_ind] = narg_derivatives + 3*indices[i] + j; nmat_ind++;
     }
   }
   unsigned base = narg_derivatives + 3*getNumberOfAtoms(); for(unsigned j=0; j<9; ++j) { matrix_indices[nmat_ind] = base + j; nmat_ind++; }
-  myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind );
+  myvals.setNumberOfMatrixRowDerivatives( nmat_ind );
 }
 
 }

diff --git a/src/clusters/ClusterDistribution.cpp b/src/clusters/ClusterDistribution.cpp
@@ -127,11 +127,9 @@ ClusterDistribution::ClusterDistribution(const ActionOptions&ao):
   }
   // Request the arguments
   requestArguments( clusters );
-  getPntrToArgument(0)->buildDataStore();
-  if( getNumberOfArguments()>1 ) getPntrToArgument(1)->buildDataStore();
   // Now create the value
   std::vector<unsigned> shape(1); shape[0]=clusters[0]->getShape()[0];
-  addValue( shape ); setNotPeriodic(); getPntrToValue()->buildDataStore();
+  addValue( shape ); setNotPeriodic();
 }
 
 unsigned ClusterDistribution::getNumberOfDerivatives() {

diff --git a/src/clusters/ClusterWeights.cpp b/src/clusters/ClusterWeights.cpp
@@ -89,7 +89,7 @@ ClusterWeights::ClusterWeights(const ActionOptions&ao):
   requestArguments( clusters );
   // Now create the value
   std::vector<unsigned> shape(1); shape[0]=clusters[0]->getShape()[0];
-  addValue( shape ); setNotPeriodic(); getPntrToComponent(0)->buildDataStore();
+  addValue( shape ); setNotPeriodic();
   // Find out which cluster we want
   parse("CLUSTER",clustr);
   if( clustr<1 ) error("cannot look for a cluster larger than the largest cluster");

diff --git a/src/clusters/ClusteringBase.cpp b/src/clusters/ClusteringBase.cpp
@@ -42,7 +42,7 @@ ClusteringBase::ClusteringBase(const ActionOptions&ao):
   // Now create a value - this holds the data on which cluster each guy is in
   std::vector<unsigned> shape(1); shape[0]=getPntrToArgument(0)->getShape()[0];
   // Build the store here to make sure that next action has all data
-  addValue( shape ); setNotPeriodic(); getPntrToValue()->buildDataStore();
+  addValue( shape ); setNotPeriodic();
   // Resize local variables
   which_cluster.resize( getPntrToArgument(0)->getShape()[0] ); cluster_sizes.resize( getPntrToArgument(0)->getShape()[0] );
 }

diff --git a/src/colvar/MultiColvarTemplate.h b/src/colvar/MultiColvarTemplate.h
@@ -44,7 +44,6 @@ class MultiColvarTemplate : public ActionWithVector {
   unsigned getNumberOfDerivatives() override ;
   void addValueWithDerivatives( const std::vector<unsigned>& shape=std::vector<unsigned>() ) override ;
   void addComponentWithDerivatives( const std::string& name, const std::vector<unsigned>& shape=std::vector<unsigned>() ) override ;
-  void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) override ;
   void performTask( const unsigned&, MultiValue& ) const override ;
   void calculate() override;
 };
@@ -114,6 +113,7 @@ unsigned MultiColvarTemplate<T>::getNumberOfDerivatives() {
 
 template <class T>
 void MultiColvarTemplate<T>::calculate() {
+  if( wholemolecules ) makeWhole();
   runAllTasks();
 }
 
@@ -127,12 +127,6 @@ void MultiColvarTemplate<T>::addComponentWithDerivatives( const std::string& nam
   std::vector<unsigned> s(1); s[0]=ablocks[0].size(); addComponent( name, s );
 }
 
-template <class T>
-void MultiColvarTemplate<T>::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) {
-  if( wholemolecules ) makeWhole();
-  ActionWithVector::setupStreamedComponents( headstr, nquants, nmat, maxcol );
-}
-
 template <class T>
 void MultiColvarTemplate<T>::performTask( const unsigned& task_index, MultiValue& myvals ) const {
   // Retrieve the positions

diff --git a/src/colvar/RMSDVector.cpp b/src/colvar/RMSDVector.cpp
@@ -94,13 +94,13 @@ RMSDVector::RMSDVector(const ActionOptions&ao):
   if( displacement && (getPntrToArgument(1)->getRank()==1 || getPntrToArgument(1)->getShape()[0]<=1) ) {
     addComponentWithDerivatives("dist"); componentIsNotPeriodic("dist");
     std::vector<unsigned> shape( 1, getPntrToArgument(0)->getNumberOfValues() );
-    addComponent( "disp", shape ); getPntrToComponent(1)->buildDataStore(); componentIsNotPeriodic("disp");
+    addComponent( "disp", shape ); componentIsNotPeriodic("disp");
   } else if( displacement ) {
     std::vector<unsigned> shape( 1, getPntrToArgument(1)->getShape()[0] );
-    addComponent( "dist", shape ); getPntrToComponent(0)->buildDataStore();
+    addComponent( "dist", shape );
     componentIsNotPeriodic("dist");
     shape.resize(2); shape[0] = getPntrToArgument(1)->getShape()[0]; shape[1] = getPntrToArgument(0)->getNumberOfValues();
-    addComponent( "disp", shape ); getPntrToComponent(1)->buildDataStore(); getPntrToComponent(1)->reshapeMatrixStore( shape[1] );
+    addComponent( "disp", shape ); getPntrToComponent(1)->reshapeMatrixStore( shape[1] );
     componentIsNotPeriodic("disp");
   } else if( (getPntrToArgument(1)->getRank()==1 || getPntrToArgument(1)->getShape()[0]==1) ) {
     addValue(); setNotPeriodic();
@@ -278,11 +278,18 @@ void RMSDVector::gatherStoredValue( const unsigned& valindex, const unsigned& co
   }
 }
 
-void RMSDVector::gatherForcesOnStoredValue( const unsigned& ival, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
-  if( getConstPntrToComponent(ival)->getRank()==1 ) { ActionWithVector::gatherForcesOnStoredValue( ival, itask, myvals, forces ); return; }
-  const std::vector<Vector>& direction( myvals.getConstFirstAtomDerivativeVector()[1] ); unsigned natoms = direction.size();
-  for(unsigned i=0; i<natoms; ++i) {
-    for(unsigned j=0; j<3; ++j ) forces[j*natoms+i] += direction[i][j];
+void RMSDVector::gatherForces( const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
+  if( checkComponentsForForce() ) {
+    for(unsigned ival=0; ival<getNumberOfComponents(); ++ival) {
+      if( getConstPntrToComponent(ival)->getRank()==1 ) {
+        gatherForcesOnStoredValue( ival, itask, myvals, forces );
+      } else {
+        const std::vector<Vector>& direction( myvals.getConstFirstAtomDerivativeVector()[1] ); unsigned natoms = direction.size();
+        for(unsigned i=0; i<natoms; ++i) {
+          for(unsigned j=0; j<3; ++j ) forces[j*natoms+i] += direction[i][j];
+        }
+      }
+    }
   }
 }
 

diff --git a/src/colvar/RMSDVector.h b/src/colvar/RMSDVector.h
@@ -46,7 +46,7 @@ class RMSDVector : public ActionWithVector {
   void performTask( const unsigned& current, MultiValue& myvals ) const override ;
   void gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
                           const unsigned& bufstart, std::vector<double>& buffer ) const override ;
-  void gatherForcesOnStoredValue( const unsigned& ival, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const override ;
+  void gatherForces( const unsigned& i, const MultiValue& myvals, std::vector<double>& forces ) const override ;
   void setReferenceConfigurations();
   void calculate() override ;
   bool checkForTaskForce( const unsigned& itask, const Value* myval ) const override ;

diff --git a/src/contour/DistanceFromContourBase.cpp b/src/contour/DistanceFromContourBase.cpp
@@ -46,10 +46,8 @@ DistanceFromContourBase::DistanceFromContourBase( const ActionOptions& ao ):
   nactive(0)
 {
   if( getNumberOfArguments()>1 ) error("should only use one argument for this action");
-  if( getNumberOfArguments()==1 ) {
-    if( getPntrToArgument(0)->getRank()!=1 ) error("ARG for distance from contour should be rank one");
-    getPntrToArgument(0)->buildDataStore();
-  }
+  if( getNumberOfArguments()==1 &&  getPntrToArgument(0)->getRank()!=1 ) error("ARG for distance from contour should be rank one");
+
   // Read in the multicolvar/atoms
   std::vector<AtomNumber> atoms; parseAtomList("POSITIONS",atoms);
   std::vector<AtomNumber> origin; parseAtomList("ATOM",origin);

diff --git a/src/contour/FindContour.cpp b/src/contour/FindContour.cpp
@@ -109,11 +109,7 @@ FindContour::FindContour(const ActionOptions&ao):
   std::vector<std::string> argn( ag->getGridCoordinateNames() );
 
   std::vector<unsigned> shape(1); shape[0]=0;
-  for(unsigned i=0; i<argn.size(); ++i ) {
-    addComponent( argn[i], shape ); componentIsNotPeriodic( argn[i] ); getPntrToComponent(i)->buildDataStore();
-  }
-  // Check for task reduction
-  updateTaskListReductionStatus();
+  for(unsigned i=0; i<argn.size(); ++i ) { addComponent( argn[i], shape ); componentIsNotPeriodic( argn[i] ); }
 }
 
 std::string FindContour::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {

diff --git a/src/contour/FindContourSurface.cpp b/src/contour/FindContourSurface.cpp
@@ -149,7 +149,6 @@ FindContourSurface::FindContourSurface(const ActionOptions&ao):
   // Now add a value
   std::vector<unsigned> shape( getInputGridObject().getDimension()-1 );
   addValueWithDerivatives( shape ); setNotPeriodic();
-  getPntrToComponent(0)->buildDataStore();
 }
 
 void FindContourSurface::setupValuesOnFirstStep() {

diff --git a/src/contour/FindSphericalContour.cpp b/src/contour/FindSphericalContour.cpp
@@ -151,7 +151,6 @@ FindSphericalContour::FindSphericalContour(const ActionOptions&ao):
   // Now create a value
   std::vector<unsigned> shape( 3 ); shape[0]=npoints; shape[1]=shape[2]=1;
   addValueWithDerivatives( shape ); setNotPeriodic();
-  getPntrToComponent(0)->buildDataStore(); checkRead();
 }
 
 unsigned FindSphericalContour::getNumberOfDerivatives() {

diff --git a/src/core/ActionToGetData.cpp b/src/core/ActionToGetData.cpp
@@ -60,7 +60,7 @@ ActionToGetData::ActionToGetData(const ActionOptions&ao):
 
   if( getNumberOfArguments()!=1 ) error("python interface works best when you ask for one argument at a time");
   if( getPntrToArgument(0)->getNumberOfValues()==0 ) error("cannot get data as shape of value " + getPntrToArgument(0)->getName() + " has not been set");
-  getPntrToArgument(0)->buildDataStore(); data.resize( getPntrToArgument(0)->getNumberOfValues() );
+  data.resize( getPntrToArgument(0)->getNumberOfValues() );
 }
 
 void ActionToGetData::get_rank( const TypesafePtr & dims ) {