Added method to calculate ensemble averages

regtest/analysis/rt-average/Makefile

+include ../../scripts/test.make

regtest/analysis/rt-average/colvar.reference

+#! FIELDS time d1 d1a t1 t1a
+#! SET min_t1 -pi
+#! SET max_t1 pi
+#! SET min_t1a -pi
+#! SET max_t1a pi
+ 0.000000   0.1522   0.0000  -3.0750   0.0000
+ 0.005000   0.1521   0.1521  -3.0833  -3.0833
+ 0.010000   0.1522   0.1522  -3.0682  -3.0757
+ 0.015000   0.1522   0.1522  -3.0406  -3.0640
+ 0.020000   0.1521   0.1522  -3.0855  -3.0694
+ 0.025000   0.1522   0.1522  -3.1023  -3.0760
+ 0.030000   0.1522   0.1522   3.0682  -3.0991
+ 0.035000   0.1522   0.1522   2.9670  -3.1300
+ 0.040000   0.1522   0.1522   3.1264  -3.1334
+ 0.045000   0.1522   0.1522  -3.1150  -3.1313
+ 0.050000   0.1522   0.1522   2.9720   3.1339
+ 0.055000   0.1522   0.1522  -3.1376   3.1350

regtest/analysis/rt-average/config

+type=driver
+arg="--plumed plumed.dat --trajectory-stride 1 --timestep 0.005 --ixyz ala12_trajectory.xyz --dump-forces forces --dump-forces-fmt=%10.6f"

regtest/analysis/rt-average/plumed.dat

+d1: DISTANCE ATOMS=2,5
+
+TORSION ATOMS=2,5,6,7  LABEL=t1
+TORSION ATOMS=5,6,7,9  LABEL=t2
+
+d1a: AVERAGE ARG=d1
+t1a: AVERAGE ARG=t1
+
+PRINT ARG=d1,d1a,t1,t1a STRIDE=1 FILE=colvar FMT=%8.4f

src/analysis/Average.cpp

+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2014,2015 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed-code.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "vesselbase/ActionWithAveraging.h"
+#include "core/ActionRegister.h"
+#include "AverageVessel.h"
+
+//+PLUMEDOC GRIDCALC AVERAGE 
+/* 
+Calculate the ensemble average of a collective variable
+
+\par Examples
+
+*/
+//+ENDPLUMEDOC
+
+namespace PLMD {
+namespace analysis {
+
+class Average : public vesselbase::ActionWithAveraging {
+private:
+  AverageVessel* myaverage;
+public:
+  static void registerKeywords( Keywords& keys );
+  explicit Average( const ActionOptions& );
+  void performOperations( const bool& from_update );
+  void finishAveraging();
+  bool isPeriodic(){ return false; } 
+  void performTask( const unsigned& , const unsigned& , MultiValue& ) const { plumed_error(); }
+};
+
+PLUMED_REGISTER_ACTION(Average,"AVERAGE")
+
+void Average::registerKeywords( Keywords& keys ){
+  vesselbase::ActionWithAveraging::registerKeywords( keys ); keys.use("ARG");
+}
+
+Average::Average( const ActionOptions& ao ):
+Action(ao),
+ActionWithAveraging(ao)
+{
+   addValue(); // Create a value so that we can output the average
+   if( getNumberOfArguments()!=1 ) error("only one quantity can be averaged at a time");
+   std::string instring; 
+   if( getPntrToArgument(0)->isPeriodic() ){
+      std::string min, max; getPntrToArgument(0)->getDomain(min,max); 
+      instring = "PERIODIC=" + min + "," + max; setPeriodic( min, max );
+   } else {
+      setNotPeriodic();
+   }
+   // Create a vessel to hold the average
+   vesselbase::VesselOptions da("myaverage","",-1,instring,this);
+   Keywords keys; AverageVessel::registerKeywords( keys );
+   vesselbase::VesselOptions dar( da, keys );
+   myaverage = new AverageVessel(dar); setAveragingAction( myaverage, false );
+}
+
+void Average::performOperations( const bool& from_update ){
+   myaverage->accumulate( cweight, getArgument(0) );
+}
+
+void Average::finishAveraging(){
+   setValue( myaverage->getAverage() );
+}
+
+}
+}

src/analysis/AverageVessel.cpp

+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2012-2016 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "AverageVessel.h"
+
+namespace PLMD {
+namespace analysis {
+
+void AverageVessel::registerKeywords( Keywords& keys ){
+  vesselbase::AveragingVessel::registerKeywords( keys );
+  keys.add("optional","PERIODIC","is the quantity being averaged periodic and what is its domain");
+}
+
+AverageVessel::AverageVessel( const vesselbase::VesselOptions& da):
+AveragingVessel(da)
+{
+  parseVector("PERIODIC",domain);
+  plumed_assert( domain.size()==2 || domain.size()==0 );
+}
+
+void AverageVessel::resize(){
+  resizeBuffer(0); 
+  if( domain.size()==2 ) setDataSize(2);
+  else setDataSize(1);
+}
+
+void AverageVessel::accumulate( const double& weight, const double& val ){
+  if( domain.size()==2 ){
+     // Average with Berry Phase 
+     double tval = 2*pi*( val - domain[0] ) / ( domain[1] - domain[0] );
+     addDataElement( 0, weight*sin(tval) ); addDataElement( 1, weight*cos(tval) ); 
+  } else addDataElement( 0, weight*val );  
+}
+
+double AverageVessel::getAverage() const {
+  if( domain.size()==2 ) return domain[0] + (( domain[1] - domain[0] )*atan2( getDataElement(0), getDataElement(1) ) / (2*pi));
+  return getDataElement(0);
+}
+
+void AverageVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
+  plumed_error();
+}
+
+}
+}

src/analysis/AverageVessel.h

+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2012-2016 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_analysis_AverageVessel_h
+#define __PLUMED_analysis_AverageVessel_h
+
+#include "vesselbase/AveragingVessel.h"
+
+namespace PLMD {
+namespace analysis {
+
+class AverageVessel : public vesselbase::AveragingVessel {
+private:
+ std::vector<double> domain;
+public:
+  /// keywords
+ static void registerKeywords( Keywords& keys );
+/// Constructor
+ explicit AverageVessel( const vesselbase::VesselOptions& );
+/// Set the size of the data vessel
+ void resize();
+/// This does nothing
+ void calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
+ std::string description(){ return ""; } 
+/// Accumulate the average
+ void accumulate( const double& weight, const double& val );
+/// Get the average value
+ double getAverage() const ;
+};
+
+}
+}
+#endif

src/gridtools/ConvertToFES.cpp

@@ -42,7 +42,6 @@ private:
 public:
   static void registerKeywords( Keywords& keys );
   explicit ConvertToFES(const ActionOptions&ao);
-  unsigned getNumberOfDerivatives(){ return 0; }
   unsigned getNumberOfQuantities() const ;
   void compute( const unsigned& current, MultiValue& myvals ) const ;
   bool isPeriodic(){ return false; }

src/gridtools/FourierTransform.cpp

@@ -80,7 +80,6 @@ public:
 #else
   void performOperations( const bool& from_update );
 #endif
-  unsigned getNumberOfDerivatives(){ return 0; }
   void compute( const unsigned& , MultiValue& ) const {}
   bool isPeriodic(){ return false; }
 };

src/gridtools/InterpolateGrid.cpp

