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 754 : Option interpretEnvString(const char* env,const char* str) {
114 754 : if(!str) {
115 : return Option::automatic;
116 : }
117 0 : if(!std::strcmp(str,"yes")) {
118 : return Option::yes;
119 : }
120 0 : 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. The code currently contain
132 : /// some heuristic to decide if the unique list should be used or not.
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 72005 : Option getenvForceUnique() {
139 : static const char* name="PLUMED_FORCE_UNIQUE";
140 72005 : static const auto opt = interpretEnvString(name,std::getenv(name));
141 72005 : return opt;
142 : }
143 :
144 : }
145 :
146 : PLUMED_REGISTER_ACTION(DomainDecomposition,"DOMAIN_DECOMPOSITION")
147 :
148 402 : void DomainDecomposition::DomainComms::enable(Communicator& c) {
149 402 : on=true;
150 402 : Set_comm(c.Get_comm());
151 402 : async=Get_size()<10;
152 402 : if(std::getenv("PLUMED_ASYNC_SHARE")) {
153 4 : std::string s(std::getenv("PLUMED_ASYNC_SHARE"));
154 4 : if(s=="yes") {
155 0 : async=true;
156 4 : } else if(s=="no") {
157 4 : async=false;
158 : } else {
159 0 : plumed_merror("PLUMED_ASYNC_SHARE variable is set to " + s + "; should be yes or no");
160 : }
161 : }
162 402 : }
163 :
164 1218 : void DomainDecomposition::registerKeywords(Keywords& keys) {
165 1218 : ActionForInterface::registerKeywords( keys );
166 2436 : keys.remove("ROLE");
167 1218 : keys.add("compulsory","NATOMS","the number of atoms across all the domains");
168 1218 : keys.add("numbered","VALUE","value to create for the domain decomposition");
169 1218 : keys.add("numbered","UNIT","unit of the ith value that is being passed through the domain decomposition");
170 1218 : keys.add("numbered","CONSTANT","is the ith value that is being passed through the domain decomposition constant");
171 1218 : 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 1218 : 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 1218 : keys.add("compulsory","PBCLABEL","Box","the label to use for the PBC action that will be created");
176 2436 : keys.setValueDescription("vector","the domain that each atom is within");
177 1218 : }
178 :
179 1216 : DomainDecomposition::DomainDecomposition(const ActionOptions&ao):
180 : Action(ao),
181 : ActionForInterface(ao),
182 1216 : ddStep(0),
183 1216 : shuffledAtoms(0),
184 1216 : asyncSent(false),
185 1216 : unique_serial(false) {
186 : // Read in the number of atoms
187 : int natoms;
188 2432 : parse("NATOMS",natoms);
189 : std::string str_natoms;
190 1216 : Tools::convert( natoms, str_natoms );
191 : // Setup the gat index
192 1216 : gatindex.resize(natoms);
193 9201219 : for(unsigned i=0; i<gatindex.size(); i++) {
194 9200003 : gatindex[i]=i;
195 : }
196 : // Now read in the values we are creating here
197 6080 : for(unsigned i=1;; ++i) {
198 : std::string valname;
199 14592 : if( !parseNumbered("VALUE",i,valname) ) {
200 : break;
201 : }
202 : std::string unit;
203 12160 : parseNumbered("UNIT",i,unit);
204 : std::string constant;
205 12160 : parseNumbered("CONSTANT",i,constant);
206 : std::string period;
207 12160 : parseNumbered("PERIODIC",i,period);
208 : std::string role;
209 12160 : parseNumbered("ROLE",i,role);
210 6080 : if( constant=="True") {
211 4864 : plumed.readInputLine( valname + ": PUT FROM_DOMAINS CONSTANT SHAPE=" + str_natoms + " UNIT=" + unit + " PERIODIC=" + period + " ROLE=" + role, !plumed.hasBeenInitialized() );
212 3648 : } else if( constant=="False") {
213 7296 : plumed.readInputLine( valname + ": PUT FROM_DOMAINS SHAPE=" + str_natoms + " UNIT=" + unit + " PERIODIC=" + period + " ROLE=" + role, !plumed.hasBeenInitialized() );
214 : } else {
215 0 : plumed_merror("missing information on whether value is constant");
216 : }
217 : // And save the list of values that are set from here
218 6080 : ActionToPutData* ap=plumed.getActionSet().selectWithLabel<ActionToPutData*>(valname);
219 6080 : ap->addDependency( this );
220 6080 : inputs.push_back( ap );
221 6080 : }
222 : std::string pbclabel;
223 1216 : parse("PBCLABEL",pbclabel);
224 1216 : plumed.readInputLine(pbclabel + ": PBC",true);
225 : // Turn on the domain decomposition
226 1216 : if( Communicator::initialized() ) {
227 401 : Set_comm(comm);
228 : }
229 1216 : }
230 :
231 402 : void DomainDecomposition::Set_comm(Communicator& comm) {
232 402 : dd.enable(comm);
233 402 : setAtomsNlocal(getNumberOfAtoms());
234 402 : setAtomsContiguous(0);
235 402 : if( dd.Get_rank()!=0 ) {
236 145 : ActionToPutData* ap=plumed.getActionSet().selectWithLabel<ActionToPutData*>("Box");
237 145 : ap->noforce=true;
238 : }
239 402 : }
240 :
241 559794 : unsigned DomainDecomposition::getNumberOfAtoms() const {
242 559794 : if( inputs.size()==0 ) {
243 : return 0;
244 : }
245 559794 : return (inputs[0]->getPntrToValue())->getShape()[0];
246 : }
247 :
248 74745 : void DomainDecomposition::resetForStepStart() {
249 448470 : for(const auto & pp : inputs) {
250 373725 : pp->resetForStepStart();
251 : }
252 74745 : }
253 :
254 365676 : void DomainDecomposition::setStart( const std::string& name, const unsigned& sss) {
255 1080761 : for(const auto & pp : inputs) {
256 1080761 : if( pp->getLabel()==name ) {
257 365676 : pp->setStart(name, sss);
258 365676 : return;
259 : }
260 : }
261 0 : plumed_error();
262 : }
263 :
264 365676 : void DomainDecomposition::setStride( const std::string& name, const unsigned& sss) {
265 1080761 : for(const auto & pp : inputs) {
266 1080761 : if( pp->getLabel()==name ) {
267 365676 : pp->setStride(name, sss);
268 365676 : return;
269 : }
270 : }
271 0 : plumed_error();
272 : }
273 :
274 369698 : bool DomainDecomposition::setValuePointer( const std::string& name, const TypesafePtr & val ) {
275 369698 : wasset=true; // Once the domain decomposition stuff is transferred moved the setting of this to where the g2l vector is setup
276 1104887 : for(const auto & pp : inputs) {
277 1100868 : if( pp->setValuePointer( name, val ) ) {
278 : return true;
279 : }
280 : }
281 : return false;
282 : }
283 :
284 224229 : bool DomainDecomposition::setForcePointer( const std::string& name, const TypesafePtr & val ) {
285 448578 : for(const auto & pp : inputs) {
286 448548 : if( pp->setForcePointer( name, val ) ) {
287 : return true;
288 : }
289 : }
290 : return false;
291 : }
292 :
293 1542 : void DomainDecomposition::setAtomsNlocal(int n) {
294 1542 : gatindex.resize(n);
295 1542 : g2l.resize(getNumberOfAtoms(),-1);
296 1542 : if(dd) {
297 : // Since these vectors are sent with MPI by using e.g.
298 : // &dd.positionsToBeSent[0]
299 : // we make sure they are non-zero-sized so as to
300 : // avoid errors when doing boundary check
301 1518 : if(n==0) {
302 8 : n++;
303 : }
304 1518 : std::size_t nvals = inputs.size(), natoms = getNumberOfAtoms();
305 1518 : dd.positionsToBeSent.resize(n*nvals,0.0);
306 1518 : dd.positionsToBeReceived.resize(natoms*nvals,0.0);
307 1518 : dd.indexToBeSent.resize(n,0);
308 1518 : dd.indexToBeReceived.resize(natoms,0);
309 : }
310 1542 : }
311 :
312 988 : void DomainDecomposition::setAtomsGatindex(const TypesafePtr & g,bool fortran) {
313 988 : plumed_massert( g || gatindex.size()==0, "NULL gatindex pointer with non-zero local atoms");
314 988 : auto gg=g.get<const int*>({gatindex.size()});
315 988 : ddStep=getStep();
316 988 : if(fortran) {
317 0 : for(unsigned i=0; i<gatindex.size(); i++) {
318 0 : gatindex[i]=gg[i]-1;
319 : }
320 : } else {
321 22140 : for(unsigned i=0; i<gatindex.size(); i++) {
322 21152 : gatindex[i]=gg[i];
323 : }
324 : }
325 57530 : for(unsigned i=0; i<g2l.size(); i++) {
326 56542 : g2l[i]=-1;
327 : }
328 988 : if( gatindex.size()==getNumberOfAtoms() ) {
329 9 : shuffledAtoms=0;
330 1005 : for(unsigned i=0; i<gatindex.size(); i++) {
331 996 : if( gatindex[i]!=i ) {
332 0 : shuffledAtoms=1;
333 0 : break;
334 : }
335 : }
336 : } else {
337 979 : shuffledAtoms=1;
338 : }
339 988 : if(dd) {
340 964 : dd.Sum(shuffledAtoms);
341 : }
342 22140 : for(unsigned i=0; i<gatindex.size(); i++) {
343 21152 : g2l[gatindex[i]]=i;
344 : }
345 : // keep in unique only those atoms that are local
346 9906 : for(unsigned i=0; i<actions.size(); i++) {
347 8918 : actions[i]->unique_local_needs_update=true;
348 : }
349 : unique.clear();
350 : forced_unique.clear();
351 988 : }
352 :
353 554 : void DomainDecomposition::setAtomsContiguous(int start) {
354 554 : ddStep=plumed.getStep();
355 218327 : for(unsigned i=0; i<gatindex.size(); i++) {
356 217773 : gatindex[i]=start+i;
357 : }
358 230639 : for(unsigned i=0; i<g2l.size(); i++) {
359 230085 : g2l[i]=-1;
360 : }
361 218327 : for(unsigned i=0; i<gatindex.size(); i++) {
362 217773 : g2l[gatindex[i]]=i;
363 : }
364 554 : if(gatindex.size()<getNumberOfAtoms()) {
365 152 : shuffledAtoms=1;
366 : }
367 : // keep in unique only those atoms that are local
368 766 : for(unsigned i=0; i<actions.size(); i++) {
369 212 : actions[i]->unique_local_needs_update=true;
370 : }
371 : unique.clear();
372 : forced_unique.clear();
373 554 : }
374 :
375 1242 : void DomainDecomposition::shareAll() {
376 1242 : unique.clear();
377 1242 : forced_unique.clear();
378 1242 : int natoms = getNumberOfAtoms();
379 1242 : if( dd && shuffledAtoms>0 ) {
380 17616 : for(int i=0; i<natoms; ++i)
381 17430 : if( g2l[i]>=0 ) {
382 8003 : unique.push_back( AtomNumber::index(i) );
383 : }
384 : } else {
385 1056 : unique.resize(natoms);
386 5143184 : for(int i=0; i<natoms; i++) {
387 5142128 : unique[i]=AtomNumber::index(i);
388 : }
389 : }
390 1242 : forced_unique.resize( unique.size() );
391 5151373 : for(unsigned i=0; i<unique.size(); ++i) {
392 5150131 : forced_unique[i] = unique[i];
393 : }
394 1242 : share(unique);
395 1226 : }
396 :
397 73133 : void DomainDecomposition::share() {
398 : // We can no longer set the pointers after the share
399 : bool atomsNeeded=false;
400 438798 : for(const auto & pp : inputs) {
401 365665 : pp->share();
402 : }
403 : // At first step I scatter all the atoms so as to store their mass and charge
404 : // Notice that this works with the assumption that charges and masses are
405 : // not changing during the simulation!
406 73133 : if( firststep ) {
407 1128 : actions = plumed.getActionSet().select<ActionAtomistic*>();
408 1128 : shareAll();
409 1112 : return;
410 : }
411 :
412 72005 : if(getenvForceUnique()==Option::automatic) {
413 : unsigned largest=0;
414 728334 : for(unsigned i=0; i<actions.size(); i++) {
415 656329 : if(actions[i]->isActive()) {
416 : auto l=actions[i]->getUnique().size();
417 451061 : if(l>largest) {
418 76450 : largest=l;
419 : }
420 : }
421 : }
422 72005 : if(largest*2<getNumberOfAtoms()) {
423 23643 : unique_serial=true;
424 : } else {
425 48362 : unique_serial=false;
426 : }
427 0 : } else if(getenvForceUnique()==Option::yes) {
428 0 : unique_serial=true;
429 0 : } else if(getenvForceUnique()==Option::no) {
430 0 : unique_serial=false;
431 : } else {
432 0 : plumed_error();
433 : }
434 :
435 72005 : if(unique_serial || !(int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0)) {
436 289079 : for(unsigned i=0; i<actions.size(); i++) {
437 265048 : if( actions[i]->unique_local_needs_update ) {
438 265778 : actions[i]->updateUniqueLocal( !(dd && shuffledAtoms>0), g2l );
439 : }
440 : }
441 : // Now reset unique for the new step
442 : gch::small_vector<const std::vector<AtomNumber>*,32> forced_vectors;
443 : gch::small_vector<const std::vector<AtomNumber>*,32> nonforced_vectors;
444 : forced_vectors.reserve(actions.size());
445 : nonforced_vectors.reserve(actions.size());
446 289079 : for(unsigned i=0; i<actions.size(); i++) {
447 265048 : if(actions[i]->isActive()) {
448 192023 : if(!actions[i]->getUnique().empty()) {
449 : // unique are the local atoms
450 106730 : if( actions[i]->actionHasForces() ) {
451 191038 : forced_vectors.push_back(&actions[i]->getUniqueLocal());
452 : } else {
453 22422 : nonforced_vectors.push_back(&actions[i]->getUniqueLocal());
454 : }
455 : }
456 : }
457 : }
458 24031 : if( !(forced_vectors.empty() && nonforced_vectors.empty()) ) {
459 : atomsNeeded=true;
460 : }
461 : // Merge the atoms from the atoms that have a force on
462 24031 : unique.clear();
463 24031 : forced_unique.clear();
464 24031 : mergeVectorTools::mergeSortedVectors(forced_vectors,forced_unique);
465 : // Merge all the atoms
466 24031 : nonforced_vectors.push_back( &forced_unique );
467 24031 : mergeVectorTools::mergeSortedVectors(nonforced_vectors,unique);
468 : } else {
469 439255 : for(unsigned i=0; i<actions.size(); i++) {
470 391281 : if(actions[i]->isActive()) {
471 259038 : if(!actions[i]->getUnique().empty()) {
472 : atomsNeeded=true;
473 : }
474 : }
475 : }
476 : }
477 :
478 : // Now we retrieve the atom numbers we need
479 72005 : if( atomsNeeded ) {
480 71300 : share( unique );
481 : }
482 : }
483 :
484 72542 : void DomainDecomposition::share(const std::vector<AtomNumber>& unique) {
485 : // This retrieves what values we need to get
486 : int ndata=0;
487 : std::vector<Value*> values_to_get;
488 72542 : if(!(int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0)) {
489 2212 : uniq_index.resize(unique.size());
490 32705 : for(unsigned i=0; i<unique.size(); i++) {
491 30493 : uniq_index[i]=g2l[unique[i].index()];
492 : }
493 13272 : for(const auto & ip : inputs) {
494 11060 : if( (!ip->fixed || firststep) && ip->wasset ) {
495 6780 : (ip->mydata)->share_data( unique, uniq_index, ip->copyOutput(0) );
496 6780 : values_to_get.push_back(ip->copyOutput(0));
497 6780 : ndata++;
498 : }
499 : }
500 70330 : } else if( unique_serial) {
501 21300 : uniq_index.resize(unique.size());
502 801379 : for(unsigned i=0; i<unique.size(); i++) {
503 780079 : uniq_index[i]=unique[i].index();
504 : }
505 127800 : for(const auto & ip : inputs) {
506 106500 : if( (!ip->fixed || firststep) && ip->wasset ) {
507 63900 : (ip->mydata)->share_data( unique, uniq_index, ip->copyOutput(0) );
508 63900 : values_to_get.push_back(ip->copyOutput(0));
509 63900 : ndata++;
510 : }
511 : }
512 : } else {
513 : // faster version, which retrieves all atoms
514 294100 : for(const auto & ip : inputs) {
515 245086 : if( (!ip->fixed || firststep) && ip->wasset ) {
516 148922 : (ip->mydata)->share_data( 0, getNumberOfAtoms(), ip->copyOutput(0) );
517 148906 : values_to_get.push_back(ip->copyOutput(0));
518 148906 : ndata++;
519 : }
520 : }
521 : }
522 :
523 72526 : if(dd && shuffledAtoms>0) {
524 2188 : if(dd.async) {
525 9598 : for(unsigned i=0; i<dd.mpi_request_positions.size(); i++) {
526 7430 : dd.mpi_request_positions[i].wait();
527 : }
528 9598 : for(unsigned i=0; i<dd.mpi_request_index.size(); i++) {
529 7430 : dd.mpi_request_index[i].wait();
530 : }
531 : }
532 :
533 2188 : int count=0;
534 32623 : for(const auto & p : unique) {
535 30435 : dd.indexToBeSent[count]=p.index();
536 : int dpoint=0;
537 126694 : for(unsigned i=0; i<values_to_get.size(); ++i) {
538 96259 : dd.positionsToBeSent[ndata*count+dpoint]=values_to_get[i]->get(p.index());
539 96259 : dpoint++;
540 : }
541 30435 : count++;
542 : }
543 :
544 2188 : if(dd.async) {
545 2168 : asyncSent=true;
546 2168 : dd.mpi_request_positions.resize(dd.Get_size());
547 2168 : dd.mpi_request_index.resize(dd.Get_size());
548 9802 : for(int i=0; i<dd.Get_size(); i++) {
549 7634 : dd.mpi_request_index[i]=dd.Isend(&dd.indexToBeSent[0],count,i,666);
550 7634 : dd.mpi_request_positions[i]=dd.Isend(&dd.positionsToBeSent[0],ndata*count,i,667);
551 : }
552 : } else {
553 20 : const int n=(dd.Get_size());
554 36 : std::vector<int> counts(n);
555 20 : std::vector<int> displ(n);
556 20 : std::vector<int> counts5(n);
557 20 : std::vector<int> displ5(n);
558 20 : dd.Allgather(count,counts);
559 20 : displ[0]=0;
560 80 : for(int i=1; i<n; ++i) {
561 60 : displ[i]=displ[i-1]+counts[i-1];
562 : }
563 100 : for(int i=0; i<n; ++i) {
564 80 : counts5[i]=counts[i]*ndata;
565 : }
566 100 : for(int i=0; i<n; ++i) {
567 80 : displ5[i]=displ[i]*ndata;
568 : }
569 20 : dd.Allgatherv(&dd.indexToBeSent[0],count,&dd.indexToBeReceived[0],&counts[0],&displ[0]);
570 20 : dd.Allgatherv(&dd.positionsToBeSent[0],ndata*count,&dd.positionsToBeReceived[0],&counts5[0],&displ5[0]);
571 20 : int tot=displ[n-1]+counts[n-1];
572 1620 : for(int i=0; i<tot; i++) {
573 : int dpoint=0;
574 7264 : for(unsigned j=0; j<values_to_get.size(); ++j) {
575 5664 : values_to_get[j]->data[dd.indexToBeReceived[i]] = dd.positionsToBeReceived[ndata*i+dpoint];
576 5664 : dpoint++;
577 : }
578 : }
579 : }
580 : }
581 72526 : }
582 :
583 72518 : void DomainDecomposition::wait() {
584 435108 : for(const auto & ip : inputs) {
585 362590 : ip->dataCanBeSet=false;
586 : }
587 :
588 72518 : if(dd && shuffledAtoms>0) {
589 : int ndata=0;
590 : std::vector<Value*> values_to_set;
591 13128 : for(const auto & ip : inputs) {
592 10940 : if( (!ip->fixed || firststep) && ip->wasset ) {
593 6708 : values_to_set.push_back(ip->copyOutput(0));
594 6708 : ndata++;
595 : }
596 : }
597 :
598 : // receive toBeReceived
599 2188 : if(asyncSent) {
600 : Communicator::Status status;
601 : std::size_t count=0;
602 9802 : for(int i=0; i<dd.Get_size(); i++) {
603 7634 : dd.Recv(&dd.indexToBeReceived[count],dd.indexToBeReceived.size()-count,i,666,status);
604 7634 : int c=status.Get_count<int>();
605 7634 : dd.Recv(&dd.positionsToBeReceived[ndata*count],dd.positionsToBeReceived.size()-ndata*count,i,667);
606 7634 : count+=c;
607 : }
608 68220 : for(int i=0; i<count; i++) {
609 : int dpoint=0;
610 276100 : for(unsigned j=0; j<values_to_set.size(); ++j) {
611 210048 : values_to_set[j]->set(dd.indexToBeReceived[i], dd.positionsToBeReceived[ndata*i+dpoint] );
612 210048 : dpoint++;
613 : }
614 : }
615 2168 : asyncSent=false;
616 : }
617 : }
618 72518 : }
619 :
620 0 : unsigned DomainDecomposition::getNumberOfForcesToRescale() const {
621 0 : return gatindex.size();
622 : }
623 :
624 72394 : void DomainDecomposition::apply() {
625 72394 : std::vector<unsigned> forced_uniq_index(forced_unique.size());
626 72394 : if(!(int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0)) {
627 24159 : for(unsigned i=0; i<forced_unique.size(); i++) {
628 22061 : forced_uniq_index[i]=g2l[forced_unique[i].index()];
629 : }
630 : } else {
631 4186780 : for(unsigned i=0; i<forced_unique.size(); i++) {
632 4116484 : forced_uniq_index[i]=forced_unique[i].index();
633 : }
634 : }
635 434354 : for(const auto & ip : inputs) {
636 361962 : if( !(ip->getPntrToValue())->forcesWereAdded() || ip->noforce ) {
637 213601 : continue;
638 148361 : } else if( ip->wasscaled || (!unique_serial && int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0) ) {
639 85772 : (ip->mydata)->add_force( gatindex, ip->getPntrToValue() );
640 : } else {
641 62589 : (ip->mydata)->add_force( forced_unique, forced_uniq_index, ip->getPntrToValue() );
642 : }
643 : }
644 72392 : }
645 :
646 72392 : void DomainDecomposition::reset() {
647 72392 : if( !unique_serial && int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0 ) {
648 : return;
649 : }
650 : // This is an optimisation to ensure that we don't call std::fill over the whole forces
651 : // array if there are a small number of atoms passed between the MD code and PLUMED
652 23388 : if( dd && shuffledAtoms>0 ) {
653 2074 : getAllActiveAtoms( unique );
654 : }
655 140328 : for(const auto & ip : inputs) {
656 116940 : (ip->copyOutput(0))->clearInputForce( unique );
657 : }
658 : }
659 :
660 114 : void DomainDecomposition::writeBinary(std::ostream&o) {
661 684 : for(const auto & ip : inputs) {
662 570 : ip->writeBinary(o);
663 : }
664 114 : }
665 :
666 114 : void DomainDecomposition::readBinary(std::istream&i) {
667 684 : for(const auto & ip : inputs) {
668 570 : ip->readBinary(i);
669 : }
670 114 : }
671 :
672 68710 : void DomainDecomposition::broadcastToDomains( Value* val ) {
673 68710 : if( dd ) {
674 20496 : dd.Bcast( val->data, 0 );
675 : }
676 68710 : }
677 :
678 3989 : void DomainDecomposition::sumOverDomains( Value* val ) {
679 3989 : if( dd && shuffledAtoms>0 ) {
680 0 : dd.Sum( val->data );
681 : }
682 3989 : }
683 :
684 900 : const long int& DomainDecomposition::getDdStep() const {
685 900 : return ddStep;
686 : }
687 :
688 459 : const std::vector<int>& DomainDecomposition::getGatindex() const {
689 459 : return gatindex;
690 : }
691 :
692 2183 : void DomainDecomposition::getAllActiveAtoms( std::vector<AtomNumber>& u ) {
693 : gch::small_vector<const std::vector<AtomNumber>*,32> vectors;
694 : vectors.reserve(actions.size());
695 21016 : for(unsigned i=0; i<actions.size(); i++) {
696 18833 : if(actions[i]->isActive()) {
697 13033 : if(!actions[i]->getUnique().empty()) {
698 : // unique are the local atoms
699 22912 : vectors.push_back(&actions[i]->getUnique());
700 : }
701 : }
702 : }
703 : u.clear();
704 2183 : mergeVectorTools::mergeSortedVectors(vectors,u);
705 2183 : }
706 :
707 116 : void DomainDecomposition::createFullList(const TypesafePtr & n) {
708 116 : if( firststep ) {
709 7 : int natoms = getNumberOfAtoms();
710 7 : n.set(int(natoms));
711 7 : fullList.resize(natoms);
712 803 : for(unsigned i=0; i<natoms; i++) {
713 796 : fullList[i]=i;
714 : }
715 : } else {
716 : // We update here the unique list defined at Atoms::unique.
717 : // This is not very clear, and probably should be coded differently.
718 : // Hopefully this fix the longstanding issue with NAMD.
719 109 : getAllActiveAtoms( unique );
720 109 : fullList.clear();
721 109 : fullList.reserve(unique.size());
722 5012 : for(const auto & p : unique) {
723 4903 : fullList.push_back(p.index());
724 : }
725 109 : n.set(int(fullList.size()));
726 : }
727 116 : }
728 :
729 116 : void DomainDecomposition::getFullList(const TypesafePtr & x) {
730 : auto xx=x.template get<const int**>();
731 116 : if(!fullList.empty()) {
732 110 : *xx=&fullList[0];
733 : } else {
734 6 : *xx=NULL;
735 : }
736 116 : }
737 :
738 116 : void DomainDecomposition::clearFullList() {
739 116 : fullList.resize(0);
740 116 : }
741 :
742 : }
|
