Skip to content
Snippets Groups Projects
Commit f6b95fb4 authored by Giovanni Bussi's avatar Giovanni Bussi
Browse files

Fix #251

Since isnan() is only available in c++11, I prefer not to backport
this so v2.3.

At variance with older commit fba71ad4
I here only perform the check at the first step
parent a702b5e1
No related branches found
No related tags found
No related merge requests found
...@@ -23,6 +23,7 @@ ...@@ -23,6 +23,7 @@
#include "ActionRegister.h" #include "ActionRegister.h"
#include "core/PlumedMain.h" #include "core/PlumedMain.h"
#include "core/Atoms.h" #include "core/Atoms.h"
#include <cmath>
using namespace std; using namespace std;
...@@ -69,6 +70,7 @@ class COM: ...@@ -69,6 +70,7 @@ class COM:
public ActionWithVirtualAtom public ActionWithVirtualAtom
{ {
bool nopbc; bool nopbc;
bool first;
public: public:
explicit COM(const ActionOptions&ao); explicit COM(const ActionOptions&ao);
void calculate(); void calculate();
...@@ -85,7 +87,8 @@ void COM::registerKeywords(Keywords& keys) { ...@@ -85,7 +87,8 @@ void COM::registerKeywords(Keywords& keys) {
COM::COM(const ActionOptions&ao): COM::COM(const ActionOptions&ao):
Action(ao), Action(ao),
ActionWithVirtualAtom(ao), ActionWithVirtualAtom(ao),
nopbc(false) nopbc(false),
first(true)
{ {
vector<AtomNumber> atoms; vector<AtomNumber> atoms;
parseAtomList("ATOMS",atoms); parseAtomList("ATOMS",atoms);
...@@ -110,6 +113,17 @@ void COM::calculate() { ...@@ -110,6 +113,17 @@ void COM::calculate() {
Vector pos; Vector pos;
if(!nopbc) makeWhole(); if(!nopbc) makeWhole();
double mass(0.0); double mass(0.0);
if( first ) {
for(unsigned i=0; i<getNumberOfAtoms(); i++) {
if(isnan(getMass(i))) {
error(
"You are trying to compute a COM but masses are not known.\n"
" If you are using plumed driver, please use the --mc option"
);
}
}
first=false;
}
vector<Tensor> deriv(getNumberOfAtoms()); vector<Tensor> deriv(getNumberOfAtoms());
for(unsigned i=0; i<getNumberOfAtoms(); i++) mass+=getMass(i); for(unsigned i=0; i<getNumberOfAtoms(); i++) mass+=getMass(i);
if( plumed.getAtoms().chargesWereSet() ) { if( plumed.getAtoms().chargesWereSet() ) {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment