# | 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(); | |
65 | AtomLocalToGlobal = info_->getGlobalAtomIndices(); | |
# | Line 141 | Line 142 | namespace OpenMD { | |
142 | } | |
143 | } | |
144 | ||
145 | < | skipsForRowAtom.clear(); |
146 | < | skipsForRowAtom.resize(nAtomsInRow_); |
145 | > | skipsForAtom.clear(); |
146 | > | skipsForAtom.resize(nAtomsInRow_); |
147 | > | toposForAtom.clear(); |
148 | > | toposForAtom.resize(nAtomsInRow_); |
149 | > | topoDist.clear(); |
150 | > | topoDist.resize(nAtomsInRow_); |
151 | for (int i = 0; i < nAtomsInRow_; i++) { | |
152 | 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 | – | } |
153 | ||
155 | – | toposForRowAtom.clear(); |
156 | – | toposForRowAtom.resize(nAtomsInRow_); |
157 | – | for (int i = 0; i < nAtomsInRow_; i++) { |
158 | – | int iglob = AtomRowToGlobal[i]; |
159 | – | int nTopos = 0; |
154 | for (int j = 0; j < nAtomsInCol_; j++) { | |
155 | < | int jglob = AtomColToGlobal[j]; |
155 | > | int jglob = AtomColToGlobal[j]; |
156 | > | |
157 | > | if (excludes.hasPair(iglob, jglob)) |
158 | > | skipsForAtom[i].push_back(j); |
159 | > | |
160 | if (oneTwo.hasPair(iglob, jglob)) { | |
161 | < | toposForRowAtom[i].push_back(j); |
162 | < | topoDistRow[i][nTopos] = 1; |
163 | < | nTopos++; |
161 | > | toposForAtom[i].push_back(j); |
162 | > | topoDist[i].push_back(1); |
163 | > | } else { |
164 | > | if (oneThree.hasPair(iglob, jglob)) { |
165 | > | toposForAtom[i].push_back(j); |
166 | > | topoDist[i].push_back(2); |
167 | > | } else { |
168 | > | if (oneFour.hasPair(iglob, jglob)) { |
169 | > | toposForAtom[i].push_back(j); |
170 | > | topoDist[i].push_back(3); |
171 | > | } |
172 | > | } |
173 | } | |
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 | – | } |
174 | } | |
175 | } | |
176 | ||
177 | #endif | |
178 | + | |
179 | groupList_.clear(); | |
180 | groupList_.resize(nGroups_); | |
181 | for (int i = 0; i < nGroups_; i++) { | |
# | Line 186 | Line 184 | namespace OpenMD { | |
184 | int aid = AtomLocalToGlobal[j]; | |
185 | if (globalGroupMembership[aid] == gid) { | |
186 | groupList_[i].push_back(j); | |
189 | – | |
187 | } | |
188 | } | |
189 | } | |
190 | ||
191 | < | skipsForLocalAtom.clear(); |
192 | < | skipsForLocalAtom.resize(nLocal_); |
191 | > | skipsForAtom.clear(); |
192 | > | skipsForAtom.resize(nLocal_); |
193 | > | toposForAtom.clear(); |
194 | > | toposForAtom.resize(nLocal_); |
195 | > | topoDist.clear(); |
196 | > | topoDist.resize(nLocal_); |
197 | ||
198 | for (int i = 0; i < nLocal_; i++) { | |
199 | int iglob = AtomLocalToGlobal[i]; | |
200 | + | |
201 | for (int j = 0; j < nLocal_; j++) { | |
202 | < | int jglob = AtomLocalToGlobal[j]; |
202 | > | int jglob = AtomLocalToGlobal[j]; |
203 | > | |
204 | if (excludes.hasPair(iglob, jglob)) | |
205 | < | skipsForLocalAtom[i].push_back(j); |
206 | < | } |
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]; |
205 | > | skipsForAtom[i].push_back(j); |
206 | > | |
207 | if (oneTwo.hasPair(iglob, jglob)) { | |
208 | < | toposForLocalAtom[i].push_back(j); |
209 | < | topoDistLocal[i][nTopos] = 1; |
210 | < | nTopos++; |
211 | < | } |
212 | < | if (oneThree.hasPair(iglob, jglob)) { |
213 | < | toposForLocalAtom[i].push_back(j); |
214 | < | topoDistLocal[i][nTopos] = 2; |
215 | < | nTopos++; |
216 | < | } |
217 | < | if (oneFour.hasPair(iglob, jglob)) { |
218 | < | toposForLocalAtom[i].push_back(j); |
219 | < | topoDistLocal[i][nTopos] = 3; |
225 | < | nTopos++; |
208 | > | toposForAtom[i].push_back(j); |
209 | > | topoDist[i].push_back(1); |
210 | > | } else { |
211 | > | if (oneThree.hasPair(iglob, jglob)) { |
212 | > | toposForAtom[i].push_back(j); |
213 | > | topoDist[i].push_back(2); |
214 | > | } else { |
215 | > | if (oneFour.hasPair(iglob, jglob)) { |
216 | > | toposForAtom[i].push_back(j); |
217 | > | topoDist[i].push_back(3); |
218 | > | } |
219 | > | } |
220 | } | |
221 | } | |
222 | < | } |
223 | < | |
222 | > | } |
223 | > | |
224 | > | createGtypeCutoffMap(); |
225 | } | |
226 | ||
227 | void ForceMatrixDecomposition::createGtypeCutoffMap() { | |
# | Line 238 | Line 233 | namespace OpenMD { | |
233 | vector<RealType> atypeCutoff; | |
234 | atypeCutoff.resize( atypes.size() ); | |
235 | ||
236 | < | for (set<AtomType*>::iterator at = atypes.begin(); at != atypes.end(); ++at){ |
236 | > | for (set<AtomType*>::iterator at = atypes.begin(); |
237 | > | at != atypes.end(); ++at){ |
238 | rc = interactionMan_->getSuggestedCutoffRadius(*at); | |
239 | atid = (*at)->getIdent(); | |
240 | atypeCutoff[atid] = rc; | |
# | Line 250 | Line 246 | namespace OpenMD { | |
246 | // largest cutoff for any atypes present in this group. | |
247 | #ifdef IS_MPI | |
248 | vector<RealType> groupCutoffRow(nGroupsInRow_, 0.0); | |
249 | + | groupRowToGtype.resize(nGroupsInRow_); |
250 | for (int cg1 = 0; cg1 < nGroupsInRow_; cg1++) { | |
251 | vector<int> atomListRow = getAtomsInGroupRow(cg1); | |
252 | for (vector<int>::iterator ia = atomListRow.begin(); | |
# | Line 275 | Line 272 | namespace OpenMD { | |
272 | ||
273 | } | |
274 | vector<RealType> groupCutoffCol(nGroupsInCol_, 0.0); | |
275 | + | groupColToGtype.resize(nGroupsInCol_); |
276 | for (int cg2 = 0; cg2 < nGroupsInCol_; cg2++) { | |
277 | vector<int> atomListCol = getAtomsInGroupColumn(cg2); | |
278 | for (vector<int>::iterator jb = atomListCol.begin(); | |
# | Line 298 | Line 296 | namespace OpenMD { | |
296 | } | |
297 | } | |
298 | #else | |
299 | + | |
300 | vector<RealType> groupCutoff(nGroups_, 0.0); | |
301 | + | groupToGtype.resize(nGroups_); |
302 | + | |
303 | + | cerr << "nGroups = " << nGroups_ << "\n"; |
304 | for (int cg1 = 0; cg1 < nGroups_; cg1++) { | |
305 | + | |
306 | groupCutoff[cg1] = 0.0; | |
307 | vector<int> atomList = getAtomsInGroupRow(cg1); | |
308 | + | |
309 | for (vector<int>::iterator ia = atomList.begin(); | |
310 | ia != atomList.end(); ++ia) { | |
311 | int atom1 = (*ia); | |
# | Line 325 | Line 329 | namespace OpenMD { | |
329 | } | |
330 | #endif | |
331 | ||
332 | + | cerr << "gTypeCutoffs.size() = " << gTypeCutoffs.size() << "\n"; |
333 | // Now we find the maximum group cutoff value present in the simulation | |
334 | ||
335 | < | vector<RealType>::iterator groupMaxLoc = max_element(gTypeCutoffs.begin(), gTypeCutoffs.end()); |
331 | < | RealType groupMax = *groupMaxLoc; |
335 | > | RealType groupMax = *max_element(gTypeCutoffs.begin(), gTypeCutoffs.end()); |
336 | ||
337 | #ifdef IS_MPI | |
338 | MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE, MPI::MAX); | |
# | Line 337 | Line 341 | namespace OpenMD { | |
341 | RealType tradRcut = groupMax; | |
342 | ||
343 | for (int i = 0; i < gTypeCutoffs.size(); i++) { | |
344 | < | for (int j = 0; j < gTypeCutoffs.size(); j++) { |
341 | < | |
344 | > | for (int j = 0; j < gTypeCutoffs.size(); j++) { |
345 | RealType thisRcut; | |
346 | switch(cutoffPolicy_) { | |
347 | case TRADITIONAL: | |
348 | thisRcut = tradRcut; | |
349 | + | break; |
350 | case MIX: | |
351 | thisRcut = 0.5 * (gTypeCutoffs[i] + gTypeCutoffs[j]); | |
352 | + | break; |
353 | case MAX: | |
354 | thisRcut = max(gTypeCutoffs[i], gTypeCutoffs[j]); | |
355 | + | break; |
356 | default: | |
357 | sprintf(painCave.errMsg, | |
358 | "ForceMatrixDecomposition::createGtypeCutoffMap " | |
359 | "hit an unknown cutoff policy!\n"); | |
360 | painCave.severity = OPENMD_ERROR; | |
361 | painCave.isFatal = 1; | |
362 | < | simError(); |
362 | > | simError(); |
363 | > | break; |
364 | } | |
365 | ||
366 | pair<int,int> key = make_pair(i,j); | |
# | Line 383 | Line 390 | namespace OpenMD { | |
390 | ||
391 | ||
392 | groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) { | |
393 | < | int i, j; |
387 | < | |
393 | > | int i, j; |
394 | #ifdef IS_MPI | |
395 | i = groupRowToGtype[cg1]; | |
396 | j = groupColToGtype[cg2]; | |
397 | #else | |
398 | i = groupToGtype[cg1]; | |
399 | j = groupToGtype[cg2]; | |
400 | < | #endif |
395 | < | |
400 | > | #endif |
401 | return gTypeCutoffMap[make_pair(i,j)]; | |
402 | } | |
403 | ||
404 | + | int ForceMatrixDecomposition::getTopologicalDistance(int atom1, int atom2) { |
405 | + | for (int j = 0; j < toposForAtom[atom1].size(); j++) { |
406 | + | if (toposForAtom[atom1][j] == atom2) |
407 | + | return topoDist[atom1][j]; |
408 | + | } |
409 | + | return 0; |
410 | + | } |
411 | ||
412 | void ForceMatrixDecomposition::zeroWorkArrays() { | |
413 | ||
# | Line 705 | Line 717 | namespace OpenMD { | |
717 | return d; | |
718 | } | |
719 | ||
720 | < | vector<int> ForceMatrixDecomposition::getSkipsForRowAtom(int atom1) { |
721 | < | #ifdef IS_MPI |
710 | < | return skipsForRowAtom[atom1]; |
711 | < | #else |
712 | < | return skipsForLocalAtom[atom1]; |
713 | < | #endif |
720 | > | vector<int> ForceMatrixDecomposition::getSkipsForAtom(int atom1) { |
721 | > | return skipsForAtom[atom1]; |
722 | } | |
723 | ||
724 | /** | |
# | Line 743 | Line 751 | namespace OpenMD { | |
751 | unique_id_2 = atom2; | |
752 | #endif | |
753 | ||
754 | < | #ifdef IS_MPI |
755 | < | for (vector<int>::iterator i = skipsForRowAtom[atom1].begin(); |
748 | < | i != skipsForRowAtom[atom1].end(); ++i) { |
754 | > | for (vector<int>::iterator i = skipsForAtom[atom1].begin(); |
755 | > | i != skipsForAtom[atom1].end(); ++i) { |
756 | if ( (*i) == unique_id_2 ) return true; | |
757 | } | |
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]; |
764 | < | } |
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 | ||
771 | – | // zero is default for unconnected (i.e. normal) pair interactions |
772 | – | return 0; |
773 | – | } |
761 | ||
762 | void ForceMatrixDecomposition::addForceToAtomRow(int atom1, Vector3d fg){ | |
763 | #ifdef IS_MPI | |
# | Line 958 | Line 945 | namespace OpenMD { | |
945 | Vector3d rs, scaled, dr; | |
946 | Vector3i whichCell; | |
947 | int cellIndex; | |
948 | + | int nCtot = nCells_.x() * nCells_.y() * nCells_.z(); |
949 | ||
950 | #ifdef IS_MPI | |
951 | + | cellListRow_.resize(nCtot); |
952 | + | cellListCol_.resize(nCtot); |
953 | + | #else |
954 | + | cellList_.resize(nCtot); |
955 | + | #endif |
956 | + | |
957 | + | #ifdef IS_MPI |
958 | for (int i = 0; i < nGroupsInRow_; i++) { | |
959 | rs = cgRowData.position[i]; | |
960 | // scaled positions relative to the box vectors |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |