# | 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 | // gather the information for atomtype IDs (atids): | |
63 | idents = info_->getIdentArray(); | |
# | Line 67 | Line 67 | namespace OpenMD { | |
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(); |
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 | ||
# | Line 143 | Line 143 | namespace OpenMD { | |
143 | } | |
144 | } | |
145 | ||
146 | < | skipsForAtom.clear(); |
147 | < | skipsForAtom.resize(nAtomsInRow_); |
146 | > | excludesForAtom.clear(); |
147 | > | excludesForAtom.resize(nAtomsInRow_); |
148 | toposForAtom.clear(); | |
149 | toposForAtom.resize(nAtomsInRow_); | |
150 | topoDist.clear(); | |
# | Line 155 | Line 155 | namespace OpenMD { | |
155 | for (int j = 0; j < nAtomsInCol_; j++) { | |
156 | int jglob = AtomColToGlobal[j]; | |
157 | ||
158 | < | if (excludes.hasPair(iglob, jglob)) |
159 | < | skipsForAtom[i].push_back(j); |
158 | > | if (excludes->hasPair(iglob, jglob)) |
159 | > | excludesForAtom[i].push_back(j); |
160 | ||
161 | < | if (oneTwo.hasPair(iglob, jglob)) { |
161 | > | if (oneTwo->hasPair(iglob, jglob)) { |
162 | toposForAtom[i].push_back(j); | |
163 | topoDist[i].push_back(1); | |
164 | } else { | |
165 | < | if (oneThree.hasPair(iglob, jglob)) { |
165 | > | if (oneThree->hasPair(iglob, jglob)) { |
166 | toposForAtom[i].push_back(j); | |
167 | topoDist[i].push_back(2); | |
168 | } else { | |
169 | < | if (oneFour.hasPair(iglob, jglob)) { |
169 | > | if (oneFour->hasPair(iglob, jglob)) { |
170 | toposForAtom[i].push_back(j); | |
171 | topoDist[i].push_back(3); | |
172 | } | |
# | Line 189 | Line 189 | namespace OpenMD { | |
189 | } | |
190 | } | |
191 | ||
192 | < | skipsForAtom.clear(); |
193 | < | skipsForAtom.resize(nLocal_); |
192 | > | excludesForAtom.clear(); |
193 | > | excludesForAtom.resize(nLocal_); |
194 | toposForAtom.clear(); | |
195 | toposForAtom.resize(nLocal_); | |
196 | topoDist.clear(); | |
# | Line 202 | Line 202 | namespace OpenMD { | |
202 | for (int j = 0; j < nLocal_; j++) { | |
203 | int jglob = AtomLocalToGlobal[j]; | |
204 | ||
205 | < | if (excludes.hasPair(iglob, jglob)) |
206 | < | skipsForAtom[i].push_back(j); |
205 | > | if (excludes->hasPair(iglob, jglob)) |
206 | > | excludesForAtom[i].push_back(j); |
207 | ||
208 | < | if (oneTwo.hasPair(iglob, jglob)) { |
208 | > | if (oneTwo->hasPair(iglob, jglob)) { |
209 | toposForAtom[i].push_back(j); | |
210 | topoDist[i].push_back(1); | |
211 | } else { | |
212 | < | if (oneThree.hasPair(iglob, jglob)) { |
212 | > | if (oneThree->hasPair(iglob, jglob)) { |
213 | toposForAtom[i].push_back(j); | |
214 | topoDist[i].push_back(2); | |
215 | } else { | |
216 | < | if (oneFour.hasPair(iglob, jglob)) { |
216 | > | if (oneFour->hasPair(iglob, jglob)) { |
217 | toposForAtom[i].push_back(j); | |
218 | topoDist[i].push_back(3); | |
219 | } | |
# | Line 223 | Line 223 | namespace OpenMD { | |
223 | } | |
224 | ||
225 | createGtypeCutoffMap(); | |
226 | + | |
227 | } | |
228 | ||
229 | void ForceMatrixDecomposition::createGtypeCutoffMap() { | |
# | Line 231 | Line 232 | namespace OpenMD { | |
232 | RealType rc; | |
233 | int atid; | |
234 | set<AtomType*> atypes = info_->getSimulatedAtomTypes(); | |
235 | < | vector<RealType> atypeCutoff; |
235 | < | atypeCutoff.resize( atypes.size() ); |
235 | > | map<int, RealType> atypeCutoff; |
236 | ||
237 | for (set<AtomType*>::iterator at = atypes.begin(); | |
238 | at != atypes.end(); ++at){ | |
239 | atid = (*at)->getIdent(); | |
240 | < | |
241 | < | if (userChoseCutoff_) |
240 | > | if (userChoseCutoff_) |
241 | atypeCutoff[atid] = userCutoff_; | |
242 | < | else |
242 | > | else |
243 | atypeCutoff[atid] = interactionMan_->getSuggestedCutoffRadius(*at); | |
244 | } | |
245 | ||
246 | vector<RealType> gTypeCutoffs; | |
248 | – | |
247 | // first we do a single loop over the cutoff groups to find the | |
248 | // largest cutoff for any atypes present in this group. | |
249 | #ifdef IS_MPI | |
# | Line 303 | Line 301 | namespace OpenMD { | |
301 | ||
302 | vector<RealType> groupCutoff(nGroups_, 0.0); | |
303 | groupToGtype.resize(nGroups_); | |
306 | – | |
304 | for (int cg1 = 0; cg1 < nGroups_; cg1++) { | |
305 | ||
306 | groupCutoff[cg1] = 0.0; | |
# | Line 455 | Line 452 | namespace OpenMD { | |
452 | } | |
453 | ||
454 | if (storageLayout_ & DataStorage::dslSkippedCharge) { | |
455 | < | fill(atomRowData.skippedCharge.begin(), atomRowData.skippedCharge.end(), 0.0); |
456 | < | fill(atomColData.skippedCharge.begin(), atomColData.skippedCharge.end(), 0.0); |
455 | > | fill(atomRowData.skippedCharge.begin(), |
456 | > | atomRowData.skippedCharge.end(), 0.0); |
457 | > | fill(atomColData.skippedCharge.begin(), |
458 | > | atomColData.skippedCharge.end(), 0.0); |
459 | } | |
460 | ||
461 | #else | |
# | Line 589 | Line 588 | namespace OpenMD { | |
588 | ||
589 | if (storageLayout_ & DataStorage::dslTorque) { | |
590 | ||
591 | < | int nt = snap_->atomData.force.size(); |
591 | > | int nt = snap_->atomData.torque.size(); |
592 | vector<Vector3d> trq_tmp(nt, V3Zero); | |
593 | ||
594 | AtomCommVectorRow->scatter(atomRowData.torque, trq_tmp); | |
595 | < | for (int i = 0; i < n; i++) { |
595 | > | for (int i = 0; i < nt; i++) { |
596 | snap_->atomData.torque[i] += trq_tmp[i]; | |
597 | trq_tmp[i] = 0.0; | |
598 | } | |
599 | ||
600 | AtomCommVectorColumn->scatter(atomColData.torque, trq_tmp); | |
601 | < | for (int i = 0; i < n; i++) |
601 | > | for (int i = 0; i < nt; i++) |
602 | snap_->atomData.torque[i] += trq_tmp[i]; | |
603 | + | } |
604 | + | |
605 | + | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
606 | + | |
607 | + | int ns = snap_->atomData.skippedCharge.size(); |
608 | + | vector<RealType> skch_tmp(ns, 0.0); |
609 | + | |
610 | + | AtomCommRealRow->scatter(atomRowData.skippedCharge, skch_tmp); |
611 | + | for (int i = 0; i < ns; i++) { |
612 | + | snap_->atomData.skippedCharge[i] = skch_tmp[i]; |
613 | + | skch_tmp[i] = 0.0; |
614 | + | } |
615 | + | |
616 | + | AtomCommRealColumn->scatter(atomColData.skippedCharge, skch_tmp); |
617 | + | for (int i = 0; i < ns; i++) |
618 | + | snap_->atomData.skippedCharge[i] += skch_tmp[i]; |
619 | } | |
620 | ||
621 | nLocal_ = snap_->getNumberOfAtoms(); | |
# | Line 724 | Line 739 | namespace OpenMD { | |
739 | return d; | |
740 | } | |
741 | ||
742 | < | vector<int> ForceMatrixDecomposition::getSkipsForAtom(int atom1) { |
743 | < | return skipsForAtom[atom1]; |
742 | > | vector<int> ForceMatrixDecomposition::getExcludesForAtom(int atom1) { |
743 | > | return excludesForAtom[atom1]; |
744 | } | |
745 | ||
746 | /** | |
747 | < | * There are a number of reasons to skip a pair or a |
733 | < | * particle. Mostly we do this to exclude atoms who are involved in |
734 | < | * short range interactions (bonds, bends, torsions), but we also |
735 | < | * need to exclude some overcounted interactions that result from |
747 | > | * We need to exclude some overcounted interactions that result from |
748 | * the parallel decomposition. | |
749 | */ | |
750 | bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) { | |
# | Line 752 | Line 764 | namespace OpenMD { | |
764 | } else { | |
765 | if ((unique_id_1 + unique_id_2) % 2 == 1) return true; | |
766 | } | |
767 | + | #endif |
768 | + | return false; |
769 | + | } |
770 | + | |
771 | + | /** |
772 | + | * We need to handle the interactions for atoms who are involved in |
773 | + | * the same rigid body as well as some short range interactions |
774 | + | * (bonds, bends, torsions) differently from other interactions. |
775 | + | * We'll still visit the pairwise routines, but with a flag that |
776 | + | * tells those routines to exclude the pair from direct long range |
777 | + | * interactions. Some indirect interactions (notably reaction |
778 | + | * field) must still be handled for these pairs. |
779 | + | */ |
780 | + | bool ForceMatrixDecomposition::excludeAtomPair(int atom1, int atom2) { |
781 | + | int unique_id_2; |
782 | + | |
783 | + | #ifdef IS_MPI |
784 | + | // in MPI, we have to look up the unique IDs for the row atom. |
785 | + | unique_id_2 = AtomColToGlobal[atom2]; |
786 | #else | |
787 | // in the normal loop, the atom numbers are unique | |
757 | – | unique_id_1 = atom1; |
788 | unique_id_2 = atom2; | |
789 | #endif | |
790 | ||
791 | < | for (vector<int>::iterator i = skipsForAtom[atom1].begin(); |
792 | < | i != skipsForAtom[atom1].end(); ++i) { |
791 | > | for (vector<int>::iterator i = excludesForAtom[atom1].begin(); |
792 | > | i != excludesForAtom[atom1].end(); ++i) { |
793 | if ( (*i) == unique_id_2 ) return true; | |
794 | } | |
795 | ||
# | Line 785 | Line 815 | namespace OpenMD { | |
815 | ||
816 | // filling interaction blocks with pointers | |
817 | void ForceMatrixDecomposition::fillInteractionData(InteractionData &idat, | |
818 | < | int atom1, int atom2) { |
818 | > | int atom1, int atom2) { |
819 | > | |
820 | > | idat.excluded = excludeAtomPair(atom1, atom2); |
821 | > | |
822 | #ifdef IS_MPI | |
823 | ||
824 | idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), | |
# | Line 826 | Line 859 | namespace OpenMD { | |
859 | idat.particlePot2 = &(atomColData.particlePot[atom2]); | |
860 | } | |
861 | ||
862 | + | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
863 | + | idat.skippedCharge1 = &(atomRowData.skippedCharge[atom1]); |
864 | + | idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]); |
865 | + | } |
866 | + | |
867 | #else | |
868 | ||
869 | idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), | |
# | Line 866 | Line 904 | namespace OpenMD { | |
904 | idat.particlePot2 = &(snap_->atomData.particlePot[atom2]); | |
905 | } | |
906 | ||
907 | + | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
908 | + | idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]); |
909 | + | idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]); |
910 | + | } |
911 | #endif | |
912 | } | |
913 | ||
# | Line 886 | Line 928 | namespace OpenMD { | |
928 | ||
929 | } | |
930 | ||
889 | – | |
890 | – | void ForceMatrixDecomposition::fillSkipData(InteractionData &idat, |
891 | – | int atom1, int atom2) { |
892 | – | #ifdef IS_MPI |
893 | – | idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), |
894 | – | ff_->getAtomType(identsCol[atom2]) ); |
895 | – | |
896 | – | if (storageLayout_ & DataStorage::dslElectroFrame) { |
897 | – | idat.eFrame1 = &(atomRowData.electroFrame[atom1]); |
898 | – | idat.eFrame2 = &(atomColData.electroFrame[atom2]); |
899 | – | } |
900 | – | |
901 | – | if (storageLayout_ & DataStorage::dslTorque) { |
902 | – | idat.t1 = &(atomRowData.torque[atom1]); |
903 | – | idat.t2 = &(atomColData.torque[atom2]); |
904 | – | } |
905 | – | |
906 | – | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
907 | – | idat.skippedCharge1 = &(atomRowData.skippedCharge[atom1]); |
908 | – | idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]); |
909 | – | } |
910 | – | #else |
911 | – | idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), |
912 | – | ff_->getAtomType(idents[atom2]) ); |
913 | – | |
914 | – | if (storageLayout_ & DataStorage::dslElectroFrame) { |
915 | – | idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]); |
916 | – | idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]); |
917 | – | } |
918 | – | |
919 | – | if (storageLayout_ & DataStorage::dslTorque) { |
920 | – | idat.t1 = &(snap_->atomData.torque[atom1]); |
921 | – | idat.t2 = &(snap_->atomData.torque[atom2]); |
922 | – | } |
923 | – | |
924 | – | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
925 | – | idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]); |
926 | – | idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]); |
927 | – | } |
928 | – | #endif |
929 | – | } |
930 | – | |
931 | – | |
932 | – | void ForceMatrixDecomposition::unpackSkipData(InteractionData &idat, int atom1, int atom2) { |
933 | – | #ifdef IS_MPI |
934 | – | pot_row[atom1] += 0.5 * *(idat.pot); |
935 | – | pot_col[atom2] += 0.5 * *(idat.pot); |
936 | – | #else |
937 | – | pairwisePot += *(idat.pot); |
938 | – | #endif |
939 | – | |
940 | – | } |
941 | – | |
942 | – | |
931 | /* | |
932 | * buildNeighborList | |
933 | * | |
# | Line 950 | Line 938 | namespace OpenMD { | |
938 | ||
939 | vector<pair<int, int> > neighborList; | |
940 | groupCutoffs cuts; | |
941 | + | bool doAllPairs = false; |
942 | + | |
943 | #ifdef IS_MPI | |
944 | cellListRow_.clear(); | |
945 | cellListCol_.clear(); | |
# | Line 969 | Line 959 | namespace OpenMD { | |
959 | nCells_.y() = (int) ( Hy.length() )/ rList_; | |
960 | nCells_.z() = (int) ( Hz.length() )/ rList_; | |
961 | ||
962 | + | // handle small boxes where the cell offsets can end up repeating cells |
963 | + | |
964 | + | if (nCells_.x() < 3) doAllPairs = true; |
965 | + | if (nCells_.y() < 3) doAllPairs = true; |
966 | + | if (nCells_.z() < 3) doAllPairs = true; |
967 | + | |
968 | Mat3x3d invHmat = snap_->getInvHmat(); | |
969 | Vector3d rs, scaled, dr; | |
970 | Vector3i whichCell; | |
# | Line 982 | Line 978 | namespace OpenMD { | |
978 | cellList_.resize(nCtot); | |
979 | #endif | |
980 | ||
981 | + | if (!doAllPairs) { |
982 | #ifdef IS_MPI | |
986 | – | for (int i = 0; i < nGroupsInRow_; i++) { |
987 | – | rs = cgRowData.position[i]; |
983 | ||
984 | < | // scaled positions relative to the box vectors |
985 | < | scaled = invHmat * rs; |
986 | < | |
987 | < | // wrap the vector back into the unit box by subtracting integer box |
988 | < | // numbers |
989 | < | for (int j = 0; j < 3; j++) { |
990 | < | scaled[j] -= roundMe(scaled[j]); |
991 | < | scaled[j] += 0.5; |
984 | > | for (int i = 0; i < nGroupsInRow_; i++) { |
985 | > | rs = cgRowData.position[i]; |
986 | > | |
987 | > | // scaled positions relative to the box vectors |
988 | > | scaled = invHmat * rs; |
989 | > | |
990 | > | // wrap the vector back into the unit box by subtracting integer box |
991 | > | // numbers |
992 | > | for (int j = 0; j < 3; j++) { |
993 | > | scaled[j] -= roundMe(scaled[j]); |
994 | > | scaled[j] += 0.5; |
995 | > | } |
996 | > | |
997 | > | // find xyz-indices of cell that cutoffGroup is in. |
998 | > | whichCell.x() = nCells_.x() * scaled.x(); |
999 | > | whichCell.y() = nCells_.y() * scaled.y(); |
1000 | > | whichCell.z() = nCells_.z() * scaled.z(); |
1001 | > | |
1002 | > | // find single index of this cell: |
1003 | > | cellIndex = Vlinear(whichCell, nCells_); |
1004 | > | |
1005 | > | // add this cutoff group to the list of groups in this cell; |
1006 | > | cellListRow_[cellIndex].push_back(i); |
1007 | } | |
1008 | < | |
1009 | < | // find xyz-indices of cell that cutoffGroup is in. |
1010 | < | whichCell.x() = nCells_.x() * scaled.x(); |
1011 | < | whichCell.y() = nCells_.y() * scaled.y(); |
1012 | < | whichCell.z() = nCells_.z() * scaled.z(); |
1013 | < | |
1014 | < | // find single index of this cell: |
1015 | < | cellIndex = Vlinear(whichCell, nCells_); |
1016 | < | |
1017 | < | // add this cutoff group to the list of groups in this cell; |
1018 | < | cellListRow_[cellIndex].push_back(i); |
1019 | < | } |
1020 | < | |
1021 | < | for (int i = 0; i < nGroupsInCol_; i++) { |
1022 | < | rs = cgColData.position[i]; |
1023 | < | |
1024 | < | // scaled positions relative to the box vectors |
1025 | < | scaled = invHmat * rs; |
1026 | < | |
1027 | < | // wrap the vector back into the unit box by subtracting integer box |
1028 | < | // numbers |
1029 | < | for (int j = 0; j < 3; j++) { |
1030 | < | scaled[j] -= roundMe(scaled[j]); |
1031 | < | scaled[j] += 0.5; |
1008 | > | |
1009 | > | for (int i = 0; i < nGroupsInCol_; i++) { |
1010 | > | rs = cgColData.position[i]; |
1011 | > | |
1012 | > | // scaled positions relative to the box vectors |
1013 | > | scaled = invHmat * rs; |
1014 | > | |
1015 | > | // wrap the vector back into the unit box by subtracting integer box |
1016 | > | // numbers |
1017 | > | for (int j = 0; j < 3; j++) { |
1018 | > | scaled[j] -= roundMe(scaled[j]); |
1019 | > | scaled[j] += 0.5; |
1020 | > | } |
1021 | > | |
1022 | > | // find xyz-indices of cell that cutoffGroup is in. |
1023 | > | whichCell.x() = nCells_.x() * scaled.x(); |
1024 | > | whichCell.y() = nCells_.y() * scaled.y(); |
1025 | > | whichCell.z() = nCells_.z() * scaled.z(); |
1026 | > | |
1027 | > | // find single index of this cell: |
1028 | > | cellIndex = Vlinear(whichCell, nCells_); |
1029 | > | |
1030 | > | // add this cutoff group to the list of groups in this cell; |
1031 | > | cellListCol_[cellIndex].push_back(i); |
1032 | } | |
1023 | – | |
1024 | – | // find xyz-indices of cell that cutoffGroup is in. |
1025 | – | whichCell.x() = nCells_.x() * scaled.x(); |
1026 | – | whichCell.y() = nCells_.y() * scaled.y(); |
1027 | – | whichCell.z() = nCells_.z() * scaled.z(); |
1028 | – | |
1029 | – | // find single index of this cell: |
1030 | – | cellIndex = Vlinear(whichCell, nCells_); |
1031 | – | |
1032 | – | // add this cutoff group to the list of groups in this cell; |
1033 | – | cellListCol_[cellIndex].push_back(i); |
1034 | – | } |
1033 | #else | |
1034 | < | for (int i = 0; i < nGroups_; i++) { |
1035 | < | rs = snap_->cgData.position[i]; |
1036 | < | |
1037 | < | // scaled positions relative to the box vectors |
1038 | < | scaled = invHmat * rs; |
1039 | < | |
1040 | < | // wrap the vector back into the unit box by subtracting integer box |
1041 | < | // numbers |
1042 | < | for (int j = 0; j < 3; j++) { |
1043 | < | scaled[j] -= roundMe(scaled[j]); |
1044 | < | scaled[j] += 0.5; |
1034 | > | for (int i = 0; i < nGroups_; i++) { |
1035 | > | rs = snap_->cgData.position[i]; |
1036 | > | |
1037 | > | // scaled positions relative to the box vectors |
1038 | > | scaled = invHmat * rs; |
1039 | > | |
1040 | > | // wrap the vector back into the unit box by subtracting integer box |
1041 | > | // numbers |
1042 | > | for (int j = 0; j < 3; j++) { |
1043 | > | scaled[j] -= roundMe(scaled[j]); |
1044 | > | scaled[j] += 0.5; |
1045 | > | } |
1046 | > | |
1047 | > | // find xyz-indices of cell that cutoffGroup is in. |
1048 | > | whichCell.x() = nCells_.x() * scaled.x(); |
1049 | > | whichCell.y() = nCells_.y() * scaled.y(); |
1050 | > | whichCell.z() = nCells_.z() * scaled.z(); |
1051 | > | |
1052 | > | // find single index of this cell: |
1053 | > | cellIndex = Vlinear(whichCell, nCells_); |
1054 | > | |
1055 | > | // add this cutoff group to the list of groups in this cell; |
1056 | > | cellList_[cellIndex].push_back(i); |
1057 | } | |
1048 | – | |
1049 | – | // find xyz-indices of cell that cutoffGroup is in. |
1050 | – | whichCell.x() = nCells_.x() * scaled.x(); |
1051 | – | whichCell.y() = nCells_.y() * scaled.y(); |
1052 | – | whichCell.z() = nCells_.z() * scaled.z(); |
1053 | – | |
1054 | – | // find single index of this cell: |
1055 | – | cellIndex = Vlinear(whichCell, nCells_); |
1056 | – | |
1057 | – | // add this cutoff group to the list of groups in this cell; |
1058 | – | cellList_[cellIndex].push_back(i); |
1059 | – | } |
1058 | #endif | |
1059 | ||
1060 | < | for (int m1z = 0; m1z < nCells_.z(); m1z++) { |
1061 | < | for (int m1y = 0; m1y < nCells_.y(); m1y++) { |
1062 | < | for (int m1x = 0; m1x < nCells_.x(); m1x++) { |
1063 | < | Vector3i m1v(m1x, m1y, m1z); |
1064 | < | int m1 = Vlinear(m1v, nCells_); |
1067 | < | |
1068 | < | for (vector<Vector3i>::iterator os = cellOffsets_.begin(); |
1069 | < | os != cellOffsets_.end(); ++os) { |
1060 | > | for (int m1z = 0; m1z < nCells_.z(); m1z++) { |
1061 | > | for (int m1y = 0; m1y < nCells_.y(); m1y++) { |
1062 | > | for (int m1x = 0; m1x < nCells_.x(); m1x++) { |
1063 | > | Vector3i m1v(m1x, m1y, m1z); |
1064 | > | int m1 = Vlinear(m1v, nCells_); |
1065 | ||
1066 | < | Vector3i m2v = m1v + (*os); |
1067 | < | |
1068 | < | if (m2v.x() >= nCells_.x()) { |
1069 | < | m2v.x() = 0; |
1070 | < | } else if (m2v.x() < 0) { |
1071 | < | m2v.x() = nCells_.x() - 1; |
1072 | < | } |
1073 | < | |
1074 | < | if (m2v.y() >= nCells_.y()) { |
1075 | < | m2v.y() = 0; |
1076 | < | } else if (m2v.y() < 0) { |
1077 | < | m2v.y() = nCells_.y() - 1; |
1078 | < | } |
1079 | < | |
1080 | < | if (m2v.z() >= nCells_.z()) { |
1081 | < | m2v.z() = 0; |
1082 | < | } else if (m2v.z() < 0) { |
1083 | < | m2v.z() = nCells_.z() - 1; |
1084 | < | } |
1085 | < | |
1086 | < | int m2 = Vlinear (m2v, nCells_); |
1087 | < | |
1066 | > | for (vector<Vector3i>::iterator os = cellOffsets_.begin(); |
1067 | > | os != cellOffsets_.end(); ++os) { |
1068 | > | |
1069 | > | Vector3i m2v = m1v + (*os); |
1070 | > | |
1071 | > | if (m2v.x() >= nCells_.x()) { |
1072 | > | m2v.x() = 0; |
1073 | > | } else if (m2v.x() < 0) { |
1074 | > | m2v.x() = nCells_.x() - 1; |
1075 | > | } |
1076 | > | |
1077 | > | if (m2v.y() >= nCells_.y()) { |
1078 | > | m2v.y() = 0; |
1079 | > | } else if (m2v.y() < 0) { |
1080 | > | m2v.y() = nCells_.y() - 1; |
1081 | > | } |
1082 | > | |
1083 | > | if (m2v.z() >= nCells_.z()) { |
1084 | > | m2v.z() = 0; |
1085 | > | } else if (m2v.z() < 0) { |
1086 | > | m2v.z() = nCells_.z() - 1; |
1087 | > | } |
1088 | > | |
1089 | > | int m2 = Vlinear (m2v, nCells_); |
1090 | > | |
1091 | #ifdef IS_MPI | |
1092 | < | for (vector<int>::iterator j1 = cellListRow_[m1].begin(); |
1093 | < | j1 != cellListRow_[m1].end(); ++j1) { |
1094 | < | for (vector<int>::iterator j2 = cellListCol_[m2].begin(); |
1095 | < | j2 != cellListCol_[m2].end(); ++j2) { |
1096 | < | |
1097 | < | // Always do this if we're in different cells or if |
1098 | < | // we're in the same cell and the global index of the |
1099 | < | // j2 cutoff group is less than the j1 cutoff group |
1100 | < | |
1101 | < | if (m2 != m1 || cgColToGlobal[(*j2)] < cgRowToGlobal[(*j1)]) { |
1102 | < | dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; |
1103 | < | snap_->wrapVector(dr); |
1104 | < | cuts = getGroupCutoffs( (*j1), (*j2) ); |
1105 | < | if (dr.lengthSquare() < cuts.third) { |
1106 | < | neighborList.push_back(make_pair((*j1), (*j2))); |
1092 | > | for (vector<int>::iterator j1 = cellListRow_[m1].begin(); |
1093 | > | j1 != cellListRow_[m1].end(); ++j1) { |
1094 | > | for (vector<int>::iterator j2 = cellListCol_[m2].begin(); |
1095 | > | j2 != cellListCol_[m2].end(); ++j2) { |
1096 | > | |
1097 | > | // Always do this if we're in different cells or if |
1098 | > | // we're in the same cell and the global index of the |
1099 | > | // j2 cutoff group is less than the j1 cutoff group |
1100 | > | |
1101 | > | if (m2 != m1 || cgColToGlobal[(*j2)] < cgRowToGlobal[(*j1)]) { |
1102 | > | dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; |
1103 | > | snap_->wrapVector(dr); |
1104 | > | cuts = getGroupCutoffs( (*j1), (*j2) ); |
1105 | > | if (dr.lengthSquare() < cuts.third) { |
1106 | > | neighborList.push_back(make_pair((*j1), (*j2))); |
1107 | > | } |
1108 | } | |
1109 | } | |
1110 | } | |
1112 | – | } |
1111 | #else | |
1112 | < | |
1113 | < | for (vector<int>::iterator j1 = cellList_[m1].begin(); |
1114 | < | j1 != cellList_[m1].end(); ++j1) { |
1115 | < | for (vector<int>::iterator j2 = cellList_[m2].begin(); |
1116 | < | j2 != cellList_[m2].end(); ++j2) { |
1117 | < | |
1118 | < | // Always do this if we're in different cells or if |
1119 | < | // we're in the same cell and the global index of the |
1120 | < | // j2 cutoff group is less than the j1 cutoff group |
1121 | < | |
1122 | < | if (m2 != m1 || (*j2) < (*j1)) { |
1123 | < | dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)]; |
1124 | < | snap_->wrapVector(dr); |
1125 | < | cuts = getGroupCutoffs( (*j1), (*j2) ); |
1126 | < | if (dr.lengthSquare() < cuts.third) { |
1127 | < | neighborList.push_back(make_pair((*j1), (*j2))); |
1112 | > | |
1113 | > | for (vector<int>::iterator j1 = cellList_[m1].begin(); |
1114 | > | j1 != cellList_[m1].end(); ++j1) { |
1115 | > | for (vector<int>::iterator j2 = cellList_[m2].begin(); |
1116 | > | j2 != cellList_[m2].end(); ++j2) { |
1117 | > | |
1118 | > | // Always do this if we're in different cells or if |
1119 | > | // we're in the same cell and the global index of the |
1120 | > | // j2 cutoff group is less than the j1 cutoff group |
1121 | > | |
1122 | > | if (m2 != m1 || (*j2) < (*j1)) { |
1123 | > | dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)]; |
1124 | > | snap_->wrapVector(dr); |
1125 | > | cuts = getGroupCutoffs( (*j1), (*j2) ); |
1126 | > | if (dr.lengthSquare() < cuts.third) { |
1127 | > | neighborList.push_back(make_pair((*j1), (*j2))); |
1128 | > | } |
1129 | } | |
1130 | } | |
1131 | } | |
1133 | – | } |
1132 | #endif | |
1133 | + | } |
1134 | } | |
1135 | } | |
1136 | } | |
1137 | + | } else { |
1138 | + | // branch to do all cutoff group pairs |
1139 | + | #ifdef IS_MPI |
1140 | + | for (int j1 = 0; j1 < nGroupsInRow_; j1++) { |
1141 | + | for (int j2 = 0; j2 < nGroupsInCol_; j2++) { |
1142 | + | dr = cgColData.position[j2] - cgRowData.position[j1]; |
1143 | + | snap_->wrapVector(dr); |
1144 | + | cuts = getGroupCutoffs( j1, j2 ); |
1145 | + | if (dr.lengthSquare() < cuts.third) { |
1146 | + | neighborList.push_back(make_pair(j1, j2)); |
1147 | + | } |
1148 | + | } |
1149 | + | } |
1150 | + | #else |
1151 | + | for (int j1 = 0; j1 < nGroups_ - 1; j1++) { |
1152 | + | for (int j2 = j1 + 1; j2 < nGroups_; j2++) { |
1153 | + | dr = snap_->cgData.position[j2] - snap_->cgData.position[j1]; |
1154 | + | snap_->wrapVector(dr); |
1155 | + | cuts = getGroupCutoffs( j1, j2 ); |
1156 | + | if (dr.lengthSquare() < cuts.third) { |
1157 | + | neighborList.push_back(make_pair(j1, j2)); |
1158 | + | } |
1159 | + | } |
1160 | + | } |
1161 | + | #endif |
1162 | } | |
1163 | < | |
1163 | > | |
1164 | // save the local cutoff group positions for the check that is | |
1165 | // done on each loop: | |
1166 | saved_CG_positions_.clear(); | |
1167 | for (int i = 0; i < nGroups_; i++) | |
1168 | saved_CG_positions_.push_back(snap_->cgData.position[i]); | |
1169 | < | |
1169 | > | |
1170 | return neighborList; | |
1171 | } | |
1172 | } //end namespace OpenMD |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |