# | Line 57 | Line 57 | namespace OpenMD { | |
---|---|---|
57 | storageLayout_ = sman_->getStorageLayout(); | |
58 | ff_ = info_->getForceField(); | |
59 | nLocal_ = snap_->getNumberOfAtoms(); | |
60 | < | |
60 | > | |
61 | nGroups_ = info_->getNLocalCutoffGroups(); | |
62 | – | cerr << "in dId, nGroups = " << nGroups_ << "\n"; |
62 | // gather the information for atomtype IDs (atids): | |
63 | < | identsLocal = info_->getIdentArray(); |
63 | > | idents = info_->getIdentArray(); |
64 | AtomLocalToGlobal = info_->getGlobalAtomIndices(); | |
65 | cgLocalToGlobal = info_->getGlobalGroupIndices(); | |
66 | vector<int> globalGroupMembership = info_->getGlobalGroupMembership(); | |
68 | – | vector<RealType> massFactorsLocal = info_->getMassFactors(); |
69 | – | PairList excludes = info_->getExcludedInteractions(); |
70 | – | PairList oneTwo = info_->getOneTwoInteractions(); |
71 | – | PairList oneThree = info_->getOneThreeInteractions(); |
72 | – | PairList oneFour = info_->getOneFourInteractions(); |
67 | ||
68 | + | massFactors = info_->getMassFactors(); |
69 | + | |
70 | + | PairList* excludes = info_->getExcludedInteractions(); |
71 | + | PairList* oneTwo = info_->getOneTwoInteractions(); |
72 | + | PairList* oneThree = info_->getOneThreeInteractions(); |
73 | + | PairList* oneFour = info_->getOneFourInteractions(); |
74 | + | |
75 | #ifdef IS_MPI | |
76 | ||
77 | AtomCommIntRow = new Communicator<Row,int>(nLocal_); | |
# | Line 108 | Line 109 | namespace OpenMD { | |
109 | identsRow.resize(nAtomsInRow_); | |
110 | identsCol.resize(nAtomsInCol_); | |
111 | ||
112 | < | AtomCommIntRow->gather(identsLocal, identsRow); |
113 | < | AtomCommIntColumn->gather(identsLocal, identsCol); |
112 | > | AtomCommIntRow->gather(idents, identsRow); |
113 | > | AtomCommIntColumn->gather(idents, identsCol); |
114 | ||
115 | + | // allocate memory for the parallel objects |
116 | + | AtomRowToGlobal.resize(nAtomsInRow_); |
117 | + | AtomColToGlobal.resize(nAtomsInCol_); |
118 | + | cgRowToGlobal.resize(nGroupsInRow_); |
119 | + | cgColToGlobal.resize(nGroupsInCol_); |
120 | + | massFactorsRow.resize(nAtomsInRow_); |
121 | + | massFactorsCol.resize(nAtomsInCol_); |
122 | + | pot_row.resize(nAtomsInRow_); |
123 | + | pot_col.resize(nAtomsInCol_); |
124 | + | |
125 | AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal); | |
126 | AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal); | |
127 | ||
128 | cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal); | |
129 | cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal); | |
130 | ||
131 | < | AtomCommRealRow->gather(massFactorsLocal, massFactorsRow); |
132 | < | AtomCommRealColumn->gather(massFactorsLocal, massFactorsCol); |
131 | > | AtomCommRealRow->gather(massFactors, massFactorsRow); |
132 | > | AtomCommRealColumn->gather(massFactors, massFactorsCol); |
133 | ||
134 | groupListRow_.clear(); | |
135 | groupListRow_.resize(nGroupsInRow_); | |
# | Line 142 | Line 153 | namespace OpenMD { | |
153 | } | |
154 | } | |
155 | ||
156 | < | skipsForAtom.clear(); |
157 | < | skipsForAtom.resize(nAtomsInRow_); |
156 | > | excludesForAtom.clear(); |
157 | > | excludesForAtom.resize(nAtomsInRow_); |
158 | toposForAtom.clear(); | |
159 | toposForAtom.resize(nAtomsInRow_); | |
160 | topoDist.clear(); | |
# | Line 154 | Line 165 | namespace OpenMD { | |
165 | for (int j = 0; j < nAtomsInCol_; j++) { | |
166 | int jglob = AtomColToGlobal[j]; | |
167 | ||
168 | < | if (excludes.hasPair(iglob, jglob)) |
169 | < | skipsForAtom[i].push_back(j); |
168 | > | if (excludes->hasPair(iglob, jglob)) |
169 | > | excludesForAtom[i].push_back(j); |
170 | ||
171 | < | if (oneTwo.hasPair(iglob, jglob)) { |
171 | > | if (oneTwo->hasPair(iglob, jglob)) { |
172 | toposForAtom[i].push_back(j); | |
173 | topoDist[i].push_back(1); | |
174 | } else { | |
175 | < | if (oneThree.hasPair(iglob, jglob)) { |
175 | > | if (oneThree->hasPair(iglob, jglob)) { |
176 | toposForAtom[i].push_back(j); | |
177 | topoDist[i].push_back(2); | |
178 | } else { | |
179 | < | if (oneFour.hasPair(iglob, jglob)) { |
179 | > | if (oneFour->hasPair(iglob, jglob)) { |
180 | toposForAtom[i].push_back(j); | |
181 | topoDist[i].push_back(3); | |
182 | } | |
# | Line 188 | Line 199 | namespace OpenMD { | |
199 | } | |
200 | } | |
201 | ||
202 | < | skipsForAtom.clear(); |
203 | < | skipsForAtom.resize(nLocal_); |
202 | > | excludesForAtom.clear(); |
203 | > | excludesForAtom.resize(nLocal_); |
204 | toposForAtom.clear(); | |
205 | toposForAtom.resize(nLocal_); | |
206 | topoDist.clear(); | |
# | Line 201 | Line 212 | namespace OpenMD { | |
212 | for (int j = 0; j < nLocal_; j++) { | |
213 | int jglob = AtomLocalToGlobal[j]; | |
214 | ||
215 | < | if (excludes.hasPair(iglob, jglob)) |
216 | < | skipsForAtom[i].push_back(j); |
215 | > | if (excludes->hasPair(iglob, jglob)) |
216 | > | excludesForAtom[i].push_back(j); |
217 | ||
218 | < | if (oneTwo.hasPair(iglob, jglob)) { |
218 | > | if (oneTwo->hasPair(iglob, jglob)) { |
219 | toposForAtom[i].push_back(j); | |
220 | topoDist[i].push_back(1); | |
221 | } else { | |
222 | < | if (oneThree.hasPair(iglob, jglob)) { |
222 | > | if (oneThree->hasPair(iglob, jglob)) { |
223 | toposForAtom[i].push_back(j); | |
224 | topoDist[i].push_back(2); | |
225 | } else { | |
226 | < | if (oneFour.hasPair(iglob, jglob)) { |
226 | > | if (oneFour->hasPair(iglob, jglob)) { |
227 | toposForAtom[i].push_back(j); | |
228 | topoDist[i].push_back(3); | |
229 | } | |
# | Line 222 | Line 233 | namespace OpenMD { | |
233 | } | |
234 | ||
235 | createGtypeCutoffMap(); | |
236 | + | |
237 | } | |
238 | ||
239 | void ForceMatrixDecomposition::createGtypeCutoffMap() { | |
240 | < | |
240 | > | |
241 | RealType tol = 1e-6; | |
242 | RealType rc; | |
243 | int atid; | |
244 | set<AtomType*> atypes = info_->getSimulatedAtomTypes(); | |
245 | < | vector<RealType> atypeCutoff; |
246 | < | atypeCutoff.resize( atypes.size() ); |
235 | < | |
245 | > | map<int, RealType> atypeCutoff; |
246 | > | |
247 | for (set<AtomType*>::iterator at = atypes.begin(); | |
248 | at != atypes.end(); ++at){ | |
238 | – | rc = interactionMan_->getSuggestedCutoffRadius(*at); |
249 | atid = (*at)->getIdent(); | |
250 | < | atypeCutoff[atid] = rc; |
250 | > | if (userChoseCutoff_) |
251 | > | atypeCutoff[atid] = userCutoff_; |
252 | > | else |
253 | > | atypeCutoff[atid] = interactionMan_->getSuggestedCutoffRadius(*at); |
254 | } | |
255 | ||
256 | vector<RealType> gTypeCutoffs; | |
244 | – | |
257 | // first we do a single loop over the cutoff groups to find the | |
258 | // largest cutoff for any atypes present in this group. | |
259 | #ifdef IS_MPI | |
# | Line 299 | Line 311 | namespace OpenMD { | |
311 | ||
312 | vector<RealType> groupCutoff(nGroups_, 0.0); | |
313 | groupToGtype.resize(nGroups_); | |
302 | – | |
303 | – | cerr << "nGroups = " << nGroups_ << "\n"; |
314 | for (int cg1 = 0; cg1 < nGroups_; cg1++) { | |
315 | ||
316 | groupCutoff[cg1] = 0.0; | |
# | Line 309 | Line 319 | namespace OpenMD { | |
319 | for (vector<int>::iterator ia = atomList.begin(); | |
320 | ia != atomList.end(); ++ia) { | |
321 | int atom1 = (*ia); | |
322 | < | atid = identsLocal[atom1]; |
322 | > | atid = idents[atom1]; |
323 | if (atypeCutoff[atid] > groupCutoff[cg1]) { | |
324 | groupCutoff[cg1] = atypeCutoff[atid]; | |
325 | } | |
# | Line 329 | Line 339 | namespace OpenMD { | |
339 | } | |
340 | #endif | |
341 | ||
332 | – | cerr << "gTypeCutoffs.size() = " << gTypeCutoffs.size() << "\n"; |
342 | // Now we find the maximum group cutoff value present in the simulation | |
343 | ||
344 | < | RealType groupMax = *max_element(gTypeCutoffs.begin(), gTypeCutoffs.end()); |
344 | > | RealType groupMax = *max_element(gTypeCutoffs.begin(), |
345 | > | gTypeCutoffs.end()); |
346 | ||
347 | #ifdef IS_MPI | |
348 | < | MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE, MPI::MAX); |
348 | > | MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE, |
349 | > | MPI::MAX); |
350 | #endif | |
351 | ||
352 | RealType tradRcut = groupMax; | |
# | Line 378 | Line 389 | namespace OpenMD { | |
389 | if (abs(gTypeCutoffMap[key].first - userCutoff_) > 0.0001) { | |
390 | sprintf(painCave.errMsg, | |
391 | "ForceMatrixDecomposition::createGtypeCutoffMap " | |
392 | < | "user-specified rCut does not match computed group Cutoff\n"); |
392 | > | "user-specified rCut (%lf) does not match computed group Cutoff\n", userCutoff_); |
393 | painCave.severity = OPENMD_ERROR; | |
394 | painCave.isFatal = 1; | |
395 | simError(); | |
# | Line 410 | Line 421 | namespace OpenMD { | |
421 | } | |
422 | ||
423 | void ForceMatrixDecomposition::zeroWorkArrays() { | |
424 | + | pairwisePot = 0.0; |
425 | + | embeddingPot = 0.0; |
426 | ||
414 | – | for (int j = 0; j < N_INTERACTION_FAMILIES; j++) { |
415 | – | longRangePot_[j] = 0.0; |
416 | – | } |
417 | – | |
427 | #ifdef IS_MPI | |
428 | if (storageLayout_ & DataStorage::dslForce) { | |
429 | fill(atomRowData.force.begin(), atomRowData.force.end(), V3Zero); | |
# | Line 430 | Line 439 | namespace OpenMD { | |
439 | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | |
440 | ||
441 | fill(pot_col.begin(), pot_col.end(), | |
442 | < | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); |
434 | < | |
435 | < | pot_local = Vector<RealType, N_INTERACTION_FAMILIES>(0.0); |
442 | > | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); |
443 | ||
444 | if (storageLayout_ & DataStorage::dslParticlePot) { | |
445 | < | fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(), 0.0); |
446 | < | fill(atomColData.particlePot.begin(), atomColData.particlePot.end(), 0.0); |
445 | > | fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(), |
446 | > | 0.0); |
447 | > | fill(atomColData.particlePot.begin(), atomColData.particlePot.end(), |
448 | > | 0.0); |
449 | } | |
450 | ||
451 | if (storageLayout_ & DataStorage::dslDensity) { | |
# | Line 445 | Line 454 | namespace OpenMD { | |
454 | } | |
455 | ||
456 | if (storageLayout_ & DataStorage::dslFunctional) { | |
457 | < | fill(atomRowData.functional.begin(), atomRowData.functional.end(), 0.0); |
458 | < | fill(atomColData.functional.begin(), atomColData.functional.end(), 0.0); |
457 | > | fill(atomRowData.functional.begin(), atomRowData.functional.end(), |
458 | > | 0.0); |
459 | > | fill(atomColData.functional.begin(), atomColData.functional.end(), |
460 | > | 0.0); |
461 | } | |
462 | ||
463 | if (storageLayout_ & DataStorage::dslFunctionalDerivative) { | |
# | Line 456 | Line 467 | namespace OpenMD { | |
467 | atomColData.functionalDerivative.end(), 0.0); | |
468 | } | |
469 | ||
470 | < | #else |
471 | < | |
470 | > | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
471 | > | fill(atomRowData.skippedCharge.begin(), |
472 | > | atomRowData.skippedCharge.end(), 0.0); |
473 | > | fill(atomColData.skippedCharge.begin(), |
474 | > | atomColData.skippedCharge.end(), 0.0); |
475 | > | } |
476 | > | |
477 | > | #endif |
478 | > | // even in parallel, we need to zero out the local arrays: |
479 | > | |
480 | if (storageLayout_ & DataStorage::dslParticlePot) { | |
481 | fill(snap_->atomData.particlePot.begin(), | |
482 | snap_->atomData.particlePot.end(), 0.0); | |
# | Line 475 | Line 494 | namespace OpenMD { | |
494 | fill(snap_->atomData.functionalDerivative.begin(), | |
495 | snap_->atomData.functionalDerivative.end(), 0.0); | |
496 | } | |
497 | < | #endif |
497 | > | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
498 | > | fill(snap_->atomData.skippedCharge.begin(), |
499 | > | snap_->atomData.skippedCharge.end(), 0.0); |
500 | > | } |
501 | ||
502 | } | |
503 | ||
# | Line 512 | Line 534 | namespace OpenMD { | |
534 | AtomCommMatrixColumn->gather(snap_->atomData.electroFrame, | |
535 | atomColData.electroFrame); | |
536 | } | |
537 | + | |
538 | #endif | |
539 | } | |
540 | ||
# | Line 578 | Line 601 | namespace OpenMD { | |
601 | AtomCommVectorColumn->scatter(atomColData.force, frc_tmp); | |
602 | for (int i = 0; i < n; i++) | |
603 | snap_->atomData.force[i] += frc_tmp[i]; | |
604 | < | |
582 | < | |
604 | > | |
605 | if (storageLayout_ & DataStorage::dslTorque) { | |
606 | ||
607 | < | int nt = snap_->atomData.force.size(); |
607 | > | int nt = snap_->atomData.torque.size(); |
608 | vector<Vector3d> trq_tmp(nt, V3Zero); | |
609 | ||
610 | AtomCommVectorRow->scatter(atomRowData.torque, trq_tmp); | |
611 | < | for (int i = 0; i < n; i++) { |
611 | > | for (int i = 0; i < nt; i++) { |
612 | snap_->atomData.torque[i] += trq_tmp[i]; | |
613 | trq_tmp[i] = 0.0; | |
614 | } | |
615 | ||
616 | AtomCommVectorColumn->scatter(atomColData.torque, trq_tmp); | |
617 | < | for (int i = 0; i < n; i++) |
617 | > | for (int i = 0; i < nt; i++) |
618 | snap_->atomData.torque[i] += trq_tmp[i]; | |
619 | + | } |
620 | + | |
621 | + | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
622 | + | |
623 | + | int ns = snap_->atomData.skippedCharge.size(); |
624 | + | vector<RealType> skch_tmp(ns, 0.0); |
625 | + | |
626 | + | AtomCommRealRow->scatter(atomRowData.skippedCharge, skch_tmp); |
627 | + | for (int i = 0; i < ns; i++) { |
628 | + | snap_->atomData.skippedCharge[i] += skch_tmp[i]; |
629 | + | skch_tmp[i] = 0.0; |
630 | + | } |
631 | + | |
632 | + | AtomCommRealColumn->scatter(atomColData.skippedCharge, skch_tmp); |
633 | + | for (int i = 0; i < ns; i++) |
634 | + | snap_->atomData.skippedCharge[i] += skch_tmp[i]; |
635 | } | |
636 | ||
637 | nLocal_ = snap_->getNumberOfAtoms(); | |
# | Line 606 | Line 644 | namespace OpenMD { | |
644 | AtomCommPotRow->scatter(pot_row, pot_temp); | |
645 | ||
646 | for (int ii = 0; ii < pot_temp.size(); ii++ ) | |
647 | < | pot_local += pot_temp[ii]; |
647 | > | pairwisePot += pot_temp[ii]; |
648 | ||
649 | fill(pot_temp.begin(), pot_temp.end(), | |
650 | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | |
# | Line 614 | Line 652 | namespace OpenMD { | |
652 | AtomCommPotColumn->scatter(pot_col, pot_temp); | |
653 | ||
654 | for (int ii = 0; ii < pot_temp.size(); ii++ ) | |
655 | < | pot_local += pot_temp[ii]; |
618 | < | |
655 | > | pairwisePot += pot_temp[ii]; |
656 | #endif | |
657 | + | |
658 | } | |
659 | ||
660 | int ForceMatrixDecomposition::getNAtomsInRow() { | |
# | Line 691 | Line 729 | namespace OpenMD { | |
729 | #ifdef IS_MPI | |
730 | return massFactorsRow[atom1]; | |
731 | #else | |
732 | < | return massFactorsLocal[atom1]; |
732 | > | return massFactors[atom1]; |
733 | #endif | |
734 | } | |
735 | ||
# | Line 699 | Line 737 | namespace OpenMD { | |
737 | #ifdef IS_MPI | |
738 | return massFactorsCol[atom2]; | |
739 | #else | |
740 | < | return massFactorsLocal[atom2]; |
740 | > | return massFactors[atom2]; |
741 | #endif | |
742 | ||
743 | } | |
# | Line 717 | Line 755 | namespace OpenMD { | |
755 | return d; | |
756 | } | |
757 | ||
758 | < | vector<int> ForceMatrixDecomposition::getSkipsForAtom(int atom1) { |
759 | < | return skipsForAtom[atom1]; |
758 | > | vector<int> ForceMatrixDecomposition::getExcludesForAtom(int atom1) { |
759 | > | return excludesForAtom[atom1]; |
760 | } | |
761 | ||
762 | /** | |
763 | < | * There are a number of reasons to skip a pair or a |
726 | < | * particle. Mostly we do this to exclude atoms who are involved in |
727 | < | * short range interactions (bonds, bends, torsions), but we also |
728 | < | * need to exclude some overcounted interactions that result from |
763 | > | * We need to exclude some overcounted interactions that result from |
764 | * the parallel decomposition. | |
765 | */ | |
766 | bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) { | |
# | Line 745 | Line 780 | namespace OpenMD { | |
780 | } else { | |
781 | if ((unique_id_1 + unique_id_2) % 2 == 1) return true; | |
782 | } | |
783 | + | #endif |
784 | + | return false; |
785 | + | } |
786 | + | |
787 | + | /** |
788 | + | * We need to handle the interactions for atoms who are involved in |
789 | + | * the same rigid body as well as some short range interactions |
790 | + | * (bonds, bends, torsions) differently from other interactions. |
791 | + | * We'll still visit the pairwise routines, but with a flag that |
792 | + | * tells those routines to exclude the pair from direct long range |
793 | + | * interactions. Some indirect interactions (notably reaction |
794 | + | * field) must still be handled for these pairs. |
795 | + | */ |
796 | + | bool ForceMatrixDecomposition::excludeAtomPair(int atom1, int atom2) { |
797 | + | int unique_id_2; |
798 | + | |
799 | + | #ifdef IS_MPI |
800 | + | // in MPI, we have to look up the unique IDs for the row atom. |
801 | + | unique_id_2 = AtomColToGlobal[atom2]; |
802 | #else | |
803 | // in the normal loop, the atom numbers are unique | |
750 | – | unique_id_1 = atom1; |
804 | unique_id_2 = atom2; | |
805 | #endif | |
806 | ||
807 | < | for (vector<int>::iterator i = skipsForAtom[atom1].begin(); |
808 | < | i != skipsForAtom[atom1].end(); ++i) { |
807 | > | for (vector<int>::iterator i = excludesForAtom[atom1].begin(); |
808 | > | i != excludesForAtom[atom1].end(); ++i) { |
809 | if ( (*i) == unique_id_2 ) return true; | |
810 | < | } |
810 | > | } |
811 | ||
812 | + | return false; |
813 | } | |
814 | ||
815 | ||
# | Line 776 | Line 830 | namespace OpenMD { | |
830 | } | |
831 | ||
832 | // filling interaction blocks with pointers | |
833 | < | InteractionData ForceMatrixDecomposition::fillInteractionData(int atom1, int atom2) { |
834 | < | InteractionData idat; |
833 | > | void ForceMatrixDecomposition::fillInteractionData(InteractionData &idat, |
834 | > | int atom1, int atom2) { |
835 | ||
836 | + | idat.excluded = excludeAtomPair(atom1, atom2); |
837 | + | |
838 | #ifdef IS_MPI | |
839 | ||
840 | idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), | |
841 | ff_->getAtomType(identsCol[atom2]) ); | |
786 | – | |
842 | ||
843 | if (storageLayout_ & DataStorage::dslAmat) { | |
844 | idat.A1 = &(atomRowData.aMat[atom1]); | |
# | Line 820 | Line 875 | namespace OpenMD { | |
875 | idat.particlePot2 = &(atomColData.particlePot[atom2]); | |
876 | } | |
877 | ||
878 | + | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
879 | + | idat.skippedCharge1 = &(atomRowData.skippedCharge[atom1]); |
880 | + | idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]); |
881 | + | } |
882 | + | |
883 | #else | |
884 | ||
885 | < | idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]), |
886 | < | ff_->getAtomType(identsLocal[atom2]) ); |
885 | > | idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), |
886 | > | ff_->getAtomType(idents[atom2]) ); |
887 | ||
888 | if (storageLayout_ & DataStorage::dslAmat) { | |
889 | idat.A1 = &(snap_->atomData.aMat[atom1]); | |
# | Line 840 | Line 900 | namespace OpenMD { | |
900 | idat.t2 = &(snap_->atomData.torque[atom2]); | |
901 | } | |
902 | ||
903 | < | if (storageLayout_ & DataStorage::dslDensity) { |
903 | > | if (storageLayout_ & DataStorage::dslDensity) { |
904 | idat.rho1 = &(snap_->atomData.density[atom1]); | |
905 | idat.rho2 = &(snap_->atomData.density[atom2]); | |
906 | } | |
# | Line 860 | Line 920 | namespace OpenMD { | |
920 | idat.particlePot2 = &(snap_->atomData.particlePot[atom2]); | |
921 | } | |
922 | ||
923 | + | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
924 | + | idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]); |
925 | + | idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]); |
926 | + | } |
927 | #endif | |
864 | – | return idat; |
928 | } | |
929 | ||
930 | ||
931 | < | void ForceMatrixDecomposition::unpackInteractionData(InteractionData idat, int atom1, int atom2) { |
931 | > | void ForceMatrixDecomposition::unpackInteractionData(InteractionData &idat, int atom1, int atom2) { |
932 | #ifdef IS_MPI | |
933 | pot_row[atom1] += 0.5 * *(idat.pot); | |
934 | pot_col[atom2] += 0.5 * *(idat.pot); | |
# | Line 873 | Line 936 | namespace OpenMD { | |
936 | atomRowData.force[atom1] += *(idat.f1); | |
937 | atomColData.force[atom2] -= *(idat.f1); | |
938 | #else | |
939 | < | longRangePot_ += *(idat.pot); |
940 | < | |
939 | > | pairwisePot += *(idat.pot); |
940 | > | |
941 | snap_->atomData.force[atom1] += *(idat.f1); | |
942 | snap_->atomData.force[atom2] -= *(idat.f1); | |
943 | #endif | |
944 | < | |
944 | > | |
945 | } | |
946 | ||
884 | – | |
885 | – | InteractionData ForceMatrixDecomposition::fillSkipData(int atom1, int atom2){ |
886 | – | |
887 | – | InteractionData idat; |
888 | – | #ifdef IS_MPI |
889 | – | idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), |
890 | – | ff_->getAtomType(identsCol[atom2]) ); |
891 | – | |
892 | – | if (storageLayout_ & DataStorage::dslElectroFrame) { |
893 | – | idat.eFrame1 = &(atomRowData.electroFrame[atom1]); |
894 | – | idat.eFrame2 = &(atomColData.electroFrame[atom2]); |
895 | – | } |
896 | – | if (storageLayout_ & DataStorage::dslTorque) { |
897 | – | idat.t1 = &(atomRowData.torque[atom1]); |
898 | – | idat.t2 = &(atomColData.torque[atom2]); |
899 | – | } |
900 | – | #else |
901 | – | idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]), |
902 | – | ff_->getAtomType(identsLocal[atom2]) ); |
903 | – | |
904 | – | if (storageLayout_ & DataStorage::dslElectroFrame) { |
905 | – | idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]); |
906 | – | idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]); |
907 | – | } |
908 | – | if (storageLayout_ & DataStorage::dslTorque) { |
909 | – | idat.t1 = &(snap_->atomData.torque[atom1]); |
910 | – | idat.t2 = &(snap_->atomData.torque[atom2]); |
911 | – | } |
912 | – | #endif |
913 | – | } |
914 | – | |
947 | /* | |
948 | * buildNeighborList | |
949 | * | |
# | Line 922 | Line 954 | namespace OpenMD { | |
954 | ||
955 | vector<pair<int, int> > neighborList; | |
956 | groupCutoffs cuts; | |
957 | + | bool doAllPairs = false; |
958 | + | |
959 | #ifdef IS_MPI | |
960 | cellListRow_.clear(); | |
961 | cellListCol_.clear(); | |
# | Line 941 | Line 975 | namespace OpenMD { | |
975 | nCells_.y() = (int) ( Hy.length() )/ rList_; | |
976 | nCells_.z() = (int) ( Hz.length() )/ rList_; | |
977 | ||
978 | + | // handle small boxes where the cell offsets can end up repeating cells |
979 | + | |
980 | + | if (nCells_.x() < 3) doAllPairs = true; |
981 | + | if (nCells_.y() < 3) doAllPairs = true; |
982 | + | if (nCells_.z() < 3) doAllPairs = true; |
983 | + | |
984 | Mat3x3d invHmat = snap_->getInvHmat(); | |
985 | Vector3d rs, scaled, dr; | |
986 | Vector3i whichCell; | |
# | Line 954 | Line 994 | namespace OpenMD { | |
994 | cellList_.resize(nCtot); | |
995 | #endif | |
996 | ||
997 | + | if (!doAllPairs) { |
998 | #ifdef IS_MPI | |
999 | < | for (int i = 0; i < nGroupsInRow_; i++) { |
1000 | < | rs = cgRowData.position[i]; |
1001 | < | // scaled positions relative to the box vectors |
1002 | < | scaled = invHmat * rs; |
1003 | < | // wrap the vector back into the unit box by subtracting integer box |
1004 | < | // numbers |
1005 | < | for (int j = 0; j < 3; j++) |
1006 | < | scaled[j] -= roundMe(scaled[j]); |
1007 | < | |
1008 | < | // find xyz-indices of cell that cutoffGroup is in. |
1009 | < | whichCell.x() = nCells_.x() * scaled.x(); |
1010 | < | whichCell.y() = nCells_.y() * scaled.y(); |
1011 | < | whichCell.z() = nCells_.z() * scaled.z(); |
1012 | < | |
1013 | < | // find single index of this cell: |
1014 | < | cellIndex = Vlinear(whichCell, nCells_); |
1015 | < | // add this cutoff group to the list of groups in this cell; |
1016 | < | cellListRow_[cellIndex].push_back(i); |
1017 | < | } |
1018 | < | |
1019 | < | for (int i = 0; i < nGroupsInCol_; i++) { |
1020 | < | rs = cgColData.position[i]; |
1021 | < | // scaled positions relative to the box vectors |
1022 | < | scaled = invHmat * rs; |
1023 | < | // wrap the vector back into the unit box by subtracting integer box |
1024 | < | // numbers |
1025 | < | for (int j = 0; j < 3; j++) |
1026 | < | scaled[j] -= roundMe(scaled[j]); |
1027 | < | |
1028 | < | // find xyz-indices of cell that cutoffGroup is in. |
1029 | < | whichCell.x() = nCells_.x() * scaled.x(); |
1030 | < | whichCell.y() = nCells_.y() * scaled.y(); |
1031 | < | whichCell.z() = nCells_.z() * scaled.z(); |
1032 | < | |
1033 | < | // find single index of this cell: |
1034 | < | cellIndex = Vlinear(whichCell, nCells_); |
1035 | < | // add this cutoff group to the list of groups in this cell; |
1036 | < | cellListCol_[cellIndex].push_back(i); |
1037 | < | } |
999 | > | |
1000 | > | for (int i = 0; i < nGroupsInRow_; i++) { |
1001 | > | rs = cgRowData.position[i]; |
1002 | > | |
1003 | > | // scaled positions relative to the box vectors |
1004 | > | scaled = invHmat * rs; |
1005 | > | |
1006 | > | // wrap the vector back into the unit box by subtracting integer box |
1007 | > | // numbers |
1008 | > | for (int j = 0; j < 3; j++) { |
1009 | > | scaled[j] -= roundMe(scaled[j]); |
1010 | > | scaled[j] += 0.5; |
1011 | > | } |
1012 | > | |
1013 | > | // find xyz-indices of cell that cutoffGroup is in. |
1014 | > | whichCell.x() = nCells_.x() * scaled.x(); |
1015 | > | whichCell.y() = nCells_.y() * scaled.y(); |
1016 | > | whichCell.z() = nCells_.z() * scaled.z(); |
1017 | > | |
1018 | > | // find single index of this cell: |
1019 | > | cellIndex = Vlinear(whichCell, nCells_); |
1020 | > | |
1021 | > | // add this cutoff group to the list of groups in this cell; |
1022 | > | cellListRow_[cellIndex].push_back(i); |
1023 | > | } |
1024 | > | |
1025 | > | for (int i = 0; i < nGroupsInCol_; i++) { |
1026 | > | rs = cgColData.position[i]; |
1027 | > | |
1028 | > | // scaled positions relative to the box vectors |
1029 | > | scaled = invHmat * rs; |
1030 | > | |
1031 | > | // wrap the vector back into the unit box by subtracting integer box |
1032 | > | // numbers |
1033 | > | for (int j = 0; j < 3; j++) { |
1034 | > | scaled[j] -= roundMe(scaled[j]); |
1035 | > | scaled[j] += 0.5; |
1036 | > | } |
1037 | > | |
1038 | > | // find xyz-indices of cell that cutoffGroup is in. |
1039 | > | whichCell.x() = nCells_.x() * scaled.x(); |
1040 | > | whichCell.y() = nCells_.y() * scaled.y(); |
1041 | > | whichCell.z() = nCells_.z() * scaled.z(); |
1042 | > | |
1043 | > | // find single index of this cell: |
1044 | > | cellIndex = Vlinear(whichCell, nCells_); |
1045 | > | |
1046 | > | // add this cutoff group to the list of groups in this cell; |
1047 | > | cellListCol_[cellIndex].push_back(i); |
1048 | > | } |
1049 | #else | |
1050 | < | for (int i = 0; i < nGroups_; i++) { |
1051 | < | rs = snap_->cgData.position[i]; |
1052 | < | // scaled positions relative to the box vectors |
1053 | < | scaled = invHmat * rs; |
1054 | < | // wrap the vector back into the unit box by subtracting integer box |
1055 | < | // numbers |
1056 | < | for (int j = 0; j < 3; j++) |
1057 | < | scaled[j] -= roundMe(scaled[j]); |
1058 | < | |
1059 | < | // find xyz-indices of cell that cutoffGroup is in. |
1060 | < | whichCell.x() = nCells_.x() * scaled.x(); |
1061 | < | whichCell.y() = nCells_.y() * scaled.y(); |
1062 | < | whichCell.z() = nCells_.z() * scaled.z(); |
1063 | < | |
1064 | < | // find single index of this cell: |
1065 | < | cellIndex = Vlinear(whichCell, nCells_); |
1066 | < | // add this cutoff group to the list of groups in this cell; |
1067 | < | cellList_[cellIndex].push_back(i); |
1068 | < | } |
1050 | > | for (int i = 0; i < nGroups_; i++) { |
1051 | > | rs = snap_->cgData.position[i]; |
1052 | > | |
1053 | > | // scaled positions relative to the box vectors |
1054 | > | scaled = invHmat * rs; |
1055 | > | |
1056 | > | // wrap the vector back into the unit box by subtracting integer box |
1057 | > | // numbers |
1058 | > | for (int j = 0; j < 3; j++) { |
1059 | > | scaled[j] -= roundMe(scaled[j]); |
1060 | > | scaled[j] += 0.5; |
1061 | > | } |
1062 | > | |
1063 | > | // find xyz-indices of cell that cutoffGroup is in. |
1064 | > | whichCell.x() = nCells_.x() * scaled.x(); |
1065 | > | whichCell.y() = nCells_.y() * scaled.y(); |
1066 | > | whichCell.z() = nCells_.z() * scaled.z(); |
1067 | > | |
1068 | > | // find single index of this cell: |
1069 | > | cellIndex = Vlinear(whichCell, nCells_); |
1070 | > | |
1071 | > | // add this cutoff group to the list of groups in this cell; |
1072 | > | cellList_[cellIndex].push_back(i); |
1073 | > | } |
1074 | #endif | |
1075 | ||
1076 | < | for (int m1z = 0; m1z < nCells_.z(); m1z++) { |
1077 | < | for (int m1y = 0; m1y < nCells_.y(); m1y++) { |
1078 | < | for (int m1x = 0; m1x < nCells_.x(); m1x++) { |
1079 | < | Vector3i m1v(m1x, m1y, m1z); |
1080 | < | int m1 = Vlinear(m1v, nCells_); |
1024 | < | |
1025 | < | for (vector<Vector3i>::iterator os = cellOffsets_.begin(); |
1026 | < | os != cellOffsets_.end(); ++os) { |
1076 | > | for (int m1z = 0; m1z < nCells_.z(); m1z++) { |
1077 | > | for (int m1y = 0; m1y < nCells_.y(); m1y++) { |
1078 | > | for (int m1x = 0; m1x < nCells_.x(); m1x++) { |
1079 | > | Vector3i m1v(m1x, m1y, m1z); |
1080 | > | int m1 = Vlinear(m1v, nCells_); |
1081 | ||
1082 | < | Vector3i m2v = m1v + (*os); |
1083 | < | |
1084 | < | if (m2v.x() >= nCells_.x()) { |
1085 | < | m2v.x() = 0; |
1086 | < | } else if (m2v.x() < 0) { |
1087 | < | m2v.x() = nCells_.x() - 1; |
1088 | < | } |
1089 | < | |
1090 | < | if (m2v.y() >= nCells_.y()) { |
1091 | < | m2v.y() = 0; |
1092 | < | } else if (m2v.y() < 0) { |
1093 | < | m2v.y() = nCells_.y() - 1; |
1094 | < | } |
1095 | < | |
1096 | < | if (m2v.z() >= nCells_.z()) { |
1097 | < | m2v.z() = 0; |
1098 | < | } else if (m2v.z() < 0) { |
1099 | < | m2v.z() = nCells_.z() - 1; |
1100 | < | } |
1101 | < | |
1102 | < | int m2 = Vlinear (m2v, nCells_); |
1103 | < | |
1082 | > | for (vector<Vector3i>::iterator os = cellOffsets_.begin(); |
1083 | > | os != cellOffsets_.end(); ++os) { |
1084 | > | |
1085 | > | Vector3i m2v = m1v + (*os); |
1086 | > | |
1087 | > | if (m2v.x() >= nCells_.x()) { |
1088 | > | m2v.x() = 0; |
1089 | > | } else if (m2v.x() < 0) { |
1090 | > | m2v.x() = nCells_.x() - 1; |
1091 | > | } |
1092 | > | |
1093 | > | if (m2v.y() >= nCells_.y()) { |
1094 | > | m2v.y() = 0; |
1095 | > | } else if (m2v.y() < 0) { |
1096 | > | m2v.y() = nCells_.y() - 1; |
1097 | > | } |
1098 | > | |
1099 | > | if (m2v.z() >= nCells_.z()) { |
1100 | > | m2v.z() = 0; |
1101 | > | } else if (m2v.z() < 0) { |
1102 | > | m2v.z() = nCells_.z() - 1; |
1103 | > | } |
1104 | > | |
1105 | > | int m2 = Vlinear (m2v, nCells_); |
1106 | > | |
1107 | #ifdef IS_MPI | |
1108 | < | for (vector<int>::iterator j1 = cellListRow_[m1].begin(); |
1109 | < | j1 != cellListRow_[m1].end(); ++j1) { |
1110 | < | for (vector<int>::iterator j2 = cellListCol_[m2].begin(); |
1111 | < | j2 != cellListCol_[m2].end(); ++j2) { |
1112 | < | |
1113 | < | // Always do this if we're in different cells or if |
1114 | < | // we're in the same cell and the global index of the |
1115 | < | // j2 cutoff group is less than the j1 cutoff group |
1116 | < | |
1117 | < | if (m2 != m1 || cgColToGlobal[(*j2)] < cgRowToGlobal[(*j1)]) { |
1118 | < | dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; |
1119 | < | snap_->wrapVector(dr); |
1120 | < | cuts = getGroupCutoffs( (*j1), (*j2) ); |
1121 | < | if (dr.lengthSquare() < cuts.third) { |
1122 | < | neighborList.push_back(make_pair((*j1), (*j2))); |
1108 | > | for (vector<int>::iterator j1 = cellListRow_[m1].begin(); |
1109 | > | j1 != cellListRow_[m1].end(); ++j1) { |
1110 | > | for (vector<int>::iterator j2 = cellListCol_[m2].begin(); |
1111 | > | j2 != cellListCol_[m2].end(); ++j2) { |
1112 | > | |
1113 | > | // Always do this if we're in different cells or if |
1114 | > | // we're in the same cell and the global index of the |
1115 | > | // j2 cutoff group is less than the j1 cutoff group |
1116 | > | |
1117 | > | if (m2 != m1 || cgColToGlobal[(*j2)] < cgRowToGlobal[(*j1)]) { |
1118 | > | dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; |
1119 | > | snap_->wrapVector(dr); |
1120 | > | cuts = getGroupCutoffs( (*j1), (*j2) ); |
1121 | > | if (dr.lengthSquare() < cuts.third) { |
1122 | > | neighborList.push_back(make_pair((*j1), (*j2))); |
1123 | > | } |
1124 | } | |
1125 | } | |
1126 | } | |
1069 | – | } |
1127 | #else | |
1128 | < | for (vector<int>::iterator j1 = cellList_[m1].begin(); |
1129 | < | j1 != cellList_[m1].end(); ++j1) { |
1130 | < | for (vector<int>::iterator j2 = cellList_[m2].begin(); |
1131 | < | j2 != cellList_[m2].end(); ++j2) { |
1132 | < | |
1133 | < | // Always do this if we're in different cells or if |
1134 | < | // we're in the same cell and the global index of the |
1135 | < | // j2 cutoff group is less than the j1 cutoff group |
1136 | < | |
1137 | < | if (m2 != m1 || (*j2) < (*j1)) { |
1138 | < | dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)]; |
1139 | < | snap_->wrapVector(dr); |
1140 | < | cuts = getGroupCutoffs( (*j1), (*j2) ); |
1141 | < | if (dr.lengthSquare() < cuts.third) { |
1142 | < | neighborList.push_back(make_pair((*j1), (*j2))); |
1128 | > | |
1129 | > | for (vector<int>::iterator j1 = cellList_[m1].begin(); |
1130 | > | j1 != cellList_[m1].end(); ++j1) { |
1131 | > | for (vector<int>::iterator j2 = cellList_[m2].begin(); |
1132 | > | j2 != cellList_[m2].end(); ++j2) { |
1133 | > | |
1134 | > | // Always do this if we're in different cells or if |
1135 | > | // we're in the same cell and the global index of the |
1136 | > | // j2 cutoff group is less than the j1 cutoff group |
1137 | > | |
1138 | > | if (m2 != m1 || (*j2) < (*j1)) { |
1139 | > | dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)]; |
1140 | > | snap_->wrapVector(dr); |
1141 | > | cuts = getGroupCutoffs( (*j1), (*j2) ); |
1142 | > | if (dr.lengthSquare() < cuts.third) { |
1143 | > | neighborList.push_back(make_pair((*j1), (*j2))); |
1144 | > | } |
1145 | } | |
1146 | } | |
1147 | } | |
1089 | – | } |
1148 | #endif | |
1149 | + | } |
1150 | } | |
1151 | } | |
1152 | } | |
1153 | + | } else { |
1154 | + | // branch to do all cutoff group pairs |
1155 | + | #ifdef IS_MPI |
1156 | + | for (int j1 = 0; j1 < nGroupsInRow_; j1++) { |
1157 | + | for (int j2 = 0; j2 < nGroupsInCol_; j2++) { |
1158 | + | dr = cgColData.position[j2] - cgRowData.position[j1]; |
1159 | + | snap_->wrapVector(dr); |
1160 | + | cuts = getGroupCutoffs( j1, j2 ); |
1161 | + | if (dr.lengthSquare() < cuts.third) { |
1162 | + | neighborList.push_back(make_pair(j1, j2)); |
1163 | + | } |
1164 | + | } |
1165 | + | } |
1166 | + | #else |
1167 | + | for (int j1 = 0; j1 < nGroups_ - 1; j1++) { |
1168 | + | for (int j2 = j1 + 1; j2 < nGroups_; j2++) { |
1169 | + | dr = snap_->cgData.position[j2] - snap_->cgData.position[j1]; |
1170 | + | snap_->wrapVector(dr); |
1171 | + | cuts = getGroupCutoffs( j1, j2 ); |
1172 | + | if (dr.lengthSquare() < cuts.third) { |
1173 | + | neighborList.push_back(make_pair(j1, j2)); |
1174 | + | } |
1175 | + | } |
1176 | + | } |
1177 | + | #endif |
1178 | } | |
1179 | < | |
1179 | > | |
1180 | // save the local cutoff group positions for the check that is | |
1181 | // done on each loop: | |
1182 | saved_CG_positions_.clear(); | |
1183 | for (int i = 0; i < nGroups_; i++) | |
1184 | saved_CG_positions_.push_back(snap_->cgData.position[i]); | |
1185 | < | |
1185 | > | |
1186 | return neighborList; | |
1187 | } | |
1188 | } //end namespace OpenMD |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |