Tentative implementation of center-of-mass and of virtual atoms.

There is a generic ActionWithVirtualAtom class which can be used
as the base for Actions able to add new virtual atoms to the
system. Virtual atoms are stored in Atoms, can be referred to using
the label of the generating action. Automatically, they are
assigned an AtomNumber which exceed the real atom numbers
(e.g. if you have 108 real atoms, virtual atoms will have
indexes 108,109,...). They are automatically understood by
the parser provided one uses the parseAtomList() function.
Each instance of this class only generate a single virtual atom.
The nice thing is that, since lists of atoms can be dynamically changed
(e.g. by neighbor lists), dependencies are handled properly
and a virtual atom is only computed of the steps on which it is needed.

I also implemented a practical case, which is center-of-mass
(GenericCOM).

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(){

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);
   }
 }
 

src/ActionWithVirtualAtom.cpp

+#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());
+}
+
+}

src/ActionWithVirtualAtom.h

+#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

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();
+}
 
 }
 

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

src/GenericCOM.cpp

+#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);
+}
+
+}

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);