# | 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 | < | #ifdef IS_MPI |
59 | < | int nLocal = snap_->getNumberOfAtoms(); |
59 | < | int nGroups = snap_->getNumberOfCutoffGroups(); |
60 | < | |
61 | < | AtomCommIntRow = new Communicator<Row,int>(nLocal); |
62 | < | AtomCommRealRow = new Communicator<Row,RealType>(nLocal); |
63 | < | AtomCommVectorRow = new Communicator<Row,Vector3d>(nLocal); |
64 | < | AtomCommMatrixRow = new Communicator<Row,Mat3x3d>(nLocal); |
58 | > | nLocal_ = snap_->getNumberOfAtoms(); |
59 | > | nGroups_ = snap_->getNumberOfCutoffGroups(); |
60 | ||
61 | < | AtomCommIntColumn = new Communicator<Column,int>(nLocal); |
62 | < | AtomCommRealColumn = new Communicator<Column,RealType>(nLocal); |
63 | < | AtomCommVectorColumn = new Communicator<Column,Vector3d>(nLocal); |
64 | < | AtomCommMatrixColumn = new Communicator<Column,Mat3x3d>(nLocal); |
61 | > | // gather the information for atomtype IDs (atids): |
62 | > | vector<int> identsLocal = info_->getIdentArray(); |
63 | > | AtomLocalToGlobal = info_->getGlobalAtomIndices(); |
64 | > | cgLocalToGlobal = info_->getGlobalGroupIndices(); |
65 | > | vector<int> globalGroupMembership = info_->getGlobalGroupMembership(); |
66 | > | vector<RealType> massFactorsLocal = info_->getMassFactors(); |
67 | > | PairList excludes = info_->getExcludedInteractions(); |
68 | > | PairList oneTwo = info_->getOneTwoInteractions(); |
69 | > | PairList oneThree = info_->getOneThreeInteractions(); |
70 | > | PairList oneFour = info_->getOneFourInteractions(); |
71 | > | vector<RealType> pot_local(N_INTERACTION_FAMILIES, 0.0); |
72 | ||
73 | < | cgCommIntRow = new Communicator<Row,int>(nGroups); |
74 | < | cgCommVectorRow = new Communicator<Row,Vector3d>(nGroups); |
75 | < | cgCommIntColumn = new Communicator<Column,int>(nGroups); |
76 | < | cgCommVectorColumn = new Communicator<Column,Vector3d>(nGroups); |
73 | > | #ifdef IS_MPI |
74 | > | |
75 | > | AtomCommIntRow = new Communicator<Row,int>(nLocal_); |
76 | > | AtomCommRealRow = new Communicator<Row,RealType>(nLocal_); |
77 | > | AtomCommVectorRow = new Communicator<Row,Vector3d>(nLocal_); |
78 | > | AtomCommMatrixRow = new Communicator<Row,Mat3x3d>(nLocal_); |
79 | ||
80 | < | int nAtomsInRow = AtomCommIntRow->getSize(); |
81 | < | int nAtomsInCol = AtomCommIntColumn->getSize(); |
82 | < | int nGroupsInRow = cgCommIntRow->getSize(); |
83 | < | int nGroupsInCol = cgCommIntColumn->getSize(); |
80 | > | AtomCommIntColumn = new Communicator<Column,int>(nLocal_); |
81 | > | AtomCommRealColumn = new Communicator<Column,RealType>(nLocal_); |
82 | > | AtomCommVectorColumn = new Communicator<Column,Vector3d>(nLocal_); |
83 | > | AtomCommMatrixColumn = new Communicator<Column,Mat3x3d>(nLocal_); |
84 | ||
85 | + | cgCommIntRow = new Communicator<Row,int>(nGroups_); |
86 | + | cgCommVectorRow = new Communicator<Row,Vector3d>(nGroups_); |
87 | + | cgCommIntColumn = new Communicator<Column,int>(nGroups_); |
88 | + | cgCommVectorColumn = new Communicator<Column,Vector3d>(nGroups_); |
89 | + | |
90 | + | nAtomsInRow_ = AtomCommIntRow->getSize(); |
91 | + | nAtomsInCol_ = AtomCommIntColumn->getSize(); |
92 | + | nGroupsInRow_ = cgCommIntRow->getSize(); |
93 | + | nGroupsInCol_ = cgCommIntColumn->getSize(); |
94 | + | |
95 | // Modify the data storage objects with the correct layouts and sizes: | |
96 | < | atomRowData.resize(nAtomsInRow); |
96 | > | atomRowData.resize(nAtomsInRow_); |
97 | atomRowData.setStorageLayout(storageLayout_); | |
98 | < | atomColData.resize(nAtomsInCol); |
98 | > | atomColData.resize(nAtomsInCol_); |
99 | atomColData.setStorageLayout(storageLayout_); | |
100 | < | cgRowData.resize(nGroupsInRow); |
100 | > | cgRowData.resize(nGroupsInRow_); |
101 | cgRowData.setStorageLayout(DataStorage::dslPosition); | |
102 | < | cgColData.resize(nGroupsInCol); |
102 | > | cgColData.resize(nGroupsInCol_); |
103 | cgColData.setStorageLayout(DataStorage::dslPosition); | |
104 | ||
105 | vector<vector<RealType> > pot_row(N_INTERACTION_FAMILIES, | |
106 | < | vector<RealType> (nAtomsInRow, 0.0)); |
106 | > | vector<RealType> (nAtomsInRow_, 0.0)); |
107 | vector<vector<RealType> > pot_col(N_INTERACTION_FAMILIES, | |
108 | < | vector<RealType> (nAtomsInCol, 0.0)); |
95 | < | |
96 | < | |
97 | < | vector<RealType> pot_local(N_INTERACTION_FAMILIES, 0.0); |
108 | > | vector<RealType> (nAtomsInCol_, 0.0)); |
109 | ||
110 | < | // gather the information for atomtype IDs (atids): |
111 | < | vector<int> identsLocal = info_->getIdentArray(); |
101 | < | identsRow.reserve(nAtomsInRow); |
102 | < | identsCol.reserve(nAtomsInCol); |
110 | > | identsRow.reserve(nAtomsInRow_); |
111 | > | identsCol.reserve(nAtomsInCol_); |
112 | ||
113 | AtomCommIntRow->gather(identsLocal, identsRow); | |
114 | AtomCommIntColumn->gather(identsLocal, identsCol); | |
115 | ||
107 | – | AtomLocalToGlobal = info_->getGlobalAtomIndices(); |
116 | AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal); | |
117 | AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal); | |
118 | ||
111 | – | cgLocalToGlobal = info_->getGlobalGroupIndices(); |
119 | cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal); | |
120 | cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal); | |
121 | ||
122 | < | // still need: |
123 | < | // topoDist |
124 | < | // exclude |
122 | > | AtomCommRealRow->gather(massFactorsLocal, massFactorsRow); |
123 | > | AtomCommRealColumn->gather(massFactorsLocal, massFactorsCol); |
124 | > | |
125 | > | groupListRow_.clear(); |
126 | > | groupListRow_.reserve(nGroupsInRow_); |
127 | > | for (int i = 0; i < nGroupsInRow_; i++) { |
128 | > | int gid = cgRowToGlobal[i]; |
129 | > | for (int j = 0; j < nAtomsInRow_; j++) { |
130 | > | int aid = AtomRowToGlobal[j]; |
131 | > | if (globalGroupMembership[aid] == gid) |
132 | > | groupListRow_[i].push_back(j); |
133 | > | } |
134 | > | } |
135 | > | |
136 | > | groupListCol_.clear(); |
137 | > | groupListCol_.reserve(nGroupsInCol_); |
138 | > | for (int i = 0; i < nGroupsInCol_; i++) { |
139 | > | int gid = cgColToGlobal[i]; |
140 | > | for (int j = 0; j < nAtomsInCol_; j++) { |
141 | > | int aid = AtomColToGlobal[j]; |
142 | > | if (globalGroupMembership[aid] == gid) |
143 | > | groupListCol_[i].push_back(j); |
144 | > | } |
145 | > | } |
146 | > | |
147 | > | skipsForRowAtom.clear(); |
148 | > | skipsForRowAtom.reserve(nAtomsInRow_); |
149 | > | for (int i = 0; i < nAtomsInRow_; i++) { |
150 | > | int iglob = AtomColToGlobal[i]; |
151 | > | for (int j = 0; j < nAtomsInCol_; j++) { |
152 | > | int jglob = AtomRowToGlobal[j]; |
153 | > | if (excludes.hasPair(iglob, jglob)) |
154 | > | skipsForRowAtom[i].push_back(j); |
155 | > | } |
156 | > | } |
157 | > | |
158 | > | toposForRowAtom.clear(); |
159 | > | toposForRowAtom.reserve(nAtomsInRow_); |
160 | > | for (int i = 0; i < nAtomsInRow_; i++) { |
161 | > | int iglob = AtomColToGlobal[i]; |
162 | > | int nTopos = 0; |
163 | > | for (int j = 0; j < nAtomsInCol_; j++) { |
164 | > | int jglob = AtomRowToGlobal[j]; |
165 | > | if (oneTwo.hasPair(iglob, jglob)) { |
166 | > | toposForRowAtom[i].push_back(j); |
167 | > | topoDistRow[i][nTopos] = 1; |
168 | > | nTopos++; |
169 | > | } |
170 | > | if (oneThree.hasPair(iglob, jglob)) { |
171 | > | toposForRowAtom[i].push_back(j); |
172 | > | topoDistRow[i][nTopos] = 2; |
173 | > | nTopos++; |
174 | > | } |
175 | > | if (oneFour.hasPair(iglob, jglob)) { |
176 | > | toposForRowAtom[i].push_back(j); |
177 | > | topoDistRow[i][nTopos] = 3; |
178 | > | nTopos++; |
179 | > | } |
180 | > | } |
181 | > | } |
182 | > | |
183 | #endif | |
119 | – | } |
120 | – | |
184 | ||
185 | + | groupList_.clear(); |
186 | + | groupList_.reserve(nGroups_); |
187 | + | for (int i = 0; i < nGroups_; i++) { |
188 | + | int gid = cgLocalToGlobal[i]; |
189 | + | for (int j = 0; j < nLocal_; j++) { |
190 | + | int aid = AtomLocalToGlobal[j]; |
191 | + | if (globalGroupMembership[aid] == gid) |
192 | + | groupList_[i].push_back(j); |
193 | + | } |
194 | + | } |
195 | ||
196 | + | skipsForLocalAtom.clear(); |
197 | + | skipsForLocalAtom.reserve(nLocal_); |
198 | + | |
199 | + | for (int i = 0; i < nLocal_; i++) { |
200 | + | int iglob = AtomLocalToGlobal[i]; |
201 | + | for (int j = 0; j < nLocal_; j++) { |
202 | + | int jglob = AtomLocalToGlobal[j]; |
203 | + | if (excludes.hasPair(iglob, jglob)) |
204 | + | skipsForLocalAtom[i].push_back(j); |
205 | + | } |
206 | + | } |
207 | + | |
208 | + | toposForLocalAtom.clear(); |
209 | + | toposForLocalAtom.reserve(nLocal_); |
210 | + | for (int i = 0; i < nLocal_; i++) { |
211 | + | int iglob = AtomLocalToGlobal[i]; |
212 | + | int nTopos = 0; |
213 | + | for (int j = 0; j < nLocal_; j++) { |
214 | + | int jglob = AtomLocalToGlobal[j]; |
215 | + | if (oneTwo.hasPair(iglob, jglob)) { |
216 | + | toposForLocalAtom[i].push_back(j); |
217 | + | topoDistLocal[i][nTopos] = 1; |
218 | + | nTopos++; |
219 | + | } |
220 | + | if (oneThree.hasPair(iglob, jglob)) { |
221 | + | toposForLocalAtom[i].push_back(j); |
222 | + | topoDistLocal[i][nTopos] = 2; |
223 | + | nTopos++; |
224 | + | } |
225 | + | if (oneFour.hasPair(iglob, jglob)) { |
226 | + | toposForLocalAtom[i].push_back(j); |
227 | + | topoDistLocal[i][nTopos] = 3; |
228 | + | nTopos++; |
229 | + | } |
230 | + | } |
231 | + | } |
232 | + | } |
233 | + | |
234 | void ForceMatrixDecomposition::distributeData() { | |
235 | snap_ = sman_->getCurrentSnapshot(); | |
236 | storageLayout_ = sman_->getStorageLayout(); | |
# | Line 229 | Line 340 | namespace OpenMD { | |
340 | snap_->atomData.torque[i] += trq_tmp[i]; | |
341 | } | |
342 | ||
343 | < | int nLocal = snap_->getNumberOfAtoms(); |
343 | > | nLocal_ = snap_->getNumberOfAtoms(); |
344 | ||
345 | vector<vector<RealType> > pot_temp(N_INTERACTION_FAMILIES, | |
346 | < | vector<RealType> (nLocal, 0.0)); |
346 | > | vector<RealType> (nLocal_, 0.0)); |
347 | ||
348 | for (int i = 0; i < N_INTERACTION_FAMILIES; i++) { | |
349 | AtomCommRealRow->scatter(pot_row[i], pot_temp[i]); | |
# | Line 243 | Line 354 | namespace OpenMD { | |
354 | #endif | |
355 | } | |
356 | ||
357 | + | int ForceMatrixDecomposition::getNAtomsInRow() { |
358 | + | #ifdef IS_MPI |
359 | + | return nAtomsInRow_; |
360 | + | #else |
361 | + | return nLocal_; |
362 | + | #endif |
363 | + | } |
364 | + | |
365 | + | /** |
366 | + | * returns the list of atoms belonging to this group. |
367 | + | */ |
368 | + | vector<int> ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){ |
369 | + | #ifdef IS_MPI |
370 | + | return groupListRow_[cg1]; |
371 | + | #else |
372 | + | return groupList_[cg1]; |
373 | + | #endif |
374 | + | } |
375 | + | |
376 | + | vector<int> ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){ |
377 | + | #ifdef IS_MPI |
378 | + | return groupListCol_[cg2]; |
379 | + | #else |
380 | + | return groupList_[cg2]; |
381 | + | #endif |
382 | + | } |
383 | ||
384 | Vector3d ForceMatrixDecomposition::getIntergroupVector(int cg1, int cg2){ | |
385 | Vector3d d; | |
# | Line 284 | Line 421 | namespace OpenMD { | |
421 | snap_->wrapVector(d); | |
422 | return d; | |
423 | } | |
424 | + | |
425 | + | RealType ForceMatrixDecomposition::getMassFactorRow(int atom1) { |
426 | + | #ifdef IS_MPI |
427 | + | return massFactorsRow[atom1]; |
428 | + | #else |
429 | + | return massFactorsLocal[atom1]; |
430 | + | #endif |
431 | + | } |
432 | + | |
433 | + | RealType ForceMatrixDecomposition::getMassFactorColumn(int atom2) { |
434 | + | #ifdef IS_MPI |
435 | + | return massFactorsCol[atom2]; |
436 | + | #else |
437 | + | return massFactorsLocal[atom2]; |
438 | + | #endif |
439 | + | |
440 | + | } |
441 | ||
442 | Vector3d ForceMatrixDecomposition::getInteratomicVector(int atom1, int atom2){ | |
443 | Vector3d d; | |
# | Line 296 | Line 450 | namespace OpenMD { | |
450 | ||
451 | snap_->wrapVector(d); | |
452 | return d; | |
453 | + | } |
454 | + | |
455 | + | vector<int> ForceMatrixDecomposition::getSkipsForRowAtom(int atom1) { |
456 | + | #ifdef IS_MPI |
457 | + | return skipsForRowAtom[atom1]; |
458 | + | #else |
459 | + | return skipsForLocalAtom[atom1]; |
460 | + | #endif |
461 | + | } |
462 | + | |
463 | + | /** |
464 | + | * there are a number of reasons to skip a pair or a particle mostly |
465 | + | * we do this to exclude atoms who are involved in short range |
466 | + | * interactions (bonds, bends, torsions), but we also need to |
467 | + | * exclude some overcounted interactions that result from the |
468 | + | * parallel decomposition. |
469 | + | */ |
470 | + | bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) { |
471 | + | int unique_id_1, unique_id_2; |
472 | + | |
473 | + | #ifdef IS_MPI |
474 | + | // in MPI, we have to look up the unique IDs for each atom |
475 | + | unique_id_1 = AtomRowToGlobal[atom1]; |
476 | + | unique_id_2 = AtomColToGlobal[atom2]; |
477 | + | |
478 | + | // this situation should only arise in MPI simulations |
479 | + | if (unique_id_1 == unique_id_2) return true; |
480 | + | |
481 | + | // this prevents us from doing the pair on multiple processors |
482 | + | if (unique_id_1 < unique_id_2) { |
483 | + | if ((unique_id_1 + unique_id_2) % 2 == 0) return true; |
484 | + | } else { |
485 | + | if ((unique_id_1 + unique_id_2) % 2 == 1) return true; |
486 | + | } |
487 | + | #else |
488 | + | // in the normal loop, the atom numbers are unique |
489 | + | unique_id_1 = atom1; |
490 | + | unique_id_2 = atom2; |
491 | + | #endif |
492 | + | |
493 | + | #ifdef IS_MPI |
494 | + | for (vector<int>::iterator i = skipsForRowAtom[atom1].begin(); |
495 | + | i != skipsForRowAtom[atom1].end(); ++i) { |
496 | + | if ( (*i) == unique_id_2 ) return true; |
497 | + | } |
498 | + | #else |
499 | + | for (vector<int>::iterator i = skipsForLocalAtom[atom1].begin(); |
500 | + | i != skipsForLocalAtom[atom1].end(); ++i) { |
501 | + | if ( (*i) == unique_id_2 ) return true; |
502 | + | } |
503 | + | #endif |
504 | + | } |
505 | + | |
506 | + | int ForceMatrixDecomposition::getTopoDistance(int atom1, int atom2) { |
507 | + | |
508 | + | #ifdef IS_MPI |
509 | + | for (int i = 0; i < toposForRowAtom[atom1].size(); i++) { |
510 | + | if ( toposForRowAtom[atom1][i] == atom2 ) return topoDistRow[atom1][i]; |
511 | + | } |
512 | + | #else |
513 | + | for (int i = 0; i < toposForLocalAtom[atom1].size(); i++) { |
514 | + | if ( toposForLocalAtom[atom1][i] == atom2 ) return topoDistLocal[atom1][i]; |
515 | + | } |
516 | + | #endif |
517 | + | |
518 | + | // zero is default for unconnected (i.e. normal) pair interactions |
519 | + | return 0; |
520 | } | |
521 | ||
522 | void ForceMatrixDecomposition::addForceToAtomRow(int atom1, Vector3d fg){ | |
# | Line 312 | Line 533 | namespace OpenMD { | |
533 | #else | |
534 | snap_->atomData.force[atom2] += fg; | |
535 | #endif | |
315 | – | |
536 | } | |
537 | ||
538 | // filling interaction blocks with pointers | |
539 | InteractionData ForceMatrixDecomposition::fillInteractionData(int atom1, int atom2) { | |
320 | – | |
540 | InteractionData idat; | |
541 | + | |
542 | #ifdef IS_MPI | |
543 | if (storageLayout_ & DataStorage::dslAmat) { | |
544 | < | idat.A1 = atomRowData.aMat[atom1]; |
545 | < | idat.A2 = atomColData.aMat[atom2]; |
544 | > | idat.A1 = &(atomRowData.aMat[atom1]); |
545 | > | idat.A2 = &(atomColData.aMat[atom2]); |
546 | } | |
547 | + | |
548 | + | if (storageLayout_ & DataStorage::dslElectroFrame) { |
549 | + | idat.eFrame1 = &(atomRowData.electroFrame[atom1]); |
550 | + | idat.eFrame2 = &(atomColData.electroFrame[atom2]); |
551 | + | } |
552 | ||
553 | + | if (storageLayout_ & DataStorage::dslTorque) { |
554 | + | idat.t1 = &(atomRowData.torque[atom1]); |
555 | + | idat.t2 = &(atomColData.torque[atom2]); |
556 | + | } |
557 | + | |
558 | + | if (storageLayout_ & DataStorage::dslDensity) { |
559 | + | idat.rho1 = &(atomRowData.density[atom1]); |
560 | + | idat.rho2 = &(atomColData.density[atom2]); |
561 | + | } |
562 | + | |
563 | + | if (storageLayout_ & DataStorage::dslFunctionalDerivative) { |
564 | + | idat.dfrho1 = &(atomRowData.functionalDerivative[atom1]); |
565 | + | idat.dfrho2 = &(atomColData.functionalDerivative[atom2]); |
566 | + | } |
567 | + | |
568 | + | #else |
569 | + | if (storageLayout_ & DataStorage::dslAmat) { |
570 | + | idat.A1 = &(snap_->atomData.aMat[atom1]); |
571 | + | idat.A2 = &(snap_->atomData.aMat[atom2]); |
572 | + | } |
573 | + | |
574 | if (storageLayout_ & DataStorage::dslElectroFrame) { | |
575 | < | idat.eFrame1 = atomRowData.electroFrame[atom1]; |
576 | < | idat.eFrame2 = atomColData.electroFrame[atom2]; |
575 | > | idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]); |
576 | > | idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]); |
577 | } | |
578 | ||
579 | if (storageLayout_ & DataStorage::dslTorque) { | |
580 | < | idat.t1 = atomRowData.torque[atom1]; |
581 | < | idat.t2 = atomColData.torque[atom2]; |
580 | > | idat.t1 = &(snap_->atomData.torque[atom1]); |
581 | > | idat.t2 = &(snap_->atomData.torque[atom2]); |
582 | } | |
583 | ||
584 | if (storageLayout_ & DataStorage::dslDensity) { | |
585 | < | idat.rho1 = atomRowData.density[atom1]; |
586 | < | idat.rho2 = atomColData.density[atom2]; |
585 | > | idat.rho1 = &(snap_->atomData.density[atom1]); |
586 | > | idat.rho2 = &(snap_->atomData.density[atom2]); |
587 | } | |
588 | ||
589 | if (storageLayout_ & DataStorage::dslFunctionalDerivative) { | |
590 | < | idat.dfrho1 = atomRowData.functionalDerivative[atom1]; |
591 | < | idat.dfrho2 = atomColData.functionalDerivative[atom2]; |
590 | > | idat.dfrho1 = &(snap_->atomData.functionalDerivative[atom1]); |
591 | > | idat.dfrho2 = &(snap_->atomData.functionalDerivative[atom2]); |
592 | } | |
593 | #endif | |
594 | < | |
594 | > | return idat; |
595 | } | |
596 | + | |
597 | InteractionData ForceMatrixDecomposition::fillSkipData(int atom1, int atom2){ | |
598 | + | |
599 | + | InteractionData idat; |
600 | + | #ifdef IS_MPI |
601 | + | if (storageLayout_ & DataStorage::dslElectroFrame) { |
602 | + | idat.eFrame1 = &(atomRowData.electroFrame[atom1]); |
603 | + | idat.eFrame2 = &(atomColData.electroFrame[atom2]); |
604 | + | } |
605 | + | if (storageLayout_ & DataStorage::dslTorque) { |
606 | + | idat.t1 = &(atomRowData.torque[atom1]); |
607 | + | idat.t2 = &(atomColData.torque[atom2]); |
608 | + | } |
609 | + | if (storageLayout_ & DataStorage::dslForce) { |
610 | + | idat.t1 = &(atomRowData.force[atom1]); |
611 | + | idat.t2 = &(atomColData.force[atom2]); |
612 | + | } |
613 | + | #else |
614 | + | if (storageLayout_ & DataStorage::dslElectroFrame) { |
615 | + | idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]); |
616 | + | idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]); |
617 | + | } |
618 | + | if (storageLayout_ & DataStorage::dslTorque) { |
619 | + | idat.t1 = &(snap_->atomData.torque[atom1]); |
620 | + | idat.t2 = &(snap_->atomData.torque[atom2]); |
621 | + | } |
622 | + | if (storageLayout_ & DataStorage::dslForce) { |
623 | + | idat.t1 = &(snap_->atomData.force[atom1]); |
624 | + | idat.t2 = &(snap_->atomData.force[atom2]); |
625 | + | } |
626 | + | #endif |
627 | + | |
628 | } | |
352 | – | SelfData ForceMatrixDecomposition::fillSelfData(int atom1) { |
353 | – | } |
629 | ||
630 | < | |
630 | > | |
631 | > | |
632 | > | |
633 | > | /* |
634 | > | * buildNeighborList |
635 | > | * |
636 | > | * first element of pair is row-indexed CutoffGroup |
637 | > | * second element of pair is column-indexed CutoffGroup |
638 | > | */ |
639 | > | vector<pair<int, int> > ForceMatrixDecomposition::buildNeighborList() { |
640 | > | |
641 | > | vector<pair<int, int> > neighborList; |
642 | > | #ifdef IS_MPI |
643 | > | cellListRow_.clear(); |
644 | > | cellListCol_.clear(); |
645 | > | #else |
646 | > | cellList_.clear(); |
647 | > | #endif |
648 | > | |
649 | > | // dangerous to not do error checking. |
650 | > | RealType rCut_; |
651 | > | |
652 | > | RealType rList_ = (rCut_ + skinThickness_); |
653 | > | RealType rl2 = rList_ * rList_; |
654 | > | Snapshot* snap_ = sman_->getCurrentSnapshot(); |
655 | > | Mat3x3d Hmat = snap_->getHmat(); |
656 | > | Vector3d Hx = Hmat.getColumn(0); |
657 | > | Vector3d Hy = Hmat.getColumn(1); |
658 | > | Vector3d Hz = Hmat.getColumn(2); |
659 | > | |
660 | > | nCells_.x() = (int) ( Hx.length() )/ rList_; |
661 | > | nCells_.y() = (int) ( Hy.length() )/ rList_; |
662 | > | nCells_.z() = (int) ( Hz.length() )/ rList_; |
663 | > | |
664 | > | Mat3x3d invHmat = snap_->getInvHmat(); |
665 | > | Vector3d rs, scaled, dr; |
666 | > | Vector3i whichCell; |
667 | > | int cellIndex; |
668 | > | |
669 | > | #ifdef IS_MPI |
670 | > | for (int i = 0; i < nGroupsInRow_; i++) { |
671 | > | rs = cgRowData.position[i]; |
672 | > | // scaled positions relative to the box vectors |
673 | > | scaled = invHmat * rs; |
674 | > | // wrap the vector back into the unit box by subtracting integer box |
675 | > | // numbers |
676 | > | for (int j = 0; j < 3; j++) |
677 | > | scaled[j] -= roundMe(scaled[j]); |
678 | > | |
679 | > | // find xyz-indices of cell that cutoffGroup is in. |
680 | > | whichCell.x() = nCells_.x() * scaled.x(); |
681 | > | whichCell.y() = nCells_.y() * scaled.y(); |
682 | > | whichCell.z() = nCells_.z() * scaled.z(); |
683 | > | |
684 | > | // find single index of this cell: |
685 | > | cellIndex = Vlinear(whichCell, nCells_); |
686 | > | // add this cutoff group to the list of groups in this cell; |
687 | > | cellListRow_[cellIndex].push_back(i); |
688 | > | } |
689 | > | |
690 | > | for (int i = 0; i < nGroupsInCol_; i++) { |
691 | > | rs = cgColData.position[i]; |
692 | > | // scaled positions relative to the box vectors |
693 | > | scaled = invHmat * rs; |
694 | > | // wrap the vector back into the unit box by subtracting integer box |
695 | > | // numbers |
696 | > | for (int j = 0; j < 3; j++) |
697 | > | scaled[j] -= roundMe(scaled[j]); |
698 | > | |
699 | > | // find xyz-indices of cell that cutoffGroup is in. |
700 | > | whichCell.x() = nCells_.x() * scaled.x(); |
701 | > | whichCell.y() = nCells_.y() * scaled.y(); |
702 | > | whichCell.z() = nCells_.z() * scaled.z(); |
703 | > | |
704 | > | // find single index of this cell: |
705 | > | cellIndex = Vlinear(whichCell, nCells_); |
706 | > | // add this cutoff group to the list of groups in this cell; |
707 | > | cellListCol_[cellIndex].push_back(i); |
708 | > | } |
709 | > | #else |
710 | > | for (int i = 0; i < nGroups_; i++) { |
711 | > | rs = snap_->cgData.position[i]; |
712 | > | // scaled positions relative to the box vectors |
713 | > | scaled = invHmat * rs; |
714 | > | // wrap the vector back into the unit box by subtracting integer box |
715 | > | // numbers |
716 | > | for (int j = 0; j < 3; j++) |
717 | > | scaled[j] -= roundMe(scaled[j]); |
718 | > | |
719 | > | // find xyz-indices of cell that cutoffGroup is in. |
720 | > | whichCell.x() = nCells_.x() * scaled.x(); |
721 | > | whichCell.y() = nCells_.y() * scaled.y(); |
722 | > | whichCell.z() = nCells_.z() * scaled.z(); |
723 | > | |
724 | > | // find single index of this cell: |
725 | > | cellIndex = Vlinear(whichCell, nCells_); |
726 | > | // add this cutoff group to the list of groups in this cell; |
727 | > | cellList_[cellIndex].push_back(i); |
728 | > | } |
729 | > | #endif |
730 | > | |
731 | > | |
732 | > | |
733 | > | for (int m1z = 0; m1z < nCells_.z(); m1z++) { |
734 | > | for (int m1y = 0; m1y < nCells_.y(); m1y++) { |
735 | > | for (int m1x = 0; m1x < nCells_.x(); m1x++) { |
736 | > | Vector3i m1v(m1x, m1y, m1z); |
737 | > | int m1 = Vlinear(m1v, nCells_); |
738 | > | |
739 | > | for (vector<Vector3i>::iterator os = cellOffsets_.begin(); |
740 | > | os != cellOffsets_.end(); ++os) { |
741 | > | |
742 | > | Vector3i m2v = m1v + (*os); |
743 | > | |
744 | > | if (m2v.x() >= nCells_.x()) { |
745 | > | m2v.x() = 0; |
746 | > | } else if (m2v.x() < 0) { |
747 | > | m2v.x() = nCells_.x() - 1; |
748 | > | } |
749 | > | |
750 | > | if (m2v.y() >= nCells_.y()) { |
751 | > | m2v.y() = 0; |
752 | > | } else if (m2v.y() < 0) { |
753 | > | m2v.y() = nCells_.y() - 1; |
754 | > | } |
755 | > | |
756 | > | if (m2v.z() >= nCells_.z()) { |
757 | > | m2v.z() = 0; |
758 | > | } else if (m2v.z() < 0) { |
759 | > | m2v.z() = nCells_.z() - 1; |
760 | > | } |
761 | > | |
762 | > | int m2 = Vlinear (m2v, nCells_); |
763 | > | |
764 | > | #ifdef IS_MPI |
765 | > | for (vector<int>::iterator j1 = cellListRow_[m1].begin(); |
766 | > | j1 != cellListRow_[m1].end(); ++j1) { |
767 | > | for (vector<int>::iterator j2 = cellListCol_[m2].begin(); |
768 | > | j2 != cellListCol_[m2].end(); ++j2) { |
769 | > | |
770 | > | // Always do this if we're in different cells or if |
771 | > | // we're in the same cell and the global index of the |
772 | > | // j2 cutoff group is less than the j1 cutoff group |
773 | > | |
774 | > | if (m2 != m1 || cgColToGlobal[(*j2)] < cgRowToGlobal[(*j1)]) { |
775 | > | dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; |
776 | > | snap_->wrapVector(dr); |
777 | > | if (dr.lengthSquare() < rl2) { |
778 | > | neighborList.push_back(make_pair((*j1), (*j2))); |
779 | > | } |
780 | > | } |
781 | > | } |
782 | > | } |
783 | > | #else |
784 | > | for (vector<int>::iterator j1 = cellList_[m1].begin(); |
785 | > | j1 != cellList_[m1].end(); ++j1) { |
786 | > | for (vector<int>::iterator j2 = cellList_[m2].begin(); |
787 | > | j2 != cellList_[m2].end(); ++j2) { |
788 | > | |
789 | > | // Always do this if we're in different cells or if |
790 | > | // we're in the same cell and the global index of the |
791 | > | // j2 cutoff group is less than the j1 cutoff group |
792 | > | |
793 | > | if (m2 != m1 || (*j2) < (*j1)) { |
794 | > | dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)]; |
795 | > | snap_->wrapVector(dr); |
796 | > | if (dr.lengthSquare() < rl2) { |
797 | > | neighborList.push_back(make_pair((*j1), (*j2))); |
798 | > | } |
799 | > | } |
800 | > | } |
801 | > | } |
802 | > | #endif |
803 | > | } |
804 | > | } |
805 | > | } |
806 | > | } |
807 | > | |
808 | > | // save the local cutoff group positions for the check that is |
809 | > | // done on each loop: |
810 | > | saved_CG_positions_.clear(); |
811 | > | for (int i = 0; i < nGroups_; i++) |
812 | > | saved_CG_positions_.push_back(snap_->cgData.position[i]); |
813 | > | |
814 | > | return neighborList; |
815 | > | } |
816 | } //end namespace OpenMD |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |