diff --git a/src/cltools/Driver.cpp b/src/cltools/Driver.cpp
index 48ddf54b129290574629eae386d1a9c1601a4afe..aa8a9deaffc4b98d096eeabf4062268577a2c372 100644
--- a/src/cltools/Driver.cpp
+++ b/src/cltools/Driver.cpp
@@ -691,7 +691,7 @@ int Driver<real>::main(FILE* in,FILE*out,Communicator& pc) {
}
}
if(intracomm.Get_rank()==0) {
- fprintf(out,"\nDRIVER: Reassigning particle decomposition\n");
+ fprintf(out,"\nDRIVER: Reassigning domain decomposition\n");
}
p.cmd("setAtomsNlocal",&dd_nlocal);
p.cmd("setAtomsGatindex",&dd_gatindex[0]);
diff --git a/src/core/ActionAtomistic.cpp b/src/core/ActionAtomistic.cpp
index 29ed09e3da598535da09afe9e7af4ac27014ce26..db3c257ff7ab8b4336258334bcded1aed1eb9b99 100644
--- a/src/core/ActionAtomistic.cpp
+++ b/src/core/ActionAtomistic.cpp
@@ -77,7 +77,8 @@ void ActionAtomistic::requestAtoms(const vector<AtomNumber> & a) {
// only real atoms are requested to lower level Atoms class
else unique.insert(indexes[i]);
}
-
+ updateUniqueLocal();
+ atoms.unique.clear();
}
Vector ActionAtomistic::pbcDistance(const Vector &v1,const Vector &v2)const {
@@ -282,4 +283,15 @@ void ActionAtomistic::makeWhole() {
}
}
+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);
+ }
+ } else {
+ unique_local.insert(unique.begin(),unique.end());
+ }
+}
+
}
diff --git a/src/core/ActionAtomistic.h b/src/core/ActionAtomistic.h
index cec10f26877788041f0e7931e5951b2e1ed437df..7269a0a1b0d2ae045fbfcfe7bbdd0fdc35501ee8 100644
--- a/src/core/ActionAtomistic.h
+++ b/src/core/ActionAtomistic.h
@@ -41,7 +41,10 @@ class ActionAtomistic :
{
std::vector<AtomNumber> indexes; // the set of needed atoms
+/// unique should be an ordered set since we later create a vector containing the corresponding indexes
std::set<AtomNumber> unique;
+/// unique_local should be an ordered set since we later create a vector containing the corresponding indexes
+ std::set<AtomNumber> unique_local;
std::vector<Vector> positions; // positions of the needed atoms
double energy;
Pbc& pbc;
@@ -146,6 +149,8 @@ public:
void makeWhole();
/// Allow calls to modifyGlobalForce()
void allowToAccessGlobalForces() {atoms.zeroallforces=true;}
+/// updates local unique atoms
+ void updateUniqueLocal();
public:
// virtual functions:
@@ -169,6 +174,7 @@ public:
void lockRequests();
void unlockRequests();
const std::set<AtomNumber> & getUnique()const;
+ const std::set<AtomNumber> & getUniqueLocal()const;
/// Read in an input file containing atom positions and calculate the action for the atomic
/// configuration therin
void readAtomsFromPDB( const PDB& pdb );
@@ -270,6 +276,11 @@ const std::set<AtomNumber> & ActionAtomistic::getUnique()const {
return unique;
}
+inline
+const std::set<AtomNumber> & ActionAtomistic::getUniqueLocal()const {
+ return unique_local;
+}
+
inline
unsigned ActionAtomistic::getTotAtoms()const {
return atoms.positions.size();
diff --git a/src/core/Atoms.cpp b/src/core/Atoms.cpp
index 674fad203228e3071294c10808081f1b9f4db125..ac71203737b445d1bc2725d2984a7301c809a015 100644
--- a/src/core/Atoms.cpp
+++ b/src/core/Atoms.cpp
@@ -108,7 +108,6 @@ void Atoms::setCharges(void*p) {
void Atoms::setVirial(void*p) {
plumed_massert( dataCanBeSet,"setVirial must be called after setStep in MD code interface");
mdatoms->setVirial(p); virialHasBeenSet=true;
-
}
void Atoms::setEnergy(void*p) {
@@ -138,7 +137,6 @@ void Atoms::setForces(void*p,int i) {
}
void Atoms::share() {
- std::set<AtomNumber> unique;
// At first step I scatter all the atoms so as to store their mass and charge
// Notice that this works with the assumption that charges and masses are
// not changing during the simulation!
@@ -146,46 +144,74 @@ void Atoms::share() {
shareAll();
return;
}
- for(unsigned i=0; i<actions.size(); i++) if(actions[i]->isActive()) {
- if(dd && shuffledAtoms>0) {
- unique.insert(actions[i]->getUnique().begin(),actions[i]->getUnique().end());
+
+ if(!(int(gatindex.size())==natoms && shuffledAtoms==0)) {
+ 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]->getUniqueLocal().begin(),actions[i]->getUniqueLocal().end());
+ }
+ }
+ }
+ } else {
+ for(unsigned i=0; i<actions.size(); i++) {
+ if(actions[i]->isActive()) {
+ if(!actions[i]->getUnique().empty()) {
+ atomsNeeded=true;
+ }
}
- if(!actions[i]->getUnique().empty()) atomsNeeded=true;
}
+
+ }
+
share(unique);
}
void Atoms::shareAll() {
- std::set<AtomNumber> unique;
- if(dd && shuffledAtoms>0)
+ 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));
+ } else {
for(int i=0; i<natoms; i++) unique.insert(AtomNumber::index(i));
+ }
atomsNeeded=true;
share(unique);
}
void Atoms::share(const std::set<AtomNumber>& unique) {
plumed_assert( positionsHaveBeenSet==3 && massesHaveBeenSet );
+
virial.zero();
if(zeroallforces || int(gatindex.size())==natoms) {
for(int i=0; i<natoms; i++) forces[i].zero();
} else {
- for(unsigned i=0; i<gatindex.size(); i++) forces[gatindex[i]].zero();
+ for(const auto & p : unique) forces[p.index()].zero();
}
for(unsigned i=getNatoms(); i<positions.size(); i++) forces[i].zero(); // virtual atoms
forceOnEnergy=0.0;
mdatoms->getBox(box);
if(!atomsNeeded) return;
-
atomsNeeded=false;
if(int(gatindex.size())==natoms && shuffledAtoms==0) {
// faster version, which retrieves all atoms
mdatoms->getPositions(0,natoms,positions);
} else {
-// version that picks only atoms that are available on this proc
- mdatoms->getPositions(gatindex,positions);
+ 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());
+ }
+ mdatoms->getPositions(unique,uniq_index,positions);
}
+
+
// how many double per atom should be scattered:
int ndata=3;
if(!massAndChargeOK) {
@@ -203,17 +229,15 @@ void Atoms::share(const std::set<AtomNumber>& unique) {
}
int count=0;
for(const auto & p : unique) {
- if(dd.g2l[p.index()]>=0) {
- dd.indexToBeSent[count]=p.index();
- dd.positionsToBeSent[ndata*count+0]=positions[p.index()][0];
- dd.positionsToBeSent[ndata*count+1]=positions[p.index()][1];
- dd.positionsToBeSent[ndata*count+2]=positions[p.index()][2];
- if(!massAndChargeOK) {
- dd.positionsToBeSent[ndata*count+3]=masses[p.index()];
- dd.positionsToBeSent[ndata*count+4]=charges[p.index()];
- }
- count++;
+ dd.indexToBeSent[count]=p.index();
+ dd.positionsToBeSent[ndata*count+0]=positions[p.index()][0];
+ dd.positionsToBeSent[ndata*count+1]=positions[p.index()][1];
+ dd.positionsToBeSent[ndata*count+2]=positions[p.index()][2];
+ if(!massAndChargeOK) {
+ dd.positionsToBeSent[ndata*count+3]=masses[p.index()];
+ dd.positionsToBeSent[ndata*count+4]=charges[p.index()];
}
+ count++;
}
if(dd.async) {
asyncSent=true;
@@ -297,8 +321,11 @@ void Atoms::updateForces() {
if(forceOnEnergy*forceOnEnergy>epsilon) {
double alpha=1.0-forceOnEnergy;
mdatoms->rescaleForces(gatindex,alpha);
+ mdatoms->updateForces(gatindex,forces);
+ } else {
+ if(int(gatindex.size())==natoms && shuffledAtoms==0) mdatoms->updateForces(gatindex,forces);
+ else mdatoms->updateForces(unique,uniq_index,forces);
}
- mdatoms->updateForces(gatindex,forces);
if( !plumed.novirial && dd.Get_rank()==0 ) {
plumed_assert( virialHasBeenSet );
mdatoms->updateVirial(virial);
@@ -316,11 +343,11 @@ void Atoms::setNatoms(int n) {
}
-void Atoms::add(const ActionAtomistic*a) {
+void Atoms::add(ActionAtomistic*a) {
actions.push_back(a);
}
-void Atoms::remove(const ActionAtomistic*a) {
+void Atoms::remove(ActionAtomistic*a) {
auto f=find(actions.begin(),actions.end(),a);
plumed_massert(f!=actions.end(),"cannot remove an action registered to atoms");
actions.erase(f);
@@ -376,6 +403,12 @@ void Atoms::setAtomsGatindex(int*g,bool fortran) {
dd.Sum(shuffledAtoms);
for(unsigned i=0; i<gatindex.size(); i++) dd.g2l[gatindex[i]]=i;
}
+
+ for(unsigned i=0; i<actions.size(); i++) {
+ // keep in unique only those atoms that are local
+ actions[i]->updateUniqueLocal();
+ }
+ unique.clear();
}
void Atoms::setAtomsContiguous(int start) {
@@ -384,6 +417,11 @@ void Atoms::setAtomsContiguous(int start) {
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;
if(gatindex.size()<natoms) shuffledAtoms=1;
+ for(unsigned i=0; i<actions.size(); i++) {
+ // keep in unique only those atoms that are local
+ actions[i]->updateUniqueLocal();
+ }
+ unique.clear();
}
void Atoms::setRealPrecision(int p) {
diff --git a/src/core/Atoms.h b/src/core/Atoms.h
index e9cccf399fe87a7ac6b33f9a9e133da4a964c86f..85d3940f91607909619dc3bb882f3c7a966071ac 100644
--- a/src/core/Atoms.h
+++ b/src/core/Atoms.h
@@ -47,6 +47,8 @@ class Atoms
friend class ActionAtomistic;
friend class ActionWithVirtualAtom;
int natoms;
+ std::set<AtomNumber> unique;
+ std::vector<unsigned> uniq_index;
std::vector<Vector> positions;
std::vector<Vector> forces;
std::vector<double> masses;
@@ -98,7 +100,7 @@ class Atoms
double kbT;
- std::vector<const ActionAtomistic*> actions;
+ std::vector<ActionAtomistic*> actions;
std::vector<int> gatindex;
bool asyncSent;
@@ -191,8 +193,8 @@ public:
void getFullList(int**);
void clearFullList();
- void add(const ActionAtomistic*);
- void remove(const ActionAtomistic*);
+ void add(ActionAtomistic*);
+ void remove(ActionAtomistic*);
double getEnergy()const {plumed_assert(collectEnergy && energyHasBeenSet); return energy;}
diff --git a/src/core/MDAtoms.cpp b/src/core/MDAtoms.cpp
index 370527d2fbccffe344174ecf586990b6fa7f1b3a..b0de91532978e27b87d1b4f5b219de8e160c1b1d 100644
--- a/src/core/MDAtoms.cpp
+++ b/src/core/MDAtoms.cpp
@@ -70,12 +70,14 @@ public:
}
void getBox(Tensor &)const;
void getPositions(const vector<int>&index,vector<Vector>&positions)const;
+ void getPositions(const std::set<AtomNumber>&index,const vector<unsigned>&i,vector<Vector>&positions)const;
void getPositions(unsigned j,unsigned k,vector<Vector>&positions)const;
void getLocalPositions(std::vector<Vector>&p)const;
void getMasses(const vector<int>&index,vector<double>&)const;
void getCharges(const vector<int>&index,vector<double>&)const;
void updateVirial(const Tensor&)const;
void updateForces(const vector<int>&index,const vector<Vector>&);
+ void updateForces(const std::set<AtomNumber>&index,const vector<unsigned>&i,const vector<Vector>&forces);
void rescaleForces(const vector<int>&index,double factor);
unsigned getRealPrecision()const;
};
@@ -112,6 +114,18 @@ void MDAtomsTyped<T>::getPositions(const vector<int>&index,vector<Vector>&positi
}
}
+template <class T>
+void MDAtomsTyped<T>::getPositions(const std::set<AtomNumber>&index,const vector<unsigned>&i, vector<Vector>&positions)const {
+// cannot be parallelized with omp because access to positions is not ordered
+ unsigned k=0;
+ for(const auto & p : index) {
+ positions[p.index()][0]=px[stride*i[k]]*scalep;
+ positions[p.index()][1]=py[stride*i[k]]*scalep;
+ positions[p.index()][2]=pz[stride*i[k]]*scalep;
+ k++;
+ }
+}
+
template <class T>
void MDAtomsTyped<T>::getPositions(unsigned j,unsigned k,vector<Vector>&positions)const {
#pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(&positions[j],(k-j)))
@@ -151,6 +165,17 @@ void MDAtomsTyped<T>::updateVirial(const Tensor&virial)const {
if(this->virial) for(int i=0; i<3; i++)for(int j=0; j<3; j++) this->virial[3*i+j]+=T(virial(i,j)*scalev);
}
+template <class T>
+void MDAtomsTyped<T>::updateForces(const std::set<AtomNumber>&index,const vector<unsigned>&i,const vector<Vector>&forces) {
+ unsigned k=0;
+ for(const auto & p : index) {
+ fx[stride*i[k]]+=scalef*T(forces[p.index()][0]);
+ fy[stride*i[k]]+=scalef*T(forces[p.index()][1]);
+ fz[stride*i[k]]+=scalef*T(forces[p.index()][2]);
+ k++;
+ }
+}
+
template <class T>
void MDAtomsTyped<T>::updateForces(const vector<int>&index,const vector<Vector>&forces) {
#pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(fx,stride*index.size()))
diff --git a/src/core/MDAtoms.h b/src/core/MDAtoms.h
index 3926d47f13d004e8edb51dd1d14e9e3f3bbe9003..03e8c29c3648cdc13099a662d2bbf8a7d34f6e9a 100644
--- a/src/core/MDAtoms.h
+++ b/src/core/MDAtoms.h
@@ -24,7 +24,9 @@
#include "tools/Tensor.h"
#include "tools/Vector.h"
+#include "tools/AtomNumber.h"
#include <vector>
+#include <set>
#include "tools/Units.h"
namespace PLMD {
@@ -85,6 +87,8 @@ public:
virtual void getPositions(const std::vector<int>&index,std::vector<Vector>&p)const=0;
/// Retrieve all atom positions from index i to index j.
virtual void getPositions(unsigned i,unsigned j,std::vector<Vector>&p)const=0;
+/// Retrieve all atom positions from atom indices and local indices.
+ virtual void getPositions(const std::set<AtomNumber>&index,const std::vector<unsigned>&i,std::vector<Vector>&p)const=0;
/// Retrieve selected masses.
/// The operation is done in such a way that m[index[i]] is equal to the mass of atom i
virtual void getMasses(const std::vector<int>&index,std::vector<double>&m)const=0;
@@ -98,6 +102,9 @@ public:
/// Increment the force on selected atoms.
/// The operation is done in such a way that f[index[i]] is added to the force on atom i
virtual void updateForces(const std::vector<int>&index,const std::vector<Vector>&f)=0;
+/// Increment the force on selected atoms.
+/// The operation is done only for local atoms used in an action
+ virtual void updateForces(const std::set<AtomNumber>&index,const std::vector<unsigned>&i,const std::vector<Vector>&forces)=0;
/// Rescale all the forces, including the virial.
/// It is applied to all atoms with local index going from 0 to index.size()-1
virtual void rescaleForces(const std::vector<int>&index,double factor)=0;
diff --git a/src/core/PlumedMain.cpp b/src/core/PlumedMain.cpp
index dd96581db44c5e97ab805e445481a82d62ac8e8b..0c5ee14150f571c0d9985529af443466603de1f4 100644
--- a/src/core/PlumedMain.cpp
+++ b/src/core/PlumedMain.cpp
@@ -595,8 +595,6 @@ void PlumedMain::prepareDependencies() {
// First switch off all actions
for(const auto & p : actionSet) {
p->deactivate();
- //I think this is already done inside deactivate
- //(*p)->clearOptions();
}
// for optimization, an "active" flag remains false if no action at all is active