diff --git a/regtest/analysis/rt-wham/Makefile b/regtest/analysis/rt-wham/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..3703b27cea227aa053fb6d1d73f861e4384dbcee
--- /dev/null
+++ b/regtest/analysis/rt-wham/Makefile
@@ -0,0 +1 @@
+include ../../scripts/test.make
diff --git a/regtest/analysis/rt-wham/alltraj.xtc b/regtest/analysis/rt-wham/alltraj.xtc
new file mode 100644
index 0000000000000000000000000000000000000000..9b30b5aa4c53cb53febd8168250f9794b161ce8b
Binary files /dev/null and b/regtest/analysis/rt-wham/alltraj.xtc differ
diff --git a/regtest/analysis/rt-wham/config b/regtest/analysis/rt-wham/config
new file mode 100644
index 0000000000000000000000000000000000000000..f5872ebf331cfde46ac28156a819f7d4f590b4b3
--- /dev/null
+++ b/regtest/analysis/rt-wham/config
@@ -0,0 +1,2 @@
+type=driver
+arg="--mf_xtc alltraj.xtc"
diff --git a/regtest/analysis/rt-wham/fes.dat.reference b/regtest/analysis/rt-wham/fes.dat.reference
new file mode 100644
index 0000000000000000000000000000000000000000..aacce303159a2f589c97cb7f9f8cee375f7c19da
--- /dev/null
+++ b/regtest/analysis/rt-wham/fes.dat.reference
@@ -0,0 +1,57 @@
+#! FIELDS ff.phi fes
+#! SET normalisation 1.0000
+#! SET min_ff.phi -pi
+#! SET max_ff.phi pi
+#! SET nbins_ff.phi 50
+#! SET periodic_ff.phi false
+ -3.1416 53.0597
+ -3.0159 47.2084
+ -2.8903 44.8862
+ -2.7646 39.8037
+ -2.6389 37.2465
+ -2.5133 34.1986
+ -2.3876 32.0535
+ -2.2619 29.4245
+ -2.1363 28.9321
+ -2.0106 25.3416
+ -1.8850 20.0828
+ -1.7593 15.1329
+ -1.6336 12.5434
+ -1.5080 9.1217
+ -1.3823 4.8916
+ -1.2566 1.0762
+ -1.1310 4.5409
+ -1.0053 11.2837
+ -0.8796 17.3947
+ -0.7540 23.3114
+ -0.6283 27.6347
+ -0.5027 35.2913
+ -0.3770 38.0352
+ -0.2513 42.3201
+ -0.1257 45.6324
+ 0.0000 44.9696
+ 0.1257 44.1181
+ 0.2513 43.5800
+ 0.3770 42.8887
+ 0.5027 37.0370
+ 0.6283 32.2551
+ 0.7540 26.1036
+ 0.8796 22.2865
+ 1.0053 20.8150
+ 1.1310 19.5912
+ 1.2566 25.6362
+ 1.3823 31.3547
+ 1.5080 40.4811
+ 1.6336 52.0285
+ 1.7593 59.3373
+ 1.8850 70.1980
+ 2.0106 77.5370
+ 2.1363 80.1242
+ 2.2619 85.1814
+ 2.3876 82.7074
+ 2.5133 83.5724
+ 2.6389 79.7823
+ 2.7646 72.3409
+ 2.8903 64.2447
+ 3.0159 59.1219
+ 3.1416 56.8935
diff --git a/regtest/analysis/rt-wham/plumed.dat b/regtest/analysis/rt-wham/plumed.dat
new file mode 100644
index 0000000000000000000000000000000000000000..d53eab06273b884a32d3d7f8ec4bfe4e9b99254a
--- /dev/null
+++ b/regtest/analysis/rt-wham/plumed.dat
@@ -0,0 +1,81 @@
+phi: TORSION ATOMS=5,7,9,15
+psi: TORSION ATOMS=7,9,15,17
+rp0: RESTRAINT ARG=phi KAPPA=200.0 AT=-3.00000000000000000000
+rp1: RESTRAINT ARG=phi KAPPA=200.0 AT=-2.80645161290322580646
+rp2: RESTRAINT ARG=phi KAPPA=200.0 AT=-2.61290322580645161292
+rp3: RESTRAINT ARG=phi KAPPA=200.0 AT=-2.41935483870967741938
+rp4: RESTRAINT ARG=phi KAPPA=200.0 AT=-2.22580645161290322584
+rp5: RESTRAINT ARG=phi KAPPA=200.0 AT=-2.03225806451612903230
+rp6: RESTRAINT ARG=phi KAPPA=200.0 AT=-1.83870967741935483876
+rp7: RESTRAINT ARG=phi KAPPA=200.0 AT=-1.64516129032258064522
+rp8: RESTRAINT ARG=phi KAPPA=200.0 AT=-1.45161290322580645168
+rp9: RESTRAINT ARG=phi KAPPA=200.0 AT=-1.25806451612903225814
+rp10: RESTRAINT ARG=phi KAPPA=200.0 AT=-1.06451612903225806460
+rp11: RESTRAINT ARG=phi KAPPA=200.0 AT=-.87096774193548387106
+rp12: RESTRAINT ARG=phi KAPPA=200.0 AT=-.67741935483870967752
+rp13: RESTRAINT ARG=phi KAPPA=200.0 AT=-.48387096774193548398
+rp14: RESTRAINT ARG=phi KAPPA=200.0 AT=-.29032258064516129044
+rp15: RESTRAINT ARG=phi KAPPA=200.0 AT=-.09677419354838709690
+rp16: RESTRAINT ARG=phi KAPPA=200.0 AT=.09677419354838709664
+rp17: RESTRAINT ARG=phi KAPPA=200.0 AT=.29032258064516129018
+rp18: RESTRAINT ARG=phi KAPPA=200.0 AT=.48387096774193548372
+rp19: RESTRAINT ARG=phi KAPPA=200.0 AT=.67741935483870967726
+rp20: RESTRAINT ARG=phi KAPPA=200.0 AT=.87096774193548387080
+rp21: RESTRAINT ARG=phi KAPPA=200.0 AT=1.06451612903225806434
+rp22: RESTRAINT ARG=phi KAPPA=200.0 AT=1.25806451612903225788
+rp23: RESTRAINT ARG=phi KAPPA=200.0 AT=1.45161290322580645142
+rp24: RESTRAINT ARG=phi KAPPA=200.0 AT=1.64516129032258064496
+rp25: RESTRAINT ARG=phi KAPPA=200.0 AT=1.83870967741935483850
+rp26: RESTRAINT ARG=phi KAPPA=200.0 AT=2.03225806451612903204
+rp27: RESTRAINT ARG=phi KAPPA=200.0 AT=2.22580645161290322558
+rp28: RESTRAINT ARG=phi KAPPA=200.0 AT=2.41935483870967741912
+rp29: RESTRAINT ARG=phi KAPPA=200.0 AT=2.61290322580645161266
+rp30: RESTRAINT ARG=phi KAPPA=200.0 AT=2.80645161290322580620
+rp31: RESTRAINT ARG=phi KAPPA=200.0 AT=2.99999999999999999974
+
+REWEIGHT_WHAM ...
+ ARG1=rp0.bias
+ ARG2=rp1.bias
+ ARG3=rp2.bias
+ ARG4=rp3.bias
+ ARG5=rp4.bias
+ ARG6=rp5.bias
+ ARG7=rp6.bias
+ ARG8=rp7.bias
+ ARG9=rp8.bias
+ ARG10=rp9.bias
+ ARG11=rp10.bias
+ ARG12=rp11.bias
+ ARG13=rp12.bias
+ ARG14=rp13.bias
+ ARG15=rp14.bias
+ ARG16=rp15.bias
+ ARG17=rp16.bias
+ ARG18=rp17.bias
+ ARG19=rp18.bias
+ ARG20=rp19.bias
+ ARG21=rp20.bias
+ ARG22=rp21.bias
+ ARG23=rp22.bias
+ ARG24=rp23.bias
+ ARG25=rp24.bias
+ ARG26=rp25.bias
+ ARG27=rp26.bias
+ ARG28=rp27.bias
+ ARG29=rp28.bias
+ ARG30=rp29.bias
+ ARG31=rp30.bias
+ ARG32=rp31.bias
+ TEMP=300
+ LABEL=ww
+... REWEIGHT_WHAM
+
+PRINT ARG=phi,psi FILE=colvar
+PRINT ARG=rp0.bias,rp1.bias,rp2.bias,rp3.bias,rp4.bias,rp5.bias,rp6.bias,rp7.bias,rp8.bias,rp9.bias,rp10.bias,rp11.bias,rp12.bias,rp13.bias,rp14.bias,rp15.bias,rp16.bias,rp17.bias,rp18.bias,rp19.bias,rp20.bias,rp21.bias,rp22.bias,rp23.bias,rp24.bias,rp25.bias,rp26.bias,rp27.bias,rp28.bias,rp29.bias,rp30.bias,rp31.bias FILE=bias
+
+ff: COLLECT_FRAMES ARG=phi STRIDE=1 LOGWEIGHTS=ww
+hh: HISTOGRAM ARG=ff.phi GRID_MIN=-pi GRID_MAX=pi GRID_BIN=50 KERNEL=DISCRETE
+
+fes: CONVERT_TO_FES GRID=hh TEMP=300
+DUMPGRID GRID=fes FILE=fes.dat FMT=%8.4f
+
diff --git a/src/analysis/DataCollectionObject.cpp b/src/analysis/DataCollectionObject.cpp
index 2b9152f2738d415258842f7fdcb847b53aae0a36..9950d7097448fd205504a1c9fb226d0490bbb6a3 100644
--- a/src/analysis/DataCollectionObject.cpp
+++ b/src/analysis/DataCollectionObject.cpp
@@ -25,8 +25,8 @@
namespace PLMD {
namespace analysis {
-void DataCollectionObject::setAtomNumbersAndArgumentNames( const std::vector<AtomNumber>& ind, const std::vector<std::string>& arg_names ){
- indices.resize( ind.size() ); positions.resize( indices.size() );
+void DataCollectionObject::setAtomNumbersAndArgumentNames( const std::string& action_label, const std::vector<AtomNumber>& ind, const std::vector<std::string>& arg_names ){
+ myaction=action_label; indices.resize( ind.size() ); positions.resize( indices.size() );
for(unsigned i=0;i<ind.size();++i) indices[i]=ind[i];
for(unsigned i=0;i<arg_names.size();++i) args.insert( std::pair<std::string,double>( arg_names[i], 0.0 ) );
}
diff --git a/src/analysis/DataCollectionObject.h b/src/analysis/DataCollectionObject.h
index df44aa3774fc1ff044310a4b263aa2f0911ab78f..9c5c723cffff88093da64b127b2e975e3fdcfd30 100644
--- a/src/analysis/DataCollectionObject.h
+++ b/src/analysis/DataCollectionObject.h
@@ -36,6 +36,8 @@ namespace analysis {
class DataCollectionObject {
friend class ReadAnalysisFrames;
private:
+/// The label of the action in which the data is stored
+ std::string myaction;
/// The list of atom numbers that are stored in the object
std::vector<AtomNumber> indices;
/// The list of atomic positions
@@ -44,7 +46,7 @@ private:
std::map<std::string,double> args;
public:
/// Set the names and atom numbers
- void setAtomNumbersAndArgumentNames( const std::vector<AtomNumber>& ind, const std::vector<std::string>& arg_names );
+ void setAtomNumbersAndArgumentNames( const std::string& action_label, const std::vector<AtomNumber>& ind, const std::vector<std::string>& arg_names );
/// Set the positions of all the atoms
void setAtomPositions( const std::vector<Vector>& pos );
/// Set the value of one of the arguments
@@ -64,6 +66,8 @@ Vector DataCollectionObject::getAtomPosition( const AtomNumber& ind ) const {
inline
double DataCollectionObject::getArgumentValue( const std::string& name ) const {
+ std::size_t dot=name.find_first_of('.'); std::string a=name.substr(0,dot);
+ if( a==myaction ) return args.find( name.substr(dot+1) )->second;
return args.find(name)->second;
}
diff --git a/src/analysis/Histogram.cpp b/src/analysis/Histogram.cpp
index 53f16ebb1b5f1d2ff130ff752b003437bb58d3ba..6b289cbc04964e40ab1a311e775476014ffb6cf0 100644
--- a/src/analysis/Histogram.cpp
+++ b/src/analysis/Histogram.cpp
@@ -28,6 +28,7 @@
#include "core/PlumedMain.h"
#include "core/ActionSet.h"
#include "AnalysisBase.h"
+#include "ReadAnalysisFrames.h"
namespace PLMD{
namespace analysis{
@@ -162,6 +163,7 @@ DUMPGRID GRID=hh FILE=histo STRIDE=100000
class Histogram : public gridtools::ActionWithGrid {
private:
double ww;
+ bool activated;
bool in_apply, mvectors;
KernelFunctions* kernel;
std::vector<double> forcesToApply, finalForces;
@@ -202,6 +204,7 @@ Histogram::Histogram(const ActionOptions&ao):
Action(ao),
ActionWithGrid(ao),
ww(0.0),
+activated(false),
in_apply(false),
mvectors(false),
kernel(NULL),
@@ -357,7 +360,10 @@ unsigned Histogram::getNumberOfQuantities() const {
bool Histogram::onStep() const {
if( !my_analysis_object ) return ActionPilot::onStep();
- return true;
+ ReadAnalysisFrames* myfram = dynamic_cast<ReadAnalysisFrames*>( my_analysis_object );
+ if( myfram && activated ) return true;
+ else if( myfram ) return ActionPilot::onStep();
+ return my_analysis_object->onStep();
}
void Histogram::prepareForAveraging(){
@@ -522,7 +528,7 @@ void Histogram::apply(){
}
void Histogram::runFinalJobs(){
- if( my_analysis_object && getStride()==0 ) update();
+ if( my_analysis_object && getStride()==0 ){ activated=true; update(); }
}
}
diff --git a/src/analysis/Makefile b/src/analysis/Makefile
index 02080ce4bf8899b32236ac292b3cb42660b01715..937de00850de8190775dec56ff505fd51612ac1a 100644
--- a/src/analysis/Makefile
+++ b/src/analysis/Makefile
@@ -1,4 +1,4 @@
-USE=core tools reference vesselbase gridtools multicolvar
+USE=core tools reference vesselbase gridtools multicolvar bias
#generic makefile
include ../maketools/make.module
diff --git a/src/analysis/ReadAnalysisFrames.cpp b/src/analysis/ReadAnalysisFrames.cpp
index 43ea1d855b77ab79fc40483b08a82c6ffe481df0..4784123465f2d6cd9138821fc7edb9470179e4db 100644
--- a/src/analysis/ReadAnalysisFrames.cpp
+++ b/src/analysis/ReadAnalysisFrames.cpp
@@ -79,7 +79,12 @@ weights_calculated(false)
}
if( wwstr.size()>0 ) log.printf("\n");
else log.printf(" weights are all equal to one\n");
+ wham_pointer = dynamic_cast<bias::ReweightWham*>( weight_vals[0]->getPntrToAction() );
+ if( wham_pointer && weight_vals.size()!=1 ) error("can only extract weights from one wham object");
requestArguments( arg );
+
+ // Now add fake components to the underlying ActionWithValue for the arguments
+ for(unsigned i=0;i<argument_names.size();++i){ addComponent( argument_names[i] ); componentIsNotPeriodic( argument_names[i] ); }
}
std::vector<Value*> ReadAnalysisFrames::getArgumentList(){
@@ -102,13 +107,18 @@ void ReadAnalysisFrames::calculateWeights(){
if( weight_vals.empty() ){
for(unsigned i=0;i<logweights.size();++i) weights[i]=1.0;
} else {
- // Find the maximum weight
- double maxweight=logweights[0];
- for(unsigned i=1;i<getNumberOfDataPoints();++i){
- if(logweights[i]>maxweight) maxweight=logweights[i];
+ if( wham_pointer ){
+ wham_pointer->calculateWeights( logweights.size() );
+ for(unsigned i=0;i<logweights.size();++i) weights[i]=wham_pointer->getWeight(i);
+ } else {
+ // Find the maximum weight
+ double maxweight=logweights[0];
+ for(unsigned i=1;i<getNumberOfDataPoints();++i){
+ if(logweights[i]>maxweight) maxweight=logweights[i];
+ }
+ // Calculate weights (no memory) -- business here with maxweight is to prevent overflows
+ for(unsigned i=0;i<logweights.size();++i) weights[i]=exp( logweights[i]-maxweight );
}
- // Calculate weights (no memory) -- business here with maxweight is to prevent overflows
- for(unsigned i=0;i<logweights.size();++i) weights[i]=exp( logweights[i]-maxweight );
}
}
@@ -118,8 +128,9 @@ void ReadAnalysisFrames::update(){
if( clearonnextstep ){
my_data_stash.clear(); my_data_stash.resize(0);
logweights.clear(); logweights.resize(0);
+ if( wham_pointer ) wham_pointer->clearData();
clearonnextstep=false;
- }
+ }
// Get the weight and store it in the weights array
double ww=0; for(unsigned i=0;i<weight_vals.size();++i) ww+=weight_vals[i]->get();
@@ -127,7 +138,7 @@ void ReadAnalysisFrames::update(){
// Now create the data collection object and push it back to be stored
unsigned index = my_data_stash.size(); my_data_stash.push_back( DataCollectionObject() );
- my_data_stash[index].setAtomNumbersAndArgumentNames( atom_numbers, argument_names );
+ my_data_stash[index].setAtomNumbersAndArgumentNames( getLabel(), atom_numbers, argument_names );
my_data_stash[index].setAtomPositions( getPositions() );
for(unsigned i=0;i<argument_names.size();++i) my_data_stash[index].setArgument( argument_names[i], getArgument(i) );
diff --git a/src/analysis/ReadAnalysisFrames.h b/src/analysis/ReadAnalysisFrames.h
index 3da79c53da99ea8d947b5e529b70a3f8f3a29545..53d766aaa683077720aea479656de89ee35e6452 100644
--- a/src/analysis/ReadAnalysisFrames.h
+++ b/src/analysis/ReadAnalysisFrames.h
@@ -23,6 +23,7 @@
#define __PLUMED_analysis_ReadAnalysisFrames_h
#include "AnalysisBase.h"
+#include "bias/ReweightWham.h"
namespace PLMD {
namespace analysis {
@@ -38,6 +39,8 @@ private:
std::vector<AtomNumber> atom_numbers;
/// The biases we are using in reweighting and the args we store them separately
std::vector<Value*> weight_vals;
+/// The object that calculates weights using WHAM
+ bias::ReweightWham* wham_pointer;
/// The weights of all the data points
bool weights_calculated;
std::vector<double> logweights, weights;
diff --git a/src/bias/ReweightBase.h b/src/bias/ReweightBase.h
index e3dfd35adb4ad69a437749287825c501d784f728..bc570c4a50c0373b8d7e01e21fb4d994fcbe86ff 100644
--- a/src/bias/ReweightBase.h
+++ b/src/bias/ReweightBase.h
@@ -40,7 +40,7 @@ public:
explicit ReweightBase(const ActionOptions&ao);
unsigned getNumberOfDerivatives(){ return 0; }
void calculate();
- virtual double getLogWeight() const = 0;
+ virtual double getLogWeight() = 0;
void apply(){}
};
diff --git a/src/bias/ReweightBias.cpp b/src/bias/ReweightBias.cpp
index d5459df962ccd55f8b302a1990ce0fbd606f2661..f3fa47fcdfcc93aa42b2fc7ad4f1aab30e820ff3 100644
--- a/src/bias/ReweightBias.cpp
+++ b/src/bias/ReweightBias.cpp
@@ -72,7 +72,7 @@ class ReweightBias : public ReweightBase {
public:
static void registerKeywords(Keywords&);
explicit ReweightBias(const ActionOptions&ao);
- double getLogWeight() const ;
+ double getLogWeight();
};
PLUMED_REGISTER_ACTION(ReweightBias,"REWEIGHT_BIAS")
@@ -88,7 +88,7 @@ ReweightBase(ao)
{
}
-double ReweightBias::getLogWeight() const {
+double ReweightBias::getLogWeight(){
// Retrieve the bias
double bias=0.0; for(unsigned i=0;i<getNumberOfArguments();++i) bias+=getArgument(i);
return bias / simtemp;
diff --git a/src/bias/ReweightMetad.cpp b/src/bias/ReweightMetad.cpp
index e53170be71c3d06c0debda8ca6b04ca704f5872c..2085e65a02032512c1d475b66c95396a98d84c9b 100644
--- a/src/bias/ReweightMetad.cpp
+++ b/src/bias/ReweightMetad.cpp
@@ -70,7 +70,7 @@ class ReweightMetad : public ReweightBase {
public:
static void registerKeywords(Keywords&);
explicit ReweightMetad(const ActionOptions&ao);
- double getLogWeight() const ;
+ double getLogWeight();
};
PLUMED_REGISTER_ACTION(ReweightMetad,"REWEIGHT_METAD")
@@ -86,7 +86,7 @@ ReweightBase(ao)
{
}
-double ReweightMetad::getLogWeight() const {
+double ReweightMetad::getLogWeight(){
// Retrieve the bias
double bias=0.0; for(unsigned i=0;i<getNumberOfArguments();++i) bias+=getArgument(i);
return bias / simtemp;
diff --git a/src/bias/ReweightTemperature.cpp b/src/bias/ReweightTemperature.cpp
index 45f1c6f2b9b9c7748cde85fd8241a0537285028d..19fb36ddbe03b323380676cf13cea84933bc1fc3 100644
--- a/src/bias/ReweightTemperature.cpp
+++ b/src/bias/ReweightTemperature.cpp
@@ -59,7 +59,7 @@ public:
static void registerKeywords(Keywords&);
explicit ReweightTemperature(const ActionOptions&ao);
void prepare();
- double getLogWeight() const ;
+ double getLogWeight();
};
PLUMED_REGISTER_ACTION(ReweightTemperature,"REWEIGHT_TEMP")
@@ -95,7 +95,7 @@ void ReweightTemperature::prepare(){
plumed.getAtoms().setCollectEnergy(true);
}
-double ReweightTemperature::getLogWeight() const {
+double ReweightTemperature::getLogWeight(){
// Retrieve the bias
double bias=0.0; for(unsigned i=0;i<biases.size();++i) bias+=biases[i]->get();
double energy=plumed.getAtoms().getEnergy()+bias;
diff --git a/src/bias/ReweightWham.cpp b/src/bias/ReweightWham.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..892c656d7f331cbbffce9094ed27c302c8338726
--- /dev/null
+++ b/src/bias/ReweightWham.cpp
@@ -0,0 +1,122 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2011-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 "core/ActionRegister.h"
+#include "ReweightWham.h"
+
+//+PLUMEDOC REWEIGHTING REWEIGHT_WHAM
+/*
+
+\par Examples
+
+*/
+//+ENDPLUMEDOC
+
+namespace PLMD {
+namespace bias {
+
+PLUMED_REGISTER_ACTION(ReweightWham,"REWEIGHT_WHAM")
+
+void ReweightWham::registerKeywords(Keywords& keys ){
+ ReweightBase::registerKeywords( keys ); keys.remove("ARG");
+ keys.add("numbered","ARG","the biases that must be taken into account when reweighting");
+ keys.add("compulsory","MAXITER","1000","maximum number of iterations for WHAM algorithm");
+ keys.add("compulsory","WHAMTOL","1e-10","threshold for convergence of WHAM algorithm");
+ keys.reset_style("ARG","compulsory"); keys.add("hidden","DATA","sneaky trick to remove error for reading numbered args");
+}
+
+ReweightWham::ReweightWham(const ActionOptions&ao):
+Action(ao),
+ReweightBase(ao)
+{
+ std::vector<Value*> targ, fagr;
+ unsigned nbias = 0; wlists.push_back( 0 );
+ for(unsigned i=1;;i++){
+ if( !parseArgumentList("ARG",i,targ ) ) break;
+ log.printf(" bias number %d involves :");
+ for(unsigned j=0;j<targ.size();++j){
+ log.printf("%s ",targ[j]->getName().c_str() );
+ fagr.push_back( targ[j] );
+ }
+ log.printf("\n"); targ.resize(0);
+ wlists.push_back( fagr.size() );
+ nbias++;
+ }
+ plumed_assert( wlists.size()==(nbias+1) ); requestArguments( fagr );
+ parse("MAXITER",maxiter); parse("WHAMTOL",thresh);
+}
+
+double ReweightWham::getLogWeight(){
+ if( getStep()==0 ) return 1.0; // This is here as first step is ignored in all analyses
+ for(unsigned i=0;i<wlists.size()-1;++i){
+ double total_bias=0;
+ for(unsigned j=wlists[i];j<wlists[i+1];++j) total_bias+=getArgument(j);
+ stored_biases.push_back( total_bias );
+ }
+ return 1.0;
+}
+
+void ReweightWham::clearData(){
+ stored_biases.resize(0);
+}
+
+void ReweightWham::calculateWeights( const unsigned& nframes ){
+ if( stored_biases.size()!=(wlists.size()-1)*nframes ) error("wrong number of weights stored");
+ // Get the minimum value of the bias
+ double minv = *min_element(std::begin(stored_biases), std::end(stored_biases));
+ // Resize final weights array
+ plumed_assert( stored_biases.size()%(wlists.size()-1)==0 );
+ final_weights.resize( stored_biases.size() / (wlists.size()-1), 1.0 );
+ // Offset and exponential of the bias
+ std::vector<double> expv( stored_biases.size() );
+ for(unsigned i=0;i<expv.size();++i) expv[i] = exp( (-stored_biases[i]+minv) / simtemp );
+ // Initialize Z
+ std::vector<double> Z( wlists.size()-1, 1.0 ), oldZ( wlists.size()-1 );
+ // Now the iterative loop to calculate the WHAM weights
+ for(unsigned iter=0;iter<maxiter;++iter){
+ // Store Z
+ for(unsigned j=0;j<Z.size();++j) oldZ[j]=Z[j];
+ // Recompute weights
+ double norm=0;
+ for(unsigned j=0;j<final_weights.size();++j){
+ double ew=0;
+ for(unsigned k=0;k<Z.size();++k) ew += expv[j*Z.size()+k] / Z[k];
+ final_weights[j] = 1.0 / ew; norm += final_weights[j];
+ }
+ // Normalize weights
+ for(unsigned j=0;j<final_weights.size();++j) final_weights[j] /= norm;
+ // Recompute Z
+ for(unsigned j=0;j<Z.size();++j) Z[j] = 0.0;
+ for(unsigned j=0;j<final_weights.size();++j){
+ for(unsigned k=0;k<Z.size();++k) Z[k] += final_weights[j]*expv[j*Z.size()+k];
+ }
+ // Normalize Z and compute change in Z
+ double change=0; norm=0; for(unsigned k=0;k<Z.size();++k) norm+=Z[k];
+ for(unsigned k=0;k<Z.size();++k){
+ Z[k] /= norm; double d = std::log( Z[k] / oldZ[k] ); change += d*d;
+ }
+ if( change<thresh ) return;
+ }
+ error("Too many iterations in WHAM" );
+}
+
+}
+}
diff --git a/src/bias/ReweightWham.h b/src/bias/ReweightWham.h
new file mode 100644
index 0000000000000000000000000000000000000000..dbf16d79edfe1aaa5a42030ff564a0509b4dcf71
--- /dev/null
+++ b/src/bias/ReweightWham.h
@@ -0,0 +1,54 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2011-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_bias_ReweightWham_h
+#define __PLUMED_bias_ReweightWham_h
+
+#include "ReweightBase.h"
+
+namespace PLMD {
+namespace bias {
+
+class ReweightWham : public ReweightBase {
+private:
+ double thresh;
+ unsigned maxiter;
+ std::vector<unsigned> wlists;
+ std::vector<double> stored_biases;
+ std::vector<double> final_weights;
+public:
+ static void registerKeywords(Keywords&);
+ explicit ReweightWham(const ActionOptions&ao);
+ void calculateWeights( const unsigned& nframes );
+ void clearData();
+ double getLogWeight();
+ double getWeight( const unsigned& iweight ) const ;
+};
+
+inline
+double ReweightWham::getWeight( const unsigned& iweight ) const {
+ plumed_dbg_assert( calculatedWeights && iweight<final_weights.size() );
+ return final_weights[iweight];
+}
+
+}
+}
+#endif
diff --git a/src/core/ActionWithArguments.cpp b/src/core/ActionWithArguments.cpp
index 9a799a2cfd03f179b31021943a8080eeda989788..955013e77ce496014d7199e552a2868f7b4f2b36 100644
--- a/src/core/ActionWithArguments.cpp
+++ b/src/core/ActionWithArguments.cpp
@@ -47,9 +47,10 @@ void ActionWithArguments::registerKeywords(Keywords& keys){
}
void ActionWithArguments::parseArgumentList(const std::string&key,std::vector<Value*>&arg){
- vector<string> c; arg.clear(); parseVector(key,c);
+ std::string def; vector<string> c; arg.clear(); parseVector(key,c);
if( c.size()==0 && (keywords.style(key,"compulsory") || keywords.style(key,"hidden")) ){
- std::string def; if( keywords.getDefaultValue(key,def) ) c.push_back( def );
+ if( keywords.getDefaultValue(key,def) ) c.push_back( def );
+ else return;
}
interpretArgumentList(c,arg);
}
diff --git a/src/gridtools/HistogramOnGrid.cpp b/src/gridtools/HistogramOnGrid.cpp
index 8adb3b13e9763db778ecbc068a9c309fa520d2b3..8b4c85fdbc2bb720c7ba2b59293732d51d5bf4ae 100644
--- a/src/gridtools/HistogramOnGrid.cpp
+++ b/src/gridtools/HistogramOnGrid.cpp
@@ -109,8 +109,8 @@ void HistogramOnGrid::calculate( const unsigned& current, MultiValue& myvals, st
for(unsigned i=0;i<dimension;++i) point[i]=myvals.get( 1+i );
// Get the kernel
- unsigned num_neigh; std::vector<unsigned> neighbors;
- std::vector<double> der( dimension );
+ unsigned num_neigh; std::vector<unsigned> neighbors(1);
+ std::vector<double> der( 0 );
KernelFunctions* kernel=getKernelAndNeighbors( point, num_neigh, neighbors );
if( !kernel && getType()=="flat" ){