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 734 : Option interpretEnvString(const char* env,const char* str) {
48 734 : if(!str) return Option::automatic;
49 0 : if(!std::strcmp(str,"yes"))return Option::yes;
50 0 : if(!std::strcmp(str,"no"))return Option::no;
51 0 : if(!std::strcmp(str,"auto"))return Option::automatic;
52 0 : plumed_error()<<"Cannot understand env var "<<env<<"\nPossible values: yes/no/auto\nActual value: "<<str;
53 : }
54 :
55 : /// Use unique list of atoms to manipulate forces and positions.
56 : /// A unique list of atoms is used to manipulate forces and positions in MPI parallel runs.
57 : /// In serial runs, this is done if convenient. The code currently contain
58 : /// some heuristic to decide if the unique list should be used or not.
59 : /// An env var can be used to override this decision.
60 : /// export PLUMED_FORCE_UNIQUE=yes # enforce using the unique list in serial runs
61 : /// export PLUMED_FORCE_UNIQUE=no # enforce not using the unique list in serial runs
62 : /// export PLUMED_FORCE_UNIQUE=auto # choose heuristically
63 : /// default: auto
64 69919 : Option getenvForceUnique() {
65 : static const char* name="PLUMED_FORCE_UNIQUE";
66 69919 : static const auto opt = interpretEnvString(name,std::getenv(name));
67 69919 : return opt;
68 : }
69 :
70 : }
71 :
72 : PLUMED_REGISTER_ACTION(DomainDecomposition,"DOMAIN_DECOMPOSITION")
73 :
74 393 : void DomainDecomposition::DomainComms::enable(Communicator& c) {
75 393 : on=true;
76 393 : Set_comm(c.Get_comm());
77 393 : async=Get_size()<10;
78 393 : if(std::getenv("PLUMED_ASYNC_SHARE")) {
79 4 : std::string s(std::getenv("PLUMED_ASYNC_SHARE"));
80 4 : if(s=="yes") async=true;
81 4 : else if(s=="no") async=false;
82 0 : else plumed_merror("PLUMED_ASYNC_SHARE variable is set to " + s + "; should be yes or no");
83 : }
84 393 : }
85 :
86 1197 : void DomainDecomposition::registerKeywords(Keywords& keys) {
87 1197 : ActionForInterface::registerKeywords( keys ); keys.remove("ROLE");
88 2394 : keys.add("compulsory","NATOMS","the number of atoms across all the domains");
89 2394 : keys.add("numbered","VALUE","value to create for the domain decomposition");
90 2394 : keys.add("numbered","UNIT","unit of the ith value that is being passed through the domain decomposition");
91 2394 : keys.add("numbered","CONSTANT","is the ith value that is being passed through the domain decomposition constant");
92 2394 : 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 "
93 : "is not periodic you must state this using PERIODIC=NO. Positions are passed with PERIODIC=NO even though special methods are used "
94 : "to deal with pbc");
95 2394 : 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");
96 2394 : keys.add("compulsory","PBCLABEL","Box","the label to use for the PBC action that will be created");
97 1197 : keys.setValueDescription("the domain that each atom is within");
98 1197 : }
99 :
100 1195 : DomainDecomposition::DomainDecomposition(const ActionOptions&ao):
101 : Action(ao),
102 : ActionForInterface(ao),
103 1195 : ddStep(0),
104 1195 : shuffledAtoms(0),
105 1195 : asyncSent(false),
106 1195 : unique_serial(false)
107 : {
108 : // Read in the number of atoms
109 2390 : int natoms; parse("NATOMS",natoms);
110 1195 : std::string str_natoms; Tools::convert( natoms, str_natoms );
111 : // Setup the gat index
112 9190732 : gatindex.resize(natoms); for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=i;
113 : // Now read in the values we are creating here
114 5975 : for(unsigned i=1;; ++i) {
115 : std::string valname;
116 14340 : if( !parseNumbered("VALUE",i,valname) ) break;
117 11950 : std::string unit; parseNumbered("UNIT",i,unit);
118 11950 : std::string constant; parseNumbered("CONSTANT",i,constant);
119 11950 : std::string period; parseNumbered("PERIODIC",i,period);
120 11950 : std::string role; parseNumbered("ROLE",i,role);
121 8365 : if( constant=="True") plumed.readInputLine( valname + ": PUT FROM_DOMAINS CONSTANT SHAPE=" + str_natoms + " UNIT=" + unit + " PERIODIC=" + period + " ROLE=" + role, true );
122 7170 : else if( constant=="False") plumed.readInputLine( valname + ": PUT FROM_DOMAINS SHAPE=" + str_natoms + " UNIT=" + unit + " PERIODIC=" + period + " ROLE=" + role, true );
123 0 : else plumed_merror("missing information on whether value is constant");
124 : // And save the list of values that are set from here
125 5975 : ActionToPutData* ap=plumed.getActionSet().selectWithLabel<ActionToPutData*>(valname); ap->addDependency( this ); inputs.push_back( ap );
126 5975 : }
127 2390 : std::string pbclabel; parse("PBCLABEL",pbclabel); plumed.readInputLine(pbclabel + ": PBC",true);
128 : // Turn on the domain decomposition
129 1195 : if( Communicator::initialized() ) Set_comm(comm);
130 1195 : }
131 :
132 393 : void DomainDecomposition::Set_comm(Communicator& comm) {
133 393 : dd.enable(comm); setAtomsNlocal(getNumberOfAtoms()); setAtomsContiguous(0);
134 393 : if( dd.Get_rank()!=0 ) {
135 284 : ActionToPutData* ap=plumed.getActionSet().selectWithLabel<ActionToPutData*>("Box"); ap->noforce=true;
136 : }
137 393 : }
138 :
139 536882 : unsigned DomainDecomposition::getNumberOfAtoms() const {
140 536882 : if( inputs.size()==0 ) return 0;
141 536882 : return (inputs[0]->getPntrToValue())->getShape()[0];
142 : }
143 :
144 72638 : void DomainDecomposition::resetForStepStart() {
145 435828 : for(const auto & pp : inputs) pp->resetForStepStart();
146 72638 : }
147 :
148 355141 : void DomainDecomposition::setStart( const std::string& name, const unsigned& sss) {
149 1049156 : for(const auto & pp : inputs) {
150 1049156 : if( pp->getLabel()==name ) { pp->setStart(name, sss); return; }
151 : }
152 0 : plumed_error();
153 : }
154 :
155 355141 : void DomainDecomposition::setStride( const std::string& name, const unsigned& sss) {
156 1049156 : for(const auto & pp : inputs) {
157 1049156 : if( pp->getLabel()==name ) { pp->setStride(name, sss); return; }
158 : }
159 0 : plumed_error();
160 : }
161 :
162 359163 : bool DomainDecomposition::setValuePointer( const std::string& name, const TypesafePtr & val ) {
163 359163 : wasset=true; // Once the domain decomposition stuff is transferred moved the setting of this to where the g2l vector is setup
164 1073282 : for(const auto & pp : inputs) {
165 1069263 : if( pp->setValuePointer( name, val ) ) return true;
166 : }
167 : return false;
168 : }
169 :
170 217908 : bool DomainDecomposition::setForcePointer( const std::string& name, const TypesafePtr & val ) {
171 435936 : for(const auto & pp : inputs) {
172 435906 : if( pp->setForcePointer( name, val ) ) return true;
173 : }
174 : return false;
175 : }
176 :
177 1533 : void DomainDecomposition::setAtomsNlocal(int n) {
178 1533 : gatindex.resize(n);
179 1533 : g2l.resize(getNumberOfAtoms(),-1);
180 1533 : if(dd) {
181 : // Since these vectors are sent with MPI by using e.g.
182 : // &dd.positionsToBeSent[0]
183 : // we make sure they are non-zero-sized so as to
184 : // avoid errors when doing boundary check
185 1509 : if(n==0) n++;
186 1509 : std::size_t nvals = inputs.size(), natoms = getNumberOfAtoms();
187 1509 : dd.positionsToBeSent.resize(n*nvals,0.0);
188 1509 : dd.positionsToBeReceived.resize(natoms*nvals,0.0);
189 1509 : dd.indexToBeSent.resize(n,0);
190 1509 : dd.indexToBeReceived.resize(natoms,0);
191 : }
192 1533 : }
193 :
194 988 : void DomainDecomposition::setAtomsGatindex(const TypesafePtr & g,bool fortran) {
195 988 : plumed_massert( g || gatindex.size()==0, "NULL gatindex pointer with non-zero local atoms");
196 988 : auto gg=g.get<const int*>({gatindex.size()});
197 988 : ddStep=getStep();
198 988 : if(fortran) {
199 0 : for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=gg[i]-1;
200 : } else {
201 22140 : for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=gg[i];
202 : }
203 57530 : for(unsigned i=0; i<g2l.size(); i++) g2l[i]=-1;
204 988 : if( gatindex.size()==getNumberOfAtoms() ) {
205 9 : shuffledAtoms=0;
206 1005 : for(unsigned i=0; i<gatindex.size(); i++) {
207 996 : if( gatindex[i]!=i ) { shuffledAtoms=1; break; }
208 : }
209 : } else {
210 979 : shuffledAtoms=1;
211 : }
212 988 : if(dd) dd.Sum(shuffledAtoms);
213 22140 : for(unsigned i=0; i<gatindex.size(); i++) g2l[gatindex[i]]=i;
214 : // keep in unique only those atoms that are local
215 9906 : for(unsigned i=0; i<actions.size(); i++) actions[i]->unique_local_needs_update=true;
216 : unique.clear(); forced_unique.clear();
217 988 : }
218 :
219 545 : void DomainDecomposition::setAtomsContiguous(int start) {
220 545 : ddStep=plumed.getStep();
221 211635 : for(unsigned i=0; i<gatindex.size(); i++) gatindex[i]=start+i;
222 223947 : for(unsigned i=0; i<g2l.size(); i++) g2l[i]=-1;
223 211635 : for(unsigned i=0; i<gatindex.size(); i++) g2l[gatindex[i]]=i;
224 545 : if(gatindex.size()<getNumberOfAtoms()) shuffledAtoms=1;
225 : // keep in unique only those atoms that are local
226 757 : for(unsigned i=0; i<actions.size(); i++) actions[i]->unique_local_needs_update=true;
227 : unique.clear(); forced_unique.clear();
228 545 : }
229 :
230 1221 : void DomainDecomposition::shareAll() {
231 1449 : unique.clear(); forced_unique.clear(); int natoms = getNumberOfAtoms();
232 1221 : if( dd && shuffledAtoms>0 ) {
233 17616 : for(int i=0; i<natoms; ++i) if( g2l[i]>=0 ) unique.push_back( AtomNumber::index(i) );
234 : } else {
235 1035 : unique.resize(natoms);
236 5132697 : for(int i=0; i<natoms; i++) unique[i]=AtomNumber::index(i);
237 : }
238 1221 : forced_unique.resize( unique.size() );
239 5140886 : for(unsigned i=0; i<unique.size(); ++i) forced_unique[i] = unique[i];
240 1221 : share(unique);
241 1205 : }
242 :
243 71026 : void DomainDecomposition::share() {
244 : // We can no longer set the pointers after the share
245 426156 : bool atomsNeeded=false; for(const auto & pp : inputs) pp->share();
246 : // At first step I scatter all the atoms so as to store their mass and charge
247 : // Notice that this works with the assumption that charges and masses are
248 : // not changing during the simulation!
249 71026 : if( firststep ) {
250 1107 : actions = plumed.getActionSet().select<ActionAtomistic*>();
251 1107 : shareAll(); return;
252 : }
253 :
254 69919 : if(getenvForceUnique()==Option::automatic) {
255 : unsigned largest=0;
256 721000 : for(unsigned i=0; i<actions.size(); i++) {
257 651081 : if(actions[i]->isActive()) {
258 : auto l=actions[i]->getUnique().size();
259 444659 : if(l>largest) largest=l;
260 : }
261 : }
262 69919 : if(largest*2<getNumberOfAtoms()) unique_serial=true;
263 46317 : else unique_serial=false;
264 0 : } else if(getenvForceUnique()==Option::yes) {
265 0 : unique_serial=true;
266 0 : } else if(getenvForceUnique()==Option::no) {
267 0 : unique_serial=false;
268 : } else {
269 0 : plumed_error();
270 : }
271 :
272 69919 : if(unique_serial || !(int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0)) {
273 288804 : for(unsigned i=0; i<actions.size(); i++) {
274 273968 : if( actions[i]->unique_local_needs_update ) actions[i]->updateUniqueLocal( !(dd && shuffledAtoms>0), g2l );
275 : }
276 : // Now reset unique for the new step
277 : gch::small_vector<const std::vector<AtomNumber>*,32> forced_vectors;
278 : gch::small_vector<const std::vector<AtomNumber>*,32> nonforced_vectors;
279 : forced_vectors.reserve(actions.size()); nonforced_vectors.reserve(actions.size());
280 288804 : for(unsigned i=0; i<actions.size(); i++) {
281 264814 : if(actions[i]->isActive()) {
282 191803 : if(!actions[i]->getUnique().empty()) {
283 : // unique are the local atoms
284 202027 : if( actions[i]->actionHasForces() ) forced_vectors.push_back(&actions[i]->getUniqueLocal());
285 22362 : else nonforced_vectors.push_back(&actions[i]->getUniqueLocal());
286 : }
287 : }
288 : }
289 23990 : if( !(forced_vectors.empty() && nonforced_vectors.empty()) ) atomsNeeded=true;
290 : // Merge the atoms from the atoms that have a force on
291 46231 : unique.clear(); forced_unique.clear();
292 23990 : mergeVectorTools::mergeSortedVectors(forced_vectors,forced_unique);
293 : // Merge all the atoms
294 23990 : nonforced_vectors.push_back( &forced_unique );
295 23990 : mergeVectorTools::mergeSortedVectors(nonforced_vectors,unique);
296 : } else {
297 432196 : for(unsigned i=0; i<actions.size(); i++) {
298 386267 : if(actions[i]->isActive()) {
299 252856 : if(!actions[i]->getUnique().empty()) { atomsNeeded=true; }
300 : }
301 : }
302 : }
303 :
304 : // Now we retrieve the atom numbers we need
305 69919 : if( atomsNeeded ) share( unique );
306 : }
307 :
308 70435 : void DomainDecomposition::share(const std::vector<AtomNumber>& unique) {
309 : // This retrieves what values we need to get
310 : int ndata=0; std::vector<Value*> values_to_get;
311 70435 : if(!(int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0)) {
312 2212 : uniq_index.resize(unique.size());
313 32705 : for(unsigned i=0; i<unique.size(); i++) uniq_index[i]=g2l[unique[i].index()];
314 13272 : for(const auto & ip : inputs) {
315 11060 : if( (!ip->fixed || firststep) && ip->wasset ) { (ip->mydata)->share_data( unique, uniq_index, ip->copyOutput(0) ); values_to_get.push_back(ip->copyOutput(0)); ndata++; }
316 : }
317 68223 : } else if( unique_serial) {
318 21259 : uniq_index.resize(unique.size());
319 789222 : for(unsigned i=0; i<unique.size(); i++) uniq_index[i]=unique[i].index();
320 127554 : for(const auto & ip : inputs) {
321 106295 : if( (!ip->fixed || firststep) && ip->wasset ) { (ip->mydata)->share_data( unique, uniq_index, ip->copyOutput(0) ); values_to_get.push_back(ip->copyOutput(0)); ndata++; }
322 : }
323 : } else {
324 : // faster version, which retrieves all atoms
325 281704 : for(const auto & ip : inputs) {
326 234756 : if( (!ip->fixed || firststep) && ip->wasset ) { (ip->mydata)->share_data( 0, getNumberOfAtoms(), ip->copyOutput(0) ); values_to_get.push_back(ip->copyOutput(0)); ndata++; }
327 : }
328 : }
329 :
330 70419 : if(dd && shuffledAtoms>0) {
331 2188 : if(dd.async) {
332 9598 : for(unsigned i=0; i<dd.mpi_request_positions.size(); i++) dd.mpi_request_positions[i].wait();
333 9598 : for(unsigned i=0; i<dd.mpi_request_index.size(); i++) dd.mpi_request_index[i].wait();
334 : }
335 :
336 2188 : int count=0;
337 32623 : for(const auto & p : unique) {
338 30435 : dd.indexToBeSent[count]=p.index(); int dpoint=0;
339 126694 : for(unsigned i=0; i<values_to_get.size(); ++i) {
340 96259 : dd.positionsToBeSent[ndata*count+dpoint]=values_to_get[i]->get(p.index());
341 96259 : dpoint++;
342 : }
343 30435 : count++;
344 : }
345 :
346 2188 : if(dd.async) {
347 2168 : asyncSent=true;
348 2168 : dd.mpi_request_positions.resize(dd.Get_size());
349 2168 : dd.mpi_request_index.resize(dd.Get_size());
350 9802 : for(int i=0; i<dd.Get_size(); i++) {
351 7634 : dd.mpi_request_index[i]=dd.Isend(&dd.indexToBeSent[0],count,i,666);
352 7634 : dd.mpi_request_positions[i]=dd.Isend(&dd.positionsToBeSent[0],ndata*count,i,667);
353 : }
354 : } else {
355 20 : const int n=(dd.Get_size());
356 36 : std::vector<int> counts(n);
357 20 : std::vector<int> displ(n);
358 20 : std::vector<int> counts5(n);
359 20 : std::vector<int> displ5(n);
360 20 : dd.Allgather(count,counts);
361 20 : displ[0]=0;
362 80 : for(int i=1; i<n; ++i) displ[i]=displ[i-1]+counts[i-1];
363 100 : for(int i=0; i<n; ++i) counts5[i]=counts[i]*ndata;
364 100 : for(int i=0; i<n; ++i) displ5[i]=displ[i]*ndata;
365 20 : dd.Allgatherv(&dd.indexToBeSent[0],count,&dd.indexToBeReceived[0],&counts[0],&displ[0]);
366 20 : dd.Allgatherv(&dd.positionsToBeSent[0],ndata*count,&dd.positionsToBeReceived[0],&counts5[0],&displ5[0]);
367 20 : int tot=displ[n-1]+counts[n-1];
368 1620 : for(int i=0; i<tot; i++) {
369 : int dpoint=0;
370 7264 : for(unsigned j=0; j<values_to_get.size(); ++j) {
371 5664 : values_to_get[j]->data[dd.indexToBeReceived[i]] = dd.positionsToBeReceived[ndata*i+dpoint]; dpoint++;
372 : }
373 : }
374 : }
375 : }
376 70419 : }
377 :
378 70411 : void DomainDecomposition::wait() {
379 422466 : for(const auto & ip : inputs) ip->dataCanBeSet=false;
380 :
381 70411 : if(dd && shuffledAtoms>0) {
382 : int ndata=0; std::vector<Value*> values_to_set;
383 13128 : for(const auto & ip : inputs) {
384 10940 : if( (!ip->fixed || firststep) && ip->wasset ) { values_to_set.push_back(ip->copyOutput(0)); ndata++; }
385 : }
386 :
387 : // receive toBeReceived
388 2188 : if(asyncSent) {
389 : Communicator::Status status;
390 : std::size_t count=0;
391 9802 : for(int i=0; i<dd.Get_size(); i++) {
392 7634 : dd.Recv(&dd.indexToBeReceived[count],dd.indexToBeReceived.size()-count,i,666,status);
393 7634 : int c=status.Get_count<int>();
394 7634 : dd.Recv(&dd.positionsToBeReceived[ndata*count],dd.positionsToBeReceived.size()-ndata*count,i,667);
395 7634 : count+=c;
396 : }
397 68220 : for(int i=0; i<count; i++) {
398 : int dpoint=0;
399 276100 : for(unsigned j=0; j<values_to_set.size(); ++j) {
400 210048 : values_to_set[j]->set(dd.indexToBeReceived[i], dd.positionsToBeReceived[ndata*i+dpoint] );
401 210048 : dpoint++;
402 : }
403 : }
404 2168 : asyncSent=false;
405 : }
406 : }
407 70411 : }
408 :
409 0 : unsigned DomainDecomposition::getNumberOfForcesToRescale() const {
410 0 : return gatindex.size();
411 : }
412 :
413 70287 : void DomainDecomposition::apply() {
414 70287 : std::vector<unsigned> forced_uniq_index(forced_unique.size());
415 70287 : if(!(int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0)) {
416 24159 : for(unsigned i=0; i<forced_unique.size(); i++) forced_uniq_index[i]=g2l[forced_unique[i].index()];
417 : } else {
418 3901524 : for(unsigned i=0; i<forced_unique.size(); i++) forced_uniq_index[i]=forced_unique[i].index();
419 : }
420 421712 : for(const auto & ip : inputs) {
421 351427 : if( !(ip->getPntrToValue())->forcesWereAdded() || ip->noforce ) {
422 209330 : continue;
423 142097 : } else if( ip->wasscaled || (!unique_serial && int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0) ) {
424 79589 : (ip->mydata)->add_force( gatindex, ip->getPntrToValue() );
425 62508 : } else { (ip->mydata)->add_force( forced_unique, forced_uniq_index, ip->getPntrToValue() ); }
426 : }
427 70285 : }
428 :
429 70285 : void DomainDecomposition::reset() {
430 70285 : if( !unique_serial && int(gatindex.size())==getNumberOfAtoms() && shuffledAtoms==0 ) return;
431 : // This is an optimisation to ensure that we don't call std::fill over the whole forces
432 : // array if there are a small number of atoms passed between the MD code and PLUMED
433 23347 : if( dd && shuffledAtoms>0 ) getAllActiveAtoms( unique );
434 140082 : for(const auto & ip : inputs) (ip->copyOutput(0))->clearInputForce( unique );
435 : }
436 :
437 114 : void DomainDecomposition::writeBinary(std::ostream&o) {
438 684 : for(const auto & ip : inputs) ip->writeBinary(o);
439 114 : }
440 :
441 114 : void DomainDecomposition::readBinary(std::istream&i) {
442 684 : for(const auto & ip : inputs) ip->readBinary(i);
443 114 : }
444 :
445 66603 : void DomainDecomposition::broadcastToDomains( Value* val ) {
446 66603 : if( dd ) dd.Bcast( val->data, 0 );
447 66603 : }
448 :
449 3989 : void DomainDecomposition::sumOverDomains( Value* val ) {
450 3989 : if( dd && shuffledAtoms>0 ) dd.Sum( val->data );
451 3989 : }
452 :
453 900 : const long int& DomainDecomposition::getDdStep() const {
454 900 : return ddStep;
455 : }
456 :
457 459 : const std::vector<int>& DomainDecomposition::getGatindex() const {
458 459 : return gatindex;
459 : }
460 :
461 2183 : void DomainDecomposition::getAllActiveAtoms( std::vector<AtomNumber>& u ) {
462 : gch::small_vector<const std::vector<AtomNumber>*,32> vectors;
463 : vectors.reserve(actions.size());
464 21016 : for(unsigned i=0; i<actions.size(); i++) {
465 37666 : if(actions[i]->isActive()) {
466 13033 : if(!actions[i]->getUnique().empty()) {
467 : // unique are the local atoms
468 22912 : vectors.push_back(&actions[i]->getUnique());
469 : }
470 : }
471 : }
472 : u.clear();
473 2183 : mergeVectorTools::mergeSortedVectors(vectors,u);
474 2183 : }
475 :
476 116 : void DomainDecomposition::createFullList(const TypesafePtr & n) {
477 116 : if( firststep ) {
478 7 : int natoms = getNumberOfAtoms();
479 7 : n.set(int(natoms)); fullList.resize(natoms);
480 803 : for(unsigned i=0; i<natoms; i++) fullList[i]=i;
481 : } else {
482 : // We update here the unique list defined at Atoms::unique.
483 : // This is not very clear, and probably should be coded differently.
484 : // Hopefully this fix the longstanding issue with NAMD.
485 109 : getAllActiveAtoms( unique );
486 109 : fullList.clear(); fullList.reserve(unique.size());
487 5012 : for(const auto & p : unique) fullList.push_back(p.index());
488 109 : n.set(int(fullList.size()));
489 : }
490 116 : }
491 :
492 116 : void DomainDecomposition::getFullList(const TypesafePtr & x) {
493 : auto xx=x.template get<const int**>();
494 116 : if(!fullList.empty()) *xx=&fullList[0];
495 6 : else *xx=NULL;
496 116 : }
497 :
498 116 : void DomainDecomposition::clearFullList() {
499 116 : fullList.resize(0);
500 116 : }
501 :
502 : }
|
