diff --git a/src/bias/MetaD.cpp b/src/bias/MetaD.cpp
index cc207b2aef1bf951cb953c383c55395f1eac48b7..27ee78e6fae808f971bd458906c33ae3b970f781 100644
--- a/src/bias/MetaD.cpp
+++ b/src/bias/MetaD.cpp
@@ -434,6 +434,7 @@ public:
void update();
static void registerKeywords(Keywords& keys);
bool checkNeedsGradients()const {if(adaptive_==FlexibleBin::geometry) {return true;} else {return false;}}
+ bool checkNeedsGradients()const;
};
PLUMED_REGISTER_ACTION(MetaD,"METAD")
@@ -1911,5 +1912,13 @@ double MetaD::getTransitionBarrierBias() {
}
}
+bool MetaD::checkNeedsGradients()const
+{
+ if(adaptive_==FlexibleBin::geometry) {
+ if(getStep()%stride_==0 && !isFirstStep) return true;
+ else return false;
+ } else return false;
+}
+
}
}
diff --git a/src/bias/PBMetaD.cpp b/src/bias/PBMetaD.cpp
index a77cb42be99de7ca69284d864e907a0987eb6085..fbfbdeb209f5210e5d98bbfd40ee382b8df3dd39 100644
--- a/src/bias/PBMetaD.cpp
+++ b/src/bias/PBMetaD.cpp
@@ -282,7 +282,7 @@ public:
void calculate();
void update();
static void registerKeywords(Keywords& keys);
- bool checkNeedsGradients()const {if(adaptive_==FlexibleBin::geometry) {return true;} else {return false;}}
+ bool checkNeedsGradients()const;
};
PLUMED_REGISTER_ACTION(PBMetaD,"PBMETAD")
@@ -1207,6 +1207,15 @@ bool PBMetaD::scanOneHill(unsigned iarg, IFile *ifile, vector<Value> &tmpvalues,
} else {
return false;
}
+
+}
+
+bool PBMetaD::checkNeedsGradients()const
+{
+ if(adaptive_==FlexibleBin::geometry) {
+ if(getStep()%stride_==0 && !isFirstStep) return true;
+ else return false;
+ } else return false;
}
}
diff --git a/src/core/FlexibleBin.cpp b/src/core/FlexibleBin.cpp
index c2f83313c310bf8598ed2697c5b7b700c1abe755..b7296658cccb0c83beef8561a42eaf83523a02b9 100644
--- a/src/core/FlexibleBin.cpp
+++ b/src/core/FlexibleBin.cpp
@@ -98,44 +98,37 @@ FlexibleBin::FlexibleBin(int type, ActionWithArguments *paction, unsigned iarg,
void FlexibleBin::update(bool nowAddAHill) {
unsigned ncv=paction->getNumberOfArguments();
unsigned dimension=ncv*(ncv+1)/2;
- // this is done all the times from scratch. It is not an accumulator
- unsigned k=0;
- unsigned i;
- vector<double> cv;
vector<double> delta;
- // if you use this below then the decay is in time units
- //double decay=paction->getTimeStep()/sigma;
- // to be consistent with the rest of the program: everything is better to be in timesteps
+ vector<double> cv;
double decay=1./sigma;
+ // this is done all the times from scratch. It is not an accumulator
// here update the flexible bin according to the needs
switch (type) {
// This should be called every time
case diffusion:
- //
+ // if you use this below then the decay is in time units
+ //double decay=paction->getTimeStep()/sigma;
+ // to be consistent with the rest of the program: everything is better to be in timesteps
// THE AVERAGE VALUE
- //
// beware: the pbc
delta.resize(ncv);
- for(i=0; i<ncv; i++)cv.push_back(paction->getArgument(i));
+ for(unsigned i=0; i<ncv; i++) cv.push_back(paction->getArgument(i));
if(average.size()==0) { // initial time: just set the initial vector
average.resize(ncv);
- for(i=0; i<ncv; i++)average[i]=cv[i];
+ for(unsigned i=0; i<ncv; i++) average[i]=cv[i];
} else { // accumulate
- for(i=0; i<ncv; i++) {
+ for(unsigned i=0; i<ncv; i++) {
delta[i]=paction->difference(i,average[i],cv[i]);
average[i]+=decay*delta[i];
average[i]=paction->bringBackInPbc(i,average[i]); // equation 8 of "Metadynamics with adaptive Gaussians"
}
-
}
- //
// THE VARIANCE
- //
if(variance.size()==0) {
variance.resize(dimension,0.); // nonredundant members dimension=ncv*(ncv+1)/2;
} else {
- k=0;
- for(i=0; i<ncv; i++) {
+ unsigned k=0;
+ for(unsigned i=0; i<ncv; i++) {
for(unsigned j=i; j<ncv; j++) { // upper diagonal loop
variance[k]+=decay*(delta[i]*delta[j]-variance[k]);
k++;
@@ -144,31 +137,25 @@ void FlexibleBin::update(bool nowAddAHill) {
}
break;
case geometry:
- //
//this calculates in variance the \nabla CV_i \dot \nabla CV_j
- //
variance.resize(dimension);
- //cerr<< "Doing geometry "<<endl;
// now the signal for retrieving the gradients should be already given by checkNeedsGradients.
// here just do the projections
// note that the call checkNeedsGradients() in BiasMetaD takes care of switching on the call to gradients
if (nowAddAHill) { // geometry is sync with hill deposition
- //cerr<< "add a hill "<<endl;
- k=0;
+ unsigned k=0;
for(unsigned i=0; i<ncv; i++) {
for(unsigned j=i; j<ncv; j++) {
// eq 12 of "Metadynamics with adaptive Gaussians"
variance[k]=sigma*sigma*(paction->getProjection(i,j));
k++;
}
- };
- };
+ }
+ }
break;
default:
- cerr<< "This flexible bin is not recognized "<<endl;
- exit(1) ;
+ plumed_merror("This flexible bin method is not recognized");
}
-
}
vector<double> FlexibleBin::getMatrix() const {
@@ -180,18 +167,16 @@ vector<double> FlexibleBin::getMatrix() const {
/// in case of gradient based: update only when you add the hill
void FlexibleBin::update(bool nowAddAHill, unsigned iarg) {
// this is done all the times from scratch. It is not an accumulator
+ // here update the flexible bin according to the needs
vector<double> cv;
vector<double> delta;
// if you use this below then the decay is in time units
// to be consistent with the rest of the program: everything is better to be in timesteps
double decay=1./sigma;
- // here update the flexible bin according to the needs
switch (type) {
// This should be called every time
case diffusion:
- //
// THE AVERAGE VALUE
- //
delta.resize(1);
cv.push_back(paction->getArgument(iarg));
if(average.size()==0) { // initial time: just set the initial vector
@@ -202,9 +187,7 @@ void FlexibleBin::update(bool nowAddAHill, unsigned iarg) {
average[0]+=decay*delta[0];
average[0]=paction->bringBackInPbc(iarg,average[0]); // equation 8 of "Metadynamics with adaptive Gaussians"
}
- //
// THE VARIANCE
- //
if(variance.size()==0) {
variance.resize(1,0.); // nonredundant members dimension=ncv*(ncv+1)/2;
} else {
@@ -212,9 +195,7 @@ void FlexibleBin::update(bool nowAddAHill, unsigned iarg) {
}
break;
case geometry:
- //
//this calculates in variance the \nabla CV_i \dot \nabla CV_j
- //
variance.resize(1);
// now the signal for retrieving the gradients should be already given by checkNeedsGradients.
// here just do the projections
@@ -225,8 +206,7 @@ void FlexibleBin::update(bool nowAddAHill, unsigned iarg) {
}
break;
default:
- cerr<< "This flexible bin is not recognized "<<endl;
- exit(1) ;
+ plumed_merror("This flexible bin is not recognized");
}
}