diff --git a/src/bias/MetaD.cpp b/src/bias/MetaD.cpp
index ed0c0170064f3632ba1d47e1b9e448fe5fe1fc8d..0fa63fcc5ceef83f639246848678a4090727437c 100644
--- a/src/bias/MetaD.cpp
+++ b/src/bias/MetaD.cpp
@@ -264,6 +264,8 @@ private:
double lowI_;
bool doInt_;
bool isFirstStep;
+/// accumulator for work
+ double work_;
void readGaussians(IFile*);
bool readChunkOfGaussians(IFile *ifile, unsigned n);
@@ -353,6 +355,7 @@ walkers_mpi(false),
acceleration(false), acc(0.0),
// Interval initialization
uppI_(-1), lowI_(-1), doInt_(false),
+work_(0.0),
isFirstStep(true)
{
// parse the flexible hills
@@ -573,6 +576,7 @@ isFirstStep(true)
}
addComponent("bias"); componentIsNotPeriodic("bias");
+ addComponent("work"); componentIsNotPeriodic("work");
if(acceleration) {
if(!welltemp_) error("The calculation of the acceleration works only if Well-Tempered Metadynamics is on");
@@ -1001,6 +1005,7 @@ void MetaD::calculate()
double mean_acc = acc/((double) getStep());
getPntrToComponent("acc")->set(mean_acc);
}
+ getPntrToComponent("work")->set(work_);
// set Forces
for(unsigned i=0;i<ncv;++i){
const double f=-der[i];
@@ -1021,6 +1026,8 @@ void MetaD::update(){
for(unsigned i=0;i<cv.size();++i){cv[i]=getArgument(i);}
+ double vbias=getBiasAndDerivatives(cv);
+
// if you use adaptive, call the FlexibleBin
if (adaptive_!=FlexibleBin::none){
flexbin->update(nowAddAHill);
@@ -1086,6 +1093,10 @@ void MetaD::update(){
writeGaussian(newhill,hillsOfile_);
}
}
+
+ double vbias1=getBiasAndDerivatives(cv);
+ work_+=vbias1-vbias;
+
// dump grid on file
if(wgridstride_>0&&getStep()%wgridstride_==0){
// in case old grids are stored, a sequence of grids should appear
diff --git a/src/bias/MovingRestraint.cpp b/src/bias/MovingRestraint.cpp
index 1f3b8f35571813028d3dc8cd8a0ab63f4c4bb042..2c7de3dc48f6b210cf733b16775a6094ed3c2c63 100644
--- a/src/bias/MovingRestraint.cpp
+++ b/src/bias/MovingRestraint.cpp
@@ -110,6 +110,7 @@ class MovingRestraint : public Bias{
std::vector<double> oldf;
std::vector<string> verse;
std::vector<double> work;
+ double tot_work;
public:
MovingRestraint(const ActionOptions&);
void calculate();
@@ -200,6 +201,8 @@ verse(getNumberOfArguments())
addComponent(comp); componentIsNotPeriodic(comp);
work.push_back(0.); // initialize the work value
}
+ addComponent("work"); componentIsNotPeriodic("work");
+ tot_work=0.0;
log<<" Bibliography ";
log<<cite("Grubmuller, Heymann, and Tavan, Science 271, 997 (1996)")<<"\n";
@@ -231,6 +234,7 @@ void MovingRestraint::calculate(){
for(unsigned j=0;j<narg;j++) kk[j]=(c1*kappa[i-1][j]+c2*kappa[i][j]);
for(unsigned j=0;j<narg;j++) aa[j]=(c1*at[i-1][j]+c2*at[i][j]);
}
+ tot_work=0.0;
for(unsigned i=0;i<narg;++i){
const double cv=difference(i,aa[i],getArgument(i)); // this gives: getArgument(i) - aa[i]
getPntrToComponent(getPntrToArgument(i)->getName()+"_cntr")->set(aa[i]);
@@ -243,10 +247,12 @@ void MovingRestraint::calculate(){
if(oldaa.size()==aa.size() && oldf.size()==f.size()) work[i]+=0.5*(oldf[i]+f[i])*(aa[i]-oldaa[i]) + 0.5*( dpotdk[i]+olddpotdk[i] )*(kk[i]-oldk[i]);
getPntrToComponent(getPntrToArgument(i)->getName()+"_work")->set(work[i]);
getPntrToComponent(getPntrToArgument(i)->getName()+"_kappa")->set(kk[i]);
+ tot_work+=work[i];
ene+=0.5*k*cv*cv;
setOutputForce(i,f[i]);
totf2+=f[i]*f[i];
};
+ getPntrToComponent("work")->set(tot_work);
oldf=f;
oldaa=aa;
oldk=kk;
diff --git a/src/core/PlumedMain.cpp b/src/core/PlumedMain.cpp
index ed827717fd6b69e69436c79556a6363f5c8a0ff3..21c72035c573bd5b7aad0d79c523f6c28a452d13 100644
--- a/src/core/PlumedMain.cpp
+++ b/src/core/PlumedMain.cpp
@@ -65,6 +65,7 @@ PlumedMain::PlumedMain():
atoms(*new Atoms(*this)),
actionSet(*new ActionSet(*this)),
bias(0.0),
+ work(0.0),
exchangePatterns(*new(ExchangePatterns)),
exchangeStep(false),
restart(false),
@@ -611,6 +612,7 @@ void PlumedMain::justCalculate(){
if(!active)return;
stopwatch.start("4 Calculating (forward loop)");
bias=0.0;
+ work=0.0;
int iaction=0;
// calculate the active actions in order (assuming *backward* dependence)
@@ -636,6 +638,7 @@ void PlumedMain::justCalculate(){
else (*p)->calculate();
// This retrieves components called bias
if(av) bias+=av->getOutputQuantity("bias");
+ if(av) work+=av->getOutputQuantity("work");
if(av)av->setGradientsIfNeeded();
ActionWithVirtualAtom*avv=dynamic_cast<ActionWithVirtualAtom*>(*p);
if(avv)avv->setGradientsIfNeeded();
@@ -736,6 +739,10 @@ double PlumedMain::getBias() const{
return bias;
}
+double PlumedMain::getWork() const{
+ return work;
+}
+
FILE* PlumedMain::fopen(const char *path, const char *mode){
std::string mmode(mode);
std::string ppath(path);
diff --git a/src/core/PlumedMain.h b/src/core/PlumedMain.h
index e8cb4fc18f5ad7645dd2b64979a3207312c657eb..659b05ecb0a35e0f0a48ea6d272244389ec83477 100644
--- a/src/core/PlumedMain.h
+++ b/src/core/PlumedMain.h
@@ -120,6 +120,10 @@ private:
/// The total bias (=total energy of the restraints)
double bias;
+/// The total work.
+/// This computed by accumulating the change in external potentials.
+ double work;
+
/// Class of possible exchange patterns, used for BIASEXCHANGE but also for future parallel tempering
ExchangePatterns& exchangePatterns;
@@ -247,6 +251,8 @@ public:
void setSuffix(const std::string&);
/// get the value of the bias
double getBias()const;
+/// get the value of the work
+ double getWork()const;
/// Opens a file.
/// Similar to plain fopen, but, if it finds an error in opening the file, it also tries with
/// path+suffix. This trick is useful for multiple replica simulations.
diff --git a/src/generic/EffectiveEnergyDrift.cpp b/src/generic/EffectiveEnergyDrift.cpp
index 51136a75bbadeedcb373ae4de2d1274b5217a445..864751e2d876c789f85ea66355302d8105112f27 100644
--- a/src/generic/EffectiveEnergyDrift.cpp
+++ b/src/generic/EffectiveEnergyDrift.cpp
@@ -211,6 +211,8 @@ void EffectiveEnergyDrift::update(){
//is not a multiple of the bias actions stride
for(int i=0;i<biases.size();i++) bias+=biases[i]->getOutputQuantity("bias");
+ bias-=plumed.getWork();
+
plumed.comm.Sum(&eedSum,1);
output.printf("%ld %f\n",plumed.getStep(),eedSum+bias);
}