Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "DomainDecomposition.h"
23 : #include "ActionToPutData.h"
24 : #include "ActionAtomistic.h"
25 : #include "ActionRegister.h"
26 : #include "PlumedMain.h"
27 : #include "ActionSet.h"
28 :
29 : #include "small_vector/small_vector.h"
30 : #include "tools/MergeVectorTools.h"
31 :
32 : //+PLUMEDOC ANALYSIS DOMAIN_DECOMPOSITION
33 : /*
34 : Pass domain decomposed properties of atoms to PLUMED
35 :
36 : Data is passed from the MD code to PLUMED by creating [PUT](PUT.md) actions.
37 : These PUT actions take the data from a void pointer that is passed to PLUMED from the MD code and transfer it into a
38 : PLMD::Value object. Passing a void pointer and using a PUT action to convert the data within it
39 : to a PLMD::Value is also used when the atomic positions, masses and charges are sent to PLUMED. However,
40 : there are some other subtleties for these quantities because MD codes use a domain decomposition and scatter the properties of atoms across
41 : multiple domains. We thus need to use the action DOMAIN_DECOMPOSITION when passing these quantities to make sense of the
42 : information in the void pointers that the MD code passes.
43 :
44 : ## Creating a DOMAIN_DECOMPOSITION action
45 :
46 : A DOMAIN_DECOMPOSITION can be created by using a call to plumed.cmd as follows:
47 :
48 : ```c++
49 : plumed.cmd("readInputLine dd: DOMAIN_DECOMPOSITION NATOMS=20 VALUE1=vv UNIT1=length PERIODIC1=NO CONSTANT1=False");
50 : ```
51 :
52 : The DOMAIN_DECOMPOSTION command above creates a PUT action with the label vv. The pointer to the data that needs to be transferred to the PLMD::Value
53 : object that is created by the PUT action is then set by using the command below:
54 :
55 : ```c++
56 : plumed.cmd("setInputValue vv, &val);
57 : ```
58 :
59 : Meanwhile, the pointer to the forces that should be modified is passed as follows:
60 :
61 : ```c++
62 : plumed.cmd("setValueForces vv", force);
63 : ```
64 :
65 : In other words, pointers to values and forces in the MD code are passed to PUT actions that are created by the DOMAIN_DECOMPOSION in
66 : [the way you pass data to other PUT actions](PUT.md).
67 :
68 : The PLMD::Value objects created by a DOMAIN_DECOMPOSITION action are always vectors with the same number of components as atoms in the system. Furthermore, you can create multiple PUT
69 : actions from a single DOMAIN_DECOMPOSITION action. To see why this is useful, consider the following DOMAIN_DECOMPOSITION command:
70 :
71 : ```plumed
72 : gromacs: DOMAIN_DECOMPOSITION ...
73 : NATOMS=2000
74 : VALUE1=myposx UNIT1=length PERIODIC1=NO CONSTANT1=False ROLE1=x
75 : VALUE2=myposy UNIT2=length PERIODIC2=NO CONSTANT2=False ROLE2=y
76 : VALUE3=myposz UNIT3=length PERIODIC3=NO CONSTANT3=False ROLE3=z
77 : VALUE4=myMasses UNIT4=mass PERIODIC4=NO CONSTANT4=True ROLE4=m
78 : VALUE5=myCharges UNIT5=charge PERIODIC5=NO CONSTANT5=True ROLE5=q
79 : PBCLABEL=mybox
80 : ...
81 : ```
82 :
83 : This action is created when you call `plumed.cmd("setNatoms",&natoms)` from gromacs. It makes 5 PLMD::Value called posx, posy, posz, Masses and Charges.
84 : These PLMD::Value objects then hold the x, y and z positions of the atoms and the masses and charges of the atoms. It is important to note that this command will
85 : also, create a PBC_ACTION to hold the cell.
86 :
87 : The ROLE keywords above are only needed because the five quantities passed by the command above play a unique role within PLUMED. If you pass
88 : some other quantities, this instruction is not required. PLMD::ActionAtomistic searches for atomic positions, masses and charges by looking for PUT Actions
89 : that have these five particular roles and for ActionWithVirtualAtom objects.
90 :
91 : ## Differences from regular PUT actions
92 :
93 : PUT actions created from a DOMAIN_DECOMPOSITION action behave differently from other PUT actions. In particular:
94 :
95 : * If a DOMAIN_DECOMPOSITION action creates a PUT action, then the PUT action depends on the DOMAIN_DECOMPOSITION action. ActionToPutData::apply thus does nothing for these PUT actions.
96 : * Similarly, when DOMAIN_DECOMPOSITION actions create PUT actions, data is transferred from the input pointer to the PLMD::Value object by the DOMAIN_DECOMPOSITION action. When ActionToPutData::wait is called for these PUT Actions `wasset=true`, ActionToPutData::wait does nothing.
97 : * Lastly, if a constant PUT action is created by DOMAIN_DECOMPOSITION, the values in the vector are set during the first step of MD rather than during initialisation.
98 :
99 : These differences are necessary because PUT actions cannot understand the information in the pointers that are passed from the MD code. These pointers are only understandable if you know
100 : which atoms are on each processor. This information is only passed to the DOMAIN_DECOMPOSITION action. DOMAIN_DECOMPOSITION must translate the information passed from the MD code before it is
101 : passed back on through PLMD::Value objects created by the PUT actions. DOMAIN_DECOMPOSITION thus keeps pointers to all the PUT actions that it creates. It sets the data in these action's PLMD::Value objects
102 : within `DomainDecomposition::share(const std::vector<AtomNumber>& unique)` and `DomainDecomposition::wait().` The forces on the PUT actions created by the DOMAIN_DECOMPOSITION action are added in `DomainDecomposition::apply()`.
103 :
104 : */
105 : //+ENDPLUMEDOC
106 :
107 : namespace PLMD {
108 :
109 : namespace {
110 :
111 : enum class Option {automatic, no, yes };
112 :
113 1039 : Option interpretEnvString(const char* env,const char* str) {
114 1039 : if(!str) {
115 : return Option::automatic;
116 : }
117 2 : if(!std::strcmp(str,"yes")) {
118 : return Option::yes;
119 : }
120 1 : if(!std::strcmp(str,"no")) {
121 : return Option::no;
122 : }
123 0 : if(!std::strcmp(str,"auto")) {
124 : return Option::automatic;
125 : }
126 0 : plumed_error()<<"Cannot understand env var "<<env<<"\nPossible values: yes/no/auto\nActual value: "<<str;
127 : }
128 :
129 : /// Use unique list of atoms to manipulate forces and positions.
130 : /// A unique list of atoms is used to manipulate forces and positions in MPI parallel runs.
131 : /// In serial runs, this is done if convenient: with auto, the unique list is used
132 : /// when the atoms needed by all the active actions are less than half of the atoms.
133 : /// An env var can be used to override this decision.
134 : /// export PLUMED_FORCE_UNIQUE=yes # enforce using the unique list in serial runs
135 : /// export PLUMED_FORCE_UNIQUE=no # enforce not using the unique list in serial runs
136 : /// export PLUMED_FORCE_UNIQUE=auto # choose heuristically
137 : /// default: auto
138 162823 : Option getenvForceUnique() {
139 : static const char* name="PLUMED_FORCE_UNIQUE";
140 162823 : static const auto opt = interpretEnvString(name,std::getenv(name));
141 162823 : return opt;
142 : }
143 :
144 : }
145 :
146 : PLUMED_REGISTER_ACTION(DomainDecomposition,"DOMAIN_DECOMPOSITION")
147 :
148 1056 : void DomainDecomposition::DomainComms::enable(Communicator& c) {
149 1056 : on=true;
150 1056 : Set_comm(c.Get_comm());
151 1056 : async=Get_size()<10;
152 1056 : if(std::getenv("PLUMED_ASYNC_SHARE")) {
153 20 : std::string s(std::getenv("PLUMED_ASYNC_SHARE"));
154 20 : if(s=="yes") {
155 8 : async=true;
156 12 : } else if(s=="no") {
157 12 : async=false;
158 : } else {
159 0 : plumed_merror("PLUMED_ASYNC_SHARE variable is set to " + s + "; should be yes or no");
160 : }
161 : }
162 1056 : }
163 :
164 2180 : void DomainDecomposition::registerKeywords(Keywords& keys) {
165 2180 : ActionForInterface::registerKeywords( keys );
166 4360 : keys.remove("ROLE");
167 2180 : keys.add("compulsory","NATOMS","the number of atoms across all the domains");
168 2180 : keys.add("numbered","VALUE","value to create for the domain decomposition");
169 2180 : keys.add("numbered","UNIT","unit of the ith value that is being passed through the domain decomposition");
170 2180 : keys.add("numbered","CONSTANT","is the ith value that is being passed through the domain decomposition constant");
171 2180 : keys.add("numbered","PERIODIC","if the value being passed to plumed is periodic then you should specify the periodicity of the function. If the value "
172 : "is not periodic you must state this using PERIODIC=NO. Positions are passed with PERIODIC=NO even though special methods are used "
173 : "to deal with pbc");
174 2180 : keys.add("numbered","ROLE","Get the role this value plays in the code can be x/y/z/m/q to signify that this is x, y, z positions of atoms or masses or charges of atoms");
175 2180 : keys.add("compulsory","PBCLABEL","Box","the label to use for the PBC action that will be created");
176 2180 : keys.remove("NUMERICAL_DERIVATIVES");
177 4360 : keys.setValueDescription("vector","the domain that each atom is within");
178 2180 : }
179 :
180 2176 : DomainDecomposition::DomainDecomposition(const ActionOptions&ao):
181 : Action(ao),
182 : ActionForInterface(ao),
183 2176 : ddStep(0),
184 2176 : shuffledAtoms(0),
185 2176 : asyncSent(false),
186 2176 : unique_serial(false),
187 2176 : g2lIsValid(false),
188 2176 : uniqueListsValid(false),
189 2176 : allActiveValid(false),
190 2176 : atomsNeeded(false) {
191 : // Read in the number of atoms
192 : int natoms;
193 4352 : parse("NATOMS",natoms);
194 : std::string str_natoms;
195 2176 : Tools::convert( natoms, str_natoms );
196 : // Setup the gat index
197 2176 : gatindex.resize(natoms);
198 51286751 : for(unsigned i=0; i<gatindex.size(); i++) {
199 51284575 : gatindex[i]=i;
200 : }
201 : // Now read in the values we are creating here
202 10880 : for(unsigned i=1;; ++i) {
203 : std::string valname;
204 26112 : if( !parseNumbered("VALUE",i,valname) ) {
205 : break;
206 : }
207 : std::string unit;
208 21760 : parseNumbered("UNIT",i,unit);
209 : std::string constant;
210 21760 : parseNumbered("CONSTANT",i,constant);
211 : std::string period;
212 21760 : parseNumbered("PERIODIC",i,period);
213 : std::string ddrole;
214 10880 : parseNumbered("ROLE",i,ddrole);
215 10880 : if( constant=="True") {
216 21760 : plumed.readInputLine( valname + ": PUT FROM_DOMAINS CONSTANT SHAPE=" + str_natoms + " UNIT=" + unit + " PERIODIC=" + period + " ROLE=" + ddrole, !plumed.hasBeenInitialized() );
217 6528 : } else if( constant=="False") {
218 32640 : plumed.readInputLine( valname + ": PUT FROM_DOMAINS SHAPE=" + str_natoms + " UNIT=" + unit + " PERIODIC=" + period + " ROLE=" + ddrole, !plumed.hasBeenInitialized() );
219 : } else {
220 0 : plumed_merror("missing information on whether value is constant");
221 : }
222 : // And save the list of values that are set from here
223 10880 : ActionToPutData* ap=plumed.getActionSet().selectWithLabel<ActionToPutData*>(valname);
224 10880 : ap->addDependency( this );
225 10880 : inputs.push_back( ap );
226 10880 : }
227 : std::string pbclabel;
228 2176 : parse("PBCLABEL",pbclabel);
229 2176 : plumed.readInputLine(pbclabel + ": PBC",true);
230 : // Turn on the domain decomposition
231 2176 : if( Communicator::initialized() ) {
232 1055 : Set_comm(comm);
233 : }
234 2176 : }
235 :
236 1056 : void DomainDecomposition::Set_comm(Communicator& newcomm) {
237 1056 : dd.enable(newcomm);
238 1056 : setAtomsNlocal(getNumberOfAtoms());
239 1056 : setAtomsContiguous(0);
240 1056 : if( dd.Get_rank()!=0 ) {
241 : // The virial computed by PLUMED is added on the first rank only. The virial of the
242 : // MD code is instead spread over all ranks, so each rank still rescales its share (ENERGY)
243 408 : ActionToPutData* ap=plumed.getActionSet().selectWithLabel<ActionToPutData*>("Box");
244 408 : ap->forces_added_elsewhere=true;
245 : }
246 1056 : }
247 :
248 992925 : unsigned DomainDecomposition::getNumberOfAtoms() const {
249 992925 : if( inputs.size()==0 ) {
250 : return 0;
251 : }
252 992925 : return (inputs[0]->getPntrToValue())->getShape()[0];
253 : }
254 :
255 166448 : void DomainDecomposition::resetForStepStart() {
256 998688 : for(const auto & pp : inputs) {
257 832240 : pp->resetForStepStart();
258 : }
259 166448 : }
260 :
261 754709 : void DomainDecomposition::setStart( const std::string& argName, const unsigned& sss) {
262 2108907 : for(const auto & pp : inputs) {
263 2108907 : if( pp->getLabel()==argName ) {
264 754709 : pp->setStart(argName, sss);
265 754709 : return;
266 : }
267 : }
268 0 : plumed_error();
269 : }
270 :
271 754709 : void DomainDecomposition::setStride( const std::string& argName, const unsigned& sss) {
272 2108907 : for(const auto & pp : inputs) {
273 2108907 : if( pp->getLabel()==argName ) {
274 754709 : pp->setStride(argName, sss);
275 754709 : return;
276 : }
277 : }
278 0 : plumed_error();
279 : }
280 :
281 758839 : bool DomainDecomposition::setValuePointer( const std::string& argName, const TypesafePtr & val ) {
282 758839 : wasset=true; // Once the domain decomposition stuff is transferred moved the setting of this to where the g2l vector is setup
283 2133681 : for(const auto & pp : inputs) {
284 2129554 : if( pp->setValuePointer( argName, val ) ) {
285 : return true;
286 : }
287 : }
288 : return false;
289 : }
290 :
291 499326 : bool DomainDecomposition::setForcePointer( const std::string& argName, const TypesafePtr & val ) {
292 998772 : for(const auto & pp : inputs) {
293 998742 : if( pp->setForcePointer( argName, val ) ) {
294 : return true;
295 : }
296 : }
297 : return false;
298 : }
299 :
300 2504 : void DomainDecomposition::setAtomsNlocal(int n) {
301 2504 : if( static_cast<std::size_t>(n)!=gatindex.size() ) {
302 : // the local atoms will be set again by setAtomsGatindex or setAtomsContiguous
303 1202 : g2lIsValid=false;
304 1202 : invalidateUniqueLists();
305 : }
306 2504 : gatindex.resize(n);
307 2504 : g2l.resize(getNumberOfAtoms(),-1);
308 2504 : if(dd) {
309 : // Since these vectors are sent with MPI by using e.g.
310 : // &dd.positionsToBeSent[0]
311 : // we make sure they are non-zero-sized so as to
312 : // avoid errors when doing boundary check
313 2332 : if(n==0) {
314 14 : n++;
315 : }
316 2332 : std::size_t nvals = inputs.size(), natoms = getNumberOfAtoms();
317 2332 : dd.positionsToBeSent.resize(n*nvals,0.0);
318 2332 : dd.positionsToBeReceived.resize(natoms*nvals,0.0);
319 2332 : dd.indexToBeSent.resize(n,0);
320 2332 : dd.indexToBeReceived.resize(natoms,0);
321 : }
322 2504 : }
323 :
324 1296 : void DomainDecomposition::setAtomsGatindex(const TypesafePtr & g,bool fortran) {
325 1296 : plumed_massert( g || gatindex.size()==0, "NULL gatindex pointer with non-zero local atoms");
326 1296 : auto gg=g.get<const int*>({gatindex.size()});
327 1296 : ddStep=getStep();
328 : // MD codes call this at every step, but the local atoms only change when the
329 : // domains are redistributed: in that case the lists of atoms are not rebuilt
330 1296 : bool localAtomsChanged=!g2lIsValid;
331 1296 : const int shift=fortran ? 1 : 0;
332 32792 : for(unsigned i=0; i<gatindex.size(); i++) {
333 31496 : if( gatindex[i]!=gg[i]-shift ) {
334 21729 : gatindex[i]=gg[i]-shift;
335 : localAtomsChanged=true;
336 : }
337 : }
338 1296 : const bool wasShuffled=shuffledAtoms>0;
339 1296 : if( gatindex.size()==getNumberOfAtoms() ) {
340 142 : shuffledAtoms=0;
341 1474 : for(unsigned i=0; i<gatindex.size(); i++) {
342 1455 : if( gatindex[i]!=static_cast<int>(i) ) {
343 123 : shuffledAtoms=1;
344 123 : break;
345 : }
346 : }
347 : } else {
348 1154 : shuffledAtoms=1;
349 : }
350 1296 : if(dd) {
351 1124 : dd.Sum(shuffledAtoms);
352 : }
353 1296 : if( !localAtomsChanged && wasShuffled==(shuffledAtoms>0) ) {
354 : return;
355 : }
356 62892 : for(unsigned i=0; i<g2l.size(); i++) {
357 61707 : g2l[i]=-1;
358 : }
359 25094 : for(unsigned i=0; i<gatindex.size(); i++) {
360 23909 : g2l[gatindex[i]]=i;
361 : }
362 1185 : g2lIsValid=true;
363 : // keep in unique only those atoms that are local
364 10482 : for(unsigned i=0; i<actions.size(); i++) {
365 9297 : actions[i]->unique_local_needs_update=true;
366 : }
367 1185 : unique.clear();
368 1185 : forced_unique.clear();
369 1185 : invalidateUniqueLists();
370 : }
371 :
372 1208 : void DomainDecomposition::setAtomsContiguous(int start) {
373 1208 : ddStep=plumed.getStep();
374 240431 : for(unsigned i=0; i<gatindex.size(); i++) {
375 239223 : gatindex[i]=start+i;
376 : }
377 252743 : for(unsigned i=0; i<g2l.size(); i++) {
378 251535 : g2l[i]=-1;
379 : }
380 240431 : for(unsigned i=0; i<gatindex.size(); i++) {
381 239223 : g2l[gatindex[i]]=i;
382 : }
383 1208 : g2lIsValid=true;
384 1208 : if(gatindex.size()<getNumberOfAtoms()) {
385 152 : shuffledAtoms=1;
386 : }
387 : // keep in unique only those atoms that are local
388 1420 : for(unsigned i=0; i<actions.size(); i++) {
389 212 : actions[i]->unique_local_needs_update=true;
390 : }
391 1208 : unique.clear();
392 1208 : forced_unique.clear();
393 1208 : invalidateUniqueLists();
394 1208 : }
395 :
396 432435 : bool DomainDecomposition::allAtomsAreLocal() const {
397 432435 : return gatindex.size()==getNumberOfAtoms() && shuffledAtoms==0;
398 : }
399 :
400 5650 : void DomainDecomposition::invalidateUniqueLists() {
401 5650 : uniqueListsValid=false;
402 5650 : allActiveValid=false;
403 5650 : }
404 :
405 5531 : void DomainDecomposition::setIndexesInMDArrays() {
406 5531 : const bool allLocal=allAtomsAreLocal();
407 5531 : uniq_index.resize(unique.size());
408 6661588 : for(unsigned i=0; i<unique.size(); i++) {
409 : // only atoms that the MD code passes to this rank can be in this list
410 : plumed_dbg_assert( allLocal || g2l[unique[i].index()]>=0 );
411 13288467 : uniq_index[i]=allLocal ? unique[i].index() : g2l[unique[i].index()];
412 : }
413 5531 : forced_uniq_index.resize(forced_unique.size());
414 6431849 : for(unsigned i=0; i<forced_unique.size(); i++) {
415 : plumed_dbg_assert( allLocal || g2l[forced_unique[i].index()]>=0 );
416 12830875 : forced_uniq_index[i]=allLocal ? forced_unique[i].index() : g2l[forced_unique[i].index()];
417 : }
418 5531 : }
419 :
420 2055 : void DomainDecomposition::shareAll() {
421 2055 : unique.clear();
422 2055 : const unsigned natoms = getNumberOfAtoms();
423 2055 : if( allAtomsAreLocal() ) {
424 1810 : unique.resize(natoms);
425 6220945 : for(unsigned i=0; i<natoms; i++) {
426 6219135 : unique[i]=AtomNumber::index(i);
427 : }
428 : } else {
429 : // only the atoms that the MD code passes to this rank can be read
430 19498 : for(unsigned i=0; i<natoms; ++i)
431 19253 : if( g2l[i]>=0 ) {
432 9302 : unique.push_back( AtomNumber::index(i) );
433 : }
434 : }
435 2055 : forced_unique=unique;
436 2055 : setIndexesInMDArrays();
437 : // these lists are not the ones requested by the actions, so they are rebuilt at the next share
438 2055 : invalidateUniqueLists();
439 2300 : shareValues( allAtomsAreLocal() && !unique_serial );
440 2036 : }
441 :
442 162816 : bool DomainDecomposition::updateUniqueLists() {
443 : // useunique is false only when the local atoms are a subset of the atoms, scattered across the ranks.
444 : // It can only change in setAtomsGatindex or setAtomsContiguous, which ask the actions to update unique_local
445 162816 : const bool useunique=!(dd && shuffledAtoms>0);
446 162816 : bool changed=!uniqueListsValid;
447 162816 : if( actionState.size()!=actions.size() ) {
448 0 : actionState.assign(actions.size(),0);
449 : changed=true;
450 : }
451 162816 : atomsNeeded=false;
452 873173 : for(unsigned i=0; i<actions.size(); i++) {
453 : // atomsNeeded is computed from the global lists of atoms, so it is the same on all ranks
454 : unsigned char state=0;
455 710357 : if( actions[i]->isActive() && !actions[i]->getUnique().empty() ) {
456 411510 : atomsNeeded=true;
457 411510 : state = actions[i]->actionHasForces() ? 2 : 1;
458 : }
459 710357 : if( state!=actionState[i] ) {
460 7309 : actionState[i]=state;
461 : changed=true;
462 : }
463 : // unique_local_needs_update is set when an action requests new atoms or when the local atoms change
464 710357 : if( state>0 && actions[i]->unique_local_needs_update ) {
465 12048 : actions[i]->updateUniqueLocal( useunique, g2l );
466 : changed=true;
467 : }
468 : }
469 162816 : if( !changed ) {
470 : return false;
471 : }
472 : gch::small_vector<const std::vector<AtomNumber>*,32> forced_vectors;
473 : gch::small_vector<const std::vector<AtomNumber>*,32> nonforced_vectors;
474 : forced_vectors.reserve(actions.size());
475 3476 : nonforced_vectors.reserve(actions.size()+1);
476 33259 : for(unsigned i=0; i<actions.size(); i++) {
477 : // unique_local are the local atoms
478 29783 : if( actionState[i]==2 ) {
479 26138 : forced_vectors.push_back(&actions[i]->getUniqueLocal());
480 16714 : } else if( actionState[i]==1 ) {
481 7300 : nonforced_vectors.push_back(&actions[i]->getUniqueLocal());
482 : }
483 : }
484 : // Merge the atoms from the atoms that have a force on
485 3476 : unique.clear();
486 3476 : forced_unique.clear();
487 3476 : mergeVectorTools::mergeSortedVectors(forced_vectors,forced_unique);
488 : // Merge all the atoms
489 3476 : nonforced_vectors.push_back( &forced_unique );
490 3476 : mergeVectorTools::mergeSortedVectors(nonforced_vectors,unique);
491 3476 : setIndexesInMDArrays();
492 3476 : allActiveValid=false;
493 3476 : uniqueListsValid=true;
494 : return true;
495 : }
496 :
497 36458 : bool DomainDecomposition::activeAtomsChangedSinceShare() const {
498 36458 : if( !uniqueListsValid ) {
499 : return true;
500 : }
501 205603 : for(unsigned i=0; i<actions.size(); i++) {
502 : // actions that were active with no atoms at share time and still have none are not considered
503 169269 : if( actions[i]->isActive() && actions[i]->unique_local_needs_update &&
504 27787 : (actionState[i]>0 || !actions[i]->getUnique().empty()) ) {
505 : return true;
506 : }
507 : }
508 : return false;
509 : }
510 :
511 164752 : void DomainDecomposition::share() {
512 : // We can no longer set the pointers after the share
513 988512 : for(const auto & pp : inputs) {
514 823760 : pp->share();
515 : }
516 : // At first step I scatter all the atoms so as to store their mass and charge
517 : // Notice that this works with the assumption that charges and masses are
518 : // not changing during the simulation!
519 164752 : if( firststep ) {
520 1929 : actions = plumed.getActionSet().select<ActionAtomistic*>();
521 1929 : actionState.assign(actions.size(),0);
522 1929 : shareAll();
523 1910 : return;
524 : }
525 :
526 162823 : const bool allLocal=allAtomsAreLocal();
527 162823 : const Option forceUnique=getenvForceUnique();
528 162823 : if( allLocal && forceUnique==Option::no ) {
529 : // all the atoms are copied at once, so the lists of atoms are not needed
530 7 : unique_serial=false;
531 7 : atomsNeeded=false;
532 7 : for(unsigned i=0; i<actions.size(); i++) {
533 7 : if( actions[i]->isActive() && !actions[i]->getUnique().empty() ) {
534 7 : atomsNeeded=true;
535 7 : break;
536 : }
537 : }
538 7 : if( atomsNeeded ) {
539 7 : shareValues( true );
540 : }
541 7 : return;
542 : }
543 :
544 : // The lists are only rebuilt when the set of active actions, the atoms they request,
545 : // or the atoms that are local to this rank change
546 162816 : updateUniqueLists();
547 :
548 162816 : if( allLocal ) {
549 160429 : if( forceUnique==Option::yes ) {
550 7 : unique_serial=true;
551 : } else {
552 : // Copying only the atoms that are needed is worthwhile when they are a small fraction of
553 : // all the atoms, as copying all the atoms reads contiguous memory and needs no indexes.
554 : // The choice is based on the union of the atoms needed by all the active actions.
555 160422 : unique_serial = unique.size()*2<getNumberOfAtoms();
556 : }
557 : }
558 :
559 : // Now we retrieve the atom numbers we need
560 162816 : if( atomsNeeded ) {
561 198449 : shareValues( allLocal && !unique_serial );
562 : }
563 : }
564 :
565 164167 : void DomainDecomposition::shareValues( const bool& allAtoms ) {
566 : // This retrieves what values we need to get. The list is stored, as wait() must receive the same values
567 164167 : values_to_get.clear();
568 984907 : for(const auto & ip : inputs) {
569 820759 : if( (!ip->fixed || firststep) && ip->wasset ) {
570 495588 : if( allAtoms ) {
571 : // faster version, which retrieves all atoms
572 385573 : (ip->mydata)->share_data( 0, getNumberOfAtoms(), ip->copyOutput(0) );
573 : } else {
574 110015 : (ip->mydata)->share_data( unique, uniq_index, gatindex.size(), ip->copyOutput(0) );
575 : }
576 495569 : values_to_get.push_back(ip->copyOutput(0));
577 : }
578 : }
579 164148 : const int ndata=values_to_get.size();
580 :
581 164148 : if(dd && shuffledAtoms>0) {
582 2460 : if(dd.async) {
583 10280 : for(unsigned i=0; i<dd.mpi_request_positions.size(); i++) {
584 7858 : dd.mpi_request_positions[i].wait();
585 : }
586 10280 : for(unsigned i=0; i<dd.mpi_request_index.size(); i++) {
587 7858 : dd.mpi_request_index[i].wait();
588 : }
589 : }
590 :
591 2460 : int count=0;
592 33839 : for(const auto & p : unique) {
593 31379 : dd.indexToBeSent[count]=p.index();
594 : int dpoint=0;
595 131518 : for(unsigned i=0; i<values_to_get.size(); ++i) {
596 100139 : dd.positionsToBeSent[ndata*count+dpoint]=values_to_get[i]->get(p.index());
597 100139 : dpoint++;
598 : }
599 31379 : count++;
600 : }
601 :
602 2460 : if(dd.async) {
603 2422 : asyncSent=true;
604 2422 : dd.mpi_request_positions.resize(dd.Get_size());
605 2422 : dd.mpi_request_index.resize(dd.Get_size());
606 10564 : for(int i=0; i<dd.Get_size(); i++) {
607 8142 : dd.mpi_request_index[i]=dd.Isend(&dd.indexToBeSent[0],count,i,666);
608 8142 : dd.mpi_request_positions[i]=dd.Isend(&dd.positionsToBeSent[0],ndata*count,i,667);
609 : }
610 : } else {
611 38 : const int n=(dd.Get_size());
612 38 : std::vector<int> counts(n);
613 38 : std::vector<int> displ(n);
614 38 : std::vector<int> counts5(n);
615 38 : std::vector<int> displ5(n);
616 38 : dd.Allgather(count,counts);
617 38 : displ[0]=0;
618 116 : for(int i=1; i<n; ++i) {
619 78 : displ[i]=displ[i-1]+counts[i-1];
620 : }
621 154 : for(int i=0; i<n; ++i) {
622 116 : counts5[i]=counts[i]*ndata;
623 : }
624 154 : for(int i=0; i<n; ++i) {
625 116 : displ5[i]=displ[i]*ndata;
626 : }
627 38 : dd.Allgatherv(&dd.indexToBeSent[0],count,&dd.indexToBeReceived[0],&counts[0],&displ[0]);
628 38 : dd.Allgatherv(&dd.positionsToBeSent[0],ndata*count,&dd.positionsToBeReceived[0],&counts5[0],&displ5[0]);
629 38 : int tot=displ[n-1]+counts[n-1];
630 1674 : for(int i=0; i<tot; i++) {
631 : int dpoint=0;
632 7440 : for(unsigned j=0; j<values_to_get.size(); ++j) {
633 5804 : values_to_get[j]->data[dd.indexToBeReceived[i]] = dd.positionsToBeReceived[ndata*i+dpoint];
634 5804 : dpoint++;
635 : }
636 : }
637 38 : }
638 : }
639 164148 : }
640 :
641 164143 : void DomainDecomposition::wait() {
642 984858 : for(const auto & ip : inputs) {
643 820715 : ip->dataCanBeSet=false;
644 : }
645 :
646 164143 : if(dd && shuffledAtoms>0) {
647 : // the values that were sent in share()
648 : const std::vector<Value*> & values_to_set=values_to_get;
649 : const std::size_t ndata=values_to_set.size();
650 :
651 : // receive toBeReceived
652 2460 : if(asyncSent) {
653 : Communicator::Status status;
654 : std::size_t count=0;
655 10564 : for(int i=0; i<dd.Get_size(); i++) {
656 8142 : plumed_assert(count*ndata<=dd.positionsToBeReceived.size() && count<=dd.indexToBeReceived.size());
657 8142 : const std::size_t remaining=dd.indexToBeReceived.size()-count;
658 8142 : const std::size_t remainingData=dd.positionsToBeReceived.size()-ndata*count;
659 : // Empty ranks still send messages after an earlier rank fills the buffer.
660 : // Avoid indexing one past the end, including for zero-sized storage.
661 8142 : int* indices=remaining ? dd.indexToBeReceived.data()+count : nullptr;
662 8142 : double* positionData=remainingData ? dd.positionsToBeReceived.data()+ndata*count : nullptr;
663 8142 : dd.Recv(indices,remaining,i,666,status);
664 8142 : const int c=status.Get_count<int>();
665 8142 : plumed_assert(c>=0 && static_cast<std::size_t>(c)<=remaining);
666 8142 : dd.Recv(positionData,remainingData,i,667);
667 8142 : count+=c;
668 : }
669 70326 : for(unsigned i=0; i<count; i++) {
670 : int dpoint=0;
671 285572 : for(unsigned j=0; j<values_to_set.size(); ++j) {
672 217668 : values_to_set[j]->set(dd.indexToBeReceived[i], dd.positionsToBeReceived[ndata*i+dpoint] );
673 217668 : dpoint++;
674 : }
675 : }
676 2422 : asyncSent=false;
677 : }
678 : }
679 164143 : }
680 :
681 150 : unsigned DomainDecomposition::getNumberOfForcesToRescale() const {
682 150 : return gatindex.size();
683 : }
684 :
685 163939 : void DomainDecomposition::apply() {
686 : // forced_uniq_index was computed together with forced_unique
687 163939 : const bool allAtoms=!unique_serial && allAtomsAreLocal();
688 983609 : for(const auto & ip : inputs) {
689 819675 : if( !(ip->getPntrToValue())->forcesWereAdded() || ip->noforce ) {
690 401578 : continue;
691 418097 : } else if( ip->wasscaled || allAtoms ) {
692 316991 : if( gatindex.size()>0 ) {
693 316991 : ip->checkForcePointer();
694 : }
695 316988 : (ip->mydata)->add_force( gatindex, ip->getPntrToValue() );
696 : } else {
697 101106 : if( forced_unique.size()>0 ) {
698 100050 : ip->checkForcePointer();
699 : }
700 101106 : (ip->mydata)->add_force( forced_unique, forced_uniq_index, gatindex.size(), ip->getPntrToValue() );
701 : }
702 : }
703 163934 : }
704 :
705 163934 : void DomainDecomposition::reset() {
706 163934 : if( !unique_serial && allAtomsAreLocal() ) {
707 : return;
708 : }
709 : // This is an optimisation to ensure that we don't call std::fill over the whole forces
710 : // array if there are a small number of atoms passed between the MD code and PLUMED.
711 : // The forces were added by the active actions on the atoms they requested
712 36458 : const std::vector<AtomNumber>* toclear=&unique;
713 36458 : if( activeAtomsChangedSinceShare() ) {
714 : // an action requested different atoms after share (e.g. a neighbor list update)
715 : // or unique was set by shareAll, so the global list is rebuilt
716 124 : if( uniqueListsValid || (dd && shuffledAtoms>0) ) {
717 116 : getAllActiveAtoms( activeAtomsBuffer );
718 : toclear=&activeAtomsBuffer;
719 : }
720 36334 : } else if( dd && shuffledAtoms>0 ) {
721 : // unique only holds the local atoms, but forces can be on all the atoms of the actions
722 2226 : if( !allActiveValid ) {
723 1226 : getAllActiveAtoms( allActiveUnique );
724 1226 : allActiveValid=true;
725 : }
726 2226 : toclear=&allActiveUnique;
727 : }
728 218748 : for(const auto & ip : inputs) {
729 182290 : (ip->copyOutput(0))->clearInputForce( *toclear );
730 : }
731 : }
732 :
733 126 : void DomainDecomposition::writeBinary(std::ostream&o) {
734 756 : for(const auto & ip : inputs) {
735 630 : ip->writeBinary(o);
736 : }
737 126 : }
738 :
739 126 : void DomainDecomposition::readBinary(std::istream&i) {
740 756 : for(const auto & ip : inputs) {
741 630 : ip->readBinary(i);
742 : }
743 126 : }
744 :
745 160048 : void DomainDecomposition::broadcastToDomains( Value* val ) {
746 160048 : if( dd ) {
747 31067 : dd.Bcast( val->data, 0 );
748 : }
749 160048 : }
750 :
751 3999 : void DomainDecomposition::sumOverDomains( Value* val ) {
752 3999 : if( dd && shuffledAtoms>0 ) {
753 0 : dd.Sum( val->data );
754 : }
755 3999 : }
756 :
757 0 : const long int& DomainDecomposition::getDdStep() const {
758 0 : return ddStep;
759 : }
760 :
761 485 : const std::vector<int>& DomainDecomposition::getGatindex() const {
762 485 : return gatindex;
763 : }
764 :
765 1472 : void DomainDecomposition::getAllActiveAtoms( std::vector<AtomNumber>& u ) {
766 : gch::small_vector<const std::vector<AtomNumber>*,32> vectors;
767 : vectors.reserve(actions.size());
768 12651 : for(unsigned i=0; i<actions.size(); i++) {
769 11179 : if(actions[i]->isActive()) {
770 7716 : if(!actions[i]->getUnique().empty()) {
771 : // unique are the local atoms
772 13608 : vectors.push_back(&actions[i]->getUnique());
773 : }
774 : }
775 : }
776 1472 : u.clear();
777 1472 : mergeVectorTools::mergeSortedVectors(vectors,u);
778 1472 : }
779 :
780 140 : void DomainDecomposition::createFullList(const TypesafePtr & n) {
781 140 : if( firststep ) {
782 10 : int natoms = getNumberOfAtoms();
783 10 : n.set(natoms);
784 10 : fullList.resize(natoms);
785 1106 : for(int i=0; i<natoms; i++) {
786 1096 : fullList[i]=i;
787 : }
788 : } else {
789 : // The list of the atoms needed by the active actions.
790 : // The MD code (e.g. NAMD) then passes only these atoms.
791 130 : getAllActiveAtoms( activeAtomsBuffer );
792 130 : fullList.clear();
793 130 : fullList.reserve(activeAtomsBuffer.size());
794 5816 : for(const auto & p : activeAtomsBuffer) {
795 5686 : fullList.push_back(p.index());
796 : }
797 130 : n.set(int(fullList.size()));
798 : }
799 140 : }
800 :
801 140 : void DomainDecomposition::getFullList(const TypesafePtr & x) {
802 : auto xx=x.template get<const int**>();
803 140 : if(!fullList.empty()) {
804 134 : *xx=&fullList[0];
805 : } else {
806 6 : *xx=NULL;
807 : }
808 140 : }
809 :
810 116 : void DomainDecomposition::clearFullList() {
811 116 : fullList.resize(0);
812 116 : }
813 :
814 : }
|
