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 : \par Examples
37 :
38 : */
39 : //+ENDPLUMEDOC
40 :
41 : namespace PLMD {
42 :
43 : namespace {
44 :
45 : enum class Option {automatic, no, yes };
46 :
47 792 : Option interpretEnvString(const char* env,const char* str) {
48 792 : if(!str) {
49 : return Option::automatic;
50 : }
51 0 : if(!std::strcmp(str,"yes")) {
52 : return Option::yes;
53 : }
54 0 : if(!std::strcmp(str,"no")) {
55 : return Option::no;
56 : }
57 0 : if(!std::strcmp(str,"auto")) {
58 : return Option::automatic;
59 : }
60 0 : plumed_error()<<"Cannot understand env var "<<env<<"\nPossible values: yes/no/auto\nActual value: "<<str;
61 : }
62 :
63 : /// Use unique list of atoms to manipulate forces and positions.
64 : /// A unique list of atoms is used to manipulate forces and positions in MPI parallel runs.
65 : /// In serial runs, this is done if convenient. The code currently contain
66 : /// some heuristic to decide if the unique list should be used or not.
67 : /// An env var can be used to override this decision.
68 : /// export PLUMED_FORCE_UNIQUE=yes # enforce using the unique list in serial runs
69 : /// export PLUMED_FORCE_UNIQUE=no # enforce not using the unique list in serial runs
70 : /// export PLUMED_FORCE_UNIQUE=auto # choose heuristically
71 : /// default: auto
72 137797 : Option getenvForceUnique() {
73 : static const char* name="PLUMED_FORCE_UNIQUE";
74 137797 : static const auto opt = interpretEnvString(name,std::getenv(name));
75 137797 : return opt;
76 : }
77 :
78 : }
79 :
80 : PLUMED_REGISTER_ACTION(DomainDecomposition,"DOMAIN_DECOMPOSITION")
81 :
82 563 : void DomainDecomposition::DomainComms::enable(Communicator& c) {
83 563 : on=true;
84 563 : Set_comm(c.Get_comm());
85 563 : async=Get_size()<10;
86 563 : if(std::getenv("PLUMED_ASYNC_SHARE")) {
87 20 : std::string s(std::getenv("PLUMED_ASYNC_SHARE"));
88 20 : if(s=="yes") {
89 8 : async=true;
90 12 : } else if(s=="no") {
91 12 : async=false;
92 : } else {
93 0 : plumed_merror("PLUMED_ASYNC_SHARE variable is set to " + s + "; should be yes or no");
94 : }
95 : }
96 563 : }
97 :
98 1538 : void DomainDecomposition::registerKeywords(Keywords& keys) {
99 1538 : ActionForInterface::registerKeywords( keys );
100 1538 : keys.remove("ROLE");
101 3076 : keys.add("compulsory","NATOMS","the number of atoms across all the domains");
102 3076 : keys.add("numbered","VALUE","value to create for the domain decomposition");
103 3076 : keys.add("numbered","UNIT","unit of the ith value that is being passed through the domain decomposition");
104 3076 : keys.add("numbered","CONSTANT","is the ith value that is being passed through the domain decomposition constant");
105 3076 : 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 "
106 : "is not periodic you must state this using PERIODIC=NO. Positions are passed with PERIODIC=NO even though special methods are used "
107 : "to deal with pbc");
108 3076 : 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");
109 3076 : keys.add("compulsory","PBCLABEL","Box","the label to use for the PBC action that will be created");
110 1538 : keys.setValueDescription("the domain that each atom is within");
111 1538 : }
112 :
113 1534 : DomainDecomposition::DomainDecomposition(const ActionOptions&ao):
114 : Action(ao),
115 : ActionForInterface(ao),
116 1534 : ddStep(0),
117 1534 : shuffledAtoms(0),
118 1534 : asyncSent(false),
119 1534 : unique_serial(false) {
120 : // Read in the number of atoms
121 : int natoms;
122 3068 : parse("NATOMS",natoms);
123 : std::string str_natoms;
124 1534 : Tools::convert( natoms, str_natoms );
125 : // Setup the gat index
126 1534 : gatindex.resize(natoms);
127 9214626 : for(unsigned i=0; i<gatindex.size(); i++) {
128 9213092 : gatindex[i]=i;
129 : }
130 : // Now read in the values we are creating here
131 7670 : for(unsigned i=1;; ++i) {
132 : std::string valname;
133 18408 : if( !parseNumbered("VALUE",i,valname) ) {
134 : break;
135 : }
136 : std::string unit;
137 15340 : parseNumbered("UNIT",i,unit);
138 : std::string constant;
139 15340 : parseNumbered("CONSTANT",i,constant);
140 : std::string period;
141 15340 : parseNumbered("PERIODIC",i,period);
142 : std::string role;
143 15340 : parseNumbered("ROLE",i,role);
144 7670 : if( constant=="True") {
145 6136 : plumed.readInputLine( valname + ": PUT FROM_DOMAINS CONSTANT SHAPE=" + str_natoms + " UNIT=" + unit + " PERIODIC=" + period + " ROLE=" + role, true );
146 4602 : } else if( constant=="False") {
147 9204 : plumed.readInputLine( valname + ": PUT FROM_DOMAINS SHAPE=" + str_natoms + " UNIT=" + unit + " PERIODIC=" + period + " ROLE=" + role, true );
148 : } else {
149 0 : plumed_merror("missing information on whether value is constant");
150 : }
151 : // And save the list of values that are set from here
152 7670 : ActionToPutData* ap=plumed.getActionSet().selectWithLabel<ActionToPutData*>(valname);
153 7670 : ap->addDependency( this );
154 7670 : inputs.push_back( ap );
155 7670 : }
156 : std::string pbclabel;
157 1534 : parse("PBCLABEL",pbclabel);
158 1534 : plumed.readInputLine(pbclabel + ": PBC",true);
159 : // Turn on the domain decomposition
160 1534 : if( Communicator::initialized() ) {
161 562 : Set_comm(comm);
162 : }
163 1534 : }
164 :
165 563 : void DomainDecomposition::Set_comm(Communicator& comm) {
166 563 : dd.enable(comm);
167 563 : setAtomsNlocal(getNumberOfAtoms());
168 563 : setAtomsContiguous(0);
169 563 : if( dd.Get_rank()!=0 ) {
170 : // The virial computed by PLUMED is added on the first rank only. The virial of the
171 : // MD code is instead spread over all ranks, so each rank still rescales its share (ENERGY)
172 224 : ActionToPutData* ap=plumed.getActionSet().selectWithLabel<ActionToPutData*>("Box");
173 224 : ap->forces_added_elsewhere=true;
174 : }
175 563 : }
176 :
177 1283474 : unsigned DomainDecomposition::getNumberOfAtoms() const {
178 1283474 : if( inputs.size()==0 ) {
179 : return 0;
180 : }
181 1283474 : return (inputs[0]->getPntrToValue())->getShape()[0];
182 : }
183 :
184 140886 : void DomainDecomposition::resetForStepStart() {
185 845316 : for(const auto & pp : inputs) {
186 704430 : pp->resetForStepStart();
187 : }
188 140886 : }
189 :
190 631254 : void DomainDecomposition::setStart( const std::string& name, const unsigned& sss) {
191 1747241 : for(const auto & pp : inputs) {
192 1747241 : if( pp->getLabel()==name ) {
193 631254 : pp->setStart(name, sss);
194 631254 : return;
195 : }
196 : }
197 0 : plumed_error();
198 : }
199 :
200 631254 : void DomainDecomposition::setStride( const std::string& name, const unsigned& sss) {
201 1747241 : for(const auto & pp : inputs) {
202 1747241 : if( pp->getLabel()==name ) {
203 631254 : pp->setStride(name, sss);
204 631254 : return;
205 : }
206 : }
207 0 : plumed_error();
208 : }
209 :
210 635284 : bool DomainDecomposition::setValuePointer( const std::string& name, const TypesafePtr & val ) {
211 635284 : wasset=true; // Once the domain decomposition stuff is transferred moved the setting of this to where the g2l vector is setup
212 1771415 : for(const auto & pp : inputs) {
213 1767388 : if( pp->setValuePointer( name, val ) ) {
214 : return true;
215 : }
216 : }
217 : return false;
218 : }
219 :
220 422652 : bool DomainDecomposition::setForcePointer( const std::string& name, const TypesafePtr & val ) {
221 845424 : for(const auto & pp : inputs) {
222 845394 : if( pp->setForcePointer( name, val ) ) {
223 : return true;
224 : }
225 : }
226 : return false;
227 : }
228 :
229 1757 : void DomainDecomposition::setAtomsNlocal(int n) {
230 1757 : gatindex.resize(n);
231 1757 : g2l.resize(getNumberOfAtoms(),-1);
232 1757 : if(dd) {
233 : // Since these vectors are sent with MPI by using e.g.
234 : // &dd.positionsToBeSent[0]
235 : // we make sure they are non-zero-sized so as to
236 : // avoid errors when doing boundary check
237 1731 : if(n==0) {
238 12 : n++;
239 : }
240 1731 : std::size_t nvals = inputs.size(), natoms = getNumberOfAtoms();
241 1731 : dd.positionsToBeSent.resize(n*nvals,0.0);
242 1731 : dd.positionsToBeReceived.resize(natoms*nvals,0.0);
243 1731 : dd.indexToBeSent.resize(n,0);
244 1731 : dd.indexToBeReceived.resize(natoms,0);
245 : }
246 1757 : }
247 :
248 1042 : void DomainDecomposition::setAtomsGatindex(const TypesafePtr & g,bool fortran) {
249 1042 : plumed_massert( g || gatindex.size()==0, "NULL gatindex pointer with non-zero local atoms");
250 1042 : auto gg=g.get<const int*>({gatindex.size()});
251 1042 : ddStep=getStep();
252 1042 : if(fortran) {
253 22 : for(unsigned i=0; i<gatindex.size(); i++) {
254 20 : gatindex[i]=gg[i]-1;
255 : }
256 : } else {
257 22352 : for(unsigned i=0; i<gatindex.size(); i++) {
258 21312 : gatindex[i]=gg[i];
259 : }
260 : }
261 57924 : for(unsigned i=0; i<g2l.size(); i++) {
262 56882 : g2l[i]=-1;
263 : }
264 1042 : if( gatindex.size()==getNumberOfAtoms() ) {
265 15 : shuffledAtoms=0;
266 1019 : for(unsigned i=0; i<gatindex.size(); i++) {
267 1006 : if( gatindex[i]!=i ) {
268 2 : shuffledAtoms=1;
269 2 : break;
270 : }
271 : }
272 : } else {
273 1027 : shuffledAtoms=1;
274 : }
275 1042 : if(dd) {
276 1016 : dd.Sum(shuffledAtoms);
277 : }
278 22374 : for(unsigned i=0; i<gatindex.size(); i++) {
279 21332 : g2l[gatindex[i]]=i;
280 : }
281 : // keep in unique only those atoms that are local
282 10048 : for(unsigned i=0; i<actions.size(); i++) {
283 9006 : actions[i]->unique_local_needs_update=true;
284 : }
285 : unique.clear();
286 : forced_unique.clear();
287 1042 : }
288 :
289 715 : void DomainDecomposition::setAtomsContiguous(int start) {
290 715 : ddStep=plumed.getStep();
291 219090 : for(unsigned i=0; i<gatindex.size(); i++) {
292 218375 : gatindex[i]=start+i;
293 : }
294 231402 : for(unsigned i=0; i<g2l.size(); i++) {
295 230687 : g2l[i]=-1;
296 : }
297 219090 : for(unsigned i=0; i<gatindex.size(); i++) {
298 218375 : g2l[gatindex[i]]=i;
299 : }
300 715 : if(gatindex.size()<getNumberOfAtoms()) {
301 152 : shuffledAtoms=1;
302 : }
303 : // keep in unique only those atoms that are local
304 927 : for(unsigned i=0; i<actions.size(); i++) {
305 212 : actions[i]->unique_local_needs_update=true;
306 : }
307 : unique.clear();
308 : forced_unique.clear();
309 715 : }
310 :
311 1546 : void DomainDecomposition::shareAll() {
312 1546 : unique.clear();
313 1546 : forced_unique.clear();
314 1546 : int natoms = getNumberOfAtoms();
315 1546 : if( dd && shuffledAtoms>0 ) {
316 17772 : for(int i=0; i<natoms; ++i)
317 17558 : if( g2l[i]>=0 ) {
318 8067 : unique.push_back( AtomNumber::index(i) );
319 : }
320 : } else {
321 1332 : unique.resize(natoms);
322 5156377 : for(int i=0; i<natoms; i++) {
323 5155045 : unique[i]=AtomNumber::index(i);
324 : }
325 : }
326 1546 : forced_unique.resize( unique.size() );
327 5164658 : for(unsigned i=0; i<unique.size(); ++i) {
328 5163112 : forced_unique[i] = unique[i];
329 : }
330 1546 : share(unique);
331 1530 : }
332 :
333 139229 : void DomainDecomposition::share() {
334 : // We can no longer set the pointers after the share
335 : bool atomsNeeded=false;
336 835374 : for(const auto & pp : inputs) {
337 696145 : pp->share();
338 : }
339 : // At first step I scatter all the atoms so as to store their mass and charge
340 : // Notice that this works with the assumption that charges and masses are
341 : // not changing during the simulation!
342 139229 : if( firststep ) {
343 1432 : actions = plumed.getActionSet().select<ActionAtomistic*>();
344 1432 : shareAll();
345 1416 : return;
346 : }
347 :
348 137797 : if(getenvForceUnique()==Option::automatic) {
349 : unsigned largest=0;
350 929455 : for(unsigned i=0; i<actions.size(); i++) {
351 791658 : if(actions[i]->isActive()) {
352 : auto l=actions[i]->getUnique().size();
353 585951 : if(l>largest) {
354 142269 : largest=l;
355 : }
356 : }
357 : }
358 137797 : if(largest*2<getNumberOfAtoms()) {
359 23970 : unique_serial=true;
360 : } else {
361 113827 : unique_serial=false;
362 : }
363 0 : } else if(getenvForceUnique()==Option::yes) {
364 0 : unique_serial=true;
365 0 : } else if(getenvForceUnique()==Option::no) {
366 0 : unique_serial=false;
367 : } else {
368 0 : plumed_error();
369 : }
370 :
371 137797 : if(unique_serial || !(int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0)) {
372 290164 : for(unsigned i=0; i<actions.size(); i++) {
373 265778 : if( actions[i]->unique_local_needs_update ) {
374 266600 : actions[i]->updateUniqueLocal( !(dd && shuffledAtoms>0), g2l );
375 : }
376 : }
377 : // Now reset unique for the new step
378 : gch::small_vector<const std::vector<AtomNumber>*,32> forced_vectors;
379 : gch::small_vector<const std::vector<AtomNumber>*,32> nonforced_vectors;
380 : forced_vectors.reserve(actions.size());
381 : nonforced_vectors.reserve(actions.size());
382 290164 : for(unsigned i=0; i<actions.size(); i++) {
383 265778 : if(actions[i]->isActive()) {
384 192570 : if(!actions[i]->getUnique().empty()) {
385 : // unique are the local atoms
386 107249 : if( actions[i]->actionHasForces() ) {
387 192002 : forced_vectors.push_back(&actions[i]->getUniqueLocal());
388 : } else {
389 22496 : nonforced_vectors.push_back(&actions[i]->getUniqueLocal());
390 : }
391 : }
392 : }
393 : }
394 24386 : if( !(forced_vectors.empty() && nonforced_vectors.empty()) ) {
395 : atomsNeeded=true;
396 : }
397 : // Merge the atoms from the atoms that have a force on
398 24386 : unique.clear();
399 24386 : forced_unique.clear();
400 24386 : mergeVectorTools::mergeSortedVectors(forced_vectors,forced_unique);
401 : // Merge all the atoms
402 24386 : nonforced_vectors.push_back( &forced_unique );
403 24386 : mergeVectorTools::mergeSortedVectors(nonforced_vectors,unique);
404 : } else {
405 639291 : for(unsigned i=0; i<actions.size(); i++) {
406 525880 : if(actions[i]->isActive()) {
407 393381 : if(!actions[i]->getUnique().empty()) {
408 : atomsNeeded=true;
409 : }
410 : }
411 : }
412 : }
413 :
414 : // Now we retrieve the atom numbers we need
415 137797 : if( atomsNeeded ) {
416 137088 : share( unique );
417 : }
418 : }
419 :
420 138634 : void DomainDecomposition::share(const std::vector<AtomNumber>& unique) {
421 : // This retrieves what values we need to get
422 : int ndata=0;
423 : std::vector<Value*> values_to_get;
424 138634 : if(!(int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0)) {
425 2286 : uniq_index.resize(unique.size());
426 32931 : for(unsigned i=0; i<unique.size(); i++) {
427 30645 : uniq_index[i]=g2l[unique[i].index()];
428 : }
429 13716 : for(const auto & ip : inputs) {
430 11430 : if( (!ip->fixed || firststep) && ip->wasset ) {
431 7060 : (ip->mydata)->share_data( unique, uniq_index, ip->copyOutput(0) );
432 7060 : values_to_get.push_back(ip->copyOutput(0));
433 7060 : ndata++;
434 : }
435 : }
436 136348 : } else if( unique_serial) {
437 21607 : uniq_index.resize(unique.size());
438 803079 : for(unsigned i=0; i<unique.size(); i++) {
439 781472 : uniq_index[i]=unique[i].index();
440 : }
441 129642 : for(const auto & ip : inputs) {
442 108035 : if( (!ip->fixed || firststep) && ip->wasset ) {
443 64821 : (ip->mydata)->share_data( unique, uniq_index, ip->copyOutput(0) );
444 64821 : values_to_get.push_back(ip->copyOutput(0));
445 64821 : ndata++;
446 : }
447 : }
448 : } else {
449 : // faster version, which retrieves all atoms
450 688366 : for(const auto & ip : inputs) {
451 573641 : if( (!ip->fixed || firststep) && ip->wasset ) {
452 346351 : (ip->mydata)->share_data( 0, getNumberOfAtoms(), ip->copyOutput(0) );
453 346335 : values_to_get.push_back(ip->copyOutput(0));
454 346335 : ndata++;
455 : }
456 : }
457 : }
458 :
459 138618 : if(dd && shuffledAtoms>0) {
460 2260 : if(dd.async) {
461 9720 : for(unsigned i=0; i<dd.mpi_request_positions.size(); i++) {
462 7498 : dd.mpi_request_positions[i].wait();
463 : }
464 9720 : for(unsigned i=0; i<dd.mpi_request_index.size(); i++) {
465 7498 : dd.mpi_request_index[i].wait();
466 : }
467 : }
468 :
469 2260 : int count=0;
470 32827 : for(const auto & p : unique) {
471 30567 : dd.indexToBeSent[count]=p.index();
472 : int dpoint=0;
473 127350 : for(unsigned i=0; i<values_to_get.size(); ++i) {
474 96783 : dd.positionsToBeSent[ndata*count+dpoint]=values_to_get[i]->get(p.index());
475 96783 : dpoint++;
476 : }
477 30567 : count++;
478 : }
479 :
480 2260 : if(dd.async) {
481 2222 : asyncSent=true;
482 2222 : dd.mpi_request_positions.resize(dd.Get_size());
483 2222 : dd.mpi_request_index.resize(dd.Get_size());
484 9964 : for(int i=0; i<dd.Get_size(); i++) {
485 7742 : dd.mpi_request_index[i]=dd.Isend(&dd.indexToBeSent[0],count,i,666);
486 7742 : dd.mpi_request_positions[i]=dd.Isend(&dd.positionsToBeSent[0],ndata*count,i,667);
487 : }
488 : } else {
489 38 : const int n=(dd.Get_size());
490 54 : std::vector<int> counts(n);
491 38 : std::vector<int> displ(n);
492 38 : std::vector<int> counts5(n);
493 38 : std::vector<int> displ5(n);
494 38 : dd.Allgather(count,counts);
495 38 : displ[0]=0;
496 116 : for(int i=1; i<n; ++i) {
497 78 : displ[i]=displ[i-1]+counts[i-1];
498 : }
499 154 : for(int i=0; i<n; ++i) {
500 116 : counts5[i]=counts[i]*ndata;
501 : }
502 154 : for(int i=0; i<n; ++i) {
503 116 : displ5[i]=displ[i]*ndata;
504 : }
505 38 : dd.Allgatherv(&dd.indexToBeSent[0],count,&dd.indexToBeReceived[0],&counts[0],&displ[0]);
506 38 : dd.Allgatherv(&dd.positionsToBeSent[0],ndata*count,&dd.positionsToBeReceived[0],&counts5[0],&displ5[0]);
507 38 : int tot=displ[n-1]+counts[n-1];
508 1674 : for(int i=0; i<tot; i++) {
509 : int dpoint=0;
510 7440 : for(unsigned j=0; j<values_to_get.size(); ++j) {
511 5804 : values_to_get[j]->data[dd.indexToBeReceived[i]] = dd.positionsToBeReceived[ndata*i+dpoint];
512 5804 : dpoint++;
513 : }
514 : }
515 : }
516 : }
517 138618 : }
518 :
519 138614 : void DomainDecomposition::wait() {
520 831684 : for(const auto & ip : inputs) {
521 693070 : ip->dataCanBeSet=false;
522 : }
523 :
524 138614 : if(dd && shuffledAtoms>0) {
525 : int ndata=0;
526 : std::vector<Value*> values_to_set;
527 13560 : for(const auto & ip : inputs) {
528 11300 : if( (!ip->fixed || firststep) && ip->wasset ) {
529 6980 : values_to_set.push_back(ip->copyOutput(0));
530 6980 : ndata++;
531 : }
532 : }
533 :
534 : // receive toBeReceived
535 2260 : if(asyncSent) {
536 : Communicator::Status status;
537 : std::size_t count=0;
538 9964 : for(int i=0; i<dd.Get_size(); i++) {
539 7742 : plumed_assert(count<=dd.positionsToBeReceived.size()/ndata);
540 7742 : const std::size_t remaining=dd.indexToBeReceived.size()-count;
541 7742 : const std::size_t remainingData=dd.positionsToBeReceived.size()-ndata*count;
542 : // Empty ranks still send messages after an earlier rank fills the buffer.
543 : // Avoid indexing one past the end, including for zero-sized storage.
544 7742 : int* indices=remaining ? dd.indexToBeReceived.data()+count : nullptr;
545 7742 : double* positionData=remainingData ? dd.positionsToBeReceived.data()+ndata*count : nullptr;
546 7742 : dd.Recv(indices,remaining,i,666,status);
547 7742 : const int c=status.Get_count<int>();
548 7742 : plumed_assert(c>=0 && static_cast<std::size_t>(c)<=remaining);
549 7742 : dd.Recv(positionData,remainingData,i,667);
550 7742 : count+=c;
551 : }
552 68502 : for(int i=0; i<count; i++) {
553 : int dpoint=0;
554 277236 : for(unsigned j=0; j<values_to_set.size(); ++j) {
555 210956 : values_to_set[j]->set(dd.indexToBeReceived[i], dd.positionsToBeReceived[ndata*i+dpoint] );
556 210956 : dpoint++;
557 : }
558 : }
559 2222 : asyncSent=false;
560 : }
561 : }
562 138614 : }
563 :
564 150 : unsigned DomainDecomposition::getNumberOfForcesToRescale() const {
565 150 : return gatindex.size();
566 : }
567 :
568 138474 : void DomainDecomposition::apply() {
569 138474 : std::vector<unsigned> forced_uniq_index(forced_unique.size());
570 138474 : if(!(int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0)) {
571 24361 : for(unsigned i=0; i<forced_unique.size(); i++) {
572 22193 : forced_uniq_index[i]=g2l[forced_unique[i].index()];
573 : }
574 : } else {
575 4474793 : for(unsigned i=0; i<forced_unique.size(); i++) {
576 4338487 : forced_uniq_index[i]=forced_unique[i].index();
577 : }
578 : }
579 830834 : for(const auto & ip : inputs) {
580 692362 : if( !(ip->getPntrToValue())->forcesWereAdded() || ip->noforce ) {
581 347513 : continue;
582 344849 : } else if( ip->wasscaled || (!unique_serial && int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0) ) {
583 281234 : (ip->mydata)->add_force( gatindex, ip->getPntrToValue() );
584 : } else {
585 63615 : (ip->mydata)->add_force( forced_unique, forced_uniq_index, ip->getPntrToValue() );
586 : }
587 : }
588 138472 : }
589 :
590 138472 : void DomainDecomposition::reset() {
591 138472 : if( !unique_serial && int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0 ) {
592 : return;
593 : }
594 : // This is an optimisation to ensure that we don't call std::fill over the whole forces
595 : // array if there are a small number of atoms passed between the MD code and PLUMED
596 23765 : if( dd && shuffledAtoms>0 ) {
597 2142 : getAllActiveAtoms( unique );
598 : }
599 142590 : for(const auto & ip : inputs) {
600 118825 : (ip->copyOutput(0))->clearInputForce( unique );
601 : }
602 : }
603 :
604 114 : void DomainDecomposition::writeBinary(std::ostream&o) {
605 684 : for(const auto & ip : inputs) {
606 570 : ip->writeBinary(o);
607 : }
608 114 : }
609 :
610 114 : void DomainDecomposition::readBinary(std::istream&i) {
611 684 : for(const auto & ip : inputs) {
612 570 : ip->readBinary(i);
613 : }
614 114 : }
615 :
616 134809 : void DomainDecomposition::broadcastToDomains( Value* val ) {
617 134809 : if( dd ) {
618 21671 : dd.Bcast( val->data, 0 );
619 : }
620 134809 : }
621 :
622 3989 : void DomainDecomposition::sumOverDomains( Value* val ) {
623 3989 : if( dd && shuffledAtoms>0 ) {
624 0 : dd.Sum( val->data );
625 : }
626 3989 : }
627 :
628 900 : const long int& DomainDecomposition::getDdStep() const {
629 900 : return ddStep;
630 : }
631 :
632 459 : const std::vector<int>& DomainDecomposition::getGatindex() const {
633 459 : return gatindex;
634 : }
635 :
636 2251 : void DomainDecomposition::getAllActiveAtoms( std::vector<AtomNumber>& u ) {
637 : gch::small_vector<const std::vector<AtomNumber>*,32> vectors;
638 : vectors.reserve(actions.size());
639 21256 : for(unsigned i=0; i<actions.size(); i++) {
640 19005 : if(actions[i]->isActive()) {
641 13133 : if(!actions[i]->getUnique().empty()) {
642 : // unique are the local atoms
643 23048 : vectors.push_back(&actions[i]->getUnique());
644 : }
645 : }
646 : }
647 : u.clear();
648 2251 : mergeVectorTools::mergeSortedVectors(vectors,u);
649 2251 : }
650 :
651 116 : void DomainDecomposition::createFullList(const TypesafePtr & n) {
652 116 : if( firststep ) {
653 7 : int natoms = getNumberOfAtoms();
654 7 : n.set(int(natoms));
655 7 : fullList.resize(natoms);
656 803 : for(unsigned i=0; i<natoms; i++) {
657 796 : fullList[i]=i;
658 : }
659 : } else {
660 : // We update here the unique list defined at Atoms::unique.
661 : // This is not very clear, and probably should be coded differently.
662 : // Hopefully this fix the longstanding issue with NAMD.
663 109 : getAllActiveAtoms( unique );
664 109 : fullList.clear();
665 109 : fullList.reserve(unique.size());
666 5012 : for(const auto & p : unique) {
667 4903 : fullList.push_back(p.index());
668 : }
669 109 : n.set(int(fullList.size()));
670 : }
671 116 : }
672 :
673 116 : void DomainDecomposition::getFullList(const TypesafePtr & x) {
674 : auto xx=x.template get<const int**>();
675 116 : if(!fullList.empty()) {
676 110 : *xx=&fullList[0];
677 : } else {
678 6 : *xx=NULL;
679 : }
680 116 : }
681 :
682 116 : void DomainDecomposition::clearFullList() {
683 116 : fullList.resize(0);
684 116 : }
685 :
686 : }
|
