diff --git a/src/ActionAtomistic.cpp b/src/ActionAtomistic.cpp
index 205bf888eeca9ecbcd8243a6c88cbcd3a6f48e9b..9f643281738a9fbdc90d5f22a137abf54a789798 100644
--- a/src/ActionAtomistic.cpp
+++ b/src/ActionAtomistic.cpp
@@ -108,6 +108,7 @@ void ActionAtomistic::retrieveAtoms(){
for(unsigned j=0;j<indexes.size();j++) positions[j]=p[indexes[j]];
for(unsigned j=0;j<indexes.size();j++) charges[j]=c[indexes[j]];
for(unsigned j=0;j<indexes.size();j++) masses[j]=m[indexes[j]];
+ energy=plumed.getAtoms().getEnergy();
}
void ActionAtomistic::applyForces(){
@@ -115,6 +116,7 @@ void ActionAtomistic::applyForces(){
Tensor & v(plumed.getAtoms().virial);
for(unsigned j=0;j<indexes.size();j++) f[indexes[j]]+=forces[j];
v+=virial;
+ plumed.getAtoms().forceOnEnergy+=forceOnEnergy;
}
diff --git a/src/ActionAtomistic.h b/src/ActionAtomistic.h
index 6250f3acc6de8ca906f39c9aee51948cdbc1b59d..296c40a00dc10efdca94a0b5a37fc0e437eec750 100644
--- a/src/ActionAtomistic.h
+++ b/src/ActionAtomistic.h
@@ -20,13 +20,16 @@ class ActionAtomistic :
std::vector<int> indexes; // the set of needed atoms
std::set<int> unique;
std::vector<Vector> positions; // positions of the needed atoms
+ double energy;
Tensor box;
Pbc pbc;
Tensor virial;
- std::vector<Vector> forces; // forces on the needed atoms
std::vector<double> masses;
std::vector<double> charges;
+ std::vector<Vector> forces; // forces on the needed atoms
+ double forceOnEnergy;
+
protected:
/// Request an array of atoms.
/// This method is used to ask for a list of atoms. Atoms
@@ -41,6 +44,8 @@ protected:
const Tensor & getBox()const;
/// Get the array of all positions
const std::vector<Vector> & getPositions()const;
+/// Get energy
+ const double & getEnergy()const;
/// Get mass of i-th atom
double getMasses(int i)const;
/// Get charge of i-th atom
@@ -49,6 +54,8 @@ protected:
std::vector<Vector> & modifyForces();
/// Get a reference to virial array
Tensor & modifyVirial();
+/// Get a reference to force on energy
+ double & modifyForceOnEnergy();
/// Get number of available atoms
int getNatoms()const{return indexes.size();};
/// Compute the pbc distance between two positions
@@ -102,6 +109,10 @@ const std::vector<Vector> & ActionAtomistic::getPositions()const{
return positions;
}
+inline
+const double & ActionAtomistic::getEnergy()const{
+ return energy;
+}
inline
const Tensor & ActionAtomistic::getBox()const{
@@ -120,11 +131,14 @@ Tensor & ActionAtomistic::modifyVirial(){
inline
void ActionAtomistic::clearOutputForces(){
- for(unsigned i=0;i<forces.size();++i){
- forces[i][0]=0.0;
- forces[i][1]=0.0;
- forces[i][2]=0.0;
- }
+ for(unsigned i=0;i<forces.size();++i)forces[i].clear();
+ forceOnEnergy=0.0;
+}
+
+
+inline
+double & ActionAtomistic::modifyForceOnEnergy(){
+ return forceOnEnergy;
}
inline
diff --git a/src/Atoms.cpp b/src/Atoms.cpp
index 455ba8f3074c9b350aced05dd95b82b0a6acf44f..afaecc1c6ddd62ab1bec65ecfbe3a16fb53cc2ed 100644
--- a/src/Atoms.cpp
+++ b/src/Atoms.cpp
@@ -132,6 +132,7 @@ void Atoms::share(){
}
virial.clear();
for(unsigned i=0;i<gatindex.size();i++) forces[gatindex[i]].clear();
+ forceOnEnergy=0.0;
}
void Atoms::wait(){
@@ -158,7 +159,6 @@ void Atoms::wait(){
}
}
if(collectEnergy) dd.Sum(&energy,1);
- forceOnEnergy=0.0;
}
}
diff --git a/src/Atoms.h b/src/Atoms.h
index c0978b2cdba85610a147e5ef708c164db368e6ae..267b5c8dfda127e9aaf422db8e16b2e6bbd34e32 100644
--- a/src/Atoms.h
+++ b/src/Atoms.h
@@ -44,7 +44,6 @@ class Atoms
public:
- void forceEnergy(double d){forceOnEnergy+=d;};
double getEnergy()const{assert(collectEnergy);return energy;};
void setCollectEnergy(bool b){collectEnergy=b;};
void setMDEnergyUnits(double d){MDEnergyUnits=d;};
diff --git a/src/Colvar.cpp b/src/Colvar.cpp
index ea462de2d6e85e1bbfb4ad7e76888ebb27092f82..66ee814801f4e666409566723479dd5024646df5 100644
--- a/src/Colvar.cpp
+++ b/src/Colvar.cpp
@@ -53,7 +53,7 @@ void Colvar::apply(){
v(2,1)+=force*derivatives[3*nat+7];
v(2,2)+=force*derivatives[3*nat+8];
}
- if(isEnergy) plumed.getAtoms().forceEnergy(getValue(0)->getForce());
+ if(isEnergy) modifyForceOnEnergy()+=getValue(0)->getForce();
}
diff --git a/src/ColvarEnergy.cpp b/src/ColvarEnergy.cpp
index 5b566b76c0ea155a74fba5dfdc3b086ed089b02a..3cf5be604f995f8cbf914715ddd3d5e17d5cbaf2 100644
--- a/src/ColvarEnergy.cpp
+++ b/src/ColvarEnergy.cpp
@@ -61,7 +61,7 @@ components(false)
// calculator
void ColvarEnergy::calculate(){
- setValue(plumed.getAtoms().getEnergy());
+ setValue(getEnergy());
}
}
diff --git a/src/PlumedMain.cpp b/src/PlumedMain.cpp
index 8634aca340b4b6b75a547122491c3105899c66b6..ceea49c00f2906c56edf0934e6037df40e51fe07 100644
--- a/src/PlumedMain.cpp
+++ b/src/PlumedMain.cpp
@@ -360,6 +360,18 @@ void PlumedMain::readInputFile(std::string str){
pilots=actionSet.select<ActionPilot*>();
}
+////////////////////////////////////////////////////////////////////////
+
+void PlumedMain::exit(int c){
+ comm.Abort(c);
+}
+
+Log& PlumedMain::getLog(){
+ return log;
+}
+
+
+
void PlumedMain::calc(){
@@ -367,21 +379,22 @@ void PlumedMain::calc(){
performCalc();
}
-void PlumedMain::prepareDependencies(){
-
- active=false;
+void PlumedMain::prepareCalc(){
+ prepareDependencies();
+ shareData();
+}
- atoms.setCollectEnergy(false);
- for(ActionSet::iterator p=actionSet.begin();p!=actionSet.end();p++){
- (*p)->deactivate();
- if(Colvar *c=dynamic_cast<Colvar*>(*p)) {
- if(c->checkIsEnergy()) atoms.setCollectEnergy(true);
- }
- }
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// here we have the main steps in "calc()"
+// they can be called individually, but the standard thing is to
+// traverse them in this order:
+void PlumedMain::prepareDependencies(){
// activate all the actions which are on step
// activation is recursive and enables also the dependencies
+// for optimization, an "active" flag remains false if no action at all is active
+ active=false;
for(unsigned i=0;i<pilots.size();++i){
if(pilots[i]->onStep()){
pilots[i]->activate();
@@ -390,36 +403,27 @@ void PlumedMain::prepareDependencies(){
};
// allow actions to update their request list before atoms are shared
- for(ActionSet::iterator p=actionSet.begin();p!=actionSet.end();p++)
- if((*p)->isActive()) (*p)->prepare();
-
-/*
-Temporary fix to allow atom requests to be update in "prepare()"
-
-Nasty bug there. Actual requests of lists of atoms were forwarded to Atoms
-object *during* dependence preparation, so that any change made immediately
-after that was not taken into account. I temporarily solved this by
-re-preparing dependencies again after calls to prepare(), but probably the
-workflow has to be adjusted to be more robust.
-*/
- for(unsigned i=0;i<pilots.size();++i){
- if(pilots[i]->onStep()){
- pilots[i]->activate();
- active=true;
- }
- };
+// also, if one of them is the total energy, tell to atoms that energy should be collected
+ bool collectEnergy=false;
+ for(ActionSet::iterator p=actionSet.begin();p!=actionSet.end();p++){
+ if((*p)->isActive()){
+ if(Colvar *c=dynamic_cast<Colvar*>(*p)) {
+ if(c->checkIsEnergy()) collectEnergy=true;
+ }
+ (*p)->prepare();
+ }
+ }
+ atoms.setCollectEnergy(collectEnergy);
+
}
void PlumedMain::shareData(){
- if(active)atoms.share();
+// atom positions are shared (but only if there is something to do)
+ if(!active)return;
+ atoms.share();
}
-void PlumedMain::prepareCalc(){
- prepareDependencies();
- shareData();
-}
-
void PlumedMain::performCalc(){
if(!active)return;
@@ -443,27 +447,22 @@ void PlumedMain::performCalc(){
}
}
-// Finally apply them in reverse order
+// apply them in reverse order
for(ActionSet::reverse_iterator p=actionSet.rbegin();p!=actionSet.rend();++p){
if((*p)->isActive()) (*p)->apply();
ActionAtomistic*a=dynamic_cast<ActionAtomistic*>(*p);
+// still ActionAtomistic has a special treatment, since they may need to add forces on atoms
if(a) if(a->isActive()) a->applyForces();
}
-// And update forces:
+// this is updating the MD copy of the forces
atoms.updateForces();
-}
-
-void PlumedMain::exit(int c){
- comm.Abort(c);
-}
+// Finally switch off all actions (they will be switched on at next step
+ for(ActionSet::iterator p=actionSet.begin();p!=actionSet.end();p++) (*p)->deactivate();
-Log& PlumedMain::getLog(){
- return log;
}
-
-
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////////