# | Line 47 | Line 47 | namespace OpenMD { | |
---|---|---|
47 | using namespace std; | |
48 | namespace OpenMD { | |
49 | ||
50 | + | ForceMatrixDecomposition::ForceMatrixDecomposition(SimInfo* info, InteractionManager* iMan) : ForceDecomposition(info, iMan) { |
51 | + | |
52 | + | // In a parallel computation, row and colum scans must visit all |
53 | + | // surrounding cells (not just the 14 upper triangular blocks that |
54 | + | // are used when the processor can see all pairs) |
55 | + | #ifdef IS_MPI |
56 | + | cellOffsets_.clear(); |
57 | + | cellOffsets_.push_back( Vector3i(-1,-1,-1) ); |
58 | + | cellOffsets_.push_back( Vector3i( 0,-1,-1) ); |
59 | + | cellOffsets_.push_back( Vector3i( 1,-1,-1) ); |
60 | + | cellOffsets_.push_back( Vector3i(-1, 0,-1) ); |
61 | + | cellOffsets_.push_back( Vector3i( 0, 0,-1) ); |
62 | + | cellOffsets_.push_back( Vector3i( 1, 0,-1) ); |
63 | + | cellOffsets_.push_back( Vector3i(-1, 1,-1) ); |
64 | + | cellOffsets_.push_back( Vector3i( 0, 1,-1) ); |
65 | + | cellOffsets_.push_back( Vector3i( 1, 1,-1) ); |
66 | + | cellOffsets_.push_back( Vector3i(-1,-1, 0) ); |
67 | + | cellOffsets_.push_back( Vector3i( 0,-1, 0) ); |
68 | + | cellOffsets_.push_back( Vector3i( 1,-1, 0) ); |
69 | + | cellOffsets_.push_back( Vector3i(-1, 0, 0) ); |
70 | + | cellOffsets_.push_back( Vector3i( 0, 0, 0) ); |
71 | + | cellOffsets_.push_back( Vector3i( 1, 0, 0) ); |
72 | + | cellOffsets_.push_back( Vector3i(-1, 1, 0) ); |
73 | + | cellOffsets_.push_back( Vector3i( 0, 1, 0) ); |
74 | + | cellOffsets_.push_back( Vector3i( 1, 1, 0) ); |
75 | + | cellOffsets_.push_back( Vector3i(-1,-1, 1) ); |
76 | + | cellOffsets_.push_back( Vector3i( 0,-1, 1) ); |
77 | + | cellOffsets_.push_back( Vector3i( 1,-1, 1) ); |
78 | + | cellOffsets_.push_back( Vector3i(-1, 0, 1) ); |
79 | + | cellOffsets_.push_back( Vector3i( 0, 0, 1) ); |
80 | + | cellOffsets_.push_back( Vector3i( 1, 0, 1) ); |
81 | + | cellOffsets_.push_back( Vector3i(-1, 1, 1) ); |
82 | + | cellOffsets_.push_back( Vector3i( 0, 1, 1) ); |
83 | + | cellOffsets_.push_back( Vector3i( 1, 1, 1) ); |
84 | + | #endif |
85 | + | } |
86 | + | |
87 | + | |
88 | /** | |
89 | * distributeInitialData is essentially a copy of the older fortran | |
90 | * SimulationSetup | |
91 | */ | |
54 | – | |
92 | void ForceMatrixDecomposition::distributeInitialData() { | |
93 | snap_ = sman_->getCurrentSnapshot(); | |
94 | storageLayout_ = sman_->getStorageLayout(); | |
# | Line 74 | Line 111 | namespace OpenMD { | |
111 | ||
112 | #ifdef IS_MPI | |
113 | ||
114 | < | AtomCommIntRow = new Communicator<Row,int>(nLocal_); |
115 | < | AtomCommRealRow = new Communicator<Row,RealType>(nLocal_); |
79 | < | AtomCommVectorRow = new Communicator<Row,Vector3d>(nLocal_); |
80 | < | AtomCommMatrixRow = new Communicator<Row,Mat3x3d>(nLocal_); |
81 | < | AtomCommPotRow = new Communicator<Row,potVec>(nLocal_); |
114 | > | MPI::Intracomm row = rowComm.getComm(); |
115 | > | MPI::Intracomm col = colComm.getComm(); |
116 | ||
117 | < | AtomCommIntColumn = new Communicator<Column,int>(nLocal_); |
118 | < | AtomCommRealColumn = new Communicator<Column,RealType>(nLocal_); |
119 | < | AtomCommVectorColumn = new Communicator<Column,Vector3d>(nLocal_); |
120 | < | AtomCommMatrixColumn = new Communicator<Column,Mat3x3d>(nLocal_); |
121 | < | AtomCommPotColumn = new Communicator<Column,potVec>(nLocal_); |
117 | > | AtomPlanIntRow = new Plan<int>(row, nLocal_); |
118 | > | AtomPlanRealRow = new Plan<RealType>(row, nLocal_); |
119 | > | AtomPlanVectorRow = new Plan<Vector3d>(row, nLocal_); |
120 | > | AtomPlanMatrixRow = new Plan<Mat3x3d>(row, nLocal_); |
121 | > | AtomPlanPotRow = new Plan<potVec>(row, nLocal_); |
122 | ||
123 | < | cgCommIntRow = new Communicator<Row,int>(nGroups_); |
124 | < | cgCommVectorRow = new Communicator<Row,Vector3d>(nGroups_); |
125 | < | cgCommIntColumn = new Communicator<Column,int>(nGroups_); |
126 | < | cgCommVectorColumn = new Communicator<Column,Vector3d>(nGroups_); |
123 | > | AtomPlanIntColumn = new Plan<int>(col, nLocal_); |
124 | > | AtomPlanRealColumn = new Plan<RealType>(col, nLocal_); |
125 | > | AtomPlanVectorColumn = new Plan<Vector3d>(col, nLocal_); |
126 | > | AtomPlanMatrixColumn = new Plan<Mat3x3d>(col, nLocal_); |
127 | > | AtomPlanPotColumn = new Plan<potVec>(col, nLocal_); |
128 | ||
129 | < | nAtomsInRow_ = AtomCommIntRow->getSize(); |
130 | < | nAtomsInCol_ = AtomCommIntColumn->getSize(); |
131 | < | nGroupsInRow_ = cgCommIntRow->getSize(); |
132 | < | nGroupsInCol_ = cgCommIntColumn->getSize(); |
129 | > | cgPlanIntRow = new Plan<int>(row, nGroups_); |
130 | > | cgPlanVectorRow = new Plan<Vector3d>(row, nGroups_); |
131 | > | cgPlanIntColumn = new Plan<int>(col, nGroups_); |
132 | > | cgPlanVectorColumn = new Plan<Vector3d>(col, nGroups_); |
133 | ||
134 | + | nAtomsInRow_ = AtomPlanIntRow->getSize(); |
135 | + | nAtomsInCol_ = AtomPlanIntColumn->getSize(); |
136 | + | nGroupsInRow_ = cgPlanIntRow->getSize(); |
137 | + | nGroupsInCol_ = cgPlanIntColumn->getSize(); |
138 | + | |
139 | // Modify the data storage objects with the correct layouts and sizes: | |
140 | atomRowData.resize(nAtomsInRow_); | |
141 | atomRowData.setStorageLayout(storageLayout_); | |
# | Line 109 | Line 149 | namespace OpenMD { | |
149 | identsRow.resize(nAtomsInRow_); | |
150 | identsCol.resize(nAtomsInCol_); | |
151 | ||
152 | < | AtomCommIntRow->gather(idents, identsRow); |
153 | < | AtomCommIntColumn->gather(idents, identsCol); |
152 | > | AtomPlanIntRow->gather(idents, identsRow); |
153 | > | AtomPlanIntColumn->gather(idents, identsCol); |
154 | ||
155 | // allocate memory for the parallel objects | |
156 | atypesRow.resize(nAtomsInRow_); | |
# | Line 126 | Line 166 | namespace OpenMD { | |
166 | ||
167 | AtomRowToGlobal.resize(nAtomsInRow_); | |
168 | AtomColToGlobal.resize(nAtomsInCol_); | |
169 | < | AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal); |
170 | < | AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal); |
171 | < | |
169 | > | AtomPlanIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal); |
170 | > | AtomPlanIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal); |
171 | > | |
172 | cgRowToGlobal.resize(nGroupsInRow_); | |
173 | cgColToGlobal.resize(nGroupsInCol_); | |
174 | < | cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal); |
175 | < | cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal); |
174 | > | cgPlanIntRow->gather(cgLocalToGlobal, cgRowToGlobal); |
175 | > | cgPlanIntColumn->gather(cgLocalToGlobal, cgColToGlobal); |
176 | ||
177 | massFactorsRow.resize(nAtomsInRow_); | |
178 | massFactorsCol.resize(nAtomsInCol_); | |
179 | < | AtomCommRealRow->gather(massFactors, massFactorsRow); |
180 | < | AtomCommRealColumn->gather(massFactors, massFactorsCol); |
179 | > | AtomPlanRealRow->gather(massFactors, massFactorsRow); |
180 | > | AtomPlanRealColumn->gather(massFactors, massFactorsCol); |
181 | ||
182 | groupListRow_.clear(); | |
183 | groupListRow_.resize(nGroupsInRow_); | |
# | Line 193 | Line 233 | namespace OpenMD { | |
233 | } | |
234 | } | |
235 | ||
236 | < | #endif |
197 | < | |
198 | < | // allocate memory for the parallel objects |
199 | < | atypesLocal.resize(nLocal_); |
200 | < | |
201 | < | for (int i = 0; i < nLocal_; i++) |
202 | < | atypesLocal[i] = ff_->getAtomType(idents[i]); |
203 | < | |
204 | < | groupList_.clear(); |
205 | < | groupList_.resize(nGroups_); |
206 | < | for (int i = 0; i < nGroups_; i++) { |
207 | < | int gid = cgLocalToGlobal[i]; |
208 | < | for (int j = 0; j < nLocal_; j++) { |
209 | < | int aid = AtomLocalToGlobal[j]; |
210 | < | if (globalGroupMembership[aid] == gid) { |
211 | < | groupList_[i].push_back(j); |
212 | < | } |
213 | < | } |
214 | < | } |
215 | < | |
236 | > | #else |
237 | excludesForAtom.clear(); | |
238 | excludesForAtom.resize(nLocal_); | |
239 | toposForAtom.clear(); | |
# | Line 226 | Line 247 | namespace OpenMD { | |
247 | for (int j = 0; j < nLocal_; j++) { | |
248 | int jglob = AtomLocalToGlobal[j]; | |
249 | ||
250 | < | if (excludes->hasPair(iglob, jglob)) |
250 | > | if (excludes->hasPair(iglob, jglob)) |
251 | excludesForAtom[i].push_back(j); | |
252 | ||
253 | + | |
254 | if (oneTwo->hasPair(iglob, jglob)) { | |
255 | toposForAtom[i].push_back(j); | |
256 | topoDist[i].push_back(1); | |
# | Line 245 | Line 267 | namespace OpenMD { | |
267 | } | |
268 | } | |
269 | } | |
270 | < | |
270 | > | #endif |
271 | > | |
272 | > | // allocate memory for the parallel objects |
273 | > | atypesLocal.resize(nLocal_); |
274 | > | |
275 | > | for (int i = 0; i < nLocal_; i++) |
276 | > | atypesLocal[i] = ff_->getAtomType(idents[i]); |
277 | > | |
278 | > | groupList_.clear(); |
279 | > | groupList_.resize(nGroups_); |
280 | > | for (int i = 0; i < nGroups_; i++) { |
281 | > | int gid = cgLocalToGlobal[i]; |
282 | > | for (int j = 0; j < nLocal_; j++) { |
283 | > | int aid = AtomLocalToGlobal[j]; |
284 | > | if (globalGroupMembership[aid] == gid) { |
285 | > | groupList_[i].push_back(j); |
286 | > | } |
287 | > | } |
288 | > | } |
289 | > | |
290 | > | |
291 | createGtypeCutoffMap(); | |
292 | ||
293 | } | |
# | Line 253 | Line 295 | namespace OpenMD { | |
295 | void ForceMatrixDecomposition::createGtypeCutoffMap() { | |
296 | ||
297 | RealType tol = 1e-6; | |
298 | + | largestRcut_ = 0.0; |
299 | RealType rc; | |
300 | int atid; | |
301 | set<AtomType*> atypes = info_->getSimulatedAtomTypes(); | |
302 | + | |
303 | map<int, RealType> atypeCutoff; | |
304 | ||
305 | for (set<AtomType*>::iterator at = atypes.begin(); | |
# | Line 263 | Line 307 | namespace OpenMD { | |
307 | atid = (*at)->getIdent(); | |
308 | if (userChoseCutoff_) | |
309 | atypeCutoff[atid] = userCutoff_; | |
310 | < | else |
310 | > | else |
311 | atypeCutoff[atid] = interactionMan_->getSuggestedCutoffRadius(*at); | |
312 | } | |
313 | < | |
313 | > | |
314 | vector<RealType> gTypeCutoffs; | |
315 | // first we do a single loop over the cutoff groups to find the | |
316 | // largest cutoff for any atypes present in this group. | |
# | Line 326 | Line 370 | namespace OpenMD { | |
370 | vector<RealType> groupCutoff(nGroups_, 0.0); | |
371 | groupToGtype.resize(nGroups_); | |
372 | for (int cg1 = 0; cg1 < nGroups_; cg1++) { | |
329 | – | |
373 | groupCutoff[cg1] = 0.0; | |
374 | vector<int> atomList = getAtomsInGroupRow(cg1); | |
332 | – | |
375 | for (vector<int>::iterator ia = atomList.begin(); | |
376 | ia != atomList.end(); ++ia) { | |
377 | int atom1 = (*ia); | |
378 | atid = idents[atom1]; | |
379 | < | if (atypeCutoff[atid] > groupCutoff[cg1]) { |
379 | > | if (atypeCutoff[atid] > groupCutoff[cg1]) |
380 | groupCutoff[cg1] = atypeCutoff[atid]; | |
339 | – | } |
381 | } | |
382 | < | |
382 | > | |
383 | bool gTypeFound = false; | |
384 | for (int gt = 0; gt < gTypeCutoffs.size(); gt++) { | |
385 | if (abs(groupCutoff[cg1] - gTypeCutoffs[gt]) < tol) { | |
# | Line 346 | Line 387 | namespace OpenMD { | |
387 | gTypeFound = true; | |
388 | } | |
389 | } | |
390 | < | if (!gTypeFound) { |
390 | > | if (!gTypeFound) { |
391 | gTypeCutoffs.push_back( groupCutoff[cg1] ); | |
392 | groupToGtype[cg1] = gTypeCutoffs.size() - 1; | |
393 | } | |
# | Line 390 | Line 431 | namespace OpenMD { | |
431 | ||
432 | pair<int,int> key = make_pair(i,j); | |
433 | gTypeCutoffMap[key].first = thisRcut; | |
393 | – | |
434 | if (thisRcut > largestRcut_) largestRcut_ = thisRcut; | |
395 | – | |
435 | gTypeCutoffMap[key].second = thisRcut*thisRcut; | |
397 | – | |
436 | gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2); | |
399 | – | |
437 | // sanity check | |
438 | ||
439 | if (userChoseCutoff_) { | |
# | Line 522 | Line 559 | namespace OpenMD { | |
559 | #ifdef IS_MPI | |
560 | ||
561 | // gather up the atomic positions | |
562 | < | AtomCommVectorRow->gather(snap_->atomData.position, |
562 | > | AtomPlanVectorRow->gather(snap_->atomData.position, |
563 | atomRowData.position); | |
564 | < | AtomCommVectorColumn->gather(snap_->atomData.position, |
564 | > | AtomPlanVectorColumn->gather(snap_->atomData.position, |
565 | atomColData.position); | |
566 | ||
567 | // gather up the cutoff group positions | |
568 | < | cgCommVectorRow->gather(snap_->cgData.position, |
568 | > | |
569 | > | cgPlanVectorRow->gather(snap_->cgData.position, |
570 | cgRowData.position); | |
571 | < | cgCommVectorColumn->gather(snap_->cgData.position, |
571 | > | |
572 | > | cgPlanVectorColumn->gather(snap_->cgData.position, |
573 | cgColData.position); | |
574 | + | |
575 | ||
576 | // if needed, gather the atomic rotation matrices | |
577 | if (storageLayout_ & DataStorage::dslAmat) { | |
578 | < | AtomCommMatrixRow->gather(snap_->atomData.aMat, |
578 | > | AtomPlanMatrixRow->gather(snap_->atomData.aMat, |
579 | atomRowData.aMat); | |
580 | < | AtomCommMatrixColumn->gather(snap_->atomData.aMat, |
580 | > | AtomPlanMatrixColumn->gather(snap_->atomData.aMat, |
581 | atomColData.aMat); | |
582 | } | |
583 | ||
584 | // if needed, gather the atomic eletrostatic frames | |
585 | if (storageLayout_ & DataStorage::dslElectroFrame) { | |
586 | < | AtomCommMatrixRow->gather(snap_->atomData.electroFrame, |
586 | > | AtomPlanMatrixRow->gather(snap_->atomData.electroFrame, |
587 | atomRowData.electroFrame); | |
588 | < | AtomCommMatrixColumn->gather(snap_->atomData.electroFrame, |
588 | > | AtomPlanMatrixColumn->gather(snap_->atomData.electroFrame, |
589 | atomColData.electroFrame); | |
590 | } | |
591 | ||
# | Line 562 | Line 602 | namespace OpenMD { | |
602 | ||
603 | if (storageLayout_ & DataStorage::dslDensity) { | |
604 | ||
605 | < | AtomCommRealRow->scatter(atomRowData.density, |
605 | > | AtomPlanRealRow->scatter(atomRowData.density, |
606 | snap_->atomData.density); | |
607 | ||
608 | int n = snap_->atomData.density.size(); | |
609 | vector<RealType> rho_tmp(n, 0.0); | |
610 | < | AtomCommRealColumn->scatter(atomColData.density, rho_tmp); |
610 | > | AtomPlanRealColumn->scatter(atomColData.density, rho_tmp); |
611 | for (int i = 0; i < n; i++) | |
612 | snap_->atomData.density[i] += rho_tmp[i]; | |
613 | } | |
# | Line 583 | Line 623 | namespace OpenMD { | |
623 | storageLayout_ = sman_->getStorageLayout(); | |
624 | #ifdef IS_MPI | |
625 | if (storageLayout_ & DataStorage::dslFunctional) { | |
626 | < | AtomCommRealRow->gather(snap_->atomData.functional, |
626 | > | AtomPlanRealRow->gather(snap_->atomData.functional, |
627 | atomRowData.functional); | |
628 | < | AtomCommRealColumn->gather(snap_->atomData.functional, |
628 | > | AtomPlanRealColumn->gather(snap_->atomData.functional, |
629 | atomColData.functional); | |
630 | } | |
631 | ||
632 | if (storageLayout_ & DataStorage::dslFunctionalDerivative) { | |
633 | < | AtomCommRealRow->gather(snap_->atomData.functionalDerivative, |
633 | > | AtomPlanRealRow->gather(snap_->atomData.functionalDerivative, |
634 | atomRowData.functionalDerivative); | |
635 | < | AtomCommRealColumn->gather(snap_->atomData.functionalDerivative, |
635 | > | AtomPlanRealColumn->gather(snap_->atomData.functionalDerivative, |
636 | atomColData.functionalDerivative); | |
637 | } | |
638 | #endif | |
# | Line 606 | Line 646 | namespace OpenMD { | |
646 | int n = snap_->atomData.force.size(); | |
647 | vector<Vector3d> frc_tmp(n, V3Zero); | |
648 | ||
649 | < | AtomCommVectorRow->scatter(atomRowData.force, frc_tmp); |
649 | > | AtomPlanVectorRow->scatter(atomRowData.force, frc_tmp); |
650 | for (int i = 0; i < n; i++) { | |
651 | snap_->atomData.force[i] += frc_tmp[i]; | |
652 | frc_tmp[i] = 0.0; | |
653 | } | |
654 | ||
655 | < | AtomCommVectorColumn->scatter(atomColData.force, frc_tmp); |
656 | < | for (int i = 0; i < n; i++) |
655 | > | AtomPlanVectorColumn->scatter(atomColData.force, frc_tmp); |
656 | > | for (int i = 0; i < n; i++) { |
657 | snap_->atomData.force[i] += frc_tmp[i]; | |
658 | + | } |
659 | ||
660 | if (storageLayout_ & DataStorage::dslTorque) { | |
661 | ||
662 | int nt = snap_->atomData.torque.size(); | |
663 | vector<Vector3d> trq_tmp(nt, V3Zero); | |
664 | ||
665 | < | AtomCommVectorRow->scatter(atomRowData.torque, trq_tmp); |
665 | > | AtomPlanVectorRow->scatter(atomRowData.torque, trq_tmp); |
666 | for (int i = 0; i < nt; i++) { | |
667 | snap_->atomData.torque[i] += trq_tmp[i]; | |
668 | trq_tmp[i] = 0.0; | |
669 | } | |
670 | ||
671 | < | AtomCommVectorColumn->scatter(atomColData.torque, trq_tmp); |
671 | > | AtomPlanVectorColumn->scatter(atomColData.torque, trq_tmp); |
672 | for (int i = 0; i < nt; i++) | |
673 | snap_->atomData.torque[i] += trq_tmp[i]; | |
674 | } | |
# | Line 637 | Line 678 | namespace OpenMD { | |
678 | int ns = snap_->atomData.skippedCharge.size(); | |
679 | vector<RealType> skch_tmp(ns, 0.0); | |
680 | ||
681 | < | AtomCommRealRow->scatter(atomRowData.skippedCharge, skch_tmp); |
681 | > | AtomPlanRealRow->scatter(atomRowData.skippedCharge, skch_tmp); |
682 | for (int i = 0; i < ns; i++) { | |
683 | snap_->atomData.skippedCharge[i] += skch_tmp[i]; | |
684 | skch_tmp[i] = 0.0; | |
685 | } | |
686 | ||
687 | < | AtomCommRealColumn->scatter(atomColData.skippedCharge, skch_tmp); |
688 | < | for (int i = 0; i < ns; i++) |
687 | > | AtomPlanRealColumn->scatter(atomColData.skippedCharge, skch_tmp); |
688 | > | for (int i = 0; i < ns; i++) |
689 | snap_->atomData.skippedCharge[i] += skch_tmp[i]; | |
690 | + | |
691 | } | |
692 | ||
693 | nLocal_ = snap_->getNumberOfAtoms(); | |
# | Line 655 | Line 697 | namespace OpenMD { | |
697 | ||
698 | // scatter/gather pot_row into the members of my column | |
699 | ||
700 | < | AtomCommPotRow->scatter(pot_row, pot_temp); |
700 | > | AtomPlanPotRow->scatter(pot_row, pot_temp); |
701 | ||
702 | for (int ii = 0; ii < pot_temp.size(); ii++ ) | |
703 | pairwisePot += pot_temp[ii]; | |
# | Line 663 | Line 705 | namespace OpenMD { | |
705 | fill(pot_temp.begin(), pot_temp.end(), | |
706 | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | |
707 | ||
708 | < | AtomCommPotColumn->scatter(pot_col, pot_temp); |
708 | > | AtomPlanPotColumn->scatter(pot_col, pot_temp); |
709 | ||
710 | for (int ii = 0; ii < pot_temp.size(); ii++ ) | |
711 | pairwisePot += pot_temp[ii]; | |
712 | + | |
713 | + | for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { |
714 | + | RealType ploc1 = pairwisePot[ii]; |
715 | + | RealType ploc2 = 0.0; |
716 | + | MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); |
717 | + | pairwisePot[ii] = ploc2; |
718 | + | } |
719 | + | |
720 | + | for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { |
721 | + | RealType ploc1 = embeddingPot[ii]; |
722 | + | RealType ploc2 = 0.0; |
723 | + | MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); |
724 | + | embeddingPot[ii] = ploc2; |
725 | + | } |
726 | + | |
727 | #endif | |
728 | ||
729 | } | |
# | Line 779 | Line 836 | namespace OpenMD { | |
836 | */ | |
837 | bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) { | |
838 | int unique_id_1, unique_id_2; | |
839 | < | |
839 | > | |
840 | #ifdef IS_MPI | |
841 | // in MPI, we have to look up the unique IDs for each atom | |
842 | unique_id_1 = AtomRowToGlobal[atom1]; | |
# | Line 808 | Line 865 | namespace OpenMD { | |
865 | * field) must still be handled for these pairs. | |
866 | */ | |
867 | bool ForceMatrixDecomposition::excludeAtomPair(int atom1, int atom2) { | |
868 | < | int unique_id_2; |
868 | > | |
869 | > | // excludesForAtom was constructed to use row/column indices in the MPI |
870 | > | // version, and to use local IDs in the non-MPI version: |
871 | ||
813 | – | #ifdef IS_MPI |
814 | – | // in MPI, we have to look up the unique IDs for the row atom. |
815 | – | unique_id_2 = AtomColToGlobal[atom2]; |
816 | – | #else |
817 | – | // in the normal loop, the atom numbers are unique |
818 | – | unique_id_2 = atom2; |
819 | – | #endif |
820 | – | |
872 | for (vector<int>::iterator i = excludesForAtom[atom1].begin(); | |
873 | i != excludesForAtom[atom1].end(); ++i) { | |
874 | < | if ( (*i) == unique_id_2 ) return true; |
874 | > | if ( (*i) == atom2 ) return true; |
875 | } | |
876 | ||
877 | return false; | |
# | Line 1036 | Line 1087 | namespace OpenMD { | |
1087 | // add this cutoff group to the list of groups in this cell; | |
1088 | cellListRow_[cellIndex].push_back(i); | |
1089 | } | |
1039 | – | |
1090 | for (int i = 0; i < nGroupsInCol_; i++) { | |
1091 | rs = cgColData.position[i]; | |
1092 | ||
# | Line 1061 | Line 1111 | namespace OpenMD { | |
1111 | // add this cutoff group to the list of groups in this cell; | |
1112 | cellListCol_[cellIndex].push_back(i); | |
1113 | } | |
1114 | + | |
1115 | #else | |
1116 | for (int i = 0; i < nGroups_; i++) { | |
1117 | rs = snap_->cgData.position[i]; | |
# | Line 1081 | Line 1132 | namespace OpenMD { | |
1132 | whichCell.z() = nCells_.z() * scaled.z(); | |
1133 | ||
1134 | // find single index of this cell: | |
1135 | < | cellIndex = Vlinear(whichCell, nCells_); |
1135 | > | cellIndex = Vlinear(whichCell, nCells_); |
1136 | ||
1137 | // add this cutoff group to the list of groups in this cell; | |
1138 | cellList_[cellIndex].push_back(i); | |
1139 | } | |
1140 | + | |
1141 | #endif | |
1142 | ||
1143 | for (int m1z = 0; m1z < nCells_.z(); m1z++) { | |
# | Line 1098 | Line 1150 | namespace OpenMD { | |
1150 | os != cellOffsets_.end(); ++os) { | |
1151 | ||
1152 | Vector3i m2v = m1v + (*os); | |
1153 | < | |
1153 | > | |
1154 | > | |
1155 | if (m2v.x() >= nCells_.x()) { | |
1156 | m2v.x() = 0; | |
1157 | } else if (m2v.x() < 0) { | |
# | Line 1116 | Line 1169 | namespace OpenMD { | |
1169 | } else if (m2v.z() < 0) { | |
1170 | m2v.z() = nCells_.z() - 1; | |
1171 | } | |
1172 | < | |
1172 | > | |
1173 | int m2 = Vlinear (m2v, nCells_); | |
1174 | ||
1175 | #ifdef IS_MPI | |
# | Line 1125 | Line 1178 | namespace OpenMD { | |
1178 | for (vector<int>::iterator j2 = cellListCol_[m2].begin(); | |
1179 | j2 != cellListCol_[m2].end(); ++j2) { | |
1180 | ||
1181 | < | // Always do this if we're in different cells or if |
1182 | < | // we're in the same cell and the global index of the |
1183 | < | // j2 cutoff group is less than the j1 cutoff group |
1184 | < | |
1185 | < | if (m2 != m1 || cgColToGlobal[(*j2)] < cgRowToGlobal[(*j1)]) { |
1186 | < | dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; |
1187 | < | snap_->wrapVector(dr); |
1188 | < | cuts = getGroupCutoffs( (*j1), (*j2) ); |
1189 | < | if (dr.lengthSquare() < cuts.third) { |
1137 | < | neighborList.push_back(make_pair((*j1), (*j2))); |
1138 | < | } |
1139 | < | } |
1181 | > | // In parallel, we need to visit *all* pairs of row |
1182 | > | // & column indicies and will divide labor in the |
1183 | > | // force evaluation later. |
1184 | > | dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; |
1185 | > | snap_->wrapVector(dr); |
1186 | > | cuts = getGroupCutoffs( (*j1), (*j2) ); |
1187 | > | if (dr.lengthSquare() < cuts.third) { |
1188 | > | neighborList.push_back(make_pair((*j1), (*j2))); |
1189 | > | } |
1190 | } | |
1191 | } | |
1192 | #else |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |