From e54c8c877b3e7a7ca667c94f782f9b89cdbb63c9 Mon Sep 17 00:00:00 2001
From: Giovanni Bussi <giovanni.bussi@gmail.com>
Date: Fri, 15 Jul 2011 19:53:38 +0200
Subject: [PATCH] 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                |  7 ++-
 src/ActionAtomistic.cpp       | 26 +++++++++--
 src/ActionWithVirtualAtom.cpp | 32 ++++++++++++++
 src/ActionWithVirtualAtom.h   | 64 +++++++++++++++++++++++++++
 src/Atoms.cpp                 | 26 ++++++++++-
 src/Atoms.h                   | 81 +++++++++++++++++++----------------
 src/GenericCOM.cpp            | 73 +++++++++++++++++++++++++++++++
 src/PlumedMain.cpp            |  6 ++-
 8 files changed, 271 insertions(+), 44 deletions(-)
 create mode 100644 src/ActionWithVirtualAtom.cpp
 create mode 100644 src/ActionWithVirtualAtom.h
 create mode 100644 src/GenericCOM.cpp

diff --git a/src/Action.cpp b/src/Action.cpp
index 3da7ff840..60a45f96f 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 9ed62297b..6a18a42ee 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 000000000..2b0a5dbf3
--- /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 000000000..f76b2fd10
--- /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 afaecc1c6..0dabd8ebc 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 267b5c8df..45afa9dcd 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 000000000..e2060f77f
--- /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 24ab1cfc7..fe1cbc28d 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);
-- 
GitLab