# | Line 59 | Line 59 | namespace OpenMD { | |
---|---|---|
59 | nLocal_ = snap_->getNumberOfAtoms(); | |
60 | ||
61 | nGroups_ = info_->getNLocalCutoffGroups(); | |
62 | + | cerr << "in dId, nGroups = " << nGroups_ << "\n"; |
63 | // gather the information for atomtype IDs (atids): | |
64 | < | identsLocal = info_->getIdentArray(); |
64 | > | idents = info_->getIdentArray(); |
65 | AtomLocalToGlobal = info_->getGlobalAtomIndices(); | |
66 | cgLocalToGlobal = info_->getGlobalGroupIndices(); | |
67 | vector<int> globalGroupMembership = info_->getGlobalGroupMembership(); | |
68 | < | vector<RealType> massFactorsLocal = info_->getMassFactors(); |
68 | > | massFactors = info_->getMassFactors(); |
69 | > | |
70 | PairList excludes = info_->getExcludedInteractions(); | |
71 | PairList oneTwo = info_->getOneTwoInteractions(); | |
72 | PairList oneThree = info_->getOneThreeInteractions(); | |
# | Line 107 | 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 | AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal); | |
116 | AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal); | |
# | Line 116 | Line 118 | namespace OpenMD { | |
118 | cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal); | |
119 | cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal); | |
120 | ||
121 | < | AtomCommRealRow->gather(massFactorsLocal, massFactorsRow); |
122 | < | AtomCommRealColumn->gather(massFactorsLocal, massFactorsCol); |
121 | > | AtomCommRealRow->gather(massFactors, massFactorsRow); |
122 | > | AtomCommRealColumn->gather(massFactors, massFactorsCol); |
123 | ||
124 | groupListRow_.clear(); | |
125 | groupListRow_.resize(nGroupsInRow_); | |
# | Line 141 | Line 143 | namespace OpenMD { | |
143 | } | |
144 | } | |
145 | ||
146 | < | skipsForRowAtom.clear(); |
147 | < | skipsForRowAtom.resize(nAtomsInRow_); |
146 | > | skipsForAtom.clear(); |
147 | > | skipsForAtom.resize(nAtomsInRow_); |
148 | > | toposForAtom.clear(); |
149 | > | toposForAtom.resize(nAtomsInRow_); |
150 | > | topoDist.clear(); |
151 | > | topoDist.resize(nAtomsInRow_); |
152 | for (int i = 0; i < nAtomsInRow_; i++) { | |
153 | int iglob = AtomRowToGlobal[i]; | |
148 | – | for (int j = 0; j < nAtomsInCol_; j++) { |
149 | – | int jglob = AtomColToGlobal[j]; |
150 | – | if (excludes.hasPair(iglob, jglob)) |
151 | – | skipsForRowAtom[i].push_back(j); |
152 | – | } |
153 | – | } |
154 | ||
155 | – | toposForRowAtom.clear(); |
156 | – | toposForRowAtom.resize(nAtomsInRow_); |
157 | – | for (int i = 0; i < nAtomsInRow_; i++) { |
158 | – | int iglob = AtomRowToGlobal[i]; |
159 | – | int nTopos = 0; |
155 | for (int j = 0; j < nAtomsInCol_; j++) { | |
156 | < | int jglob = AtomColToGlobal[j]; |
156 | > | int jglob = AtomColToGlobal[j]; |
157 | > | |
158 | > | if (excludes.hasPair(iglob, jglob)) |
159 | > | skipsForAtom[i].push_back(j); |
160 | > | |
161 | if (oneTwo.hasPair(iglob, jglob)) { | |
162 | < | toposForRowAtom[i].push_back(j); |
163 | < | topoDistRow[i][nTopos] = 1; |
164 | < | nTopos++; |
162 | > | toposForAtom[i].push_back(j); |
163 | > | topoDist[i].push_back(1); |
164 | > | } else { |
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)) { |
170 | > | toposForAtom[i].push_back(j); |
171 | > | topoDist[i].push_back(3); |
172 | > | } |
173 | > | } |
174 | } | |
167 | – | if (oneThree.hasPair(iglob, jglob)) { |
168 | – | toposForRowAtom[i].push_back(j); |
169 | – | topoDistRow[i][nTopos] = 2; |
170 | – | nTopos++; |
171 | – | } |
172 | – | if (oneFour.hasPair(iglob, jglob)) { |
173 | – | toposForRowAtom[i].push_back(j); |
174 | – | topoDistRow[i][nTopos] = 3; |
175 | – | nTopos++; |
176 | – | } |
175 | } | |
176 | } | |
177 | ||
178 | #endif | |
179 | + | |
180 | groupList_.clear(); | |
181 | groupList_.resize(nGroups_); | |
182 | for (int i = 0; i < nGroups_; i++) { | |
# | Line 186 | Line 185 | namespace OpenMD { | |
185 | int aid = AtomLocalToGlobal[j]; | |
186 | if (globalGroupMembership[aid] == gid) { | |
187 | groupList_[i].push_back(j); | |
189 | – | |
188 | } | |
189 | } | |
190 | } | |
191 | ||
192 | < | skipsForLocalAtom.clear(); |
193 | < | skipsForLocalAtom.resize(nLocal_); |
192 | > | skipsForAtom.clear(); |
193 | > | skipsForAtom.resize(nLocal_); |
194 | > | toposForAtom.clear(); |
195 | > | toposForAtom.resize(nLocal_); |
196 | > | topoDist.clear(); |
197 | > | topoDist.resize(nLocal_); |
198 | ||
199 | for (int i = 0; i < nLocal_; i++) { | |
200 | int iglob = AtomLocalToGlobal[i]; | |
201 | + | |
202 | for (int j = 0; j < nLocal_; j++) { | |
203 | < | int jglob = AtomLocalToGlobal[j]; |
203 | > | int jglob = AtomLocalToGlobal[j]; |
204 | > | |
205 | if (excludes.hasPair(iglob, jglob)) | |
206 | < | skipsForLocalAtom[i].push_back(j); |
207 | < | } |
204 | < | } |
205 | < | toposForLocalAtom.clear(); |
206 | < | toposForLocalAtom.resize(nLocal_); |
207 | < | for (int i = 0; i < nLocal_; i++) { |
208 | < | int iglob = AtomLocalToGlobal[i]; |
209 | < | int nTopos = 0; |
210 | < | for (int j = 0; j < nLocal_; j++) { |
211 | < | int jglob = AtomLocalToGlobal[j]; |
206 | > | skipsForAtom[i].push_back(j); |
207 | > | |
208 | if (oneTwo.hasPair(iglob, jglob)) { | |
209 | < | toposForLocalAtom[i].push_back(j); |
210 | < | topoDistLocal[i][nTopos] = 1; |
211 | < | nTopos++; |
209 | > | toposForAtom[i].push_back(j); |
210 | > | topoDist[i].push_back(1); |
211 | > | } else { |
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)) { |
217 | > | toposForAtom[i].push_back(j); |
218 | > | topoDist[i].push_back(3); |
219 | > | } |
220 | > | } |
221 | } | |
217 | – | if (oneThree.hasPair(iglob, jglob)) { |
218 | – | toposForLocalAtom[i].push_back(j); |
219 | – | topoDistLocal[i][nTopos] = 2; |
220 | – | nTopos++; |
221 | – | } |
222 | – | if (oneFour.hasPair(iglob, jglob)) { |
223 | – | toposForLocalAtom[i].push_back(j); |
224 | – | topoDistLocal[i][nTopos] = 3; |
225 | – | nTopos++; |
226 | – | } |
222 | } | |
223 | < | } |
224 | < | |
223 | > | } |
224 | > | |
225 | > | createGtypeCutoffMap(); |
226 | } | |
227 | ||
228 | void ForceMatrixDecomposition::createGtypeCutoffMap() { | |
# | Line 237 | Line 233 | namespace OpenMD { | |
233 | set<AtomType*> atypes = info_->getSimulatedAtomTypes(); | |
234 | vector<RealType> atypeCutoff; | |
235 | atypeCutoff.resize( atypes.size() ); | |
236 | < | |
237 | < | for (set<AtomType*>::iterator at = atypes.begin(); at != atypes.end(); ++at){ |
238 | < | rc = interactionMan_->getSuggestedCutoffRadius(*at); |
236 | > | |
237 | > | for (set<AtomType*>::iterator at = atypes.begin(); |
238 | > | at != atypes.end(); ++at){ |
239 | atid = (*at)->getIdent(); | |
240 | < | atypeCutoff[atid] = rc; |
240 | > | |
241 | > | if (userChoseCutoff_) |
242 | > | atypeCutoff[atid] = userCutoff_; |
243 | > | else |
244 | > | atypeCutoff[atid] = interactionMan_->getSuggestedCutoffRadius(*at); |
245 | } | |
246 | ||
247 | vector<RealType> gTypeCutoffs; | |
# | Line 250 | Line 250 | namespace OpenMD { | |
250 | // largest cutoff for any atypes present in this group. | |
251 | #ifdef IS_MPI | |
252 | vector<RealType> groupCutoffRow(nGroupsInRow_, 0.0); | |
253 | + | groupRowToGtype.resize(nGroupsInRow_); |
254 | for (int cg1 = 0; cg1 < nGroupsInRow_; cg1++) { | |
255 | vector<int> atomListRow = getAtomsInGroupRow(cg1); | |
256 | for (vector<int>::iterator ia = atomListRow.begin(); | |
# | Line 275 | Line 276 | namespace OpenMD { | |
276 | ||
277 | } | |
278 | vector<RealType> groupCutoffCol(nGroupsInCol_, 0.0); | |
279 | + | groupColToGtype.resize(nGroupsInCol_); |
280 | for (int cg2 = 0; cg2 < nGroupsInCol_; cg2++) { | |
281 | vector<int> atomListCol = getAtomsInGroupColumn(cg2); | |
282 | for (vector<int>::iterator jb = atomListCol.begin(); | |
# | Line 298 | Line 300 | namespace OpenMD { | |
300 | } | |
301 | } | |
302 | #else | |
303 | + | |
304 | vector<RealType> groupCutoff(nGroups_, 0.0); | |
305 | + | groupToGtype.resize(nGroups_); |
306 | + | |
307 | + | cerr << "nGroups = " << nGroups_ << "\n"; |
308 | for (int cg1 = 0; cg1 < nGroups_; cg1++) { | |
309 | + | |
310 | groupCutoff[cg1] = 0.0; | |
311 | vector<int> atomList = getAtomsInGroupRow(cg1); | |
312 | + | |
313 | for (vector<int>::iterator ia = atomList.begin(); | |
314 | ia != atomList.end(); ++ia) { | |
315 | int atom1 = (*ia); | |
316 | < | atid = identsLocal[atom1]; |
316 | > | atid = idents[atom1]; |
317 | if (atypeCutoff[atid] > groupCutoff[cg1]) { | |
318 | groupCutoff[cg1] = atypeCutoff[atid]; | |
319 | } | |
# | Line 325 | Line 333 | namespace OpenMD { | |
333 | } | |
334 | #endif | |
335 | ||
336 | + | cerr << "gTypeCutoffs.size() = " << gTypeCutoffs.size() << "\n"; |
337 | // Now we find the maximum group cutoff value present in the simulation | |
338 | ||
339 | < | vector<RealType>::iterator groupMaxLoc = max_element(gTypeCutoffs.begin(), gTypeCutoffs.end()); |
331 | < | RealType groupMax = *groupMaxLoc; |
339 | > | RealType groupMax = *max_element(gTypeCutoffs.begin(), gTypeCutoffs.end()); |
340 | ||
341 | #ifdef IS_MPI | |
342 | MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE, MPI::MAX); | |
# | Line 337 | Line 345 | namespace OpenMD { | |
345 | RealType tradRcut = groupMax; | |
346 | ||
347 | for (int i = 0; i < gTypeCutoffs.size(); i++) { | |
348 | < | for (int j = 0; j < gTypeCutoffs.size(); j++) { |
341 | < | |
348 | > | for (int j = 0; j < gTypeCutoffs.size(); j++) { |
349 | RealType thisRcut; | |
350 | switch(cutoffPolicy_) { | |
351 | case TRADITIONAL: | |
352 | thisRcut = tradRcut; | |
353 | + | break; |
354 | case MIX: | |
355 | thisRcut = 0.5 * (gTypeCutoffs[i] + gTypeCutoffs[j]); | |
356 | + | break; |
357 | case MAX: | |
358 | thisRcut = max(gTypeCutoffs[i], gTypeCutoffs[j]); | |
359 | + | break; |
360 | default: | |
361 | sprintf(painCave.errMsg, | |
362 | "ForceMatrixDecomposition::createGtypeCutoffMap " | |
363 | "hit an unknown cutoff policy!\n"); | |
364 | painCave.severity = OPENMD_ERROR; | |
365 | painCave.isFatal = 1; | |
366 | < | simError(); |
366 | > | simError(); |
367 | > | break; |
368 | } | |
369 | ||
370 | pair<int,int> key = make_pair(i,j); | |
# | Line 371 | Line 382 | namespace OpenMD { | |
382 | if (abs(gTypeCutoffMap[key].first - userCutoff_) > 0.0001) { | |
383 | sprintf(painCave.errMsg, | |
384 | "ForceMatrixDecomposition::createGtypeCutoffMap " | |
385 | < | "user-specified rCut does not match computed group Cutoff\n"); |
385 | > | "user-specified rCut (%lf) does not match computed group Cutoff\n", userCutoff_); |
386 | painCave.severity = OPENMD_ERROR; | |
387 | painCave.isFatal = 1; | |
388 | simError(); | |
# | Line 383 | Line 394 | namespace OpenMD { | |
394 | ||
395 | ||
396 | groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) { | |
397 | < | int i, j; |
387 | < | |
397 | > | int i, j; |
398 | #ifdef IS_MPI | |
399 | i = groupRowToGtype[cg1]; | |
400 | j = groupColToGtype[cg2]; | |
401 | #else | |
402 | i = groupToGtype[cg1]; | |
403 | j = groupToGtype[cg2]; | |
404 | < | #endif |
395 | < | |
404 | > | #endif |
405 | return gTypeCutoffMap[make_pair(i,j)]; | |
406 | } | |
407 | ||
408 | + | int ForceMatrixDecomposition::getTopologicalDistance(int atom1, int atom2) { |
409 | + | for (int j = 0; j < toposForAtom[atom1].size(); j++) { |
410 | + | if (toposForAtom[atom1][j] == atom2) |
411 | + | return topoDist[atom1][j]; |
412 | + | } |
413 | + | return 0; |
414 | + | } |
415 | ||
416 | void ForceMatrixDecomposition::zeroWorkArrays() { | |
417 | + | pairwisePot = 0.0; |
418 | + | embeddingPot = 0.0; |
419 | ||
402 | – | for (int j = 0; j < N_INTERACTION_FAMILIES; j++) { |
403 | – | longRangePot_[j] = 0.0; |
404 | – | } |
405 | – | |
420 | #ifdef IS_MPI | |
421 | if (storageLayout_ & DataStorage::dslForce) { | |
422 | fill(atomRowData.force.begin(), atomRowData.force.end(), V3Zero); | |
# | Line 418 | Line 432 | namespace OpenMD { | |
432 | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | |
433 | ||
434 | fill(pot_col.begin(), pot_col.end(), | |
435 | < | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); |
422 | < | |
423 | < | pot_local = Vector<RealType, N_INTERACTION_FAMILIES>(0.0); |
435 | > | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); |
436 | ||
437 | if (storageLayout_ & DataStorage::dslParticlePot) { | |
438 | fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(), 0.0); | |
# | Line 594 | Line 606 | namespace OpenMD { | |
606 | AtomCommPotRow->scatter(pot_row, pot_temp); | |
607 | ||
608 | for (int ii = 0; ii < pot_temp.size(); ii++ ) | |
609 | < | pot_local += pot_temp[ii]; |
609 | > | pairwisePot += pot_temp[ii]; |
610 | ||
611 | fill(pot_temp.begin(), pot_temp.end(), | |
612 | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | |
# | Line 602 | Line 614 | namespace OpenMD { | |
614 | AtomCommPotColumn->scatter(pot_col, pot_temp); | |
615 | ||
616 | for (int ii = 0; ii < pot_temp.size(); ii++ ) | |
617 | < | pot_local += pot_temp[ii]; |
606 | < | |
617 | > | pairwisePot += pot_temp[ii]; |
618 | #endif | |
619 | + | |
620 | } | |
621 | ||
622 | int ForceMatrixDecomposition::getNAtomsInRow() { | |
# | Line 679 | Line 691 | namespace OpenMD { | |
691 | #ifdef IS_MPI | |
692 | return massFactorsRow[atom1]; | |
693 | #else | |
694 | < | return massFactorsLocal[atom1]; |
694 | > | cerr << "mfs = " << massFactors.size() << " atom1 = " << atom1 << "\n"; |
695 | > | return massFactors[atom1]; |
696 | #endif | |
697 | } | |
698 | ||
# | Line 687 | Line 700 | namespace OpenMD { | |
700 | #ifdef IS_MPI | |
701 | return massFactorsCol[atom2]; | |
702 | #else | |
703 | < | return massFactorsLocal[atom2]; |
703 | > | return massFactors[atom2]; |
704 | #endif | |
705 | ||
706 | } | |
# | Line 705 | Line 718 | namespace OpenMD { | |
718 | return d; | |
719 | } | |
720 | ||
721 | < | vector<int> ForceMatrixDecomposition::getSkipsForRowAtom(int atom1) { |
722 | < | #ifdef IS_MPI |
710 | < | return skipsForRowAtom[atom1]; |
711 | < | #else |
712 | < | return skipsForLocalAtom[atom1]; |
713 | < | #endif |
721 | > | vector<int> ForceMatrixDecomposition::getSkipsForAtom(int atom1) { |
722 | > | return skipsForAtom[atom1]; |
723 | } | |
724 | ||
725 | /** | |
# | Line 743 | Line 752 | namespace OpenMD { | |
752 | unique_id_2 = atom2; | |
753 | #endif | |
754 | ||
755 | < | #ifdef IS_MPI |
756 | < | for (vector<int>::iterator i = skipsForRowAtom[atom1].begin(); |
748 | < | i != skipsForRowAtom[atom1].end(); ++i) { |
755 | > | for (vector<int>::iterator i = skipsForAtom[atom1].begin(); |
756 | > | i != skipsForAtom[atom1].end(); ++i) { |
757 | if ( (*i) == unique_id_2 ) return true; | |
750 | – | } |
751 | – | #else |
752 | – | for (vector<int>::iterator i = skipsForLocalAtom[atom1].begin(); |
753 | – | i != skipsForLocalAtom[atom1].end(); ++i) { |
754 | – | if ( (*i) == unique_id_2 ) return true; |
755 | – | } |
756 | – | #endif |
757 | – | } |
758 | – | |
759 | – | int ForceMatrixDecomposition::getTopoDistance(int atom1, int atom2) { |
760 | – | |
761 | – | #ifdef IS_MPI |
762 | – | for (int i = 0; i < toposForRowAtom[atom1].size(); i++) { |
763 | – | if ( toposForRowAtom[atom1][i] == atom2 ) return topoDistRow[atom1][i]; |
758 | } | |
765 | – | #else |
766 | – | for (int i = 0; i < toposForLocalAtom[atom1].size(); i++) { |
767 | – | if ( toposForLocalAtom[atom1][i] == atom2 ) return topoDistLocal[atom1][i]; |
768 | – | } |
769 | – | #endif |
759 | ||
760 | < | // zero is default for unconnected (i.e. normal) pair interactions |
772 | < | return 0; |
760 | > | return false; |
761 | } | |
762 | ||
763 | + | |
764 | void ForceMatrixDecomposition::addForceToAtomRow(int atom1, Vector3d fg){ | |
765 | #ifdef IS_MPI | |
766 | atomRowData.force[atom1] += fg; | |
# | Line 789 | Line 778 | namespace OpenMD { | |
778 | } | |
779 | ||
780 | // filling interaction blocks with pointers | |
781 | < | InteractionData ForceMatrixDecomposition::fillInteractionData(int atom1, int atom2) { |
782 | < | InteractionData idat; |
794 | < | |
781 | > | void ForceMatrixDecomposition::fillInteractionData(InteractionData &idat, |
782 | > | int atom1, int atom2) { |
783 | #ifdef IS_MPI | |
784 | ||
785 | idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), | |
786 | ff_->getAtomType(identsCol[atom2]) ); | |
799 | – | |
787 | ||
788 | if (storageLayout_ & DataStorage::dslAmat) { | |
789 | idat.A1 = &(atomRowData.aMat[atom1]); | |
# | Line 835 | Line 822 | namespace OpenMD { | |
822 | ||
823 | #else | |
824 | ||
825 | < | idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]), |
826 | < | ff_->getAtomType(identsLocal[atom2]) ); |
825 | > | idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), |
826 | > | ff_->getAtomType(idents[atom2]) ); |
827 | ||
828 | if (storageLayout_ & DataStorage::dslAmat) { | |
829 | idat.A1 = &(snap_->atomData.aMat[atom1]); | |
# | Line 853 | Line 840 | namespace OpenMD { | |
840 | idat.t2 = &(snap_->atomData.torque[atom2]); | |
841 | } | |
842 | ||
843 | < | if (storageLayout_ & DataStorage::dslDensity) { |
843 | > | if (storageLayout_ & DataStorage::dslDensity) { |
844 | idat.rho1 = &(snap_->atomData.density[atom1]); | |
845 | idat.rho2 = &(snap_->atomData.density[atom2]); | |
846 | } | |
# | Line 874 | Line 861 | namespace OpenMD { | |
861 | } | |
862 | ||
863 | #endif | |
877 | – | return idat; |
864 | } | |
865 | ||
866 | ||
867 | < | void ForceMatrixDecomposition::unpackInteractionData(InteractionData idat, int atom1, int atom2) { |
867 | > | void ForceMatrixDecomposition::unpackInteractionData(InteractionData &idat, int atom1, int atom2) { |
868 | #ifdef IS_MPI | |
869 | pot_row[atom1] += 0.5 * *(idat.pot); | |
870 | pot_col[atom2] += 0.5 * *(idat.pot); | |
# | Line 886 | Line 872 | namespace OpenMD { | |
872 | atomRowData.force[atom1] += *(idat.f1); | |
873 | atomColData.force[atom2] -= *(idat.f1); | |
874 | #else | |
875 | < | longRangePot_ += *(idat.pot); |
876 | < | |
875 | > | pairwisePot += *(idat.pot); |
876 | > | |
877 | snap_->atomData.force[atom1] += *(idat.f1); | |
878 | snap_->atomData.force[atom2] -= *(idat.f1); | |
879 | #endif | |
# | Line 895 | Line 881 | namespace OpenMD { | |
881 | } | |
882 | ||
883 | ||
884 | < | InteractionData ForceMatrixDecomposition::fillSkipData(int atom1, int atom2){ |
885 | < | |
900 | < | InteractionData idat; |
884 | > | void ForceMatrixDecomposition::fillSkipData(InteractionData &idat, |
885 | > | int atom1, int atom2) { |
886 | #ifdef IS_MPI | |
887 | idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), | |
888 | ff_->getAtomType(identsCol[atom2]) ); | |
# | Line 906 | Line 891 | namespace OpenMD { | |
891 | idat.eFrame1 = &(atomRowData.electroFrame[atom1]); | |
892 | idat.eFrame2 = &(atomColData.electroFrame[atom2]); | |
893 | } | |
894 | + | |
895 | if (storageLayout_ & DataStorage::dslTorque) { | |
896 | idat.t1 = &(atomRowData.torque[atom1]); | |
897 | idat.t2 = &(atomColData.torque[atom2]); | |
898 | } | |
899 | + | |
900 | + | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
901 | + | idat.skippedCharge1 = &(atomRowData.skippedCharge[atom1]); |
902 | + | idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]); |
903 | + | } |
904 | #else | |
905 | < | idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]), |
906 | < | ff_->getAtomType(identsLocal[atom2]) ); |
905 | > | idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), |
906 | > | ff_->getAtomType(idents[atom2]) ); |
907 | ||
908 | if (storageLayout_ & DataStorage::dslElectroFrame) { | |
909 | idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]); | |
910 | idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]); | |
911 | } | |
912 | + | |
913 | if (storageLayout_ & DataStorage::dslTorque) { | |
914 | idat.t1 = &(snap_->atomData.torque[atom1]); | |
915 | idat.t2 = &(snap_->atomData.torque[atom2]); | |
916 | } | |
917 | + | |
918 | + | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
919 | + | idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]); |
920 | + | idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]); |
921 | + | } |
922 | #endif | |
923 | } | |
924 | ||
925 | + | |
926 | + | void ForceMatrixDecomposition::unpackSkipData(InteractionData &idat, int atom1, int atom2) { |
927 | + | #ifdef IS_MPI |
928 | + | pot_row[atom1] += 0.5 * *(idat.pot); |
929 | + | pot_col[atom2] += 0.5 * *(idat.pot); |
930 | + | #else |
931 | + | pairwisePot += *(idat.pot); |
932 | + | #endif |
933 | + | |
934 | + | } |
935 | + | |
936 | + | |
937 | /* | |
938 | * buildNeighborList | |
939 | * | |
# | Line 958 | Line 967 | namespace OpenMD { | |
967 | Vector3d rs, scaled, dr; | |
968 | Vector3i whichCell; | |
969 | int cellIndex; | |
970 | + | int nCtot = nCells_.x() * nCells_.y() * nCells_.z(); |
971 | ||
972 | #ifdef IS_MPI | |
973 | + | cellListRow_.resize(nCtot); |
974 | + | cellListCol_.resize(nCtot); |
975 | + | #else |
976 | + | cellList_.resize(nCtot); |
977 | + | #endif |
978 | + | |
979 | + | #ifdef IS_MPI |
980 | for (int i = 0; i < nGroupsInRow_; i++) { | |
981 | rs = cgRowData.position[i]; | |
982 | + | |
983 | // scaled positions relative to the box vectors | |
984 | scaled = invHmat * rs; | |
985 | + | |
986 | // wrap the vector back into the unit box by subtracting integer box | |
987 | // numbers | |
988 | < | for (int j = 0; j < 3; j++) |
988 | > | for (int j = 0; j < 3; j++) { |
989 | scaled[j] -= roundMe(scaled[j]); | |
990 | + | scaled[j] += 0.5; |
991 | + | } |
992 | ||
993 | // find xyz-indices of cell that cutoffGroup is in. | |
994 | whichCell.x() = nCells_.x() * scaled.x(); | |
# | Line 976 | Line 997 | namespace OpenMD { | |
997 | ||
998 | // find single index of this cell: | |
999 | cellIndex = Vlinear(whichCell, nCells_); | |
1000 | + | |
1001 | // add this cutoff group to the list of groups in this cell; | |
1002 | cellListRow_[cellIndex].push_back(i); | |
1003 | } | |
1004 | ||
1005 | for (int i = 0; i < nGroupsInCol_; i++) { | |
1006 | rs = cgColData.position[i]; | |
1007 | + | |
1008 | // scaled positions relative to the box vectors | |
1009 | scaled = invHmat * rs; | |
1010 | + | |
1011 | // wrap the vector back into the unit box by subtracting integer box | |
1012 | // numbers | |
1013 | < | for (int j = 0; j < 3; j++) |
1013 | > | for (int j = 0; j < 3; j++) { |
1014 | scaled[j] -= roundMe(scaled[j]); | |
1015 | + | scaled[j] += 0.5; |
1016 | + | } |
1017 | ||
1018 | // find xyz-indices of cell that cutoffGroup is in. | |
1019 | whichCell.x() = nCells_.x() * scaled.x(); | |
# | Line 996 | Line 1022 | namespace OpenMD { | |
1022 | ||
1023 | // find single index of this cell: | |
1024 | cellIndex = Vlinear(whichCell, nCells_); | |
1025 | + | |
1026 | // add this cutoff group to the list of groups in this cell; | |
1027 | cellListCol_[cellIndex].push_back(i); | |
1028 | } | |
1029 | #else | |
1030 | for (int i = 0; i < nGroups_; i++) { | |
1031 | rs = snap_->cgData.position[i]; | |
1032 | + | |
1033 | // scaled positions relative to the box vectors | |
1034 | scaled = invHmat * rs; | |
1035 | + | |
1036 | // wrap the vector back into the unit box by subtracting integer box | |
1037 | // numbers | |
1038 | < | for (int j = 0; j < 3; j++) |
1038 | > | for (int j = 0; j < 3; j++) { |
1039 | scaled[j] -= roundMe(scaled[j]); | |
1040 | + | scaled[j] += 0.5; |
1041 | + | } |
1042 | ||
1043 | // find xyz-indices of cell that cutoffGroup is in. | |
1044 | whichCell.x() = nCells_.x() * scaled.x(); | |
# | Line 1015 | Line 1046 | namespace OpenMD { | |
1046 | whichCell.z() = nCells_.z() * scaled.z(); | |
1047 | ||
1048 | // find single index of this cell: | |
1049 | < | cellIndex = Vlinear(whichCell, nCells_); |
1049 | > | cellIndex = Vlinear(whichCell, nCells_); |
1050 | > | |
1051 | // add this cutoff group to the list of groups in this cell; | |
1052 | cellList_[cellIndex].push_back(i); | |
1053 | } | |
# | Line 1073 | Line 1105 | namespace OpenMD { | |
1105 | } | |
1106 | } | |
1107 | #else | |
1108 | + | |
1109 | for (vector<int>::iterator j1 = cellList_[m1].begin(); | |
1110 | j1 != cellList_[m1].end(); ++j1) { | |
1111 | for (vector<int>::iterator j2 = cellList_[m2].begin(); | |
1112 | j2 != cellList_[m2].end(); ++j2) { | |
1113 | < | |
1113 | > | |
1114 | // Always do this if we're in different cells or if | |
1115 | // we're in the same cell and the global index of the | |
1116 | // j2 cutoff group is less than the j1 cutoff group | |
# | Line 1097 | Line 1130 | namespace OpenMD { | |
1130 | } | |
1131 | } | |
1132 | } | |
1133 | < | |
1133 | > | |
1134 | // save the local cutoff group positions for the check that is | |
1135 | // done on each loop: | |
1136 | saved_CG_positions_.clear(); | |
1137 | for (int i = 0; i < nGroups_; i++) | |
1138 | saved_CG_positions_.push_back(snap_->cgData.position[i]); | |
1139 | < | |
1139 | > | |
1140 | return neighborList; | |
1141 | } | |
1142 | } //end namespace OpenMD |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |