Fix on NAMD-like interface.

- when using NAMD-like interface, also update the unique list of atoms.
- moved g2l out of dd, since it is also needed when not using MPI
  (in NAMD-like interface)

Partly implemented by @hanatok, see #359

src/core/ActionAtomistic.cpp

@@ -287,7 +287,7 @@ void ActionAtomistic::updateUniqueLocal() {
   unique_local.clear();
   if(atoms.dd && atoms.shuffledAtoms>0) {
     for(auto pp=unique.begin(); pp!=unique.end(); ++pp) {
-      if(atoms.dd.g2l[pp->index()]>=0) unique_local.insert(*pp);
+      if(atoms.g2l[pp->index()]>=0) unique_local.insert(*pp);
     }
   } else {
     unique_local.insert(unique.begin(),unique.end());

src/core/Atoms.cpp

@@ -175,7 +175,7 @@ void Atoms::shareAll() {
   unique.clear();
   // keep in unique only those atoms that are local
   if(dd && shuffledAtoms>0) {
-    for(int i=0; i<natoms; i++) if(dd.g2l[i]>=0) unique.insert(AtomNumber::index(i));
+    for(int i=0; i<natoms; i++) if(g2l[i]>=0) unique.insert(AtomNumber::index(i));
   } else {
     for(int i=0; i<natoms; i++) unique.insert(AtomNumber::index(i));
   }
@@ -205,10 +205,8 @@ void Atoms::share(const std::set<AtomNumber>& unique) {
   } else {
     uniq_index.clear();
     uniq_index.reserve(unique.size());
-    if(dd && shuffledAtoms>0) {
-      for(const auto & p : unique) uniq_index.push_back(dd.g2l[p.index()]);
-    } else {
-      for(const auto & p : unique) uniq_index.push_back(p.index());
+    if(shuffledAtoms>0) {
+      for(const auto & p : unique) uniq_index.push_back(g2l[p.index()]);
     }
     mdatoms->getPositions(unique,uniq_index,positions);
   }
@@ -370,8 +368,8 @@ void Atoms::DomainDecomposition::enable(Communicator& c) {
 
 void Atoms::setAtomsNlocal(int n) {
   gatindex.resize(n);
+  g2l.resize(natoms,-1);
   if(dd) {
-    dd.g2l.resize(natoms,-1);
 // Since these vectors are sent with MPI by using e.g.
 // &dd.positionsToBeSent[0]
 // we make sure they are non-zero-sized so as to
@@ -381,7 +379,7 @@ void Atoms::setAtomsNlocal(int n) {
     dd.positionsToBeReceived.resize(natoms*5,0.0);
     dd.indexToBeSent.resize(n,0);
     dd.indexToBeReceived.resize(natoms,0);
-  };
+  }
 }
 
 void Atoms::setAtomsGatindex(int*g,bool fortran) {
@@ -392,7 +390,7 @@ void Atoms::setAtomsGatindex(int*g,bool fortran) {
   } else {
     for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=g[i];
   }
-  for(unsigned i=0; i<dd.g2l.size(); i++) dd.g2l[i]=-1;
+  for(unsigned i=0; i<g2l.size(); i++) g2l[i]=-1;
   if( gatindex.size()==natoms ) {
     shuffledAtoms=0;
     for(unsigned i=0; i<gatindex.size(); i++) {
@@ -403,8 +401,8 @@ void Atoms::setAtomsGatindex(int*g,bool fortran) {
   }
   if(dd) {
     dd.Sum(shuffledAtoms);
-    for(unsigned i=0; i<gatindex.size(); i++) dd.g2l[gatindex[i]]=i;
   }
+  for(unsigned i=0; i<gatindex.size(); i++) g2l[gatindex[i]]=i;
 
   for(unsigned i=0; i<actions.size(); i++) {
     // keep in unique only those atoms that are local
@@ -416,8 +414,8 @@ void Atoms::setAtomsGatindex(int*g,bool fortran) {
 void Atoms::setAtomsContiguous(int start) {
   ddStep=plumed.getStep();
   for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=start+i;
-  for(unsigned i=0; i<dd.g2l.size(); i++) dd.g2l[i]=-1;
-  if(dd) for(unsigned i=0; i<gatindex.size(); i++) dd.g2l[gatindex[i]]=i;
+  for(unsigned i=0; i<g2l.size(); i++) g2l[i]=-1;
+  for(unsigned i=0; i<gatindex.size(); i++) g2l[gatindex[i]]=i;
   if(gatindex.size()<natoms) shuffledAtoms=1;
   for(unsigned i=0; i<actions.size(); i++) {
     // keep in unique only those atoms that are local
@@ -465,14 +463,29 @@ double Atoms::getKbT()const {
 
 
 void Atoms::createFullList(int*n) {
-  vector<AtomNumber> fullListTmp;
-  for(unsigned i=0; i<actions.size(); i++) if(actions[i]->isActive())
-      fullListTmp.insert(fullListTmp.end(),actions[i]->getUnique().begin(),actions[i]->getUnique().end());
-  std::sort(fullListTmp.begin(),fullListTmp.end());
-  int nn=std::unique(fullListTmp.begin(),fullListTmp.end())-fullListTmp.begin();
-  fullList.resize(nn);
-  for(int i=0; i<nn; ++i) fullList[i]=fullListTmp[i].index();
-  *n=nn;
+  if(!massAndChargeOK && shareMassAndChargeOnlyAtFirstStep) {
+    *n=natoms;
+    fullList.resize(natoms);
+    for(unsigned i=0; i<natoms; i++) fullList[i]=i;
+  } else {
+// We update here the unique list defined at Atoms::unique.
+// This is not very clear, and probably should be coded differently.
+// Hopefully this fix the longstanding issue with NAMD.
+    unique.clear();
+    for(unsigned i=0; i<actions.size(); i++) {
+      if(actions[i]->isActive()) {
+        if(!actions[i]->getUnique().empty()) {
+          atomsNeeded=true;
+          // unique are the local atoms
+          unique.insert(actions[i]->getUnique().begin(),actions[i]->getUnique().end());
+        }
+      }
+    }
+    fullList.resize(0);
+    fullList.reserve(unique.size());
+    for(const auto & p : unique) fullList.push_back(p.index());
+    *n=fullList.size();
+  }
 }
 
 void Atoms::getFullList(int**x) {

src/core/Atoms.h

@@ -49,6 +49,7 @@ class Atoms
   int natoms;
   std::set<AtomNumber> unique;
   std::vector<unsigned> uniq_index;
+  std::vector<int> g2l;
   std::vector<Vector> positions;
   std::vector<Vector> forces;
   std::vector<double> masses;
@@ -112,7 +113,6 @@ class Atoms
   public:
     bool on;
     bool async;
-    std::vector<int>    g2l;
 
     std::vector<Communicator::Request> mpi_request_positions;
     std::vector<Communicator::Request> mpi_request_index;