diff --git a/src/Action.cpp b/src/Action.cpp
index 3da7ff8405417f9c24603cd391f8fed0393451cc..60a45f96f92b9f0e3da2b572a243f0918969c836 100644
--- a/src/Action.cpp
+++ b/src/Action.cpp
@@ -62,8 +62,13 @@ void Action::addDependency(Action*action){
}
void Action::activate(){
- active=true;
+// preparation step is called only the first time an Action is activated.
+// since it could change its dependences (e.g. in an ActionAtomistic which is
+// accessing to a virtual atom), this is done just before dependencies are
+// activated
+ if(!active) prepare();
for(Dependencies::iterator p=after.begin();p!=after.end();p++) (*p)->activate();
+ active=true;
}
void Action::clearDependencies(){
diff --git a/src/ActionAtomistic.cpp b/src/ActionAtomistic.cpp
index 9ed62297b1a7ae1cb20fec787735547f019269c2..6a18a42eea046e00c0d8da4ec5dd9c7ea0b69212 100644
--- a/src/ActionAtomistic.cpp
+++ b/src/ActionAtomistic.cpp
@@ -5,6 +5,7 @@
#include <cassert>
#include "ActionWithValue.h"
#include "Colvar.h"
+#include "ActionWithVirtualAtom.h"
using namespace std;
using namespace PLMD;
@@ -21,6 +22,7 @@ Action(ao)
}
void ActionAtomistic::requestAtoms(const vector<AtomNumber> & a){
+ Atoms&atoms(plumed.getAtoms());
int nat=a.size();
indexes.resize(nat);
for(int i=0;i<nat;i++) indexes[i]=a[i].index();
@@ -28,8 +30,12 @@ void ActionAtomistic::requestAtoms(const vector<AtomNumber> & a){
forces.resize(nat);
masses.resize(nat);
charges.resize(nat);
- unsigned n=plumed.getAtoms().natoms;
- for(unsigned i=0;i<indexes.size();i++) assert(indexes[i]<n);
+ int n=atoms.positions.size();
+ clearDependencies();
+ for(unsigned i=0;i<indexes.size();i++){
+ assert(indexes[i]<n);
+ if(indexes[i]>=atoms.getNatoms()) addDependency(atoms.virtualAtomsActions[indexes[i]-atoms.getNatoms()]);
+ }
unique.clear();
unique.insert(indexes.begin(),indexes.end());
@@ -96,7 +102,21 @@ void ActionAtomistic::parseAtomList(const std::string&key,std::vector<AtomNumber
Tools::interpretRanges(strings);
t.resize(strings.size());
for(unsigned i=0;i<t.size();++i){
- Tools::convert(strings[i],t[i]); // this is converting strings to AtomNumbers
+ bool ok=false;
+ ok=Tools::convert(strings[i],t[i]); // this is converting strings to AtomNumbers
+// here we check if the atom name is the name of an added virtual atom
+ if(!ok){
+ const ActionSet&actionSet(plumed.getActionSet());
+ for(ActionSet::const_iterator a=actionSet.begin();a!=actionSet.end();++a){
+ ActionWithVirtualAtom* c=dynamic_cast<ActionWithVirtualAtom*>(*a);
+ if(c) if(c->getLabel()==strings[i]){
+ ok=true;
+ t[i]=c->getIndex();
+ break;
+ }
+ }
+ }
+ assert(ok);
}
}
diff --git a/src/ActionWithVirtualAtom.cpp b/src/ActionWithVirtualAtom.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2b0a5dbf3a580a037699010f6561022b0a3e171d
--- /dev/null
+++ b/src/ActionWithVirtualAtom.cpp
@@ -0,0 +1,32 @@
+#include "ActionWithVirtualAtom.h"
+#include "Atoms.h"
+#include "PlumedMain.h"
+
+using namespace std;
+
+namespace PLMD{
+
+ActionWithVirtualAtom::ActionWithVirtualAtom(const ActionOptions&ao):
+ Action(ao),
+ ActionAtomistic(ao)
+{
+ index=plumed.getAtoms().addVirtualAtom(this);
+ AtomNumber a=AtomNumber::index(index);
+ log.printf(" serial associated to this virtual atom is %d\n",a.serial());
+}
+
+ActionWithVirtualAtom::~ActionWithVirtualAtom(){
+ plumed.getAtoms().removeVirtualAtom(this);
+}
+
+void ActionWithVirtualAtom::apply(){
+ const Vector & f(plumed.getAtoms().forces[index]);
+ for(int i=0;i<getNatoms();i++) modifyForces()[i]=matmul(derivatives[i],f);
+}
+
+void ActionWithVirtualAtom::requestAtoms(const std::vector<AtomNumber> & a){
+ ActionAtomistic::requestAtoms(a);
+ derivatives.resize(a.size());
+}
+
+}
diff --git a/src/ActionWithVirtualAtom.h b/src/ActionWithVirtualAtom.h
new file mode 100644
index 0000000000000000000000000000000000000000..f76b2fd10ef550c7c56d5c46d44c8d79c9ab84e3
--- /dev/null
+++ b/src/ActionWithVirtualAtom.h
@@ -0,0 +1,64 @@
+#ifndef __PLUMED_ActionWithVirtualAtom_h
+#define __PLUMED_ActionWithVirtualAtom_h
+
+#include "ActionAtomistic.h"
+#include "AtomNumber.h"
+#include "Vector.h"
+#include "Tensor.h"
+
+namespace PLMD{
+
+/// Class to add a single virtual atom to the system.
+/// (it might be extended to add multiple virtual atoms).
+class ActionWithVirtualAtom:
+ public ActionAtomistic
+{
+ unsigned index;
+ std::vector<Tensor> derivatives;
+ void apply();
+protected:
+/// Set position of the virtual atom
+ void setPosition(const Vector &);
+/// Set its mass
+ void setMass(double);
+/// Set its charge
+ void setCharge(double);
+/// Request atoms on which the calculation depends
+ void requestAtoms(const std::vector<AtomNumber> & a);
+/// Set the derivatives of virtual atom coordinate wrt atoms on which it dependes
+ void setAtomsDerivatives(const std::vector<Tensor> &d);
+public:
+/// Return the atom id of the corresponding virtual atom
+ AtomNumber getIndex()const;
+ ActionWithVirtualAtom(const ActionOptions&ao);
+ ~ActionWithVirtualAtom();
+};
+
+inline
+AtomNumber ActionWithVirtualAtom::getIndex()const{
+ return AtomNumber::index(index);
+}
+
+inline
+void ActionWithVirtualAtom::setPosition(const Vector & pos){
+ plumed.getAtoms().positions[index]=pos;
+}
+
+inline
+void ActionWithVirtualAtom::setMass(double m){
+ plumed.getAtoms().masses[index]=m;
+}
+
+inline
+void ActionWithVirtualAtom::setCharge(double c){
+ plumed.getAtoms().charges[index]=c;
+}
+
+inline
+void ActionWithVirtualAtom::setAtomsDerivatives(const std::vector<Tensor> &d){
+ derivatives=d;
+}
+
+}
+
+#endif
diff --git a/src/Atoms.cpp b/src/Atoms.cpp
index afaecc1c6ddd62ab1bec65ecfbe3a16fb53cc2ed..0dabd8ebc20587da9fe940575abe950df4d75927 100644
--- a/src/Atoms.cpp
+++ b/src/Atoms.cpp
@@ -32,7 +32,7 @@ Atoms::Atoms(PlumedMain&plumed):
Atoms::~Atoms(){
assert(actions.size()==0);
- if(mdatoms) delete mdatoms;
+ delete mdatoms;
}
void Atoms::setBox(void*p){
@@ -132,6 +132,7 @@ void Atoms::share(){
}
virial.clear();
for(unsigned i=0;i<gatindex.size();i++) forces[gatindex[i]].clear();
+ for(unsigned i=getNatoms();i<positions.size();i++) forces[i].clear(); // virtual atoms
forceOnEnergy=0.0;
}
@@ -275,7 +276,30 @@ void Atoms::init(){
}
}
+void Atoms::setDomainDecomposition(PlumedCommunicator& comm){
+ dd.enable(comm);
+}
+void Atoms::resizeVectors(unsigned n){
+ positions.resize(n);
+ forces.resize(n);
+ masses.resize(n);
+ charges.resize(n);
+}
+
+unsigned Atoms::addVirtualAtom(ActionWithVirtualAtom*a){
+ unsigned n=positions.size();
+ resizeVectors(n+1);
+ virtualAtomsActions.push_back(a);
+ return n;
+}
+
+void Atoms::removeVirtualAtom(ActionWithVirtualAtom*a){
+ unsigned n=positions.size();
+ assert(a==virtualAtomsActions[n-1]);
+ resizeVectors(n-1);
+ virtualAtomsActions.pop_back();
+}
}
diff --git a/src/Atoms.h b/src/Atoms.h
index 267b5c8dfda127e9aaf422db8e16b2e6bbd34e32..45afa9dcdb1581c9230d77f5acc21bd6b67f27af 100644
--- a/src/Atoms.h
+++ b/src/Atoms.h
@@ -19,16 +19,20 @@ class Atoms
{
friend class ActionAtomistic;
friend class GenericWholeMolecules;
+ friend class ActionWithVirtualAtom;
int natoms;
std::vector<Vector> positions;
std::vector<Vector> forces;
std::vector<double> masses;
std::vector<double> charges;
+ std::vector<ActionWithVirtualAtom*> virtualAtomsActions;
Tensor box;
Tensor virial;
double energy;
bool collectEnergy;
+ void resizeVectors(unsigned);
+
std::vector<int> fullList;
MDAtomsBase* mdatoms;
@@ -41,29 +45,6 @@ class Atoms
double timestep;
double forceOnEnergy;
-
-public:
-
- double getEnergy()const{assert(collectEnergy);return energy;};
- void setCollectEnergy(bool b){collectEnergy=b;};
- void setMDEnergyUnits(double d){MDEnergyUnits=d;};
- void setMDLengthUnits(double d){MDLengthUnits=d;};
- void setMDTimeUnits(double d){MDTimeUnits=d;};
- double getMDEnergyUnits()const{return MDEnergyUnits;};
- double getMDLengthUnits()const{return MDLengthUnits;};
- double getMDTimeUnits()const{return MDTimeUnits;};
- void setInternalEnergyUnits(double d){internalEnergyUnits=d;};
- void setInternalLengthUnits(double d){internalLengthUnits=d;};
- void setInternalTimeUnits(double d){internalTimeUnits=d;};
- double getInternalEnergyUnits()const{return internalEnergyUnits;};
- double getInternalLengthUnits()const{return internalLengthUnits;};
- double getInternalTimeUnits()const{return internalTimeUnits;};
- void updateUnits();
-
- void setTimeStep(void*);
- double getTimeStep()const;
-
-private:
std::vector<const ActionAtomistic*> actions;
std::vector<int> gatindex;
@@ -91,38 +72,69 @@ private:
DomainDecomposition dd;
public:
+
Atoms(PlumedMain&plumed);
~Atoms();
+
void init();
+
+ void share();
+ void wait();
+ void updateForces();
+
+ void setRealPrecision(int);
+ int getRealPrecision()const;
+
+ void setTimeStep(void*);
+ double getTimeStep()const;
+
void setNatoms(int);
const int & getNatoms()const;
- void updateForces();
+
+ void setCollectEnergy(bool b){collectEnergy=b;};
+
void setDomainDecomposition(PlumedCommunicator&);
void setAtomsGatindex(int*);
void setAtomsContiguous(int);
void setAtomsNlocal(int);
+ void setEnergy(void*);
void setBox(void*);
+ void setVirial(void*);
void setPositions(void*);
void setPositions(void*,int);
- void setMasses(void*);
- void setCharges(void*);
- void setVirial(void*);
void setForces(void*);
void setForces(void*,int);
- void setEnergy(void*);
+ void setMasses(void*);
+ void setCharges(void*);
- void share();
- void wait();
void MD2double(const void*m,double&d)const;
void double2MD(const double&d,void*m)const;
- void setRealPrecision(int);
- int getRealPrecision()const;
+
void createFullList(int*);
void getFullList(int**);
void clearFullList();
+
void add(const ActionAtomistic*);
void remove(const ActionAtomistic*);
+
+ double getEnergy()const{assert(collectEnergy);return energy;};
+
+ void setMDEnergyUnits(double d){MDEnergyUnits=d;};
+ void setMDLengthUnits(double d){MDLengthUnits=d;};
+ void setMDTimeUnits(double d){MDTimeUnits=d;};
+ double getMDEnergyUnits()const{return MDEnergyUnits;};
+ double getMDLengthUnits()const{return MDLengthUnits;};
+ double getMDTimeUnits()const{return MDTimeUnits;};
+ void setInternalEnergyUnits(double d){internalEnergyUnits=d;};
+ void setInternalLengthUnits(double d){internalLengthUnits=d;};
+ void setInternalTimeUnits(double d){internalTimeUnits=d;};
+ double getInternalEnergyUnits()const{return internalEnergyUnits;};
+ double getInternalLengthUnits()const{return internalLengthUnits;};
+ double getInternalTimeUnits()const{return internalTimeUnits;};
+ void updateUnits();
+ unsigned int addVirtualAtom(ActionWithVirtualAtom*);
+ void removeVirtualAtom(ActionWithVirtualAtom*);
};
inline
@@ -130,11 +142,6 @@ const int & Atoms::getNatoms()const{
return natoms;
}
-inline
-void Atoms::setDomainDecomposition(PlumedCommunicator& comm){
- dd.enable(comm);
-}
-
}
#endif
diff --git a/src/GenericCOM.cpp b/src/GenericCOM.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e2060f77f3421ace01d4c073d66f3b60dafd387c
--- /dev/null
+++ b/src/GenericCOM.cpp
@@ -0,0 +1,73 @@
+#include "ActionWithVirtualAtom.h"
+#include "ActionRegister.h"
+
+using namespace std;
+
+namespace PLMD{
+
+//+PLUMEDOC GENERIC COM
+/**
+Calculate the center of mass of a group of atoms
+
+\par Syntax
+\verbatim
+COM LABEL=label ATOMS=x,y,z,...
+\endverbatim
+The center of mass of atoms x,y,z,... is computed and stored in a virtual
+atom which can be accessed through the label "label"
+
+\par Example
+The following input is printing the distance between the
+center of mass of atoms 1,2,3,4,5,6,7 and that of atoms 15,20:
+\verbatim
+COM ATOMS=1-7 LABEL=c1
+COM ATOMS=15,20 LABEL=c2
+DISTANCE ATOMS=c1,c2 LABEL=d1
+PRINT ARG=d1
+\endverbatim
+(See also \ref DISTANCE and \ref PRINT).
+
+*/
+//+ENDPLUMEDOC
+
+
+class GenericCOM:
+ public ActionWithVirtualAtom
+{
+public:
+ GenericCOM(const ActionOptions&ao);
+ void calculate();
+};
+
+PLUMED_REGISTER_ACTION(GenericCOM,"COM")
+
+GenericCOM::GenericCOM(const ActionOptions&ao):
+ Action(ao),
+ ActionWithVirtualAtom(ao)
+{
+ vector<AtomNumber> atoms;
+ parseAtomList("ATOMS",atoms);
+ checkRead();
+ log.printf(" of atoms");
+ for(unsigned i=0;i<atoms.size();++i) log.printf(" %d",atoms[i].serial());
+ log.printf("\n");
+ requestAtoms(atoms);
+}
+
+void GenericCOM::calculate(){
+ Vector pos;
+ double mass(0.0),charge(0.0);
+ vector<Tensor> deriv(getNatoms());
+ for(int i=0;i<getNatoms();i++) mass+=getMasses(i);
+ for(int i=0;i<getNatoms();i++) charge+=getCharges(i);
+ for(int i=0;i<getNatoms();i++){
+ pos+=(getMasses(i)/mass)*getPositions(i);
+ deriv[i]=(getMasses(i)/mass)*Tensor::identity();
+ }
+ setPosition(pos);
+ setMass(mass);
+ setCharge(charge);
+ setAtomsDerivatives(deriv);
+}
+
+}
diff --git a/src/PlumedMain.cpp b/src/PlumedMain.cpp
index 24ab1cfc7daa2f140e739b56a88eeaf6c5632e35..fe1cbc28d78392d797047e2673e82e4e1619f3f3 100644
--- a/src/PlumedMain.cpp
+++ b/src/PlumedMain.cpp
@@ -394,6 +394,10 @@ void PlumedMain::prepareDependencies(){
// activate all the actions which are on step
// activation is recursive and enables also the dependencies
+// before doing that, the prepare() method is called to see if there is some
+// new/changed dependency (up to now, only useful for dependences on virtual atoms,
+// which can be dynamically changed).
+//
// for optimization, an "active" flag remains false if no action at all is active
active=false;
for(unsigned i=0;i<pilots.size();++i){
@@ -403,7 +407,6 @@ void PlumedMain::prepareDependencies(){
}
};
-// allow actions to update their request list before atoms are shared
// 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++){
@@ -411,7 +414,6 @@ void PlumedMain::prepareDependencies(){
if(Colvar *c=dynamic_cast<Colvar*>(*p)) {
if(c->checkIsEnergy()) collectEnergy=true;
}
- (*p)->prepare();
}
}
atoms.setCollectEnergy(collectEnergy);