From 5b21dacdb53355f735282b83e02a9a3676da04e9 Mon Sep 17 00:00:00 2001
From: Giovanni Bussi <giovanni.bussi@gmail.com>
Date: Fri, 15 Jul 2011 11:55:30 +0200
Subject: [PATCH] Further cleanup
I moved around stuff related to the total energy, in such a way that
energy CV is less different from the others. I also clarified the
main flux in PlumedMain with several comments
---
src/ActionAtomistic.cpp | 2 +
src/ActionAtomistic.h | 26 +++++++++---
src/Atoms.cpp | 2 +-
src/Atoms.h | 1 -
src/Colvar.cpp | 2 +-
src/ColvarEnergy.cpp | 2 +-
src/PlumedMain.cpp | 89 ++++++++++++++++++++---------------------
7 files changed, 69 insertions(+), 55 deletions(-)
diff --git a/src/ActionAtomistic.cpp b/src/ActionAtomistic.cpp
index 205bf888e..9f6432817 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 6250f3acc..296c40a00 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 455ba8f30..afaecc1c6 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 c0978b2cd..267b5c8df 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 ea462de2d..66ee81480 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 5b566b76c..3cf5be604 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 8634aca34..ceea49c00 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;
}
-
-
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
--
GitLab