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 756 : Option interpretEnvString(const char* env,const char* str) {
48 756 : 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 72010 : Option getenvForceUnique() {
73 : static const char* name="PLUMED_FORCE_UNIQUE";
74 72010 : static const auto opt = interpretEnvString(name,std::getenv(name));
75 72010 : return opt;
76 : }
77 :
78 : }
79 :
80 : PLUMED_REGISTER_ACTION(DomainDecomposition,"DOMAIN_DECOMPOSITION")
81 :
82 402 : void DomainDecomposition::DomainComms::enable(Communicator& c) {
83 402 : on=true;
84 402 : Set_comm(c.Get_comm());
85 402 : async=Get_size()<10;
86 402 : if(std::getenv("PLUMED_ASYNC_SHARE")) {
87 4 : std::string s(std::getenv("PLUMED_ASYNC_SHARE"));
88 4 : if(s=="yes") {
89 0 : async=true;
90 4 : } else if(s=="no") {
91 4 : async=false;
92 : } else {
93 0 : plumed_merror("PLUMED_ASYNC_SHARE variable is set to " + s + "; should be yes or no");
94 : }
95 : }
96 402 : }
97 :
98 1226 : void DomainDecomposition::registerKeywords(Keywords& keys) {
99 1226 : ActionForInterface::registerKeywords( keys );
100 1226 : keys.remove("ROLE");
101 2452 : keys.add("compulsory","NATOMS","the number of atoms across all the domains");
102 2452 : keys.add("numbered","VALUE","value to create for the domain decomposition");
103 2452 : keys.add("numbered","UNIT","unit of the ith value that is being passed through the domain decomposition");
104 2452 : keys.add("numbered","CONSTANT","is the ith value that is being passed through the domain decomposition constant");
105 2452 : 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 2452 : 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 2452 : keys.add("compulsory","PBCLABEL","Box","the label to use for the PBC action that will be created");
110 1226 : keys.setValueDescription("the domain that each atom is within");
111 1226 : }
112 :
113 1222 : DomainDecomposition::DomainDecomposition(const ActionOptions&ao):
114 : Action(ao),
115 : ActionForInterface(ao),
116 1222 : ddStep(0),
117 1222 : shuffledAtoms(0),
118 1222 : asyncSent(false),
119 1222 : unique_serial(false) {
120 : // Read in the number of atoms
121 : int natoms;
122 2444 : parse("NATOMS",natoms);
123 : std::string str_natoms;
124 1222 : Tools::convert( natoms, str_natoms );
125 : // Setup the gat index
126 1222 : gatindex.resize(natoms);
127 9202518 : for(unsigned i=0; i<gatindex.size(); i++) {
128 9201296 : gatindex[i]=i;
129 : }
130 : // Now read in the values we are creating here
131 6110 : for(unsigned i=1;; ++i) {
132 : std::string valname;
133 14664 : if( !parseNumbered("VALUE",i,valname) ) {
134 : break;
135 : }
136 : std::string unit;
137 12220 : parseNumbered("UNIT",i,unit);
138 : std::string constant;
139 12220 : parseNumbered("CONSTANT",i,constant);
140 : std::string period;
141 12220 : parseNumbered("PERIODIC",i,period);
142 : std::string role;
143 12220 : parseNumbered("ROLE",i,role);
144 6110 : if( constant=="True") {
145 4888 : plumed.readInputLine( valname + ": PUT FROM_DOMAINS CONSTANT SHAPE=" + str_natoms + " UNIT=" + unit + " PERIODIC=" + period + " ROLE=" + role, true );
146 3666 : } else if( constant=="False") {
147 7332 : 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 6110 : ActionToPutData* ap=plumed.getActionSet().selectWithLabel<ActionToPutData*>(valname);
153 6110 : ap->addDependency( this );
154 6110 : inputs.push_back( ap );
155 6110 : }
156 : std::string pbclabel;
157 1222 : parse("PBCLABEL",pbclabel);
158 1222 : plumed.readInputLine(pbclabel + ": PBC",true);
159 : // Turn on the domain decomposition
160 1222 : if( Communicator::initialized() ) {
161 401 : Set_comm(comm);
162 : }
163 1222 : }
164 :
165 402 : void DomainDecomposition::Set_comm(Communicator& comm) {
166 402 : dd.enable(comm);
167 402 : setAtomsNlocal(getNumberOfAtoms());
168 402 : setAtomsContiguous(0);
169 402 : if( dd.Get_rank()!=0 ) {
170 145 : ActionToPutData* ap=plumed.getActionSet().selectWithLabel<ActionToPutData*>("Box");
171 145 : ap->noforce=true;
172 : }
173 402 : }
174 :
175 559880 : unsigned DomainDecomposition::getNumberOfAtoms() const {
176 559880 : if( inputs.size()==0 ) {
177 : return 0;
178 : }
179 559880 : return (inputs[0]->getPntrToValue())->getShape()[0];
180 : }
181 :
182 74756 : void DomainDecomposition::resetForStepStart() {
183 448536 : for(const auto & pp : inputs) {
184 373780 : pp->resetForStepStart();
185 : }
186 74756 : }
187 :
188 365729 : void DomainDecomposition::setStart( const std::string& name, const unsigned& sss) {
189 1080916 : for(const auto & pp : inputs) {
190 1080916 : if( pp->getLabel()==name ) {
191 365729 : pp->setStart(name, sss);
192 365729 : return;
193 : }
194 : }
195 0 : plumed_error();
196 : }
197 :
198 365729 : void DomainDecomposition::setStride( const std::string& name, const unsigned& sss) {
199 1080916 : for(const auto & pp : inputs) {
200 1080916 : if( pp->getLabel()==name ) {
201 365729 : pp->setStride(name, sss);
202 365729 : return;
203 : }
204 : }
205 0 : plumed_error();
206 : }
207 :
208 369751 : bool DomainDecomposition::setValuePointer( const std::string& name, const TypesafePtr & val ) {
209 369751 : wasset=true; // Once the domain decomposition stuff is transferred moved the setting of this to where the g2l vector is setup
210 1105042 : for(const auto & pp : inputs) {
211 1101023 : if( pp->setValuePointer( name, val ) ) {
212 : return true;
213 : }
214 : }
215 : return false;
216 : }
217 :
218 224262 : bool DomainDecomposition::setForcePointer( const std::string& name, const TypesafePtr & val ) {
219 448644 : for(const auto & pp : inputs) {
220 448614 : if( pp->setForcePointer( name, val ) ) {
221 : return true;
222 : }
223 : }
224 : return false;
225 : }
226 :
227 1544 : void DomainDecomposition::setAtomsNlocal(int n) {
228 1544 : gatindex.resize(n);
229 1544 : g2l.resize(getNumberOfAtoms(),-1);
230 1544 : if(dd) {
231 : // Since these vectors are sent with MPI by using e.g.
232 : // &dd.positionsToBeSent[0]
233 : // we make sure they are non-zero-sized so as to
234 : // avoid errors when doing boundary check
235 1518 : if(n==0) {
236 8 : n++;
237 : }
238 1518 : std::size_t nvals = inputs.size(), natoms = getNumberOfAtoms();
239 1518 : dd.positionsToBeSent.resize(n*nvals,0.0);
240 1518 : dd.positionsToBeReceived.resize(natoms*nvals,0.0);
241 1518 : dd.indexToBeSent.resize(n,0);
242 1518 : dd.indexToBeReceived.resize(natoms,0);
243 : }
244 1544 : }
245 :
246 990 : void DomainDecomposition::setAtomsGatindex(const TypesafePtr & g,bool fortran) {
247 990 : plumed_massert( g || gatindex.size()==0, "NULL gatindex pointer with non-zero local atoms");
248 990 : auto gg=g.get<const int*>({gatindex.size()});
249 990 : ddStep=getStep();
250 990 : if(fortran) {
251 22 : for(unsigned i=0; i<gatindex.size(); i++) {
252 20 : gatindex[i]=gg[i]-1;
253 : }
254 : } else {
255 22140 : for(unsigned i=0; i<gatindex.size(); i++) {
256 21152 : gatindex[i]=gg[i];
257 : }
258 : }
259 57552 : for(unsigned i=0; i<g2l.size(); i++) {
260 56562 : g2l[i]=-1;
261 : }
262 990 : if( gatindex.size()==getNumberOfAtoms() ) {
263 11 : shuffledAtoms=0;
264 1007 : for(unsigned i=0; i<gatindex.size(); i++) {
265 998 : if( gatindex[i]!=i ) {
266 2 : shuffledAtoms=1;
267 2 : break;
268 : }
269 : }
270 : } else {
271 979 : shuffledAtoms=1;
272 : }
273 990 : if(dd) {
274 964 : dd.Sum(shuffledAtoms);
275 : }
276 22162 : for(unsigned i=0; i<gatindex.size(); i++) {
277 21172 : g2l[gatindex[i]]=i;
278 : }
279 : // keep in unique only those atoms that are local
280 9908 : for(unsigned i=0; i<actions.size(); i++) {
281 8918 : actions[i]->unique_local_needs_update=true;
282 : }
283 : unique.clear();
284 : forced_unique.clear();
285 990 : }
286 :
287 554 : void DomainDecomposition::setAtomsContiguous(int start) {
288 554 : ddStep=plumed.getStep();
289 218327 : for(unsigned i=0; i<gatindex.size(); i++) {
290 217773 : gatindex[i]=start+i;
291 : }
292 230639 : for(unsigned i=0; i<g2l.size(); i++) {
293 230085 : g2l[i]=-1;
294 : }
295 218327 : for(unsigned i=0; i<gatindex.size(); i++) {
296 217773 : g2l[gatindex[i]]=i;
297 : }
298 554 : if(gatindex.size()<getNumberOfAtoms()) {
299 152 : shuffledAtoms=1;
300 : }
301 : // keep in unique only those atoms that are local
302 766 : for(unsigned i=0; i<actions.size(); i++) {
303 212 : actions[i]->unique_local_needs_update=true;
304 : }
305 : unique.clear();
306 : forced_unique.clear();
307 554 : }
308 :
309 1248 : void DomainDecomposition::shareAll() {
310 1248 : unique.clear();
311 1248 : forced_unique.clear();
312 1248 : int natoms = getNumberOfAtoms();
313 1248 : if( dd && shuffledAtoms>0 ) {
314 17616 : for(int i=0; i<natoms; ++i)
315 17430 : if( g2l[i]>=0 ) {
316 8003 : unique.push_back( AtomNumber::index(i) );
317 : }
318 : } else {
319 1062 : unique.resize(natoms);
320 5144483 : for(int i=0; i<natoms; i++) {
321 5143421 : unique[i]=AtomNumber::index(i);
322 : }
323 : }
324 1248 : forced_unique.resize( unique.size() );
325 5152672 : for(unsigned i=0; i<unique.size(); ++i) {
326 5151424 : forced_unique[i] = unique[i];
327 : }
328 1248 : share(unique);
329 1232 : }
330 :
331 73144 : void DomainDecomposition::share() {
332 : // We can no longer set the pointers after the share
333 : bool atomsNeeded=false;
334 438864 : for(const auto & pp : inputs) {
335 365720 : pp->share();
336 : }
337 : // At first step I scatter all the atoms so as to store their mass and charge
338 : // Notice that this works with the assumption that charges and masses are
339 : // not changing during the simulation!
340 73144 : if( firststep ) {
341 1134 : actions = plumed.getActionSet().select<ActionAtomistic*>();
342 1134 : shareAll();
343 1118 : return;
344 : }
345 :
346 72010 : if(getenvForceUnique()==Option::automatic) {
347 : unsigned largest=0;
348 728346 : for(unsigned i=0; i<actions.size(); i++) {
349 656336 : if(actions[i]->isActive()) {
350 : auto l=actions[i]->getUnique().size();
351 451068 : if(l>largest) {
352 76455 : largest=l;
353 : }
354 : }
355 : }
356 72010 : if(largest*2<getNumberOfAtoms()) {
357 23647 : unique_serial=true;
358 : } else {
359 48363 : unique_serial=false;
360 : }
361 0 : } else if(getenvForceUnique()==Option::yes) {
362 0 : unique_serial=true;
363 0 : } else if(getenvForceUnique()==Option::no) {
364 0 : unique_serial=false;
365 : } else {
366 0 : plumed_error();
367 : }
368 :
369 72010 : if(unique_serial || !(int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0)) {
370 289091 : for(unsigned i=0; i<actions.size(); i++) {
371 265056 : if( actions[i]->unique_local_needs_update ) {
372 265786 : actions[i]->updateUniqueLocal( !(dd && shuffledAtoms>0), g2l );
373 : }
374 : }
375 : // Now reset unique for the new step
376 : gch::small_vector<const std::vector<AtomNumber>*,32> forced_vectors;
377 : gch::small_vector<const std::vector<AtomNumber>*,32> nonforced_vectors;
378 : forced_vectors.reserve(actions.size());
379 : nonforced_vectors.reserve(actions.size());
380 289091 : for(unsigned i=0; i<actions.size(); i++) {
381 265056 : if(actions[i]->isActive()) {
382 192031 : if(!actions[i]->getUnique().empty()) {
383 : // unique are the local atoms
384 106738 : if( actions[i]->actionHasForces() ) {
385 191054 : forced_vectors.push_back(&actions[i]->getUniqueLocal());
386 : } else {
387 22422 : nonforced_vectors.push_back(&actions[i]->getUniqueLocal());
388 : }
389 : }
390 : }
391 : }
392 24035 : if( !(forced_vectors.empty() && nonforced_vectors.empty()) ) {
393 : atomsNeeded=true;
394 : }
395 : // Merge the atoms from the atoms that have a force on
396 24035 : unique.clear();
397 24035 : forced_unique.clear();
398 24035 : mergeVectorTools::mergeSortedVectors(forced_vectors,forced_unique);
399 : // Merge all the atoms
400 24035 : nonforced_vectors.push_back( &forced_unique );
401 24035 : mergeVectorTools::mergeSortedVectors(nonforced_vectors,unique);
402 : } else {
403 439255 : for(unsigned i=0; i<actions.size(); i++) {
404 391280 : if(actions[i]->isActive()) {
405 259037 : if(!actions[i]->getUnique().empty()) {
406 : atomsNeeded=true;
407 : }
408 : }
409 : }
410 : }
411 :
412 : // Now we retrieve the atom numbers we need
413 72010 : if( atomsNeeded ) {
414 71305 : share( unique );
415 : }
416 : }
417 :
418 72553 : void DomainDecomposition::share(const std::vector<AtomNumber>& unique) {
419 : // This retrieves what values we need to get
420 : int ndata=0;
421 : std::vector<Value*> values_to_get;
422 72553 : if(!(int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0)) {
423 2214 : uniq_index.resize(unique.size());
424 32727 : for(unsigned i=0; i<unique.size(); i++) {
425 30513 : uniq_index[i]=g2l[unique[i].index()];
426 : }
427 13284 : for(const auto & ip : inputs) {
428 11070 : if( (!ip->fixed || firststep) && ip->wasset ) {
429 6788 : (ip->mydata)->share_data( unique, uniq_index, ip->copyOutput(0) );
430 6788 : values_to_get.push_back(ip->copyOutput(0));
431 6788 : ndata++;
432 : }
433 : }
434 70339 : } else if( unique_serial) {
435 21304 : uniq_index.resize(unique.size());
436 801395 : for(unsigned i=0; i<unique.size(); i++) {
437 780091 : uniq_index[i]=unique[i].index();
438 : }
439 127824 : for(const auto & ip : inputs) {
440 106520 : if( (!ip->fixed || firststep) && ip->wasset ) {
441 63912 : (ip->mydata)->share_data( unique, uniq_index, ip->copyOutput(0) );
442 63912 : values_to_get.push_back(ip->copyOutput(0));
443 63912 : ndata++;
444 : }
445 : }
446 : } else {
447 : // faster version, which retrieves all atoms
448 294130 : for(const auto & ip : inputs) {
449 245111 : if( (!ip->fixed || firststep) && ip->wasset ) {
450 148945 : (ip->mydata)->share_data( 0, getNumberOfAtoms(), ip->copyOutput(0) );
451 148929 : values_to_get.push_back(ip->copyOutput(0));
452 148929 : ndata++;
453 : }
454 : }
455 : }
456 :
457 72537 : if(dd && shuffledAtoms>0) {
458 2188 : if(dd.async) {
459 9598 : for(unsigned i=0; i<dd.mpi_request_positions.size(); i++) {
460 7430 : dd.mpi_request_positions[i].wait();
461 : }
462 9598 : for(unsigned i=0; i<dd.mpi_request_index.size(); i++) {
463 7430 : dd.mpi_request_index[i].wait();
464 : }
465 : }
466 :
467 2188 : int count=0;
468 32623 : for(const auto & p : unique) {
469 30435 : dd.indexToBeSent[count]=p.index();
470 : int dpoint=0;
471 126694 : for(unsigned i=0; i<values_to_get.size(); ++i) {
472 96259 : dd.positionsToBeSent[ndata*count+dpoint]=values_to_get[i]->get(p.index());
473 96259 : dpoint++;
474 : }
475 30435 : count++;
476 : }
477 :
478 2188 : if(dd.async) {
479 2168 : asyncSent=true;
480 2168 : dd.mpi_request_positions.resize(dd.Get_size());
481 2168 : dd.mpi_request_index.resize(dd.Get_size());
482 9802 : for(int i=0; i<dd.Get_size(); i++) {
483 7634 : dd.mpi_request_index[i]=dd.Isend(&dd.indexToBeSent[0],count,i,666);
484 7634 : dd.mpi_request_positions[i]=dd.Isend(&dd.positionsToBeSent[0],ndata*count,i,667);
485 : }
486 : } else {
487 20 : const int n=(dd.Get_size());
488 36 : std::vector<int> counts(n);
489 20 : std::vector<int> displ(n);
490 20 : std::vector<int> counts5(n);
491 20 : std::vector<int> displ5(n);
492 20 : dd.Allgather(count,counts);
493 20 : displ[0]=0;
494 80 : for(int i=1; i<n; ++i) {
495 60 : displ[i]=displ[i-1]+counts[i-1];
496 : }
497 100 : for(int i=0; i<n; ++i) {
498 80 : counts5[i]=counts[i]*ndata;
499 : }
500 100 : for(int i=0; i<n; ++i) {
501 80 : displ5[i]=displ[i]*ndata;
502 : }
503 20 : dd.Allgatherv(&dd.indexToBeSent[0],count,&dd.indexToBeReceived[0],&counts[0],&displ[0]);
504 20 : dd.Allgatherv(&dd.positionsToBeSent[0],ndata*count,&dd.positionsToBeReceived[0],&counts5[0],&displ5[0]);
505 20 : int tot=displ[n-1]+counts[n-1];
506 1620 : for(int i=0; i<tot; i++) {
507 : int dpoint=0;
508 7264 : for(unsigned j=0; j<values_to_get.size(); ++j) {
509 5664 : values_to_get[j]->data[dd.indexToBeReceived[i]] = dd.positionsToBeReceived[ndata*i+dpoint];
510 5664 : dpoint++;
511 : }
512 : }
513 : }
514 : }
515 72537 : }
516 :
517 72529 : void DomainDecomposition::wait() {
518 435174 : for(const auto & ip : inputs) {
519 362645 : ip->dataCanBeSet=false;
520 : }
521 :
522 72529 : if(dd && shuffledAtoms>0) {
523 : int ndata=0;
524 : std::vector<Value*> values_to_set;
525 13128 : for(const auto & ip : inputs) {
526 10940 : if( (!ip->fixed || firststep) && ip->wasset ) {
527 6708 : values_to_set.push_back(ip->copyOutput(0));
528 6708 : ndata++;
529 : }
530 : }
531 :
532 : // receive toBeReceived
533 2188 : if(asyncSent) {
534 : Communicator::Status status;
535 : std::size_t count=0;
536 9802 : for(int i=0; i<dd.Get_size(); i++) {
537 7634 : dd.Recv(&dd.indexToBeReceived[count],dd.indexToBeReceived.size()-count,i,666,status);
538 7634 : int c=status.Get_count<int>();
539 7634 : dd.Recv(&dd.positionsToBeReceived[ndata*count],dd.positionsToBeReceived.size()-ndata*count,i,667);
540 7634 : count+=c;
541 : }
542 68220 : for(int i=0; i<count; i++) {
543 : int dpoint=0;
544 276100 : for(unsigned j=0; j<values_to_set.size(); ++j) {
545 210048 : values_to_set[j]->set(dd.indexToBeReceived[i], dd.positionsToBeReceived[ndata*i+dpoint] );
546 210048 : dpoint++;
547 : }
548 : }
549 2168 : asyncSent=false;
550 : }
551 : }
552 72529 : }
553 :
554 0 : unsigned DomainDecomposition::getNumberOfForcesToRescale() const {
555 0 : return gatindex.size();
556 : }
557 :
558 72405 : void DomainDecomposition::apply() {
559 72405 : std::vector<unsigned> forced_uniq_index(forced_unique.size());
560 72405 : if(!(int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0)) {
561 24181 : for(unsigned i=0; i<forced_unique.size(); i++) {
562 22081 : forced_uniq_index[i]=g2l[forced_unique[i].index()];
563 : }
564 : } else {
565 4189154 : for(unsigned i=0; i<forced_unique.size(); i++) {
566 4118849 : forced_uniq_index[i]=forced_unique[i].index();
567 : }
568 : }
569 434420 : for(const auto & ip : inputs) {
570 362017 : if( !(ip->getPntrToValue())->forcesWereAdded() || ip->noforce ) {
571 213644 : continue;
572 148373 : } else if( ip->wasscaled || (!unique_serial && int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0) ) {
573 85784 : (ip->mydata)->add_force( gatindex, ip->getPntrToValue() );
574 : } else {
575 62589 : (ip->mydata)->add_force( forced_unique, forced_uniq_index, ip->getPntrToValue() );
576 : }
577 : }
578 72403 : }
579 :
580 72403 : void DomainDecomposition::reset() {
581 72403 : if( !unique_serial && int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0 ) {
582 : return;
583 : }
584 : // This is an optimisation to ensure that we don't call std::fill over the whole forces
585 : // array if there are a small number of atoms passed between the MD code and PLUMED
586 23394 : if( dd && shuffledAtoms>0 ) {
587 2074 : getAllActiveAtoms( unique );
588 : }
589 140364 : for(const auto & ip : inputs) {
590 116970 : (ip->copyOutput(0))->clearInputForce( unique );
591 : }
592 : }
593 :
594 114 : void DomainDecomposition::writeBinary(std::ostream&o) {
595 684 : for(const auto & ip : inputs) {
596 570 : ip->writeBinary(o);
597 : }
598 114 : }
599 :
600 114 : void DomainDecomposition::readBinary(std::istream&i) {
601 684 : for(const auto & ip : inputs) {
602 570 : ip->readBinary(i);
603 : }
604 114 : }
605 :
606 68721 : void DomainDecomposition::broadcastToDomains( Value* val ) {
607 68721 : if( dd ) {
608 20496 : dd.Bcast( val->data, 0 );
609 : }
610 68721 : }
611 :
612 3989 : void DomainDecomposition::sumOverDomains( Value* val ) {
613 3989 : if( dd && shuffledAtoms>0 ) {
614 0 : dd.Sum( val->data );
615 : }
616 3989 : }
617 :
618 900 : const long int& DomainDecomposition::getDdStep() const {
619 900 : return ddStep;
620 : }
621 :
622 459 : const std::vector<int>& DomainDecomposition::getGatindex() const {
623 459 : return gatindex;
624 : }
625 :
626 2183 : void DomainDecomposition::getAllActiveAtoms( std::vector<AtomNumber>& u ) {
627 : gch::small_vector<const std::vector<AtomNumber>*,32> vectors;
628 : vectors.reserve(actions.size());
629 21016 : for(unsigned i=0; i<actions.size(); i++) {
630 18833 : if(actions[i]->isActive()) {
631 13033 : if(!actions[i]->getUnique().empty()) {
632 : // unique are the local atoms
633 22912 : vectors.push_back(&actions[i]->getUnique());
634 : }
635 : }
636 : }
637 : u.clear();
638 2183 : mergeVectorTools::mergeSortedVectors(vectors,u);
639 2183 : }
640 :
641 116 : void DomainDecomposition::createFullList(const TypesafePtr & n) {
642 116 : if( firststep ) {
643 7 : int natoms = getNumberOfAtoms();
644 7 : n.set(int(natoms));
645 7 : fullList.resize(natoms);
646 803 : for(unsigned i=0; i<natoms; i++) {
647 796 : fullList[i]=i;
648 : }
649 : } else {
650 : // We update here the unique list defined at Atoms::unique.
651 : // This is not very clear, and probably should be coded differently.
652 : // Hopefully this fix the longstanding issue with NAMD.
653 109 : getAllActiveAtoms( unique );
654 109 : fullList.clear();
655 109 : fullList.reserve(unique.size());
656 5012 : for(const auto & p : unique) {
657 4903 : fullList.push_back(p.index());
658 : }
659 109 : n.set(int(fullList.size()));
660 : }
661 116 : }
662 :
663 116 : void DomainDecomposition::getFullList(const TypesafePtr & x) {
664 : auto xx=x.template get<const int**>();
665 116 : if(!fullList.empty()) {
666 110 : *xx=&fullList[0];
667 : } else {
668 6 : *xx=NULL;
669 : }
670 116 : }
671 :
672 116 : void DomainDecomposition::clearFullList() {
673 116 : fullList.resize(0);
674 116 : }
675 :
676 : }
|
