diff --git a/src/cltools/Driver.cpp b/src/cltools/Driver.cpp
index a9bd0ecbe408c921bb869a02b1222f53e886cb66..ea62a50f7258c0ea012503d69650a3c1b4dd7e36 100644
--- a/src/cltools/Driver.cpp
+++ b/src/cltools/Driver.cpp
@@ -380,6 +380,7 @@ int Driver<real>::main(FILE* in,FILE*out,Communicator& pc){
}
}
+ p.cmd("setStep",&step);
if(!noatoms){
bool ok=Tools::getline(fp,line);
plumed_massert(ok,"premature end of file");
@@ -437,7 +438,6 @@ int Driver<real>::main(FILE* in,FILE*out,Communicator& pc){
p.cmd("setBox",&cell[0]);
p.cmd("setVirial",&virial[0]);
}
- p.cmd("setStep",&step);
p.cmd("calc");
// this is necessary as only processor zero is adding to the virial:
diff --git a/src/cltools/SimpleMD.cpp b/src/cltools/SimpleMD.cpp
index 1294a025fee807f575bb95b43d92bda666f35fa5..6409cf9d318d1e4a0810c1c5db0fdd0263ac7aae 100644
--- a/src/cltools/SimpleMD.cpp
+++ b/src/cltools/SimpleMD.cpp
@@ -466,7 +466,7 @@ int main(FILE* in,FILE*out,PLMD::Communicator& pc){
list.resize(listsize);
// masses are hard-coded to 1
- for(unsigned i=0;i<natoms;++i) masses[i]=1.0;
+ for(unsigned i=0;i<natoms;++i) masses[i]=1.0;
// energy integral initialized to 0
engint=0.0;
@@ -478,6 +478,7 @@ int main(FILE* in,FILE*out,PLMD::Communicator& pc){
randomize_velocities(natoms,ndim,temperature,masses,velocities,random);
if(plumed){
+ plumed->cmd("setNoVirial");
plumed->cmd("setNatoms",&natoms);
plumed->cmd("setMDEngine","simpleMD");
plumed->cmd("setTimestep",&tstep);
@@ -533,9 +534,9 @@ int main(FILE* in,FILE*out,PLMD::Communicator& pc){
int istepplusone=istep+1;
for(int i=0;i<3;i++)for(int k=0;k<3;k++) cell9[i][k]=0.0;
for(int i=0;i<3;i++) cell9[i][i]=cell[i];
+ plumed->cmd("setStep",&istepplusone);
plumed->cmd("setMasses",&masses[0]);
plumed->cmd("setForces",&forces[0]);
- plumed->cmd("setStep",&istepplusone);
plumed->cmd("setEnergy",&engconf);
plumed->cmd("setPositions",&positions[0]);
plumed->cmd("setBox",cell9);
diff --git a/src/core/ActionAtomistic.cpp b/src/core/ActionAtomistic.cpp
index c612cc0cfdd26e13e3719e606bd306d01a36fe23..0621442b7919ee9580114fe7c3ca2a445ed97cc4 100644
--- a/src/core/ActionAtomistic.cpp
+++ b/src/core/ActionAtomistic.cpp
@@ -182,6 +182,7 @@ void ActionAtomistic::parseAtomList(const std::string&key,const int num, std::ve
void ActionAtomistic::retrieveAtoms(){
box=atoms.box;
pbc.setBox(box);
+ chargesWereSet=atoms.chargesWereSet();
const vector<Vector> & p(atoms.positions);
const vector<double> & c(atoms.charges);
const vector<double> & m(atoms.masses);
diff --git a/src/core/ActionAtomistic.h b/src/core/ActionAtomistic.h
index 48960bacf1f49ea87735ab694ae2fff804d577aa..6ed09140c85a603628f2645debf50edc36e33faa 100644
--- a/src/core/ActionAtomistic.h
+++ b/src/core/ActionAtomistic.h
@@ -47,6 +47,7 @@ class ActionAtomistic :
Pbc& pbc;
Tensor virial;
std::vector<double> masses;
+ bool chargesWereSet;
std::vector<double> charges;
std::vector<Vector> forces; // forces on the needed atoms
@@ -132,6 +133,7 @@ double ActionAtomistic::getMass(int i)const{
inline
double ActionAtomistic::getCharge(int i)const{
+ plumed_assert( chargesWereSet );
return charges[i];
}
diff --git a/src/core/Atoms.cpp b/src/core/Atoms.cpp
index 5d9a5ff27526df4f95a13bd49fac04f54e7eebea..aaaf932c24e943f641e55382b4ce6028a38ad38a 100644
--- a/src/core/Atoms.cpp
+++ b/src/core/Atoms.cpp
@@ -36,8 +36,15 @@ class PlumedMain;
Atoms::Atoms(PlumedMain&plumed):
natoms(0),
energy(0.0),
+ dataCanBeSet(false),
collectEnergy(0.0),
energyHasBeenSet(false),
+ positionsHaveBeenSet(0),
+ massesHaveBeenSet(false),
+ chargesHaveBeenSet(false),
+ boxHasBeenSet(false),
+ forcesHaveBeenSet(0),
+ virialHasBeenSet(false),
plumed(plumed),
naturalUnits(false),
timestep(0.0),
@@ -53,43 +60,60 @@ Atoms::~Atoms(){
delete mdatoms;
}
+void Atoms::startStep(){
+ collectEnergy=false; energyHasBeenSet=false; positionsHaveBeenSet=0;
+ massesHaveBeenSet=false; chargesHaveBeenSet=false; boxHasBeenSet=false;
+ forcesHaveBeenSet=0; virialHasBeenSet=false; dataCanBeSet=true;
+}
+
void Atoms::setBox(void*p){
mdatoms->setBox(p);
- Tensor b; mdatoms->getBox(b);
+ Tensor b; mdatoms->getBox(b); boxHasBeenSet=true;
}
void Atoms::setPositions(void*p){
- mdatoms->setp(p);
+ plumed_massert( dataCanBeSet ,"setPositions must be called after setStep in MD code interface");
+ mdatoms->setp(p); positionsHaveBeenSet=2;
}
void Atoms::setMasses(void*p){
- mdatoms->setm(p);
+ plumed_massert( dataCanBeSet ,"setMasses must be called after setStep in MD code interface");
+ mdatoms->setm(p); massesHaveBeenSet=true;
+
}
void Atoms::setCharges(void*p){
- mdatoms->setc(p);
+ plumed_massert( dataCanBeSet, "setCharges must be called after setStep in MD code interface");
+ mdatoms->setc(p); chargesHaveBeenSet=true;
}
void Atoms::setVirial(void*p){
- mdatoms->setVirial(p);
+ plumed_massert( dataCanBeSet ,"setVirial must be called after setStep in MD code interface");
+ mdatoms->setVirial(p); virialHasBeenSet=true;
+
}
void Atoms::setEnergy(void*p){
+ plumed_massert( dataCanBeSet ,"setEnergy must be called after setStep in MD code interface");
MD2double(p,energy);
energy*=MDUnits.getEnergy()/units.getEnergy();
energyHasBeenSet=true;
}
void Atoms::setForces(void*p){
+ plumed_massert( dataCanBeSet ,"setForces must be called after setStep in MD code interface");
+ forcesHaveBeenSet=2;
mdatoms->setf(p);
}
void Atoms::setPositions(void*p,int i){
- mdatoms->setp(p,i);
+ plumed_massert( dataCanBeSet ,"setPositions must be called after setStep in MD code interface");
+ mdatoms->setp(p,i); positionsHaveBeenSet++;
}
void Atoms::setForces(void*p,int i){
- mdatoms->setf(p,i);
+ plumed_massert( dataCanBeSet ,"setForces must be called after setStep in MD code interface");
+ mdatoms->setf(p,i); forcesHaveBeenSet++;
}
void Atoms::share(){
@@ -110,6 +134,7 @@ void Atoms::shareAll(){
}
void Atoms::share(const std::set<AtomNumber>& unique){
+ plumed_assert( positionsHaveBeenSet==2 && massesHaveBeenSet );
mdatoms->getBox(box);
mdatoms->getMasses(gatindex,masses);
mdatoms->getCharges(gatindex,charges);
@@ -168,6 +193,7 @@ void Atoms::share(const std::set<AtomNumber>& unique){
}
void Atoms::wait(){
+ dataCanBeSet=false; // Everything should be set by this stage
if(dd){
dd.Bcast(&box[0][0],9,0);
@@ -197,14 +223,16 @@ void Atoms::wait(){
}
void Atoms::updateForces(){
+ plumed_assert( forcesHaveBeenSet==2 );
if(forceOnEnergy*forceOnEnergy>epsilon){
double alpha=1.0-forceOnEnergy;
mdatoms->rescaleForces(gatindex,alpha);
}
- energyHasBeenSet=false;
mdatoms->updateForces(gatindex,forces);
- if(!plumed.novirial && dd.Get_rank()==0) mdatoms->updateVirial(virial);
-
+ if( !plumed.novirial && dd.Get_rank()==0 ){
+ plumed_assert( virialHasBeenSet );
+ mdatoms->updateVirial(virial);
+ }
}
void Atoms::setNatoms(int n){
diff --git a/src/core/Atoms.h b/src/core/Atoms.h
index ebff0c099585fb11480889d2a84995bd14d01776..50f044c043f2c9521873f3d61f11032d001bf1fd 100644
--- a/src/core/Atoms.h
+++ b/src/core/Atoms.h
@@ -54,8 +54,16 @@ class Atoms
Tensor box;
Tensor virial;
double energy;
+
+ bool dataCanBeSet;
bool collectEnergy;
bool energyHasBeenSet;
+ unsigned positionsHaveBeenSet;
+ bool massesHaveBeenSet;
+ bool chargesHaveBeenSet;
+ bool boxHasBeenSet;
+ unsigned forcesHaveBeenSet;
+ bool virialHasBeenSet;
std::map<std::string,std::vector<AtomNumber> > groups;
@@ -133,6 +141,7 @@ public:
void setAtomsContiguous(int);
void setAtomsNlocal(int);
+ void startStep();
void setEnergy(void*);
void setBox(void*);
void setVirial(void*);
@@ -142,6 +151,8 @@ public:
void setForces(void*,int);
void setMasses(void*);
void setCharges(void*);
+ bool chargesWereSet() const ;
+ bool boxWasSet() const ;
void MD2double(const void*m,double&d)const;
void double2MD(const double&d,void*m)const;
@@ -199,6 +210,15 @@ bool Atoms::usingNaturalUnits() const {
return naturalUnits;
}
+inline
+bool Atoms::chargesWereSet() const {
+ return chargesHaveBeenSet;
+}
+
+inline
+bool Atoms::boxWasSet() const {
+ return boxHasBeenSet;
+}
}
diff --git a/src/core/PlumedMain.cpp b/src/core/PlumedMain.cpp
index a888f1c816cad1c74b01d79f7bdae4eb35b6f78a..9c2975bf0dc107b48df0edda7585dabd3fa32be8 100644
--- a/src/core/PlumedMain.cpp
+++ b/src/core/PlumedMain.cpp
@@ -177,6 +177,7 @@ void PlumedMain::cmd(const std::string & word,void*val){
CHECK_INIT(initialized,word);
CHECK_NULL(val,word);
step=(*static_cast<int*>(val));
+ atoms.startStep();
// words used less frequently:
} else if(word=="setAtomsNlocal"){
CHECK_INIT(initialized,word);
@@ -244,6 +245,9 @@ void PlumedMain::cmd(const std::string & word,void*val){
CHECK_NOTINIT(initialized,word);
CHECK_NULL(val,word);
atoms.setMDNaturalUnits(true);
+ } else if(word=="setNoVirial"){
+ CHECK_NOTINIT(initialized,word);
+ novirial=true;
} else if(word=="setPlumedDat"){
CHECK_NOTINIT(initialized,word);
CHECK_NULL(val,word);
@@ -476,7 +480,7 @@ void PlumedMain::shareData(){
// atom positions are shared (but only if there is something to do)
if(!active)return;
stopwatch.start("2 Sharing data");
- atoms.share();
+ if(atoms.getNatoms()>0) atoms.share();
stopwatch.stop("2 Sharing data");
}
@@ -489,7 +493,7 @@ void PlumedMain::performCalc(){
void PlumedMain::waitData(){
if(!active)return;
stopwatch.start("3 Waiting for data");
- atoms.wait();
+ if(atoms.getNatoms()>0) atoms.wait();
stopwatch.stop("3 Waiting for data");
}
@@ -537,7 +541,7 @@ void PlumedMain::justApply(){
}
// this is updating the MD copy of the forces
- atoms.updateForces();
+ if(atoms.getNatoms()>0) atoms.updateForces();
// update step (for statistics, etc)
for(ActionSet::iterator p=actionSet.begin();p!=actionSet.end();++p){
diff --git a/src/core/PlumedMain.h b/src/core/PlumedMain.h
index a498189880ff5cbea8614e25520076ed0a646274..b95bffcb7bfa3fc7de393a009cad9227c7edf426 100644
--- a/src/core/PlumedMain.h
+++ b/src/core/PlumedMain.h
@@ -132,7 +132,7 @@ private:
bool stopNow;
public:
-/// Flag to switch off virial calculation (for debug)
+/// Flag to switch off virial calculation (for debug and MD codes with no barostat)
bool novirial;
diff --git a/src/vatom/COM.cpp b/src/vatom/COM.cpp
index faebd8c072002543d528e4844e1b1d5be1d20e2c..c8cda37ed0e796b3503529acd400757f4bfd82e5 100644
--- a/src/vatom/COM.cpp
+++ b/src/vatom/COM.cpp
@@ -21,6 +21,8 @@
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
#include "ActionWithVirtualAtom.h"
#include "ActionRegister.h"
+#include "core/PlumedMain.h"
+#include "core/Atoms.h"
using namespace std;
@@ -78,17 +80,22 @@ COM::COM(const ActionOptions&ao):
void COM::calculate(){
Vector pos;
- double mass(0.0),charge(0.0);
+ double mass(0.0);
vector<Tensor> deriv(getNumberOfAtoms());
for(unsigned i=0;i<getNumberOfAtoms();i++) mass+=getMass(i);
- for(unsigned i=0;i<getNumberOfAtoms();i++) charge+=getCharge(i);
+ if( plumed.getAtoms().chargesWereSet() ){
+ double charge(0.0);
+ for(unsigned i=0;i<getNumberOfAtoms();i++) charge+=getCharge(i);
+ setCharge(charge);
+ } else {
+ setCharge(0.0);
+ }
for(unsigned i=0;i<getNumberOfAtoms();i++){
pos+=(getMass(i)/mass)*getPosition(i);
deriv[i]=(getMass(i)/mass)*Tensor::identity();
}
setPosition(pos);
setMass(mass);
- setCharge(charge);
setAtomsDerivatives(deriv);
}