# | 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 | idents = info_->getIdentArray(); | |
64 | AtomLocalToGlobal = info_->getGlobalAtomIndices(); | |
65 | cgLocalToGlobal = info_->getGlobalGroupIndices(); | |
66 | vector<int> globalGroupMembership = info_->getGlobalGroupMembership(); | |
67 | + | |
68 | massFactors = info_->getMassFactors(); | |
69 | – | PairList excludes = info_->getExcludedInteractions(); |
70 | – | PairList oneTwo = info_->getOneTwoInteractions(); |
71 | – | PairList oneThree = info_->getOneThreeInteractions(); |
72 | – | PairList oneFour = info_->getOneFourInteractions(); |
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 111 | Line 112 | namespace OpenMD { | |
112 | AtomCommIntRow->gather(idents, identsRow); | |
113 | AtomCommIntColumn->gather(idents, identsCol); | |
114 | ||
115 | + | vector<int>::iterator it; |
116 | + | for (it = AtomLocalToGlobal.begin(); it != AtomLocalToGlobal.end(); ++it) { |
117 | + | cerr << "my AtomLocalToGlobal = " << (*it) << "\n"; |
118 | + | } |
119 | AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal); | |
120 | AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal); | |
121 | ||
# | Line 142 | Line 147 | namespace OpenMD { | |
147 | } | |
148 | } | |
149 | ||
150 | < | skipsForAtom.clear(); |
151 | < | skipsForAtom.resize(nAtomsInRow_); |
150 | > | excludesForAtom.clear(); |
151 | > | excludesForAtom.resize(nAtomsInRow_); |
152 | toposForAtom.clear(); | |
153 | toposForAtom.resize(nAtomsInRow_); | |
154 | topoDist.clear(); | |
# | Line 154 | Line 159 | namespace OpenMD { | |
159 | for (int j = 0; j < nAtomsInCol_; j++) { | |
160 | int jglob = AtomColToGlobal[j]; | |
161 | ||
162 | < | if (excludes.hasPair(iglob, jglob)) |
163 | < | skipsForAtom[i].push_back(j); |
162 | > | if (excludes->hasPair(iglob, jglob)) |
163 | > | excludesForAtom[i].push_back(j); |
164 | ||
165 | < | if (oneTwo.hasPair(iglob, jglob)) { |
165 | > | if (oneTwo->hasPair(iglob, jglob)) { |
166 | toposForAtom[i].push_back(j); | |
167 | topoDist[i].push_back(1); | |
168 | } else { | |
169 | < | if (oneThree.hasPair(iglob, jglob)) { |
169 | > | if (oneThree->hasPair(iglob, jglob)) { |
170 | toposForAtom[i].push_back(j); | |
171 | topoDist[i].push_back(2); | |
172 | } else { | |
173 | < | if (oneFour.hasPair(iglob, jglob)) { |
173 | > | if (oneFour->hasPair(iglob, jglob)) { |
174 | toposForAtom[i].push_back(j); | |
175 | topoDist[i].push_back(3); | |
176 | } | |
# | Line 188 | Line 193 | namespace OpenMD { | |
193 | } | |
194 | } | |
195 | ||
196 | < | skipsForAtom.clear(); |
197 | < | skipsForAtom.resize(nLocal_); |
196 | > | excludesForAtom.clear(); |
197 | > | excludesForAtom.resize(nLocal_); |
198 | toposForAtom.clear(); | |
199 | toposForAtom.resize(nLocal_); | |
200 | topoDist.clear(); | |
# | Line 201 | Line 206 | namespace OpenMD { | |
206 | for (int j = 0; j < nLocal_; j++) { | |
207 | int jglob = AtomLocalToGlobal[j]; | |
208 | ||
209 | < | if (excludes.hasPair(iglob, jglob)) |
210 | < | skipsForAtom[i].push_back(j); |
209 | > | if (excludes->hasPair(iglob, jglob)) |
210 | > | excludesForAtom[i].push_back(j); |
211 | ||
212 | < | if (oneTwo.hasPair(iglob, jglob)) { |
212 | > | if (oneTwo->hasPair(iglob, jglob)) { |
213 | toposForAtom[i].push_back(j); | |
214 | topoDist[i].push_back(1); | |
215 | } else { | |
216 | < | if (oneThree.hasPair(iglob, jglob)) { |
216 | > | if (oneThree->hasPair(iglob, jglob)) { |
217 | toposForAtom[i].push_back(j); | |
218 | topoDist[i].push_back(2); | |
219 | } else { | |
220 | < | if (oneFour.hasPair(iglob, jglob)) { |
220 | > | if (oneFour->hasPair(iglob, jglob)) { |
221 | toposForAtom[i].push_back(j); | |
222 | topoDist[i].push_back(3); | |
223 | } | |
# | Line 222 | Line 227 | namespace OpenMD { | |
227 | } | |
228 | ||
229 | createGtypeCutoffMap(); | |
230 | + | |
231 | } | |
232 | ||
233 | void ForceMatrixDecomposition::createGtypeCutoffMap() { | |
234 | < | |
234 | > | |
235 | RealType tol = 1e-6; | |
236 | RealType rc; | |
237 | int atid; | |
238 | set<AtomType*> atypes = info_->getSimulatedAtomTypes(); | |
239 | < | vector<RealType> atypeCutoff; |
234 | < | atypeCutoff.resize( atypes.size() ); |
239 | > | map<int, RealType> atypeCutoff; |
240 | ||
241 | for (set<AtomType*>::iterator at = atypes.begin(); | |
242 | at != atypes.end(); ++at){ | |
243 | atid = (*at)->getIdent(); | |
244 | < | |
240 | < | if (userChoseCutoff_) |
244 | > | if (userChoseCutoff_) |
245 | atypeCutoff[atid] = userCutoff_; | |
246 | < | else |
246 | > | else |
247 | atypeCutoff[atid] = interactionMan_->getSuggestedCutoffRadius(*at); | |
248 | } | |
249 | ||
250 | vector<RealType> gTypeCutoffs; | |
247 | – | |
251 | // first we do a single loop over the cutoff groups to find the | |
252 | // largest cutoff for any atypes present in this group. | |
253 | #ifdef IS_MPI | |
# | Line 302 | Line 305 | namespace OpenMD { | |
305 | ||
306 | vector<RealType> groupCutoff(nGroups_, 0.0); | |
307 | groupToGtype.resize(nGroups_); | |
305 | – | |
306 | – | cerr << "nGroups = " << nGroups_ << "\n"; |
308 | for (int cg1 = 0; cg1 < nGroups_; cg1++) { | |
309 | ||
310 | groupCutoff[cg1] = 0.0; | |
# | Line 332 | Line 333 | namespace OpenMD { | |
333 | } | |
334 | #endif | |
335 | ||
335 | – | cerr << "gTypeCutoffs.size() = " << gTypeCutoffs.size() << "\n"; |
336 | // Now we find the maximum group cutoff value present in the simulation | |
337 | ||
338 | RealType groupMax = *max_element(gTypeCutoffs.begin(), gTypeCutoffs.end()); | |
# | Line 453 | Line 453 | namespace OpenMD { | |
453 | atomRowData.functionalDerivative.end(), 0.0); | |
454 | fill(atomColData.functionalDerivative.begin(), | |
455 | atomColData.functionalDerivative.end(), 0.0); | |
456 | + | } |
457 | + | |
458 | + | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
459 | + | fill(atomRowData.skippedCharge.begin(), |
460 | + | atomRowData.skippedCharge.end(), 0.0); |
461 | + | fill(atomColData.skippedCharge.begin(), |
462 | + | atomColData.skippedCharge.end(), 0.0); |
463 | } | |
464 | ||
465 | #else | |
# | Line 474 | Line 481 | namespace OpenMD { | |
481 | fill(snap_->atomData.functionalDerivative.begin(), | |
482 | snap_->atomData.functionalDerivative.end(), 0.0); | |
483 | } | |
484 | + | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
485 | + | fill(snap_->atomData.skippedCharge.begin(), |
486 | + | snap_->atomData.skippedCharge.end(), 0.0); |
487 | + | } |
488 | #endif | |
489 | ||
490 | } | |
# | Line 581 | Line 592 | namespace OpenMD { | |
592 | ||
593 | if (storageLayout_ & DataStorage::dslTorque) { | |
594 | ||
595 | < | int nt = snap_->atomData.force.size(); |
595 | > | int nt = snap_->atomData.torque.size(); |
596 | vector<Vector3d> trq_tmp(nt, V3Zero); | |
597 | ||
598 | AtomCommVectorRow->scatter(atomRowData.torque, trq_tmp); | |
599 | < | for (int i = 0; i < n; i++) { |
599 | > | for (int i = 0; i < nt; i++) { |
600 | snap_->atomData.torque[i] += trq_tmp[i]; | |
601 | trq_tmp[i] = 0.0; | |
602 | } | |
603 | ||
604 | AtomCommVectorColumn->scatter(atomColData.torque, trq_tmp); | |
605 | < | for (int i = 0; i < n; i++) |
605 | > | for (int i = 0; i < nt; i++) |
606 | snap_->atomData.torque[i] += trq_tmp[i]; | |
607 | + | } |
608 | + | |
609 | + | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
610 | + | |
611 | + | int ns = snap_->atomData.skippedCharge.size(); |
612 | + | vector<RealType> skch_tmp(ns, 0.0); |
613 | + | |
614 | + | AtomCommRealRow->scatter(atomRowData.skippedCharge, skch_tmp); |
615 | + | for (int i = 0; i < ns; i++) { |
616 | + | snap_->atomData.skippedCharge[i] = skch_tmp[i]; |
617 | + | skch_tmp[i] = 0.0; |
618 | + | } |
619 | + | |
620 | + | AtomCommRealColumn->scatter(atomColData.skippedCharge, skch_tmp); |
621 | + | for (int i = 0; i < ns; i++) |
622 | + | snap_->atomData.skippedCharge[i] += skch_tmp[i]; |
623 | } | |
624 | ||
625 | nLocal_ = snap_->getNumberOfAtoms(); | |
# | Line 716 | Line 743 | namespace OpenMD { | |
743 | return d; | |
744 | } | |
745 | ||
746 | < | vector<int> ForceMatrixDecomposition::getSkipsForAtom(int atom1) { |
747 | < | return skipsForAtom[atom1]; |
746 | > | vector<int> ForceMatrixDecomposition::getExcludesForAtom(int atom1) { |
747 | > | return excludesForAtom[atom1]; |
748 | } | |
749 | ||
750 | /** | |
751 | < | * There are a number of reasons to skip a pair or a |
725 | < | * particle. Mostly we do this to exclude atoms who are involved in |
726 | < | * short range interactions (bonds, bends, torsions), but we also |
727 | < | * need to exclude some overcounted interactions that result from |
751 | > | * We need to exclude some overcounted interactions that result from |
752 | * the parallel decomposition. | |
753 | */ | |
754 | bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) { | |
# | Line 744 | Line 768 | namespace OpenMD { | |
768 | } else { | |
769 | if ((unique_id_1 + unique_id_2) % 2 == 1) return true; | |
770 | } | |
771 | + | #endif |
772 | + | return false; |
773 | + | } |
774 | + | |
775 | + | /** |
776 | + | * We need to handle the interactions for atoms who are involved in |
777 | + | * the same rigid body as well as some short range interactions |
778 | + | * (bonds, bends, torsions) differently from other interactions. |
779 | + | * We'll still visit the pairwise routines, but with a flag that |
780 | + | * tells those routines to exclude the pair from direct long range |
781 | + | * interactions. Some indirect interactions (notably reaction |
782 | + | * field) must still be handled for these pairs. |
783 | + | */ |
784 | + | bool ForceMatrixDecomposition::excludeAtomPair(int atom1, int atom2) { |
785 | + | int unique_id_2; |
786 | + | |
787 | + | #ifdef IS_MPI |
788 | + | // in MPI, we have to look up the unique IDs for the row atom. |
789 | + | unique_id_2 = AtomColToGlobal[atom2]; |
790 | #else | |
791 | // in the normal loop, the atom numbers are unique | |
749 | – | unique_id_1 = atom1; |
792 | unique_id_2 = atom2; | |
793 | #endif | |
794 | ||
795 | < | for (vector<int>::iterator i = skipsForAtom[atom1].begin(); |
796 | < | i != skipsForAtom[atom1].end(); ++i) { |
795 | > | for (vector<int>::iterator i = excludesForAtom[atom1].begin(); |
796 | > | i != excludesForAtom[atom1].end(); ++i) { |
797 | if ( (*i) == unique_id_2 ) return true; | |
798 | } | |
799 | ||
# | Line 777 | Line 819 | namespace OpenMD { | |
819 | ||
820 | // filling interaction blocks with pointers | |
821 | void ForceMatrixDecomposition::fillInteractionData(InteractionData &idat, | |
822 | < | int atom1, int atom2) { |
822 | > | int atom1, int atom2) { |
823 | > | |
824 | > | idat.excluded = excludeAtomPair(atom1, atom2); |
825 | > | |
826 | #ifdef IS_MPI | |
827 | ||
828 | idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), | |
# | Line 818 | Line 863 | namespace OpenMD { | |
863 | idat.particlePot2 = &(atomColData.particlePot[atom2]); | |
864 | } | |
865 | ||
866 | + | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
867 | + | idat.skippedCharge1 = &(atomRowData.skippedCharge[atom1]); |
868 | + | idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]); |
869 | + | } |
870 | + | |
871 | #else | |
872 | ||
873 | idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), | |
# | Line 858 | Line 908 | namespace OpenMD { | |
908 | idat.particlePot2 = &(snap_->atomData.particlePot[atom2]); | |
909 | } | |
910 | ||
911 | + | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
912 | + | idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]); |
913 | + | idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]); |
914 | + | } |
915 | #endif | |
916 | } | |
917 | ||
# | Line 875 | Line 929 | namespace OpenMD { | |
929 | snap_->atomData.force[atom1] += *(idat.f1); | |
930 | snap_->atomData.force[atom2] -= *(idat.f1); | |
931 | #endif | |
932 | < | |
932 | > | |
933 | } | |
934 | ||
881 | – | |
882 | – | void ForceMatrixDecomposition::fillSkipData(InteractionData &idat, |
883 | – | int atom1, int atom2) { |
884 | – | // Still Missing:: skippedCharge fill must be added to DataStorage |
885 | – | #ifdef IS_MPI |
886 | – | idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), |
887 | – | ff_->getAtomType(identsCol[atom2]) ); |
888 | – | |
889 | – | if (storageLayout_ & DataStorage::dslElectroFrame) { |
890 | – | idat.eFrame1 = &(atomRowData.electroFrame[atom1]); |
891 | – | idat.eFrame2 = &(atomColData.electroFrame[atom2]); |
892 | – | } |
893 | – | if (storageLayout_ & DataStorage::dslTorque) { |
894 | – | idat.t1 = &(atomRowData.torque[atom1]); |
895 | – | idat.t2 = &(atomColData.torque[atom2]); |
896 | – | } |
897 | – | #else |
898 | – | idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), |
899 | – | ff_->getAtomType(idents[atom2]) ); |
900 | – | |
901 | – | if (storageLayout_ & DataStorage::dslElectroFrame) { |
902 | – | idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]); |
903 | – | idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]); |
904 | – | } |
905 | – | if (storageLayout_ & DataStorage::dslTorque) { |
906 | – | idat.t1 = &(snap_->atomData.torque[atom1]); |
907 | – | idat.t2 = &(snap_->atomData.torque[atom2]); |
908 | – | } |
909 | – | #endif |
910 | – | } |
911 | – | |
912 | – | |
913 | – | void ForceMatrixDecomposition::unpackSkipData(InteractionData &idat, int atom1, int atom2) { |
914 | – | #ifdef IS_MPI |
915 | – | pot_row[atom1] += 0.5 * *(idat.pot); |
916 | – | pot_col[atom2] += 0.5 * *(idat.pot); |
917 | – | #else |
918 | – | pairwisePot += *(idat.pot); |
919 | – | #endif |
920 | – | |
921 | – | } |
922 | – | |
923 | – | |
935 | /* | |
936 | * buildNeighborList | |
937 | * | |
# | Line 931 | Line 942 | namespace OpenMD { | |
942 | ||
943 | vector<pair<int, int> > neighborList; | |
944 | groupCutoffs cuts; | |
945 | + | bool doAllPairs = false; |
946 | + | |
947 | #ifdef IS_MPI | |
948 | cellListRow_.clear(); | |
949 | cellListCol_.clear(); | |
# | Line 950 | Line 963 | namespace OpenMD { | |
963 | nCells_.y() = (int) ( Hy.length() )/ rList_; | |
964 | nCells_.z() = (int) ( Hz.length() )/ rList_; | |
965 | ||
966 | + | // handle small boxes where the cell offsets can end up repeating cells |
967 | + | |
968 | + | if (nCells_.x() < 3) doAllPairs = true; |
969 | + | if (nCells_.y() < 3) doAllPairs = true; |
970 | + | if (nCells_.z() < 3) doAllPairs = true; |
971 | + | |
972 | Mat3x3d invHmat = snap_->getInvHmat(); | |
973 | Vector3d rs, scaled, dr; | |
974 | Vector3i whichCell; | |
# | Line 963 | Line 982 | namespace OpenMD { | |
982 | cellList_.resize(nCtot); | |
983 | #endif | |
984 | ||
985 | + | if (!doAllPairs) { |
986 | #ifdef IS_MPI | |
967 | – | for (int i = 0; i < nGroupsInRow_; i++) { |
968 | – | rs = cgRowData.position[i]; |
987 | ||
988 | < | // scaled positions relative to the box vectors |
989 | < | scaled = invHmat * rs; |
990 | < | |
991 | < | // wrap the vector back into the unit box by subtracting integer box |
992 | < | // numbers |
993 | < | for (int j = 0; j < 3; j++) { |
994 | < | scaled[j] -= roundMe(scaled[j]); |
995 | < | scaled[j] += 0.5; |
996 | < | } |
997 | < | |
998 | < | // find xyz-indices of cell that cutoffGroup is in. |
999 | < | whichCell.x() = nCells_.x() * scaled.x(); |
1000 | < | whichCell.y() = nCells_.y() * scaled.y(); |
1001 | < | whichCell.z() = nCells_.z() * scaled.z(); |
1002 | < | |
1003 | < | // find single index of this cell: |
1004 | < | cellIndex = Vlinear(whichCell, nCells_); |
1005 | < | |
1006 | < | // add this cutoff group to the list of groups in this cell; |
1007 | < | cellListRow_[cellIndex].push_back(i); |
1008 | < | } |
1009 | < | |
1010 | < | for (int i = 0; i < nGroupsInCol_; i++) { |
1011 | < | rs = cgColData.position[i]; |
1012 | < | |
1013 | < | // scaled positions relative to the box vectors |
1014 | < | scaled = invHmat * rs; |
1015 | < | |
1016 | < | // wrap the vector back into the unit box by subtracting integer box |
1017 | < | // numbers |
1018 | < | for (int j = 0; j < 3; j++) { |
1019 | < | scaled[j] -= roundMe(scaled[j]); |
1020 | < | scaled[j] += 0.5; |
1021 | < | } |
1022 | < | |
1023 | < | // find xyz-indices of cell that cutoffGroup is in. |
1024 | < | whichCell.x() = nCells_.x() * scaled.x(); |
1025 | < | whichCell.y() = nCells_.y() * scaled.y(); |
1026 | < | whichCell.z() = nCells_.z() * scaled.z(); |
1027 | < | |
1028 | < | // find single index of this cell: |
1029 | < | cellIndex = Vlinear(whichCell, nCells_); |
1030 | < | |
1031 | < | // add this cutoff group to the list of groups in this cell; |
1032 | < | cellListCol_[cellIndex].push_back(i); |
1033 | < | } |
988 | > | for (int i = 0; i < nGroupsInRow_; i++) { |
989 | > | rs = cgRowData.position[i]; |
990 | > | |
991 | > | // scaled positions relative to the box vectors |
992 | > | scaled = invHmat * rs; |
993 | > | |
994 | > | // wrap the vector back into the unit box by subtracting integer box |
995 | > | // numbers |
996 | > | for (int j = 0; j < 3; j++) { |
997 | > | scaled[j] -= roundMe(scaled[j]); |
998 | > | scaled[j] += 0.5; |
999 | > | } |
1000 | > | |
1001 | > | // find xyz-indices of cell that cutoffGroup is in. |
1002 | > | whichCell.x() = nCells_.x() * scaled.x(); |
1003 | > | whichCell.y() = nCells_.y() * scaled.y(); |
1004 | > | whichCell.z() = nCells_.z() * scaled.z(); |
1005 | > | |
1006 | > | // find single index of this cell: |
1007 | > | cellIndex = Vlinear(whichCell, nCells_); |
1008 | > | |
1009 | > | // add this cutoff group to the list of groups in this cell; |
1010 | > | cellListRow_[cellIndex].push_back(i); |
1011 | > | } |
1012 | > | |
1013 | > | for (int i = 0; i < nGroupsInCol_; i++) { |
1014 | > | rs = cgColData.position[i]; |
1015 | > | |
1016 | > | // scaled positions relative to the box vectors |
1017 | > | scaled = invHmat * rs; |
1018 | > | |
1019 | > | // wrap the vector back into the unit box by subtracting integer box |
1020 | > | // numbers |
1021 | > | for (int j = 0; j < 3; j++) { |
1022 | > | scaled[j] -= roundMe(scaled[j]); |
1023 | > | scaled[j] += 0.5; |
1024 | > | } |
1025 | > | |
1026 | > | // find xyz-indices of cell that cutoffGroup is in. |
1027 | > | whichCell.x() = nCells_.x() * scaled.x(); |
1028 | > | whichCell.y() = nCells_.y() * scaled.y(); |
1029 | > | whichCell.z() = nCells_.z() * scaled.z(); |
1030 | > | |
1031 | > | // find single index of this cell: |
1032 | > | cellIndex = Vlinear(whichCell, nCells_); |
1033 | > | |
1034 | > | // add this cutoff group to the list of groups in this cell; |
1035 | > | cellListCol_[cellIndex].push_back(i); |
1036 | > | } |
1037 | #else | |
1038 | < | for (int i = 0; i < nGroups_; i++) { |
1039 | < | rs = snap_->cgData.position[i]; |
1040 | < | |
1041 | < | // scaled positions relative to the box vectors |
1042 | < | scaled = invHmat * rs; |
1043 | < | |
1044 | < | // wrap the vector back into the unit box by subtracting integer box |
1045 | < | // numbers |
1046 | < | for (int j = 0; j < 3; j++) { |
1047 | < | scaled[j] -= roundMe(scaled[j]); |
1048 | < | scaled[j] += 0.5; |
1038 | > | for (int i = 0; i < nGroups_; i++) { |
1039 | > | rs = snap_->cgData.position[i]; |
1040 | > | |
1041 | > | // scaled positions relative to the box vectors |
1042 | > | scaled = invHmat * rs; |
1043 | > | |
1044 | > | // wrap the vector back into the unit box by subtracting integer box |
1045 | > | // numbers |
1046 | > | for (int j = 0; j < 3; j++) { |
1047 | > | scaled[j] -= roundMe(scaled[j]); |
1048 | > | scaled[j] += 0.5; |
1049 | > | } |
1050 | > | |
1051 | > | // find xyz-indices of cell that cutoffGroup is in. |
1052 | > | whichCell.x() = nCells_.x() * scaled.x(); |
1053 | > | whichCell.y() = nCells_.y() * scaled.y(); |
1054 | > | whichCell.z() = nCells_.z() * scaled.z(); |
1055 | > | |
1056 | > | // find single index of this cell: |
1057 | > | cellIndex = Vlinear(whichCell, nCells_); |
1058 | > | |
1059 | > | // add this cutoff group to the list of groups in this cell; |
1060 | > | cellList_[cellIndex].push_back(i); |
1061 | } | |
1029 | – | |
1030 | – | // find xyz-indices of cell that cutoffGroup is in. |
1031 | – | whichCell.x() = nCells_.x() * scaled.x(); |
1032 | – | whichCell.y() = nCells_.y() * scaled.y(); |
1033 | – | whichCell.z() = nCells_.z() * scaled.z(); |
1034 | – | |
1035 | – | // find single index of this cell: |
1036 | – | cellIndex = Vlinear(whichCell, nCells_); |
1037 | – | |
1038 | – | // add this cutoff group to the list of groups in this cell; |
1039 | – | cellList_[cellIndex].push_back(i); |
1040 | – | } |
1062 | #endif | |
1063 | ||
1064 | < | for (int m1z = 0; m1z < nCells_.z(); m1z++) { |
1065 | < | for (int m1y = 0; m1y < nCells_.y(); m1y++) { |
1066 | < | for (int m1x = 0; m1x < nCells_.x(); m1x++) { |
1067 | < | Vector3i m1v(m1x, m1y, m1z); |
1068 | < | int m1 = Vlinear(m1v, nCells_); |
1048 | < | |
1049 | < | for (vector<Vector3i>::iterator os = cellOffsets_.begin(); |
1050 | < | os != cellOffsets_.end(); ++os) { |
1064 | > | for (int m1z = 0; m1z < nCells_.z(); m1z++) { |
1065 | > | for (int m1y = 0; m1y < nCells_.y(); m1y++) { |
1066 | > | for (int m1x = 0; m1x < nCells_.x(); m1x++) { |
1067 | > | Vector3i m1v(m1x, m1y, m1z); |
1068 | > | int m1 = Vlinear(m1v, nCells_); |
1069 | ||
1070 | < | Vector3i m2v = m1v + (*os); |
1071 | < | |
1072 | < | if (m2v.x() >= nCells_.x()) { |
1073 | < | m2v.x() = 0; |
1074 | < | } else if (m2v.x() < 0) { |
1075 | < | m2v.x() = nCells_.x() - 1; |
1076 | < | } |
1077 | < | |
1078 | < | if (m2v.y() >= nCells_.y()) { |
1079 | < | m2v.y() = 0; |
1080 | < | } else if (m2v.y() < 0) { |
1081 | < | m2v.y() = nCells_.y() - 1; |
1082 | < | } |
1083 | < | |
1084 | < | if (m2v.z() >= nCells_.z()) { |
1085 | < | m2v.z() = 0; |
1086 | < | } else if (m2v.z() < 0) { |
1087 | < | m2v.z() = nCells_.z() - 1; |
1088 | < | } |
1089 | < | |
1090 | < | int m2 = Vlinear (m2v, nCells_); |
1091 | < | |
1070 | > | for (vector<Vector3i>::iterator os = cellOffsets_.begin(); |
1071 | > | os != cellOffsets_.end(); ++os) { |
1072 | > | |
1073 | > | Vector3i m2v = m1v + (*os); |
1074 | > | |
1075 | > | if (m2v.x() >= nCells_.x()) { |
1076 | > | m2v.x() = 0; |
1077 | > | } else if (m2v.x() < 0) { |
1078 | > | m2v.x() = nCells_.x() - 1; |
1079 | > | } |
1080 | > | |
1081 | > | if (m2v.y() >= nCells_.y()) { |
1082 | > | m2v.y() = 0; |
1083 | > | } else if (m2v.y() < 0) { |
1084 | > | m2v.y() = nCells_.y() - 1; |
1085 | > | } |
1086 | > | |
1087 | > | if (m2v.z() >= nCells_.z()) { |
1088 | > | m2v.z() = 0; |
1089 | > | } else if (m2v.z() < 0) { |
1090 | > | m2v.z() = nCells_.z() - 1; |
1091 | > | } |
1092 | > | |
1093 | > | int m2 = Vlinear (m2v, nCells_); |
1094 | > | |
1095 | #ifdef IS_MPI | |
1096 | < | for (vector<int>::iterator j1 = cellListRow_[m1].begin(); |
1097 | < | j1 != cellListRow_[m1].end(); ++j1) { |
1098 | < | for (vector<int>::iterator j2 = cellListCol_[m2].begin(); |
1099 | < | j2 != cellListCol_[m2].end(); ++j2) { |
1100 | < | |
1101 | < | // Always do this if we're in different cells or if |
1102 | < | // we're in the same cell and the global index of the |
1103 | < | // j2 cutoff group is less than the j1 cutoff group |
1104 | < | |
1105 | < | if (m2 != m1 || cgColToGlobal[(*j2)] < cgRowToGlobal[(*j1)]) { |
1106 | < | dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; |
1107 | < | snap_->wrapVector(dr); |
1108 | < | cuts = getGroupCutoffs( (*j1), (*j2) ); |
1109 | < | if (dr.lengthSquare() < cuts.third) { |
1110 | < | neighborList.push_back(make_pair((*j1), (*j2))); |
1096 | > | for (vector<int>::iterator j1 = cellListRow_[m1].begin(); |
1097 | > | j1 != cellListRow_[m1].end(); ++j1) { |
1098 | > | for (vector<int>::iterator j2 = cellListCol_[m2].begin(); |
1099 | > | j2 != cellListCol_[m2].end(); ++j2) { |
1100 | > | |
1101 | > | // Always do this if we're in different cells or if |
1102 | > | // we're in the same cell and the global index of the |
1103 | > | // j2 cutoff group is less than the j1 cutoff group |
1104 | > | |
1105 | > | if (m2 != m1 || cgColToGlobal[(*j2)] < cgRowToGlobal[(*j1)]) { |
1106 | > | dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; |
1107 | > | snap_->wrapVector(dr); |
1108 | > | cuts = getGroupCutoffs( (*j1), (*j2) ); |
1109 | > | if (dr.lengthSquare() < cuts.third) { |
1110 | > | neighborList.push_back(make_pair((*j1), (*j2))); |
1111 | > | } |
1112 | } | |
1113 | } | |
1114 | } | |
1093 | – | } |
1115 | #else | |
1116 | < | |
1117 | < | for (vector<int>::iterator j1 = cellList_[m1].begin(); |
1118 | < | j1 != cellList_[m1].end(); ++j1) { |
1119 | < | for (vector<int>::iterator j2 = cellList_[m2].begin(); |
1120 | < | j2 != cellList_[m2].end(); ++j2) { |
1121 | < | |
1122 | < | // Always do this if we're in different cells or if |
1123 | < | // we're in the same cell and the global index of the |
1124 | < | // j2 cutoff group is less than the j1 cutoff group |
1125 | < | |
1126 | < | if (m2 != m1 || (*j2) < (*j1)) { |
1127 | < | dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)]; |
1128 | < | snap_->wrapVector(dr); |
1129 | < | cuts = getGroupCutoffs( (*j1), (*j2) ); |
1130 | < | if (dr.lengthSquare() < cuts.third) { |
1131 | < | neighborList.push_back(make_pair((*j1), (*j2))); |
1116 | > | |
1117 | > | for (vector<int>::iterator j1 = cellList_[m1].begin(); |
1118 | > | j1 != cellList_[m1].end(); ++j1) { |
1119 | > | for (vector<int>::iterator j2 = cellList_[m2].begin(); |
1120 | > | j2 != cellList_[m2].end(); ++j2) { |
1121 | > | |
1122 | > | // Always do this if we're in different cells or if |
1123 | > | // we're in the same cell and the global index of the |
1124 | > | // j2 cutoff group is less than the j1 cutoff group |
1125 | > | |
1126 | > | if (m2 != m1 || (*j2) < (*j1)) { |
1127 | > | dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)]; |
1128 | > | snap_->wrapVector(dr); |
1129 | > | cuts = getGroupCutoffs( (*j1), (*j2) ); |
1130 | > | if (dr.lengthSquare() < cuts.third) { |
1131 | > | neighborList.push_back(make_pair((*j1), (*j2))); |
1132 | > | } |
1133 | } | |
1134 | } | |
1135 | } | |
1114 | – | } |
1136 | #endif | |
1137 | + | } |
1138 | } | |
1139 | } | |
1140 | } | |
1141 | + | } else { |
1142 | + | // branch to do all cutoff group pairs |
1143 | + | #ifdef IS_MPI |
1144 | + | for (int j1 = 0; j1 < nGroupsInRow_; j1++) { |
1145 | + | for (int j2 = 0; j2 < nGroupsInCol_; j2++) { |
1146 | + | dr = cgColData.position[j2] - cgRowData.position[j1]; |
1147 | + | snap_->wrapVector(dr); |
1148 | + | cuts = getGroupCutoffs( j1, j2 ); |
1149 | + | if (dr.lengthSquare() < cuts.third) { |
1150 | + | neighborList.push_back(make_pair(j1, j2)); |
1151 | + | } |
1152 | + | } |
1153 | + | } |
1154 | + | #else |
1155 | + | for (int j1 = 0; j1 < nGroups_ - 1; j1++) { |
1156 | + | for (int j2 = j1 + 1; j2 < nGroups_; j2++) { |
1157 | + | dr = snap_->cgData.position[j2] - snap_->cgData.position[j1]; |
1158 | + | snap_->wrapVector(dr); |
1159 | + | cuts = getGroupCutoffs( j1, j2 ); |
1160 | + | if (dr.lengthSquare() < cuts.third) { |
1161 | + | neighborList.push_back(make_pair(j1, j2)); |
1162 | + | } |
1163 | + | } |
1164 | + | } |
1165 | + | #endif |
1166 | } | |
1167 | < | |
1167 | > | |
1168 | // save the local cutoff group positions for the check that is | |
1169 | // done on each loop: | |
1170 | saved_CG_positions_.clear(); | |
1171 | for (int i = 0; i < nGroups_; i++) | |
1172 | saved_CG_positions_.push_back(snap_->cgData.position[i]); | |
1173 | < | |
1173 | > | |
1174 | return neighborList; | |
1175 | } | |
1176 | } //end namespace OpenMD |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |