# | Line 225 | Line 225 | namespace OpenMD { | |
---|---|---|
225 | nTopos++; | |
226 | } | |
227 | } | |
228 | < | } |
228 | > | } |
229 | > | |
230 | } | |
231 | ||
232 | + | void ForceMatrixDecomposition::createGtypeCutoffMap() { |
233 | + | |
234 | + | RealType tol = 1e-6; |
235 | + | RealType rc; |
236 | + | int atid; |
237 | + | set<AtomType*> atypes = info_->getSimulatedAtomTypes(); |
238 | + | vector<RealType> atypeCutoff; |
239 | + | atypeCutoff.reserve( atypes.size() ); |
240 | + | |
241 | + | for (set<AtomType*>::iterator at = atypes.begin(); at != atypes.end(); ++at){ |
242 | + | rc = interactionMan_->getSuggestedCutoffRadius(*at); |
243 | + | atid = (*at)->getIdent(); |
244 | + | atypeCutoff[atid] = rc; |
245 | + | } |
246 | + | |
247 | + | vector<RealType> gTypeCutoffs; |
248 | + | |
249 | + | // first we do a single loop over the cutoff groups to find the |
250 | + | // largest cutoff for any atypes present in this group. |
251 | + | #ifdef IS_MPI |
252 | + | vector<RealType> groupCutoffRow(nGroupsInRow_, 0.0); |
253 | + | for (int cg1 = 0; cg1 < nGroupsInRow_; cg1++) { |
254 | + | vector<int> atomListRow = getAtomsInGroupRow(cg1); |
255 | + | for (vector<int>::iterator ia = atomListRow.begin(); |
256 | + | ia != atomListRow.end(); ++ia) { |
257 | + | int atom1 = (*ia); |
258 | + | atid = identsRow[atom1]; |
259 | + | if (atypeCutoff[atid] > groupCutoffRow[cg1]) { |
260 | + | groupCutoffRow[cg1] = atypeCutoff[atid]; |
261 | + | } |
262 | + | } |
263 | + | |
264 | + | bool gTypeFound = false; |
265 | + | for (int gt = 0; gt < gTypeCutoffs.size(); gt++) { |
266 | + | if (abs(groupCutoffRow[cg1] - gTypeCutoffs[gt]) < tol) { |
267 | + | groupRowToGtype[cg1] = gt; |
268 | + | gTypeFound = true; |
269 | + | } |
270 | + | } |
271 | + | if (!gTypeFound) { |
272 | + | gTypeCutoffs.push_back( groupCutoffRow[cg1] ); |
273 | + | groupRowToGtype[cg1] = gTypeCutoffs.size() - 1; |
274 | + | } |
275 | + | |
276 | + | } |
277 | + | vector<RealType> groupCutoffCol(nGroupsInCol_, 0.0); |
278 | + | for (int cg2 = 0; cg2 < nGroupsInCol_; cg2++) { |
279 | + | vector<int> atomListCol = getAtomsInGroupColumn(cg2); |
280 | + | for (vector<int>::iterator jb = atomListCol.begin(); |
281 | + | jb != atomListCol.end(); ++jb) { |
282 | + | int atom2 = (*jb); |
283 | + | atid = identsCol[atom2]; |
284 | + | if (atypeCutoff[atid] > groupCutoffCol[cg2]) { |
285 | + | groupCutoffCol[cg2] = atypeCutoff[atid]; |
286 | + | } |
287 | + | } |
288 | + | bool gTypeFound = false; |
289 | + | for (int gt = 0; gt < gTypeCutoffs.size(); gt++) { |
290 | + | if (abs(groupCutoffCol[cg2] - gTypeCutoffs[gt]) < tol) { |
291 | + | groupColToGtype[cg2] = gt; |
292 | + | gTypeFound = true; |
293 | + | } |
294 | + | } |
295 | + | if (!gTypeFound) { |
296 | + | gTypeCutoffs.push_back( groupCutoffCol[cg2] ); |
297 | + | groupColToGtype[cg2] = gTypeCutoffs.size() - 1; |
298 | + | } |
299 | + | } |
300 | + | #else |
301 | + | vector<RealType> groupCutoff(nGroups_, 0.0); |
302 | + | for (int cg1 = 0; cg1 < nGroups_; cg1++) { |
303 | + | groupCutoff[cg1] = 0.0; |
304 | + | vector<int> atomList = getAtomsInGroupRow(cg1); |
305 | + | for (vector<int>::iterator ia = atomList.begin(); |
306 | + | ia != atomList.end(); ++ia) { |
307 | + | int atom1 = (*ia); |
308 | + | atid = identsLocal[atom1]; |
309 | + | if (atypeCutoff[atid] > groupCutoff[cg1]) { |
310 | + | groupCutoff[cg1] = atypeCutoff[atid]; |
311 | + | } |
312 | + | } |
313 | + | |
314 | + | bool gTypeFound = false; |
315 | + | for (int gt = 0; gt < gTypeCutoffs.size(); gt++) { |
316 | + | if (abs(groupCutoff[cg1] - gTypeCutoffs[gt]) < tol) { |
317 | + | groupToGtype[cg1] = gt; |
318 | + | gTypeFound = true; |
319 | + | } |
320 | + | } |
321 | + | if (!gTypeFound) { |
322 | + | gTypeCutoffs.push_back( groupCutoff[cg1] ); |
323 | + | groupToGtype[cg1] = gTypeCutoffs.size() - 1; |
324 | + | } |
325 | + | } |
326 | + | #endif |
327 | + | |
328 | + | // Now we find the maximum group cutoff value present in the simulation |
329 | + | |
330 | + | vector<RealType>::iterator groupMaxLoc = max_element(gTypeCutoffs.begin(), gTypeCutoffs.end()); |
331 | + | RealType groupMax = *groupMaxLoc; |
332 | + | |
333 | + | #ifdef IS_MPI |
334 | + | MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE, MPI::MAX); |
335 | + | #endif |
336 | + | |
337 | + | RealType tradRcut = groupMax; |
338 | + | |
339 | + | for (int i = 0; i < gTypeCutoffs.size(); i++) { |
340 | + | for (int j = 0; j < gTypeCutoffs.size(); j++) { |
341 | + | |
342 | + | RealType thisRcut; |
343 | + | switch(cutoffPolicy_) { |
344 | + | case TRADITIONAL: |
345 | + | thisRcut = tradRcut; |
346 | + | case MIX: |
347 | + | thisRcut = 0.5 * (gTypeCutoffs[i] + gTypeCutoffs[j]); |
348 | + | case MAX: |
349 | + | thisRcut = max(gTypeCutoffs[i], gTypeCutoffs[j]); |
350 | + | default: |
351 | + | sprintf(painCave.errMsg, |
352 | + | "ForceMatrixDecomposition::createGtypeCutoffMap " |
353 | + | "hit an unknown cutoff policy!\n"); |
354 | + | painCave.severity = OPENMD_ERROR; |
355 | + | painCave.isFatal = 1; |
356 | + | simError(); |
357 | + | } |
358 | + | |
359 | + | pair<int,int> key = make_pair(i,j); |
360 | + | gTypeCutoffMap[key].first = thisRcut; |
361 | + | |
362 | + | if (thisRcut > largestRcut_) largestRcut_ = thisRcut; |
363 | + | |
364 | + | gTypeCutoffMap[key].second = thisRcut*thisRcut; |
365 | + | |
366 | + | gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2); |
367 | + | |
368 | + | // sanity check |
369 | + | |
370 | + | if (userChoseCutoff_) { |
371 | + | if (abs(gTypeCutoffMap[key].first - userCutoff_) > 0.0001) { |
372 | + | sprintf(painCave.errMsg, |
373 | + | "ForceMatrixDecomposition::createGtypeCutoffMap " |
374 | + | "user-specified rCut does not match computed group Cutoff\n"); |
375 | + | painCave.severity = OPENMD_ERROR; |
376 | + | painCave.isFatal = 1; |
377 | + | simError(); |
378 | + | } |
379 | + | } |
380 | + | } |
381 | + | } |
382 | + | } |
383 | + | |
384 | + | |
385 | + | groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) { |
386 | + | int i, j; |
387 | + | |
388 | + | #ifdef IS_MPI |
389 | + | i = groupRowToGtype[cg1]; |
390 | + | j = groupColToGtype[cg2]; |
391 | + | #else |
392 | + | i = groupToGtype[cg1]; |
393 | + | j = groupToGtype[cg2]; |
394 | + | #endif |
395 | + | |
396 | + | return gTypeCutoffMap[make_pair(i,j)]; |
397 | + | } |
398 | + | |
399 | + | |
400 | void ForceMatrixDecomposition::zeroWorkArrays() { | |
401 | ||
402 | for (int j = 0; j < N_INTERACTION_FAMILIES; j++) { | |
# | Line 765 | Line 934 | namespace OpenMD { | |
934 | vector<pair<int, int> > ForceMatrixDecomposition::buildNeighborList() { | |
935 | ||
936 | vector<pair<int, int> > neighborList; | |
937 | + | groupCutoffs cuts; |
938 | #ifdef IS_MPI | |
939 | cellListRow_.clear(); | |
940 | cellListCol_.clear(); | |
# | Line 772 | Line 942 | namespace OpenMD { | |
942 | cellList_.clear(); | |
943 | #endif | |
944 | ||
945 | < | // dangerous to not do error checking. |
776 | < | RealType rCut_; |
777 | < | |
778 | < | RealType rList_ = (rCut_ + skinThickness_); |
945 | > | RealType rList_ = (largestRcut_ + skinThickness_); |
946 | RealType rl2 = rList_ * rList_; | |
947 | Snapshot* snap_ = sman_->getCurrentSnapshot(); | |
948 | Mat3x3d Hmat = snap_->getHmat(); | |
# | Line 898 | Line 1065 | namespace OpenMD { | |
1065 | if (m2 != m1 || cgColToGlobal[(*j2)] < cgRowToGlobal[(*j1)]) { | |
1066 | dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; | |
1067 | snap_->wrapVector(dr); | |
1068 | < | if (dr.lengthSquare() < rl2) { |
1068 | > | cuts = getGroupCutoffs( (*j1), (*j2) ); |
1069 | > | if (dr.lengthSquare() < cuts.third) { |
1070 | neighborList.push_back(make_pair((*j1), (*j2))); | |
1071 | } | |
1072 | } | |
# | Line 917 | Line 1085 | namespace OpenMD { | |
1085 | if (m2 != m1 || (*j2) < (*j1)) { | |
1086 | dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)]; | |
1087 | snap_->wrapVector(dr); | |
1088 | < | if (dr.lengthSquare() < rl2) { |
1088 | > | cuts = getGroupCutoffs( (*j1), (*j2) ); |
1089 | > | if (dr.lengthSquare() < cuts.third) { |
1090 | neighborList.push_back(make_pair((*j1), (*j2))); | |
1091 | } | |
1092 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |