# | Line 42 | Line 42 | |
---|---|---|
42 | #include "math/SquareMatrix3.hpp" | |
43 | #include "nonbonded/NonBondedInteraction.hpp" | |
44 | #include "brains/SnapshotManager.hpp" | |
45 | + | #include "brains/PairList.hpp" |
46 | ||
47 | using namespace std; | |
48 | namespace OpenMD { | |
# | Line 54 | Line 55 | namespace OpenMD { | |
55 | void ForceMatrixDecomposition::distributeInitialData() { | |
56 | snap_ = sman_->getCurrentSnapshot(); | |
57 | storageLayout_ = sman_->getStorageLayout(); | |
58 | + | ff_ = info_->getForceField(); |
59 | nLocal_ = snap_->getNumberOfAtoms(); | |
60 | nGroups_ = snap_->getNumberOfCutoffGroups(); | |
61 | ||
62 | + | // gather the information for atomtype IDs (atids): |
63 | + | identsLocal = info_->getIdentArray(); |
64 | + | AtomLocalToGlobal = info_->getGlobalAtomIndices(); |
65 | + | cgLocalToGlobal = info_->getGlobalGroupIndices(); |
66 | + | vector<int> globalGroupMembership = info_->getGlobalGroupMembership(); |
67 | + | vector<RealType> massFactorsLocal = info_->getMassFactors(); |
68 | + | PairList excludes = info_->getExcludedInteractions(); |
69 | + | PairList oneTwo = info_->getOneTwoInteractions(); |
70 | + | PairList oneThree = info_->getOneThreeInteractions(); |
71 | + | PairList oneFour = info_->getOneFourInteractions(); |
72 | + | vector<RealType> pot_local(N_INTERACTION_FAMILIES, 0.0); |
73 | + | |
74 | #ifdef IS_MPI | |
75 | ||
76 | AtomCommIntRow = new Communicator<Row,int>(nLocal_); | |
# | Line 93 | Line 107 | namespace OpenMD { | |
107 | vector<RealType> (nAtomsInRow_, 0.0)); | |
108 | vector<vector<RealType> > pot_col(N_INTERACTION_FAMILIES, | |
109 | vector<RealType> (nAtomsInCol_, 0.0)); | |
96 | – | |
97 | – | |
98 | – | vector<RealType> pot_local(N_INTERACTION_FAMILIES, 0.0); |
110 | ||
100 | – | // gather the information for atomtype IDs (atids): |
101 | – | vector<int> identsLocal = info_->getIdentArray(); |
111 | identsRow.reserve(nAtomsInRow_); | |
112 | identsCol.reserve(nAtomsInCol_); | |
113 | ||
114 | AtomCommIntRow->gather(identsLocal, identsRow); | |
115 | AtomCommIntColumn->gather(identsLocal, identsCol); | |
116 | ||
108 | – | AtomLocalToGlobal = info_->getGlobalAtomIndices(); |
117 | AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal); | |
118 | AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal); | |
119 | ||
112 | – | cgLocalToGlobal = info_->getGlobalGroupIndices(); |
120 | cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal); | |
121 | cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal); | |
122 | ||
123 | < | // still need: |
124 | < | // topoDist |
125 | < | // exclude |
123 | > | AtomCommRealRow->gather(massFactorsLocal, massFactorsRow); |
124 | > | AtomCommRealColumn->gather(massFactorsLocal, massFactorsCol); |
125 | > | |
126 | > | groupListRow_.clear(); |
127 | > | groupListRow_.reserve(nGroupsInRow_); |
128 | > | for (int i = 0; i < nGroupsInRow_; i++) { |
129 | > | int gid = cgRowToGlobal[i]; |
130 | > | for (int j = 0; j < nAtomsInRow_; j++) { |
131 | > | int aid = AtomRowToGlobal[j]; |
132 | > | if (globalGroupMembership[aid] == gid) |
133 | > | groupListRow_[i].push_back(j); |
134 | > | } |
135 | > | } |
136 | > | |
137 | > | groupListCol_.clear(); |
138 | > | groupListCol_.reserve(nGroupsInCol_); |
139 | > | for (int i = 0; i < nGroupsInCol_; i++) { |
140 | > | int gid = cgColToGlobal[i]; |
141 | > | for (int j = 0; j < nAtomsInCol_; j++) { |
142 | > | int aid = AtomColToGlobal[j]; |
143 | > | if (globalGroupMembership[aid] == gid) |
144 | > | groupListCol_[i].push_back(j); |
145 | > | } |
146 | > | } |
147 | > | |
148 | > | skipsForRowAtom.clear(); |
149 | > | skipsForRowAtom.reserve(nAtomsInRow_); |
150 | > | for (int i = 0; i < nAtomsInRow_; i++) { |
151 | > | int iglob = AtomRowToGlobal[i]; |
152 | > | for (int j = 0; j < nAtomsInCol_; j++) { |
153 | > | int jglob = AtomColToGlobal[j]; |
154 | > | if (excludes.hasPair(iglob, jglob)) |
155 | > | skipsForRowAtom[i].push_back(j); |
156 | > | } |
157 | > | } |
158 | > | |
159 | > | toposForRowAtom.clear(); |
160 | > | toposForRowAtom.reserve(nAtomsInRow_); |
161 | > | for (int i = 0; i < nAtomsInRow_; i++) { |
162 | > | int iglob = AtomRowToGlobal[i]; |
163 | > | int nTopos = 0; |
164 | > | for (int j = 0; j < nAtomsInCol_; j++) { |
165 | > | int jglob = AtomColToGlobal[j]; |
166 | > | if (oneTwo.hasPair(iglob, jglob)) { |
167 | > | toposForRowAtom[i].push_back(j); |
168 | > | topoDistRow[i][nTopos] = 1; |
169 | > | nTopos++; |
170 | > | } |
171 | > | if (oneThree.hasPair(iglob, jglob)) { |
172 | > | toposForRowAtom[i].push_back(j); |
173 | > | topoDistRow[i][nTopos] = 2; |
174 | > | nTopos++; |
175 | > | } |
176 | > | if (oneFour.hasPair(iglob, jglob)) { |
177 | > | toposForRowAtom[i].push_back(j); |
178 | > | topoDistRow[i][nTopos] = 3; |
179 | > | nTopos++; |
180 | > | } |
181 | > | } |
182 | > | } |
183 | > | |
184 | #endif | |
120 | – | } |
121 | – | |
185 | ||
186 | + | groupList_.clear(); |
187 | + | groupList_.reserve(nGroups_); |
188 | + | for (int i = 0; i < nGroups_; i++) { |
189 | + | int gid = cgLocalToGlobal[i]; |
190 | + | for (int j = 0; j < nLocal_; j++) { |
191 | + | int aid = AtomLocalToGlobal[j]; |
192 | + | if (globalGroupMembership[aid] == gid) |
193 | + | groupList_[i].push_back(j); |
194 | + | } |
195 | + | } |
196 | ||
197 | + | skipsForLocalAtom.clear(); |
198 | + | skipsForLocalAtom.reserve(nLocal_); |
199 | + | |
200 | + | for (int i = 0; i < nLocal_; i++) { |
201 | + | int iglob = AtomLocalToGlobal[i]; |
202 | + | for (int j = 0; j < nLocal_; j++) { |
203 | + | int jglob = AtomLocalToGlobal[j]; |
204 | + | if (excludes.hasPair(iglob, jglob)) |
205 | + | skipsForLocalAtom[i].push_back(j); |
206 | + | } |
207 | + | } |
208 | + | |
209 | + | toposForLocalAtom.clear(); |
210 | + | toposForLocalAtom.reserve(nLocal_); |
211 | + | for (int i = 0; i < nLocal_; i++) { |
212 | + | int iglob = AtomLocalToGlobal[i]; |
213 | + | int nTopos = 0; |
214 | + | for (int j = 0; j < nLocal_; j++) { |
215 | + | int jglob = AtomLocalToGlobal[j]; |
216 | + | if (oneTwo.hasPair(iglob, jglob)) { |
217 | + | toposForLocalAtom[i].push_back(j); |
218 | + | topoDistLocal[i][nTopos] = 1; |
219 | + | nTopos++; |
220 | + | } |
221 | + | if (oneThree.hasPair(iglob, jglob)) { |
222 | + | toposForLocalAtom[i].push_back(j); |
223 | + | topoDistLocal[i][nTopos] = 2; |
224 | + | nTopos++; |
225 | + | } |
226 | + | if (oneFour.hasPair(iglob, jglob)) { |
227 | + | toposForLocalAtom[i].push_back(j); |
228 | + | topoDistLocal[i][nTopos] = 3; |
229 | + | nTopos++; |
230 | + | } |
231 | + | } |
232 | + | } |
233 | + | } |
234 | + | |
235 | void ForceMatrixDecomposition::distributeData() { | |
236 | snap_ = sman_->getCurrentSnapshot(); | |
237 | storageLayout_ = sman_->getStorageLayout(); | |
# | Line 244 | Line 355 | namespace OpenMD { | |
355 | #endif | |
356 | } | |
357 | ||
358 | + | int ForceMatrixDecomposition::getNAtomsInRow() { |
359 | + | #ifdef IS_MPI |
360 | + | return nAtomsInRow_; |
361 | + | #else |
362 | + | return nLocal_; |
363 | + | #endif |
364 | + | } |
365 | + | |
366 | + | /** |
367 | + | * returns the list of atoms belonging to this group. |
368 | + | */ |
369 | + | vector<int> ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){ |
370 | + | #ifdef IS_MPI |
371 | + | return groupListRow_[cg1]; |
372 | + | #else |
373 | + | return groupList_[cg1]; |
374 | + | #endif |
375 | + | } |
376 | + | |
377 | + | vector<int> ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){ |
378 | + | #ifdef IS_MPI |
379 | + | return groupListCol_[cg2]; |
380 | + | #else |
381 | + | return groupList_[cg2]; |
382 | + | #endif |
383 | + | } |
384 | ||
385 | Vector3d ForceMatrixDecomposition::getIntergroupVector(int cg1, int cg2){ | |
386 | Vector3d d; | |
# | Line 285 | Line 422 | namespace OpenMD { | |
422 | snap_->wrapVector(d); | |
423 | return d; | |
424 | } | |
425 | + | |
426 | + | RealType ForceMatrixDecomposition::getMassFactorRow(int atom1) { |
427 | + | #ifdef IS_MPI |
428 | + | return massFactorsRow[atom1]; |
429 | + | #else |
430 | + | return massFactorsLocal[atom1]; |
431 | + | #endif |
432 | + | } |
433 | + | |
434 | + | RealType ForceMatrixDecomposition::getMassFactorColumn(int atom2) { |
435 | + | #ifdef IS_MPI |
436 | + | return massFactorsCol[atom2]; |
437 | + | #else |
438 | + | return massFactorsLocal[atom2]; |
439 | + | #endif |
440 | + | |
441 | + | } |
442 | ||
443 | Vector3d ForceMatrixDecomposition::getInteratomicVector(int atom1, int atom2){ | |
444 | Vector3d d; | |
# | Line 297 | Line 451 | namespace OpenMD { | |
451 | ||
452 | snap_->wrapVector(d); | |
453 | return d; | |
454 | + | } |
455 | + | |
456 | + | vector<int> ForceMatrixDecomposition::getSkipsForRowAtom(int atom1) { |
457 | + | #ifdef IS_MPI |
458 | + | return skipsForRowAtom[atom1]; |
459 | + | #else |
460 | + | return skipsForLocalAtom[atom1]; |
461 | + | #endif |
462 | + | } |
463 | + | |
464 | + | /** |
465 | + | * there are a number of reasons to skip a pair or a particle mostly |
466 | + | * we do this to exclude atoms who are involved in short range |
467 | + | * interactions (bonds, bends, torsions), but we also need to |
468 | + | * exclude some overcounted interactions that result from the |
469 | + | * parallel decomposition. |
470 | + | */ |
471 | + | bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) { |
472 | + | int unique_id_1, unique_id_2; |
473 | + | |
474 | + | #ifdef IS_MPI |
475 | + | // in MPI, we have to look up the unique IDs for each atom |
476 | + | unique_id_1 = AtomRowToGlobal[atom1]; |
477 | + | unique_id_2 = AtomColToGlobal[atom2]; |
478 | + | |
479 | + | // this situation should only arise in MPI simulations |
480 | + | if (unique_id_1 == unique_id_2) return true; |
481 | + | |
482 | + | // this prevents us from doing the pair on multiple processors |
483 | + | if (unique_id_1 < unique_id_2) { |
484 | + | if ((unique_id_1 + unique_id_2) % 2 == 0) return true; |
485 | + | } else { |
486 | + | if ((unique_id_1 + unique_id_2) % 2 == 1) return true; |
487 | + | } |
488 | + | #else |
489 | + | // in the normal loop, the atom numbers are unique |
490 | + | unique_id_1 = atom1; |
491 | + | unique_id_2 = atom2; |
492 | + | #endif |
493 | + | |
494 | + | #ifdef IS_MPI |
495 | + | for (vector<int>::iterator i = skipsForRowAtom[atom1].begin(); |
496 | + | i != skipsForRowAtom[atom1].end(); ++i) { |
497 | + | if ( (*i) == unique_id_2 ) return true; |
498 | + | } |
499 | + | #else |
500 | + | for (vector<int>::iterator i = skipsForLocalAtom[atom1].begin(); |
501 | + | i != skipsForLocalAtom[atom1].end(); ++i) { |
502 | + | if ( (*i) == unique_id_2 ) return true; |
503 | + | } |
504 | + | #endif |
505 | } | |
506 | ||
507 | + | int ForceMatrixDecomposition::getTopoDistance(int atom1, int atom2) { |
508 | + | |
509 | + | #ifdef IS_MPI |
510 | + | for (int i = 0; i < toposForRowAtom[atom1].size(); i++) { |
511 | + | if ( toposForRowAtom[atom1][i] == atom2 ) return topoDistRow[atom1][i]; |
512 | + | } |
513 | + | #else |
514 | + | for (int i = 0; i < toposForLocalAtom[atom1].size(); i++) { |
515 | + | if ( toposForLocalAtom[atom1][i] == atom2 ) return topoDistLocal[atom1][i]; |
516 | + | } |
517 | + | #endif |
518 | + | |
519 | + | // zero is default for unconnected (i.e. normal) pair interactions |
520 | + | return 0; |
521 | + | } |
522 | + | |
523 | void ForceMatrixDecomposition::addForceToAtomRow(int atom1, Vector3d fg){ | |
524 | #ifdef IS_MPI | |
525 | atomRowData.force[atom1] += fg; | |
# | Line 320 | Line 541 | namespace OpenMD { | |
541 | InteractionData idat; | |
542 | ||
543 | #ifdef IS_MPI | |
544 | + | |
545 | + | idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), |
546 | + | ff_->getAtomType(identsCol[atom2]) ); |
547 | + | |
548 | if (storageLayout_ & DataStorage::dslAmat) { | |
549 | idat.A1 = &(atomRowData.aMat[atom1]); | |
550 | idat.A2 = &(atomColData.aMat[atom2]); | |
# | Line 344 | Line 569 | namespace OpenMD { | |
569 | idat.dfrho1 = &(atomRowData.functionalDerivative[atom1]); | |
570 | idat.dfrho2 = &(atomColData.functionalDerivative[atom2]); | |
571 | } | |
572 | + | |
573 | #else | |
574 | + | |
575 | + | idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]), |
576 | + | ff_->getAtomType(identsLocal[atom2]) ); |
577 | + | |
578 | if (storageLayout_ & DataStorage::dslAmat) { | |
579 | idat.A1 = &(snap_->atomData.aMat[atom1]); | |
580 | idat.A2 = &(snap_->atomData.aMat[atom2]); | |
# | Line 377 | Line 607 | namespace OpenMD { | |
607 | ||
608 | InteractionData idat; | |
609 | #ifdef IS_MPI | |
610 | + | idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), |
611 | + | ff_->getAtomType(identsCol[atom2]) ); |
612 | + | |
613 | if (storageLayout_ & DataStorage::dslElectroFrame) { | |
614 | idat.eFrame1 = &(atomRowData.electroFrame[atom1]); | |
615 | idat.eFrame2 = &(atomColData.electroFrame[atom2]); | |
# | Line 390 | Line 623 | namespace OpenMD { | |
623 | idat.t2 = &(atomColData.force[atom2]); | |
624 | } | |
625 | #else | |
626 | + | idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]), |
627 | + | ff_->getAtomType(identsLocal[atom2]) ); |
628 | + | |
629 | if (storageLayout_ & DataStorage::dslElectroFrame) { | |
630 | idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]); | |
631 | idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]); | |
# | Line 402 | Line 638 | namespace OpenMD { | |
638 | idat.t1 = &(snap_->atomData.force[atom1]); | |
639 | idat.t2 = &(snap_->atomData.force[atom2]); | |
640 | } | |
641 | < | #endif |
406 | < | |
641 | > | #endif |
642 | } | |
643 | ||
409 | – | |
410 | – | |
411 | – | |
644 | /* | |
645 | * buildNeighborList | |
646 | * |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |