diff --git a/src/gridtools/ActionWithInputGrid.cpp b/src/gridtools/ActionWithInputGrid.cpp
index 3ec6ab5d8533084684228ac3aa2f4665f346d53f..f0dfe559e4b8ad68c44656e4d4d179ffd4a19b16 100644
--- a/src/gridtools/ActionWithInputGrid.cpp
+++ b/src/gridtools/ActionWithInputGrid.cpp
@@ -54,7 +54,11 @@ 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( unormalised ){
+ log.printf(" working with unormalised grid \n");
+ mygrid->switchOffNormalisation();
+ }
+
if( keywords.exists("USE_ALL_DATA") ){
parseFlag("USE_ALL_DATA",single_run);
if( !single_run ) log.printf(" outputting grid every %u steps \n", getStride() );
diff --git a/src/gridtools/GridVessel.cpp b/src/gridtools/GridVessel.cpp
index a6f110e79419ae23c42aeb705b2d6dee5b1f4056..9b424a516da1e82347ad48a2f26315afa565c67f 100644
--- a/src/gridtools/GridVessel.cpp
+++ b/src/gridtools/GridVessel.cpp
@@ -46,7 +46,6 @@ cube_units(1.0)
dimension=coordnames.size();
std::vector<std::string> spbc( dimension ); parseVector("PBC",spbc);
str_min.resize( dimension); str_max.resize( dimension );
- nbin.resize( dimension ); parseVector("NBINS",nbin);
parseFlag("NOMEMORY",nomemory);
tmp_indices.resize( str_min.size() );
@@ -60,30 +59,43 @@ cube_units(1.0)
for(unsigned j=0;j<coordnames.size();++j){ arg_names[n] = "d" + compnames[i] + "_" + coordnames[j]; n++; }
}
- npoints=1; dx.resize( dimension ); pbc.resize( dimension );
- min.resize( dimension ); stride.resize( dimension ); max.resize( dimension );
+ pbc.resize( dimension );
for(unsigned i=0;i<dimension;++i){
- if( spbc[i]=="F" ){ pbc[i]=false; nbin[i]+=1; }
+ if( spbc[i]=="F" ) pbc[i]=false;
else if( spbc[i]=="T" ) pbc[i]=true;
else plumed_error();
- stride[i]=npoints;
- npoints*=nbin[i];
}
}
-void GridVessel::setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax ){
- bounds_set=true; plumed_dbg_assert( smin.size()==dimension && smax.size()==dimension );
+void GridVessel::setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax,
+ const std::vector<unsigned>& binsin, const std::vector<double>& spacing ){
+ plumed_dbg_assert( smin.size()==dimension && smax.size()==dimension );
+ plumed_assert( spacing.size()==dimension || binsin.size()==dimension );
+
+ npoints=1; bounds_set=true;
+ stride.resize( dimension ); max.resize( dimension );
+ dx.resize( dimension ); nbin.resize( dimension ); min.resize( dimension );
for(unsigned i=0;i<dimension;++i){
str_min[i]=smin[i]; str_max[i]=smax[i];
Tools::convert( str_min[i], min[i] );
Tools::convert( str_max[i], max[i] );
- if( !pbc[i] ){ nbin[i]-=1; }
+ if( spacing.size()==dimension && binsin.size()==dimension ){
+ unsigned spc = std::floor(( max[i] - min[i] ) / spacing[i]) + 1;
+ if( spc>binsin[i] ) nbin[i]=spc;
+ else nbin[i]=binsin[i];
+ } else if( binsin.size()==dimension ) nbin[i]=binsin[i];
+ else if( spacing.size()==dimension ) nbin[i] = std::floor(( max[i] - min[i] ) / spacing[i]) + 1;
+ else plumed_error();
dx[i] = ( max[i] - min[i] ) / static_cast<double>( nbin[i] );
if( !pbc[i] ){ max[i] +=dx[i]; nbin[i]+=1; }
+ stride[i]=npoints;
+ npoints*=nbin[i];
}
}
-std::string GridVessel::getGridDescription() const {
+std::string GridVessel::description(){
+ if( !bounds_set ) return "";
+
std::string des="grid of "; std::string num;
for(unsigned i=0;i<dimension-1;++i){
Tools::convert( nbin[i], num );
diff --git a/src/gridtools/GridVessel.h b/src/gridtools/GridVessel.h
index c0df8d1734cd8be95c2c4ef137b6c9b04fe804b3..9dfdd8952a17fe0cf1e447fb70687669dc698cea 100644
--- a/src/gridtools/GridVessel.h
+++ b/src/gridtools/GridVessel.h
@@ -81,9 +81,9 @@ public:
/// Constructor
explicit GridVessel( const vesselbase::VesselOptions& );
/// Set the minimum and maximum of the grid
- virtual void setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax );
+ virtual void setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax, const std::vector<unsigned>& nbins, const std::vector<double>& spacing );
/// Get a description of the grid to output to the log
- std::string getGridDescription() const ;
+ std::string description();
/// Convert an index into indices
void convertIndexToIndices( const unsigned& index, const std::vector<unsigned>& nnbin, std::vector<unsigned>& indices ) const ;
/// Get the indices fof a point
@@ -139,6 +139,7 @@ public:
double getCubeUnits() const ;
/// This gives the normalisation of histograms
virtual double getNorm() const ;
+ virtual void switchOffNormalisation(){}
};
inline
diff --git a/src/gridtools/HistogramOnGrid.cpp b/src/gridtools/HistogramOnGrid.cpp
index 050810e1f4a2eabe0effe9cc4a803aea1bade3a2..6ca321c07133fa1be643a8d953e61322d99afe3a 100644
--- a/src/gridtools/HistogramOnGrid.cpp
+++ b/src/gridtools/HistogramOnGrid.cpp
@@ -29,18 +29,25 @@ void HistogramOnGrid::registerKeywords( Keywords& keys ){
GridVessel::registerKeywords( keys );
keys.add("compulsory","KERNEL","the type of kernel to use");
keys.add("compulsory","BANDWIDTH","the bandwidths");
+ keys.addFlag("STORE_NORMED",false,"are we to store the data normalised");
}
HistogramOnGrid::HistogramOnGrid( const vesselbase::VesselOptions& da ):
GridVessel(da),
norm(0),
+store_normed(false),
bandwidths(dimension)
{
- parse("KERNEL",kerneltype); parseVector("BANDWIDTH",bandwidths);
+ parseFlag("STORE_NORMED",store_normed); parse("KERNEL",kerneltype); parseVector("BANDWIDTH",bandwidths);
}
-void HistogramOnGrid::setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax ){
- GridVessel::setBounds( smin, smax );
+void HistogramOnGrid::switchOffNormalisation(){
+ store_normed=false;
+}
+
+void HistogramOnGrid::setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax,
+ const std::vector<unsigned>& nbins, const std::vector<double>& spacing ){
+ GridVessel::setBounds( smin, smax, nbins, spacing );
std::vector<double> point(dimension,0);
KernelFunctions kernel( point, bandwidths, kerneltype, false, 1.0, true );
nneigh=kernel.getSupport( dx );
@@ -80,8 +87,13 @@ void HistogramOnGrid::finish( const std::vector<double>& buffer ){
}
void HistogramOnGrid::clear(){
- if( !nomemory ) return ;
- norm = 0.; GridVessel::clear();
+ if( !nomemory && !store_normed ) return ;
+ if( nomemory ){
+ norm = 0.; GridVessel::clear(); return;
+ }
+ if( store_normed ){
+ for(unsigned i=0;i<data.size();++i) data[i] /= norm;
+ }
}
}
diff --git a/src/gridtools/HistogramOnGrid.h b/src/gridtools/HistogramOnGrid.h
index 0bdae980e082f0154202ec0ce02bbaf249cdcf95..cea0b9c02ccb96b3d51307c788f8c50d8856fe63 100644
--- a/src/gridtools/HistogramOnGrid.h
+++ b/src/gridtools/HistogramOnGrid.h
@@ -30,20 +30,22 @@ namespace gridtools {
class HistogramOnGrid : public GridVessel {
private:
double norm;
+ bool store_normed;
std::string kerneltype;
std::vector<double> bandwidths;
std::vector<unsigned> nneigh;
public:
static void registerKeywords( Keywords& keys );
explicit HistogramOnGrid( const vesselbase::VesselOptions& da );
- void setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax );
- std::string description(){ return ""; }
+ void setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax,
+ const std::vector<unsigned>& nbins, const std::vector<double>& spacing );
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 switchOffNormalisation();
void clear();
};
diff --git a/src/multicolvar/MultiColvarDensity.cpp b/src/multicolvar/MultiColvarDensity.cpp
index 230e32746ec76656d6b7255d66657c3c695b6ab0..697bc822b0db55a31c0f2ba1ebd175a6d050075f 100644
--- a/src/multicolvar/MultiColvarDensity.cpp
+++ b/src/multicolvar/MultiColvarDensity.cpp
@@ -63,6 +63,8 @@ class MultiColvarDensity :
bool fractional;
unsigned rstride;
MultiColvarBase* mycolv;
+ std::vector<unsigned> nbins;
+ std::vector<double> gspacing;
vesselbase::StoreDataVessel* stash;
gridtools::HistogramOnGrid* mygrid;
Vector origin;
@@ -91,7 +93,8 @@ void MultiColvarDensity::registerKeywords( Keywords& keys ){
keys.add("compulsory","STRIDE","1","the frequency with which the data should be collected and added to the grid");
keys.add("atoms","ORIGIN","we will use the position of this atom as the origin");
keys.add("compulsory","DIR","the direction in which to calculate the density profile");
- keys.add("compulsory","NBINS","the number of bins to use to represent the density profile");
+ keys.add("optional","NBINS","the number of bins to use to represent the density profile");
+ keys.add("optional","SPACING","the approximate grid spacing (to be used as an alternative or together with NBINS)");
keys.add("compulsory","BANDWIDTH","the bandwidths for kernel density esimtation");
keys.add("compulsory","KERNEL","gaussian","the kernel function you are using. More details on the kernels available "
"in plumed plumed can be found in \\ref kernelfunctions.");
@@ -148,8 +151,10 @@ MultiColvarDensity::MultiColvarDensity(const ActionOptions&ao):
error( direction + " is not valid gradient direction");
}
log.printf(" for colvars calculated by action %s \n",mycolv->getLabel().c_str() );
+ parseVector("NBINS",nbins); parseVector("SPACING",gspacing);
+ if( nbins.size()!=directions.size() && gspacing.size()!=directions.size() ) error("NBINS or SPACING must be set");
- std::string vstring = getKeyword("NBINS") + " " + getKeyword("KERNEL") + " " + getKeyword("BANDWIDTH");
+ std::string vstring = getKeyword("KERNEL") + " " + getKeyword("BANDWIDTH");
vstring += " PBC=T"; for(unsigned i=1;i<directions.size();++i) vstring+=",T";
vstring +=" COMPONENTS=" + getPntrToArgument()->getLabel() + ".dens";
vstring +=" COORDINATES=";
@@ -169,7 +174,6 @@ MultiColvarDensity::MultiColvarDensity(const ActionOptions&ao):
Keywords keys; gridtools::HistogramOnGrid::registerKeywords( keys );
vesselbase::VesselOptions dar( da, keys );
mygrid = new gridtools::HistogramOnGrid(dar); addVessel( mygrid );
- resizeFunctions();
// Enusre units for cube files are set correctly
if( !fractional ){
@@ -204,7 +208,7 @@ void MultiColvarDensity::update(){
if( plumed.getRestart() ){
error("restarting of MultiColvarDensity is not yet implemented");
} else {
- mygrid->setBounds( gmin, gmax );
+ mygrid->setBounds( gmin, gmax, nbins, gspacing ); resizeFunctions();
}
firststep=false; // We only have the first step once
} else {