@@ -40,7 +40,6 @@ class InterpolateGrid : public ActionWithInputGrid {
 public:
   static void registerKeywords( Keywords& keys );
   explicit InterpolateGrid(const ActionOptions&ao);
-  unsigned getNumberOfDerivatives(){ return 0; }
   unsigned getNumberOfQuantities() const ;
   void compute( const unsigned& current, MultiValue& myvals ) const ;
   bool isPeriodic(){ return false; }

src/multicolvar/MultiColvarDensity.cpp

@@ -100,7 +100,6 @@ public:
   static void registerKeywords( Keywords& keys );
   unsigned getNumberOfQuantities() const ;
   bool isPeriodic(){ return false; }
-  unsigned getNumberOfDerivatives(){ return 0; }
   void clearAverage();
   void prepareForAveraging();
   void compute( const unsigned& , MultiValue& ) const ;

src/vesselbase/ActionWithAveraging.cpp

@@ -27,13 +27,13 @@ namespace PLMD {
 namespace vesselbase {
 
 void ActionWithAveraging::registerKeywords( Keywords& keys ){
-  Action::registerKeywords( keys ); ActionPilot::registerKeywords( keys );
-  ActionAtomistic::registerKeywords( keys ); ActionWithArguments::registerKeywords( keys ); ActionWithVessel::registerKeywords( keys );
+  Action::registerKeywords( keys ); ActionPilot::registerKeywords( keys ); ActionAtomistic::registerKeywords( keys ); 
+  ActionWithArguments::registerKeywords( keys ); ActionWithValue::registerKeywords( keys ); ActionWithVessel::registerKeywords( keys );
   keys.add("compulsory","STRIDE","1","the frequency with which the data should be collected and added to the quantity being averaged");
   keys.add("compulsory","CLEAR","0","the frequency with which to clear all the accumulated data.  The default value "
                                     "of 0 implies that all the data will be used and that the grid will never be cleared");
   keys.add("optional","LOGWEIGHTS","list of actions that calculates log weights that should be used to weight configurations when calculating averages");
-  keys.addFlag("UNORMALIZED",false,"output the unaveraged quantity/quantities.");
+  keys.addFlag("UNORMALIZED",false,"output the unaveraged quantity/quantities."); keys.remove("NUMERICAL_DERIVATIVES");
 }
 
 ActionWithAveraging::ActionWithAveraging( const ActionOptions& ao ):
@@ -41,6 +41,7 @@ Action(ao),
 ActionPilot(ao),
 ActionAtomistic(ao),
 ActionWithArguments(ao),
+ActionWithValue(ao),
 ActionWithVessel(ao),
 myaverage(NULL),
 useRunAllTasks(false),
@@ -55,7 +56,7 @@ clearstride(0)
   }
   if( keywords.exists("LOGWEIGHTS") ){
       std::vector<std::string> wwstr; parseVector("LOGWEIGHTS",wwstr);
-      log.printf("  reweighting using weights from ");
+      if( wwstr.size()>0 ) log.printf("  reweighting using weights from ");
       std::vector<Value*> arg( getArguments() );
       for(unsigned i=0;i<wwstr.size();++i){
           ActionWithValue* val = plumed.getActionSet().selectWithLabel<ActionWithValue*>(wwstr[i]);
@@ -64,7 +65,9 @@ clearstride(0)
           arg.push_back( val->copyOutput(val->getLabel()) );
           log.printf("%s ",wwstr[i].c_str() );
       }
-     log.printf("\n"); requestArguments( arg );
+      if( wwstr.size()>0 ) log.printf("\n"); 
+      else log.printf("  weights are all equal to one\n");
+      requestArguments( arg );
   }
   if( keywords.exists("UNORMALIZED") ) parseFlag("UNORMALIZED",unormalised);
 }

src/vesselbase/ActionWithAveraging.h

@@ -25,6 +25,7 @@
 #include "core/ActionPilot.h"
 #include "core/ActionWithValue.h"
 #include "core/ActionAtomistic.h"
+#include "core/ActionWithValue.h"
 #include "core/ActionWithArguments.h"
 #include "ActionWithVessel.h"
 #include "AveragingVessel.h"
@@ -45,6 +46,7 @@ class ActionWithAveraging :
   public ActionPilot,
   public ActionAtomistic,
   public ActionWithArguments,
+  public ActionWithValue,
   public ActionWithVessel
 {
 friend class AveragingVessel;
@@ -70,6 +72,7 @@ public:
   void lockRequests();
   void unlockRequests();
   void calculateNumericalDerivatives(PLMD::ActionWithValue*);
+  virtual unsigned getNumberOfDerivatives(){ return 0; }
   unsigned getNumberOfArguments() const ;
   void calculate(){}
   void apply(){}