Began to re-add normalisation stuff

src/gridtools/ActionWithInputGrid.cpp

@@ -32,11 +32,13 @@ void ActionWithInputGrid::registerKeywords( Keywords& keys ){
   keys.add("compulsory","GRID","the action that creates the input grid you would like to use");
   keys.add("compulsory","STRIDE","the frequency with which to output the grid");
   keys.addFlag("USE_ALL_DATA",false,"use the data from the entire trajectory to perform the analysis");
+  keys.addFlag("UNORMALISED",false,"output an unormalised grid"); 
 }
 
 ActionWithInputGrid::ActionWithInputGrid(const ActionOptions&ao):
 Action(ao),
 ActionPilot(ao),
+norm(0.0),
 mygrid(NULL)
 {
   std::string mlab; parse("GRID",mlab);
@@ -51,6 +53,8 @@ mygrid(NULL)
   }
   if( !mygrid ) error("input action does not calculate a grid");
 
+  parseFlag("UNORMALISED",unormalised);
+  if( unormalised ) log.printf("  working with unormalised grid \n");
   if( keywords.exists("USE_ALL_DATA") ){
      parseFlag("USE_ALL_DATA",single_run);
      if( !single_run ) log.printf("  outputting grid every %u steps \n", getStride() );
@@ -59,10 +63,14 @@ mygrid(NULL)
 }
 
 void ActionWithInputGrid::update(){
+  if( unormalised ) norm = 1.0;
+  else norm=1.0/mygrid->getNorm();
   performOperationsWithGrid( true );
 }
 
 void ActionWithInputGrid::runFinalJobs(){
+  if( unormalised ) norm = 1.0;
+  else norm=1.0/mygrid->getNorm();
   performOperationsWithGrid( false );
 }
 

src/gridtools/ActionWithInputGrid.h

@@ -30,9 +30,14 @@ namespace gridtools {
 
 class ActionWithInputGrid : 
 public ActionPilot {
+private:
+  double norm;
+  bool unormalised;
 protected:
   bool single_run;
   GridVessel* mygrid;
+  double getGridElement( const std::vector<unsigned>& pp, const unsigned& ind ) const ;
+  double getGridElement( const unsigned& ind, const unsigned& jind ) const ;
 public:
   static void registerKeywords( Keywords& keys );
   explicit ActionWithInputGrid(const ActionOptions&ao);
@@ -43,6 +48,16 @@ public:
   virtual void performOperationsWithGrid( const bool& from_update )=0;
 };
 
+inline
+double ActionWithInputGrid::getGridElement( const std::vector<unsigned>& pp, const unsigned& ind ) const {
+  return norm*mygrid->getGridElement( pp, ind );
+}
+
+inline
+double ActionWithInputGrid::getGridElement( const unsigned& ind, const unsigned& jind ) const {
+  return norm*mygrid->getGridElement( ind, jind );
+}
+
 }
 }
 #endif

src/gridtools/GridVessel.h

@@ -32,8 +32,6 @@ namespace gridtools {
 
 class GridVessel : public vesselbase::Vessel {
 private:
-/// Are we deleting the data after print
- bool nomemory;
 /// Have the minimum and maximum for the grid been set
  bool bounds_set;
 /// These two variables are used to 
@@ -60,6 +58,8 @@ private:
 /// The grid point that was requested last by getGridPointCoordinates
  unsigned currentGridPoint;
 protected:
+/// Are we deleting the data after print
+ bool nomemory;
 /// The number of pieces of information we are storing for each 
 /// point in the grid
  unsigned nper;
@@ -132,11 +132,13 @@ public:
 /// Get the extent of the grid in one of the axis
  double getGridExtent( const unsigned& i ) const ;
 /// Clear all the data stored on the grid
- void clear();
+ virtual void clear();
 /// This ensures that Gaussian cube fies are in correct units
  void setCubeUnits( const double& units );
 /// This ensures that Gaussian cube files are in correct units
  double getCubeUnits() const ;
+/// This gives the normalisation of histograms
+ virtual double getNorm() const ;
 };
 
 inline
@@ -185,6 +187,11 @@ double GridVessel::getGridExtent( const unsigned& i ) const {
   return max[i] - min[i];
 }
 
+inline
+double GridVessel::getNorm() const {
+  return 1.0;
+}
+
 }
 }
 #endif

src/gridtools/HistogramOnGrid.cpp

@@ -33,6 +33,7 @@ void HistogramOnGrid::registerKeywords( Keywords& keys ){
 
 HistogramOnGrid::HistogramOnGrid( const vesselbase::VesselOptions& da ):
 GridVessel(da),
+norm(0),
 bandwidths(dimension)
 {
   parse("KERNEL",kerneltype); parseVector("BANDWIDTH",bandwidths);
@@ -78,5 +79,10 @@ void HistogramOnGrid::finish( const std::vector<double>& buffer ){
   for(unsigned i=0;i<data.size();++i) data[i]+=buffer[bufstart + i];
 }
 
+void HistogramOnGrid::clear(){
+  if( !nomemory ) return ;
+  norm = 0.; GridVessel::clear();
+}
+
 }
 }

src/gridtools/HistogramOnGrid.h

@@ -29,6 +29,7 @@ namespace gridtools {
 
 class HistogramOnGrid : public GridVessel {
 private:
+  double norm;
   std::string kerneltype;
   std::vector<double> bandwidths;
   std::vector<unsigned> nneigh;
@@ -40,8 +41,27 @@ public:
   bool calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
   void finish( const std::vector<double>& );
   bool applyForce(  std::vector<double>& forces ){ return false; }
+  void addToNorm( const double& anorm );
+  void setNorm( const double& snorm );
+  double getNorm() const ;
+  void clear();
 };
 
+inline
+void HistogramOnGrid::addToNorm( const double& anorm ){
+  norm+=anorm;
+}
+
+inline
+void HistogramOnGrid::setNorm( const double& snorm ){
+  norm=snorm;
+}
+
+inline
+double HistogramOnGrid::getNorm() const {
+  return norm;
+}
+
 }
 }
 #endif

src/gridtools/PrintCube.cpp

@@ -78,7 +78,7 @@ void PrintCube::performOperationsWithGrid( const bool& from_update ){
   for(pp[0]=0;pp[0]<nbin[0];++pp[0]){
       for(pp[1]=0;pp[1]<nbin[1];++pp[1]){
           for(pp[2]=0;pp[2]<nbin[2];++pp[2]){
-              ofile.printf("%f ",mygrid->getGridElement( pp, 0 ) );
+              ofile.printf("%f ",getGridElement( pp, 0 ) );
               if(pp[2]%6==5) ofile.printf("\n");
           }
           ofile.printf("\n");

src/gridtools/PrintGrid.cpp

@@ -85,7 +85,7 @@ void PrintGrid::performOperationsWithGrid( const bool& from_update ){
      mygrid->getGridPointCoordinates(i, xx ); 
      for(unsigned j=0;j<mygrid->getDimension();++j){ ofile.fmtField(fmt); ofile.printField(mygrid->getComponentName(j),xx[j]); }
      for(unsigned j=0;j<mygrid->getNumberOfQuantities();++j){ 
-        ofile.fmtField(fmt); ofile.printField(mygrid->getComponentName(mygrid->getDimension()+j), mygrid->getGridElement( i, j ) ); 
+        ofile.fmtField(fmt); ofile.printField(mygrid->getComponentName(mygrid->getDimension()+j), getGridElement( i, j ) ); 
      }
      ofile.printField();
   }

src/multicolvar/MultiColvarDensity.cpp

@@ -59,7 +59,6 @@ class MultiColvarDensity :
   public vesselbase::ActionWithInputVessel
 {
   std::string kerneltype;
-  double norm;
   bool firststep;
   bool fractional;
   unsigned rstride;
@@ -106,7 +105,6 @@ MultiColvarDensity::MultiColvarDensity(const ActionOptions&ao):
   ActionAtomistic(ao),
   ActionWithVessel(ao),
   ActionWithInputVessel(ao),
-  norm(0),
   firststep(true)
 {
 
@@ -217,31 +215,7 @@ void MultiColvarDensity::update(){
   }
   // Now perform All Tasks
   origin = getPosition(0);
-  runAllTasks(); norm += 1.0;
-
-//  // Output and reset the counter if required
-//  if( getStep()%rstride==0 ){  // && getStep()>0 ){
-//      // Normalise prior to output
-//      gg->scaleAllValuesAndDerivatives( 1.0 / norm );
-//
-//      OFile gridfile; gridfile.link(*this); gridfile.setBackupString("analysis");
-//      gridfile.open( filename ); 
-//      if( cube ){
-//         // Cube files are in au so I convert from "Angstrom" to AU so that when
-//         // VMD converts this number back to Angstroms (from AU) it looks right
-//         if( plumed.getAtoms().usingNaturalUnits() ) gg->writeCubeFile( gridfile, 1.0/0.5292 );  
-//         else gg->writeCubeFile( gridfile, plumed.getAtoms().getUnits().getLength()/.05929 );
-//      } else gg->writeToFile( gridfile ); 
-//      gridfile.close();
-//
-//      if( nomemory ){ 
-//        gg->clear(); norm=0.0; 
-//      } else {
-//        // Unormalise after output
-//        gg->scaleAllValuesAndDerivatives( norm );
-//      }
-//  }
-
+  runAllTasks(); mygrid->addToNorm( 1.0 ); 
 }
 
 void MultiColvarDensity::performTask( const unsigned& tindex, const unsigned& current, MultiValue& myvals ) const {