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