From 12dd30350d43947d1306ffb4a0474764dd41a976 Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Mon, 2 May 2016 20:40:04 +0100
Subject: [PATCH] Using new weight actions in analysis

---
 .../rt-mds/analysis.0.list_embed.reference    |   2 +-
 regtest/analysis/rt-mds/list_embed.reference  |   2 +-
 regtest/analysis/rt-mds/plumed.dat            |   5 +-
 src/analysis/Analysis.cpp                     | 182 ++++--------------
 src/analysis/Analysis.h                       |  74 +------
 src/bias/ReweightTemperature.cpp              |   6 +-
 src/vesselbase/ActionWithAveraging.cpp        |   5 +-
 src/vesselbase/ActionWithAveraging.h          |  13 +-
 src/vesselbase/ActionWithVessel.h             |   2 -
 9 files changed, 71 insertions(+), 220 deletions(-)

diff --git a/regtest/analysis/rt-mds/analysis.0.list_embed.reference b/regtest/analysis/rt-mds/analysis.0.list_embed.reference
index a5fdcd999..a0932b381 100644
--- a/regtest/analysis/rt-mds/analysis.0.list_embed.reference
+++ b/regtest/analysis/rt-mds/analysis.0.list_embed.reference
@@ -1,4 +1,4 @@
-#! FIELDS @17.1 @17.2
+#! FIELDS @18.1 @18.2
   0.0873   0.0013 
   0.0535   0.0044 
   0.0150   0.0044 
diff --git a/regtest/analysis/rt-mds/list_embed.reference b/regtest/analysis/rt-mds/list_embed.reference
index 6f40ac6e8..e653af32a 100644
--- a/regtest/analysis/rt-mds/list_embed.reference
+++ b/regtest/analysis/rt-mds/list_embed.reference
@@ -1,4 +1,4 @@
-#! FIELDS @17.1 @17.2
+#! FIELDS @18.1 @18.2
   0.0594   0.0053 
   0.0683   0.0034 
   0.0705   0.0006 
diff --git a/regtest/analysis/rt-mds/plumed.dat b/regtest/analysis/rt-mds/plumed.dat
index 385718a58..7fe2ce0ee 100755
--- a/regtest/analysis/rt-mds/plumed.dat
+++ b/regtest/analysis/rt-mds/plumed.dat
@@ -16,13 +16,16 @@ DISTANCE ATOMS=7,com LABEL=d7
 UPPER_WALLS ARG=d7 AT=2.0 KAPPA=100.
 
 COORDINATIONNUMBER SPECIES=1-7 MOMENTS=2-3 SWITCH={RATIONAL R_0=1.5 NN=8} LABEL=c1
+
+ww: REWEIGHT_TEMP REWEIGHT_TEMP=0.1 TEMP=0.2
+
 CLASSICAL_MDS ... 
   ARG=c1.moment-2,c1.moment-3 
-  REWEIGHT_TEMP=0.1 TEMP=0.2
   STRIDE=10 
   RUN=1000
   NLOW_DIM=2
   FMT=%8.4f
   OUTPUT_FILE=list_embed
   EMBEDDING_OFILE=embed
+  LOGWEIGHTS=ww
 ... CLASSICAL_MDS
diff --git a/src/analysis/Analysis.cpp b/src/analysis/Analysis.cpp
index a05f92f63..8dd4a7866 100644
--- a/src/analysis/Analysis.cpp
+++ b/src/analysis/Analysis.cpp
@@ -57,47 +57,28 @@ P(s',t) = \frac{ \sum_{t'}^t \delta( s(x) - s' ) \exp\left( +\frac{V(x,t')}{k_B
 //+ENDPLUMEDOC
 
 void Analysis::registerKeywords( Keywords& keys ){
-  Action::registerKeywords( keys );
-  ActionPilot::registerKeywords( keys );
-  ActionAtomistic::registerKeywords( keys );
-  ActionWithArguments::registerKeywords( keys );
+  vesselbase::ActionWithAveraging::registerKeywords( keys );
   keys.use("ARG"); keys.reset_style("ARG","optional");
   keys.add("atoms","ATOMS","the atoms whose positions we are tracking for the purpose of analysing the data");
   keys.add("compulsory","METRIC","EUCLIDEAN","how are we measuring the distances between configurations");
-  keys.add("compulsory","STRIDE","1","the frequency with which data should be stored for analysis");
-  keys.addFlag("USE_ALL_DATA",false,"use the data from the entire trajectory to perform the analysis");
-  keys.add("compulsory","RUN","the frequency with which to run the analysis algorithm. This is not required if you specify USE_ALL_DATA");
+  keys.add("compulsory","RUN","0","the frequency with which to run the analysis algorithm. The default value of zero assumes you want to analyse the whole trajectory");
   keys.add("optional","FMT","the format that should be used in analysis output files");
-  keys.addFlag("REWEIGHT_BIAS",false,"reweight the data using all the biases acting on the dynamics. For more information see \\ref reweighting.");
-  keys.add("optional","TEMP","the system temperature.  This is required if you are reweighting or doing free energies.");
-  keys.add("optional","REWEIGHT_TEMP","reweight data from a trajectory at one temperature and output the probability "
-                                      "distribution at a second temperature. For more information see \\ref reweighting. "
-                                      "This is not possible during postprocessing.");
   keys.addFlag("WRITE_CHECKPOINT",false,"write out a checkpoint so that the analysis can be restarted in a later run");
   keys.add("hidden","REUSE_DATA_FROM","eventually this will allow you to analyse the same set of data multiple times");
   keys.add("hidden","IGNORE_REWEIGHTING","this allows you to ignore any reweighting factors");
-  keys.use("RESTART");
-  keys.use("UPDATE_FROM");
-  keys.use("UPDATE_UNTIL");
-  ActionWithVessel::registerKeywords( keys ); keys.remove("TOL"); 
+  keys.use("RESTART"); keys.use("UPDATE_FROM"); keys.use("UPDATE_UNTIL"); keys.remove("TOL"); 
 }
 
 Analysis::Analysis(const ActionOptions&ao):
 Action(ao),
-ActionPilot(ao),
-ActionAtomistic(ao),
-ActionWithArguments(ao),
-ActionWithVessel(ao),
-single_run(true),
+ActionWithAveraging(ao),
 nomemory(true),
 write_chq(false),
 reusing_data(false),
 ignore_reweight(false),
-freq(0),
-needeng(false),
 idata(0),
-firstAnalysisDone(false),
-old_norm(0.0),
+//firstAnalysisDone(false),
+//old_norm(0.0),
 ofmt("%f"),
 current_args(getNumberOfArguments()),
 argument_names(getNumberOfArguments())
@@ -132,59 +113,18 @@ argument_names(getNumberOfArguments())
       if( ignore_reweight ) log.printf("  reusing data stored by %s but ignoring all reweighting\n",prev_analysis.c_str() );
       else log.printf("  reusing data stored by %s\n",prev_analysis.c_str() ); 
   } else { 
-      if( keywords.exists("REWEIGHT_BIAS") ){
-         bool dobias; parseFlag("REWEIGHT_BIAS",dobias);
-         if( dobias ){
-             std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>();
-             if( all.empty() ) error("your input file is not telling plumed to calculate anything");
-             std::vector<Value*> arg( getArguments() );
-             log.printf("  reweigting using the following biases ");
-             for(unsigned j=0;j<all.size();j++){
-                 std::string flab; flab=all[j]->getLabel() + ".rbias";
-                 if( all[j]->exists(flab) ){ 
-                    biases.push_back( all[j]->copyOutput(flab) ); 
-                    arg.push_back( all[j]->copyOutput(flab) ); 
-                    log.printf(" %s",flab.c_str()); 
-                 } else {
-                     std::string flab2; flab2=all[j]->getLabel() + ".bias";
-                     if( all[j]->exists(flab2) ){
-                        biases.push_back( all[j]->copyOutput(flab2) );
-                        arg.push_back( all[j]->copyOutput(flab2) );
-                        log.printf(" %s",flab2.c_str());
-                     }
-                 } 
-             }
-             log.printf("\n");
-             if( biases.empty() ) error("you are asking to reweight bias but there does not appear to be a bias acting on your system");
-             requestArguments( arg ); 
-         }
-      }
-
-      rtemp=0;      
-      if( keywords.exists("REWEIGHT_TEMP") ) parse("REWEIGHT_TEMP",rtemp);
-      if( rtemp!=0 ){
-         needeng=true;
-         log.printf("  reweighting simulation to probabilities at temperature %f\n",rtemp);
-         rtemp*=plumed.getAtoms().getKBoltzmann(); 
+      parse("RUN",freq); 
+      if( freq==0 ){
+          log.printf("  analyzing all data in trajectory\n");
+      } else {
+          if( freq%getStride()!=0 ) error("frequncy of running is not a multiple of the stride");
+          log.printf("  running analysis every %u steps\n",freq);
+          ndata=freq/getStride(); data.resize( ndata ); logweights.resize( ndata );
+          for(unsigned i=0;i<ndata;++i){
+             data[i]=metricRegister().create<ReferenceConfiguration>( metricname );
+             data[i]->setNamesAndAtomNumbers( getAbsoluteIndexes(), argument_names );
+          } 
       } 
-      simtemp=0.; parse("TEMP",simtemp);
-      if(simtemp>0) simtemp*=plumed.getAtoms().getKBoltzmann();
-      else simtemp=plumed.getAtoms().getKbT();
-
-      if( rtemp>0 || !biases.empty() ){
-         if(simtemp==0) error("The MD engine does not pass the temperature to plumed so you have to specify it using TEMP");
-      }
-
-      if( keywords.exists("USE_ALL_DATA") ){
-          parseFlag("USE_ALL_DATA",single_run); 
-          if( !single_run ){
-              unsigned astride; parse("RUN",astride);
-              log.printf("  running analysis every %u steps\n",astride);
-              setAnalysisStride( false, astride );
-          } else {       
-              log.printf("  analyzing all data in trajectory\n");
-          }
-      }
       parseFlag("WRITE_CHECKPOINT",write_chq);
       if( write_chq && single_run ){
           write_chq=false;
@@ -208,26 +148,12 @@ argument_names(getNumberOfArguments())
           rfile.open( filename.c_str() );  // In overwrite mode automatically because there is no restart
       }
       if( write_chq ){
-         rfile.addConstantField("old_normalization");
+         //rfile.addConstantField("old_normalization");
          for(unsigned i=0;i<getNumberOfArguments();++i) rfile.setupPrintValue( getPntrToArgument(i) );
       }
   }
 }
 
-void Analysis::setAnalysisStride( const bool& use_all, const unsigned& astride ){
-  if( freq>0 && astride!=freq ) error("each histogram can only be output with one stride"); 
-  else if( use_all ) return;  
-
-  freq=astride; single_run=false;
-  if( astride%getStride()!= 0 ) error("Frequncy of running is not a multiple of the stride");
-  ndata=freq/getStride(); data.resize( ndata );
-  for(unsigned i=0;i<ndata;++i){
-     data[i]=metricRegister().create<ReferenceConfiguration>( metricname );
-     data[i]->setNamesAndAtomNumbers( getAbsoluteIndexes(), argument_names );
-  }
-  logweights.resize( ndata );
-}
-
 void Analysis::readDataFromFile( const std::string& filename ){
   FILE* fp=fopen(filename.c_str(),"r");
   if(fp!=NULL){
@@ -243,7 +169,7 @@ void Analysis::readDataFromFile( const std::string& filename ){
               error("frequency of data storage in " + filename + " is not equal to frequency of data storage plumed.dat file");
            }
            data[idata]->parse("LOG_WEIGHT",logweights[idata]);
-           data[idata]->parse("OLD_NORM",old_norm);
+           //data[idata]->parse("OLD_NORM",old_norm);
            data[idata]->checkRead();
            idata++; first=false; oldtstep=tstep;
         } else{
@@ -252,7 +178,7 @@ void Analysis::readDataFromFile( const std::string& filename ){
      }
     fclose(fp);
   }
-  if(old_norm>0) firstAnalysisDone=true;
+  // if(old_norm>0) firstAnalysisDone=true;
 }
 
 void Analysis::parseOutputFile( const std::string& key, std::string& filename ){
@@ -266,32 +192,11 @@ void Analysis::parseOutputFile( const std::string& key, std::string& filename ){
   } 
 }
 
-void Analysis::prepare(){
-  if(needeng) plumed.getAtoms().setCollectEnergy(true);
-}
-
-void Analysis::calculate(){
-// do nothing
-}
-
 void Analysis::accumulate(){
   // Don't store the first step (also don't store if we are getting data from elsewhere)
   if( (!single_run && getStep()==0) || reusing_data ) return;
   // This is used when we have a full quota of data from the first run
   if( !single_run && idata==logweights.size() ) return; 
-
-  // Retrieve the bias
-  double bias=0.0; for(unsigned i=0;i<biases.size();++i) bias+=biases[i]->get();
-
-  double ww=0;
-  if(needeng){
-     double energy=plumed.getAtoms().getEnergy()+bias;
-     // Reweighting because of temperature difference
-     ww=-( (1.0/rtemp) - (1.0/simtemp) )*(energy+bias);
-  }
-  // Reweighting because of biases
-  if( !biases.empty() ) ww += bias/ simtemp;
-
   // Get the arguments ready to transfer to reference configuration
   for(unsigned i=0;i<getNumberOfArguments();++i) current_args[i]=getArgument(i);
 
@@ -300,17 +205,17 @@ void Analysis::accumulate(){
      plumed_dbg_assert( data.size()==idata+1 );
      data[idata]->setNamesAndAtomNumbers( getAbsoluteIndexes(), argument_names );
      data[idata]->setReferenceConfig( getPositions(), current_args, getMetric() );
-     logweights.push_back(ww);
+     logweights.push_back(lweight);
   } else {
      // Get the arguments and store them in a vector of vectors
      data[idata]->setReferenceConfig( getPositions(), current_args, getMetric() );
-     logweights[idata] = ww; 
+     logweights[idata] = lweight; 
   }
 
   // Write data to checkpoint file
   if( write_chq ){
      rfile.rewind();
-     data[idata]->print( rfile, getTime(), logweights[idata], atoms.getUnits().getLength()/0.1, old_norm );
+     data[idata]->print( rfile, getTime(), logweights[idata], atoms.getUnits().getLength()/0.1, 1.0 ); //old_norm );
      rfile.flush();
   }
   // Increment data counter
@@ -344,7 +249,7 @@ void Analysis::finalizeWeights( const bool& ignore_weights ){
   // Check that we have the correct ammount of data
   if( !reusing_data && idata!=logweights.size() ) error("something has gone wrong.  Am trying to run analysis but I don't have sufficient data");
 
-  norm=0;  // Reset normalization constant
+  double norm=0;  // Reset normalization constant
   if( ignore_weights ){
       for(unsigned i=0;i<logweights.size();++i){
           data[i]->setWeight(1.0); norm+=1.0;
@@ -365,18 +270,20 @@ void Analysis::finalizeWeights( const bool& ignore_weights ){
       }
   // Calculate normalized weights (with memory)
   } else {
-      if( !firstAnalysisDone ) finalizeWeightsNoLogSums( 1.0 );
-      else finalizeWeightsNoLogSums( old_norm );
+      plumed_merror("analysis can now only support block averages");
+      // if( !firstAnalysisDone ) 
+      // finalizeWeightsNoLogSums( 1.0 );
+      // else finalizeWeightsNoLogSums( old_norm );
   }
 }
 
-void Analysis::finalizeWeightsNoLogSums( const double& onorm ){
-  if( !reusing_data && idata!=logweights.size() ) error("something has gone wrong.  Am trying to run analysis but I don't have sufficient data");
-  // Calculate normalization constant
-  norm=0; for(unsigned i=0;i<logweights.size();++i) norm+=exp( logweights[i] );
-  // Calculate weights (with memory)
-  for(unsigned i=0;i<logweights.size();++i) data[i]->setWeight( exp( logweights[i] ) / onorm );
-}
+// void Analysis::finalizeWeightsNoLogSums( const double& onorm ){
+//   if( !reusing_data && idata!=logweights.size() ) error("something has gone wrong.  Am trying to run analysis but I don't have sufficient data");
+//   // Calculate normalization constant
+//   double norm=0; for(unsigned i=0;i<logweights.size();++i) norm+=exp( logweights[i] );
+//   // Calculate weights (with memory)
+//   for(unsigned i=0;i<logweights.size();++i) data[i]->setWeight( exp( logweights[i] ) / norm );
+// }
 
 void Analysis::getDataPoint( const unsigned& idata, std::vector<double>& point, double& weight ) const {
   plumed_dbg_assert( getNumberOfAtoms()==0 );
@@ -399,27 +306,18 @@ void Analysis::runAnalysis(){
      finalizeWeights( ignore_reweight ); 
   } else {
      mydatastash->finalizeWeights( ignore_reweight );
-     norm=mydatastash->retrieveNorm();
+     // norm=mydatastash->retrieveNorm();
   }
   // This ensures everything is set up to run the calculation
-  if( single_run ) setAnalysisStride( single_run, freq );
+  // if( single_run ) setAnalysisStride( single_run, freq );
   // And run the analysis
   performAnalysis(); idata=0;
   // Update total normalization constant
-  old_norm+=norm; firstAnalysisDone=true;
-
-}
-
-double Analysis::getNormalization() const {
-  if( nomemory || !firstAnalysisDone ) return norm;
-  return ( 1. + norm/old_norm );
-}
+  // old_norm+=norm; firstAnalysisDone=true;
 
-double Analysis::getTemp() const {
-  return simtemp;
 }
 
-void Analysis::update(){
+void Analysis::performOperations( const bool& from_update ){
   accumulate();
   if( !single_run ){
     if( getStep()>0 && getStep()%freq==0 ) runAnalysis(); 
@@ -434,7 +332,7 @@ bool Analysis::getPeriodicityInformation(const unsigned& i, std::string& dmin, s
 }
 
 void Analysis::runFinalJobs() {
-  if( !single_run ) return;
+  if( freq>0 ) return;
   if( getNumberOfDataPoints()==0 ) error("no data is available for analysis");
   runAnalysis(); 
 }
diff --git a/src/analysis/Analysis.h b/src/analysis/Analysis.h
index ac716c064..0a59bb770 100644
--- a/src/analysis/Analysis.h
+++ b/src/analysis/Analysis.h
@@ -22,10 +22,7 @@
 #ifndef __PLUMED_analysis_Analysis_h
 #define __PLUMED_analysis_Analysis_h
 
-#include "core/ActionPilot.h"
-#include "core/ActionAtomistic.h"
-#include "core/ActionWithArguments.h"
-#include "vesselbase/ActionWithVessel.h"
+#include "vesselbase/ActionWithAveraging.h"
 
 #define PLUMED_ANALYSIS_INIT(ao) Action(ao),Analysis(ao)
 
@@ -42,12 +39,7 @@ is information as to how to go about implementing a new analysis method.
 
 */
 
-class Analysis :
-  public ActionPilot,
-  public ActionAtomistic,
-  public ActionWithArguments,
-  public vesselbase::ActionWithVessel
-  {
+class Analysis : public vesselbase::ActionWithAveraging {
 private:
 /// Are we running only once for the whole trajectory
   bool single_run;
@@ -71,16 +63,14 @@ private:
   double rtemp;
 /// Do we need the energy (are we reweighting at a different temperature)
   bool needeng;
-/// The biases we are using in reweighting and the args we store them separately
-  std::vector<Value*> biases;
 /// The piece of data we are inserting
   unsigned idata;
 /// The weights of all the data points
   std::vector<double> logweights;
 /// Have we analyzed the data for the first time
-  bool firstAnalysisDone;
+//  bool firstAnalysisDone;
 /// The value of the old normalization constant
-  double norm, old_norm;
+//  double norm, old_norm;
 /// The format to use in output files
   std::string ofmt;
 /// Tempory vector to store values of arguments
@@ -91,11 +81,6 @@ private:
   OFile rfile;
 /// Read in data from a file
   void readDataFromFile( const std::string& filename );
-/// This retrieves the value of norm from the analysis action.
-/// It is only used to transfer data from one analysis action to
-/// another. You should never need to use it.  If you think you need it
-/// you probably need getNormalization()
-  double retrieveNorm() const ;
 /// Get the metric if we are using malonobius distance and flexible hill
   std::vector<double> getMetric() const ;
 protected:
@@ -109,8 +94,6 @@ protected:
   std::vector<ReferenceConfiguration*> data;
 /// Get the name of the metric we are using to measure distances
   std::string getMetricName() const ;
-/// Return the number of arguments (this overwrites the one in ActionWithArguments)
-  unsigned getNumberOfArguments() const;
 /// Return the number of data points
   unsigned getNumberOfDataPoints() const;
 /// Return the weight of the ith point
@@ -119,37 +102,22 @@ protected:
   void getDataPoint( const unsigned& idata, std::vector<double>& point, double& weight ) const ;
 /// Returns true if argument i is periodic together with the domain 
   bool getPeriodicityInformation(const unsigned& i, std::string& dmin, std::string& dmax);
-/// Return the normalization constant
-  double getNormalization() const;
-/// Return the set temperature (N.B. k_B T is what is returned by this function)
-  double getTemp () const;
 /// Are we analyzing each data block separately (if we are not this also returns the old normalization )
   bool usingMemory() const; 
-/// Overwrite ActionWithArguments getArguments() so that we don't return
-/// the bias
-  std::vector<Value*> getArguments();
 /// Return the format to use for numbers in output files
   std::string getOutputFormat() const ;
 /// Finalize the weights without using the log sums
-  void finalizeWeightsNoLogSums( const double& onorm );
+//  void finalizeWeightsNoLogSums( const double& onorm );
 public:
   static void registerKeywords( Keywords& keys );
   explicit Analysis(const ActionOptions&);
   ~Analysis();
-  void prepare();
-  void calculate();
-  void update();
   void accumulate();
+  void performOperations( const bool& from_update );
   virtual void performAnalysis()=0;
-  void apply(){}
   void runFinalJobs();
   void runAnalysis();
-  void lockRequests();
-  void unlockRequests();
-  void calculateNumericalDerivatives( ActionWithValue* a=NULL ){ plumed_error(); }
   bool isPeriodic(){ plumed_error(); return false; }
-  unsigned getNumberOfDerivatives(){ return 0; }
-  virtual void setAnalysisStride( const bool& use_all, const unsigned& astride );
   /// Convert the stored log weights to proper weights
   virtual void finalizeWeights( const bool& ignore_weights );
 };
@@ -159,18 +127,6 @@ std::string Analysis::getMetricName() const {
   return metricname;
 }
 
-inline 
-void Analysis::lockRequests(){
-  ActionAtomistic::lockRequests();
-  ActionWithArguments::lockRequests();
-} 
-
-inline
-void Analysis::unlockRequests(){ 
-  ActionAtomistic::unlockRequests();
-  ActionWithArguments::unlockRequests();
-}
-
 inline
 unsigned Analysis::getNumberOfDataPoints() const {
   if( !reusing_data ){
@@ -190,24 +146,6 @@ bool Analysis::usingMemory() const {
   }
 }
 
-inline
-unsigned Analysis::getNumberOfArguments() const {
-  unsigned nargs=ActionWithArguments::getNumberOfArguments();
-  return nargs - biases.size(); 
-}
-
-inline
-double Analysis::retrieveNorm() const {
-  return norm;
-}
-
-inline
-std::vector<Value*> Analysis::getArguments(){
-  std::vector<Value*> arg_vals( ActionWithArguments::getArguments() );
-  for(unsigned i=0;i<biases.size();++i) arg_vals.erase(arg_vals.end()-1);
-  return arg_vals;
-}
-
 inline
 std::string Analysis::getOutputFormat() const {
   return ofmt;
diff --git a/src/bias/ReweightTemperature.cpp b/src/bias/ReweightTemperature.cpp
index e503ebfd7..2d25876c7 100644
--- a/src/bias/ReweightTemperature.cpp
+++ b/src/bias/ReweightTemperature.cpp
@@ -44,6 +44,7 @@ private:
 public:
   static void registerKeywords(Keywords&);
   explicit ReweightTemperature(const ActionOptions&ao);
+  void prepare();
   double getLogWeight() const ;
 };
 
@@ -59,7 +60,6 @@ ReweightTemperature::ReweightTemperature(const ActionOptions&ao):
 Action(ao),
 ReweightBase(ao)
 {
-   plumed.getAtoms().setCollectEnergy(true);
    parse("REWEIGHT_TEMP",rtemp);
    log.printf("  reweighting simulation to probabilities at temperature %f\n",rtemp);
    rtemp*=plumed.getAtoms().getKBoltzmann();
@@ -67,6 +67,10 @@ ReweightBase(ao)
    retrieveAllBiases( "bias", biases );
 }
 
+void ReweightTemperature::prepare(){
+   plumed.getAtoms().setCollectEnergy(true);
+}
+
 double ReweightTemperature::getLogWeight() const {
    // Retrieve the bias
    double bias=0.0; for(unsigned i=0;i<biases.size();++i) bias+=biases[i]->get(); 
diff --git a/src/vesselbase/ActionWithAveraging.cpp b/src/vesselbase/ActionWithAveraging.cpp
index abe50da5c..372bb272c 100644
--- a/src/vesselbase/ActionWithAveraging.cpp
+++ b/src/vesselbase/ActionWithAveraging.cpp
@@ -43,6 +43,7 @@ ActionAtomistic(ao),
 ActionWithArguments(ao),
 ActionWithValue(ao),
 ActionWithVessel(ao),
+lweight(0),cweight(0),
 myaverage(NULL),
 useRunAllTasks(false),
 clearstride(0)
@@ -100,9 +101,9 @@ void ActionWithAveraging::update(){
   // Calculate the weight for all reweighting
   if ( weights.size()>0 ){
      double sum=0; for(unsigned i=0;i<weights.size();++i) sum+=weights[i]->get();
-     cweight = exp( sum );
+     lweight=sum; cweight = exp( sum );
   } else {
-     cweight=1.0;
+     lweight=0; cweight=1.0;
   }
   // Prepare to do the averaging
   prepareForAveraging();
diff --git a/src/vesselbase/ActionWithAveraging.h b/src/vesselbase/ActionWithAveraging.h
index 8a91b2657..295085591 100644
--- a/src/vesselbase/ActionWithAveraging.h
+++ b/src/vesselbase/ActionWithAveraging.h
@@ -62,8 +62,8 @@ private:
 /// Are we accumulated the unormalized quantity
   bool unormalised;
 protected:
-/// The current weight
-  double cweight;
+/// The current weight and its logarithm
+  double lweight, cweight;
 /// Set the averaging action
   void setAveragingAction( AveragingVessel* av_vessel, const bool& usetasks );
 public:
@@ -74,6 +74,8 @@ public:
   void calculateNumericalDerivatives(PLMD::ActionWithValue*);
   virtual unsigned getNumberOfDerivatives(){ return 0; }
   unsigned getNumberOfArguments() const ;
+/// Overwrite ActionWithArguments getArguments() so that we don't return the bias
+  std::vector<Value*> getArguments();  
   void calculate(){}
   void apply(){}
   void update();
@@ -92,6 +94,13 @@ unsigned ActionWithAveraging::getNumberOfArguments() const {
   return ActionWithArguments::getNumberOfArguments() - weights.size();
 }
 
+inline
+std::vector<Value*> ActionWithAveraging::getArguments(){
+  std::vector<Value*> arg_vals( ActionWithArguments::getArguments() );
+  for(unsigned i=0;i<weights.size();++i) arg_vals.erase(arg_vals.end()-1);
+  return arg_vals;
+}
+
 }
 }
 #endif
diff --git a/src/vesselbase/ActionWithVessel.h b/src/vesselbase/ActionWithVessel.h
index 2abc5d441..6225bfe93 100644
--- a/src/vesselbase/ActionWithVessel.h
+++ b/src/vesselbase/ActionWithVessel.h
@@ -188,8 +188,6 @@ public:
   Vessel* getVesselWithName( const std::string& mynam );
 /// Does the weight have derivatives
   bool weightWithDerivatives() const ;
-/// This is used to set the stride for the output of histograms
-  virtual void setAnalysisStride( const bool& use_all, const unsigned& astride ){}
 /// Return the position in the current task list
   unsigned getPositionInCurrentTaskList( const unsigned& myind ) const ;
 };
-- 
GitLab