From b46193c835a869ee2330bc8274154232898e856a Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Sun, 27 Dec 2015 09:54:05 +0000
Subject: [PATCH] Began to re-add normalisation stuff
---
src/gridtools/ActionWithInputGrid.cpp | 8 ++++++++
src/gridtools/ActionWithInputGrid.h | 15 ++++++++++++++
src/gridtools/GridVessel.h | 13 +++++++++---
src/gridtools/HistogramOnGrid.cpp | 6 ++++++
src/gridtools/HistogramOnGrid.h | 20 ++++++++++++++++++
src/gridtools/PrintCube.cpp | 2 +-
src/gridtools/PrintGrid.cpp | 2 +-
src/multicolvar/MultiColvarDensity.cpp | 28 +-------------------------
8 files changed, 62 insertions(+), 32 deletions(-)
diff --git a/src/gridtools/ActionWithInputGrid.cpp b/src/gridtools/ActionWithInputGrid.cpp
index 633174c8e..3ec6ab5d8 100644
--- a/src/gridtools/ActionWithInputGrid.cpp
+++ b/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 );
}
diff --git a/src/gridtools/ActionWithInputGrid.h b/src/gridtools/ActionWithInputGrid.h
index 4902b07df..3016bee2c 100644
--- a/src/gridtools/ActionWithInputGrid.h
+++ b/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
diff --git a/src/gridtools/GridVessel.h b/src/gridtools/GridVessel.h
index 95239ef84..c0df8d173 100644
--- a/src/gridtools/GridVessel.h
+++ b/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
diff --git a/src/gridtools/HistogramOnGrid.cpp b/src/gridtools/HistogramOnGrid.cpp
index 6e2bcb479..050810e1f 100644
--- a/src/gridtools/HistogramOnGrid.cpp
+++ b/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();
+}
+
}
}
diff --git a/src/gridtools/HistogramOnGrid.h b/src/gridtools/HistogramOnGrid.h
index 5ff4d6a12..0bdae980e 100644
--- a/src/gridtools/HistogramOnGrid.h
+++ b/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
diff --git a/src/gridtools/PrintCube.cpp b/src/gridtools/PrintCube.cpp
index 913a0da9e..993bad2c3 100644
--- a/src/gridtools/PrintCube.cpp
+++ b/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");
diff --git a/src/gridtools/PrintGrid.cpp b/src/gridtools/PrintGrid.cpp
index eab6b8c7c..8ae5f491c 100644
--- a/src/gridtools/PrintGrid.cpp
+++ b/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();
}
diff --git a/src/multicolvar/MultiColvarDensity.cpp b/src/multicolvar/MultiColvarDensity.cpp
index 7614bf173..230e32746 100644
--- a/src/multicolvar/MultiColvarDensity.cpp
+++ b/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 {
--
GitLab