43 |
|
|
44 |
|
namespace OpenMD { |
45 |
|
|
46 |
< |
InteractionManager* InteractionManager::_instance = NULL; |
47 |
< |
SimInfo* InteractionManager::info_ = NULL; |
48 |
< |
bool InteractionManager::initialized_ = false; |
46 |
> |
InteractionManager::InteractionManager() { |
47 |
|
|
48 |
< |
RealType InteractionManager::rCut_ = 0.0; |
49 |
< |
RealType InteractionManager::rSwitch_ = 0.0; |
50 |
< |
CutoffMethod InteractionManager::cutoffMethod_ = SHIFTED_FORCE; |
51 |
< |
SwitchingFunctionType InteractionManager::sft_ = cubic; |
52 |
< |
RealType InteractionManager::vdwScale_[4] = {1.0, 0.0, 0.0, 0.0}; |
53 |
< |
RealType InteractionManager::electrostaticScale_[4] = {1.0, 0.0, 0.0, 0.0}; |
54 |
< |
|
55 |
< |
map<int, AtomType*> InteractionManager::typeMap_; |
56 |
< |
map<pair<AtomType*, AtomType*>, set<NonBondedInteraction*> > InteractionManager::interactions_; |
57 |
< |
|
60 |
< |
LJ* InteractionManager::lj_ = new LJ(); |
61 |
< |
GB* InteractionManager::gb_ = new GB(); |
62 |
< |
Sticky* InteractionManager::sticky_ = new Sticky(); |
63 |
< |
Morse* InteractionManager::morse_ = new Morse(); |
64 |
< |
EAM* InteractionManager::eam_ = new EAM(); |
65 |
< |
SC* InteractionManager::sc_ = new SC(); |
66 |
< |
Electrostatic* InteractionManager::electrostatic_ = new Electrostatic(); |
67 |
< |
MAW* InteractionManager::maw_ = new MAW(); |
68 |
< |
SwitchingFunction* InteractionManager::switcher_ = new SwitchingFunction(); |
69 |
< |
|
70 |
< |
InteractionManager* InteractionManager::Instance() { |
71 |
< |
if (!_instance) { |
72 |
< |
_instance = new InteractionManager(); |
73 |
< |
} |
74 |
< |
return _instance; |
48 |
> |
initialized_ = false; |
49 |
> |
|
50 |
> |
lj_ = new LJ(); |
51 |
> |
gb_ = new GB(); |
52 |
> |
sticky_ = new Sticky(); |
53 |
> |
morse_ = new Morse(); |
54 |
> |
eam_ = new EAM(); |
55 |
> |
sc_ = new SC(); |
56 |
> |
electrostatic_ = new Electrostatic(); |
57 |
> |
maw_ = new MAW(); |
58 |
|
} |
59 |
|
|
60 |
|
void InteractionManager::initialize() { |
67 |
|
eam_->setForceField(forceField_); |
68 |
|
sc_->setForceField(forceField_); |
69 |
|
morse_->setForceField(forceField_); |
70 |
+ |
electrostatic_->setSimInfo(info_); |
71 |
|
electrostatic_->setForceField(forceField_); |
72 |
|
maw_->setForceField(forceField_); |
73 |
|
|
90 |
– |
ForceFieldOptions& fopts = forceField_->getForceFieldOptions(); |
91 |
– |
|
92 |
– |
// Force fields can set options on how to scale van der Waals and electrostatic |
93 |
– |
// interactions for atoms connected via bonds, bends and torsions |
94 |
– |
// in this case the topological distance between atoms is: |
95 |
– |
// 0 = topologically unconnected |
96 |
– |
// 1 = bonded together |
97 |
– |
// 2 = connected via a bend |
98 |
– |
// 3 = connected via a torsion |
99 |
– |
|
100 |
– |
vdwScale_[0] = 1.0; |
101 |
– |
vdwScale_[1] = fopts.getvdw12scale(); |
102 |
– |
vdwScale_[2] = fopts.getvdw13scale(); |
103 |
– |
vdwScale_[3] = fopts.getvdw14scale(); |
104 |
– |
|
105 |
– |
electrostaticScale_[0] = 1.0; |
106 |
– |
electrostaticScale_[1] = fopts.getelectrostatic12scale(); |
107 |
– |
electrostaticScale_[2] = fopts.getelectrostatic13scale(); |
108 |
– |
electrostaticScale_[3] = fopts.getelectrostatic14scale(); |
109 |
– |
|
74 |
|
ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes(); |
75 |
|
ForceField::AtomTypeContainer::MapTypeIterator i1, i2; |
76 |
|
AtomType* atype1; |
242 |
|
} |
243 |
|
|
244 |
|
|
245 |
< |
// make sure every pair of atom types in this simulation has a |
246 |
< |
// non-bonded interaction: |
245 |
> |
// Make sure every pair of atom types in this simulation has a |
246 |
> |
// non-bonded interaction. If not, just inform the user. |
247 |
|
|
248 |
|
set<AtomType*> simTypes = info_->getSimulatedAtomTypes(); |
249 |
|
set<AtomType*>::iterator it, jt; |
250 |
+ |
|
251 |
|
for (it = simTypes.begin(); it != simTypes.end(); ++it) { |
252 |
|
atype1 = (*it); |
253 |
< |
for (jt = simTypes.begin(); jt != simTypes.end(); ++jt) { |
253 |
> |
for (jt = it; jt != simTypes.end(); ++jt) { |
254 |
|
atype2 = (*jt); |
255 |
|
key = make_pair(atype1, atype2); |
256 |
|
|
257 |
|
if (interactions_[key].size() == 0) { |
258 |
|
sprintf( painCave.errMsg, |
259 |
< |
"InteractionManager unable to find an appropriate non-bonded\n" |
260 |
< |
"\tinteraction for atom types %s - %s\n", |
259 |
> |
"InteractionManager could not find a matching non-bonded\n" |
260 |
> |
"\tinteraction for atom types %s - %s\n" |
261 |
> |
"\tProceeding without this interaction.\n", |
262 |
|
atype1->getName().c_str(), atype2->getName().c_str()); |
263 |
|
painCave.severity = OPENMD_INFO; |
264 |
< |
painCave.isFatal = 1; |
264 |
> |
painCave.isFatal = 0; |
265 |
|
simError(); |
266 |
|
} |
267 |
|
} |
268 |
|
} |
269 |
|
|
270 |
< |
setupCutoffs(); |
271 |
< |
setupSwitching(); |
270 |
> |
initialized_ = true; |
271 |
> |
} |
272 |
|
|
273 |
< |
//int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; |
274 |
< |
//int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; |
275 |
< |
//notifyFortranCutoffs(&rCut_, &rSwitch_, &ljsp, &ljsf); |
273 |
> |
void InteractionManager::setCutoffRadius(RealType rcut) { |
274 |
> |
electrostatic_->setCutoffRadius(rcut); |
275 |
> |
eam_->setCutoffRadius(rcut); |
276 |
> |
} |
277 |
|
|
278 |
< |
initialized_ = true; |
278 |
> |
void InteractionManager::setSwitchingRadius(RealType rswitch) { |
279 |
> |
electrostatic_->setSwitchingRadius(rswitch); |
280 |
|
} |
281 |
|
|
314 |
– |
/** |
315 |
– |
* setupCutoffs |
316 |
– |
* |
317 |
– |
* Sets the values of cutoffRadius and cutoffMethod |
318 |
– |
* |
319 |
– |
* cutoffRadius : realType |
320 |
– |
* If the cutoffRadius was explicitly set, use that value. |
321 |
– |
* If the cutoffRadius was not explicitly set: |
322 |
– |
* Are there electrostatic atoms? Use 12.0 Angstroms. |
323 |
– |
* No electrostatic atoms? Poll the atom types present in the |
324 |
– |
* simulation for suggested cutoff values (e.g. 2.5 * sigma). |
325 |
– |
* Use the maximum suggested value that was found. |
326 |
– |
* |
327 |
– |
* cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL) |
328 |
– |
* If cutoffMethod was explicitly set, use that choice. |
329 |
– |
* If cutoffMethod was not explicitly set, use SHIFTED_FORCE |
330 |
– |
*/ |
331 |
– |
void InteractionManager::setupCutoffs() { |
332 |
– |
|
333 |
– |
Globals* simParams_ = info_->getSimParams(); |
334 |
– |
|
335 |
– |
if (simParams_->haveCutoffRadius()) { |
336 |
– |
rCut_ = simParams_->getCutoffRadius(); |
337 |
– |
} else { |
338 |
– |
if (info_->usesElectrostaticAtoms()) { |
339 |
– |
sprintf(painCave.errMsg, |
340 |
– |
"InteractionManager::setupCutoffs: No value was set for the cutoffRadius.\n" |
341 |
– |
"\tOpenMD will use a default value of 12.0 angstroms" |
342 |
– |
"\tfor the cutoffRadius.\n"); |
343 |
– |
painCave.isFatal = 0; |
344 |
– |
painCave.severity = OPENMD_INFO; |
345 |
– |
simError(); |
346 |
– |
rCut_ = 12.0; |
347 |
– |
} else { |
348 |
– |
RealType thisCut; |
349 |
– |
set<AtomType*>::iterator i; |
350 |
– |
set<AtomType*> atomTypes; |
351 |
– |
atomTypes = info_->getSimulatedAtomTypes(); |
352 |
– |
for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
353 |
– |
thisCut = getSuggestedCutoffRadius((*i)); |
354 |
– |
rCut_ = max(thisCut, rCut_); |
355 |
– |
} |
356 |
– |
sprintf(painCave.errMsg, |
357 |
– |
"InteractionManager::setupCutoffs: No value was set for the cutoffRadius.\n" |
358 |
– |
"\tOpenMD will use %lf angstroms.\n", |
359 |
– |
rCut_); |
360 |
– |
painCave.isFatal = 0; |
361 |
– |
painCave.severity = OPENMD_INFO; |
362 |
– |
simError(); |
363 |
– |
} |
364 |
– |
} |
365 |
– |
|
366 |
– |
map<string, CutoffMethod> stringToCutoffMethod; |
367 |
– |
stringToCutoffMethod["HARD"] = HARD; |
368 |
– |
stringToCutoffMethod["SWITCHED"] = SWITCHED; |
369 |
– |
stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL; |
370 |
– |
stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE; |
371 |
– |
|
372 |
– |
if (simParams_->haveCutoffMethod()) { |
373 |
– |
string cutMeth = toUpperCopy(simParams_->getCutoffMethod()); |
374 |
– |
map<string, CutoffMethod>::iterator i; |
375 |
– |
i = stringToCutoffMethod.find(cutMeth); |
376 |
– |
if (i == stringToCutoffMethod.end()) { |
377 |
– |
sprintf(painCave.errMsg, |
378 |
– |
"InteractionManager::setupCutoffs: Could not find chosen cutoffMethod %s\n" |
379 |
– |
"\tShould be one of: " |
380 |
– |
"HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n", |
381 |
– |
cutMeth.c_str()); |
382 |
– |
painCave.isFatal = 1; |
383 |
– |
painCave.severity = OPENMD_ERROR; |
384 |
– |
simError(); |
385 |
– |
} else { |
386 |
– |
cutoffMethod_ = i->second; |
387 |
– |
} |
388 |
– |
} else { |
389 |
– |
sprintf(painCave.errMsg, |
390 |
– |
"InteractionManager::setupCutoffs: No value was set for the cutoffMethod.\n" |
391 |
– |
"\tOpenMD will use SHIFTED_FORCE.\n"); |
392 |
– |
painCave.isFatal = 0; |
393 |
– |
painCave.severity = OPENMD_INFO; |
394 |
– |
simError(); |
395 |
– |
cutoffMethod_ = SHIFTED_FORCE; |
396 |
– |
} |
397 |
– |
} |
398 |
– |
|
399 |
– |
|
400 |
– |
/** |
401 |
– |
* setupSwitching |
402 |
– |
* |
403 |
– |
* Sets the values of switchingRadius and |
404 |
– |
* If the switchingRadius was explicitly set, use that value (but check it) |
405 |
– |
* If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_ |
406 |
– |
*/ |
407 |
– |
void InteractionManager::setupSwitching() { |
408 |
– |
Globals* simParams_ = info_->getSimParams(); |
409 |
– |
|
410 |
– |
if (simParams_->haveSwitchingRadius()) { |
411 |
– |
rSwitch_ = simParams_->getSwitchingRadius(); |
412 |
– |
if (rSwitch_ > rCut_) { |
413 |
– |
sprintf(painCave.errMsg, |
414 |
– |
"InteractionManager::setupSwitching: switchingRadius (%f) is larger than cutoffRadius(%f)\n", |
415 |
– |
rSwitch_, rCut_); |
416 |
– |
painCave.isFatal = 1; |
417 |
– |
painCave.severity = OPENMD_ERROR; |
418 |
– |
simError(); |
419 |
– |
} |
420 |
– |
} else { |
421 |
– |
rSwitch_ = 0.85 * rCut_; |
422 |
– |
sprintf(painCave.errMsg, |
423 |
– |
"InteractionManager::setupSwitching: No value was set for the switchingRadius.\n" |
424 |
– |
"\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" |
425 |
– |
"\tswitchingRadius = %f. for this simulation\n", rSwitch_); |
426 |
– |
painCave.isFatal = 0; |
427 |
– |
painCave.severity = OPENMD_WARNING; |
428 |
– |
simError(); |
429 |
– |
} |
430 |
– |
|
431 |
– |
if (simParams_->haveSwitchingFunctionType()) { |
432 |
– |
string funcType = simParams_->getSwitchingFunctionType(); |
433 |
– |
toUpper(funcType); |
434 |
– |
if (funcType == "CUBIC") { |
435 |
– |
sft_ = cubic; |
436 |
– |
} else { |
437 |
– |
if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
438 |
– |
sft_ = fifth_order_poly; |
439 |
– |
} else { |
440 |
– |
// throw error |
441 |
– |
sprintf( painCave.errMsg, |
442 |
– |
"InteractionManager::setupSwitching : Unknown switchingFunctionType. (Input file specified %s .)\n" |
443 |
– |
"\tswitchingFunctionType must be one of: " |
444 |
– |
"\"cubic\" or \"fifth_order_polynomial\".", |
445 |
– |
funcType.c_str() ); |
446 |
– |
painCave.isFatal = 1; |
447 |
– |
painCave.severity = OPENMD_ERROR; |
448 |
– |
simError(); |
449 |
– |
} |
450 |
– |
} |
451 |
– |
} |
452 |
– |
|
453 |
– |
switcher_->setSwitchType(sft_); |
454 |
– |
switcher_->setSwitch(rSwitch_, rCut_); |
455 |
– |
} |
456 |
– |
|
282 |
|
void InteractionManager::doPrePair(InteractionData idat){ |
283 |
|
|
284 |
|
if (!initialized_) initialize(); |
285 |
|
|
286 |
|
set<NonBondedInteraction*>::iterator it; |
287 |
|
|
288 |
< |
for (it = interactions_[ *(idat.atypes) ].begin(); |
289 |
< |
it != interactions_[ *(idat.atypes) ].end(); ++it){ |
288 |
> |
for (it = interactions_[ idat.atypes ].begin(); |
289 |
> |
it != interactions_[ idat.atypes ].end(); ++it){ |
290 |
|
if ((*it)->getFamily() == METALLIC_FAMILY) { |
291 |
|
dynamic_cast<MetallicInteraction*>(*it)->calcDensity(idat); |
292 |
|
} |
317 |
|
|
318 |
|
set<NonBondedInteraction*>::iterator it; |
319 |
|
|
320 |
< |
for (it = interactions_[ *(idat.atypes) ].begin(); |
321 |
< |
it != interactions_[ *(idat.atypes) ].end(); ++it) |
320 |
> |
for (it = interactions_[ idat.atypes ].begin(); |
321 |
> |
it != interactions_[ idat.atypes ].end(); ++it) |
322 |
|
(*it)->calcForce(idat); |
323 |
|
|
324 |
|
return; |
327 |
|
void InteractionManager::doSkipCorrection(InteractionData idat){ |
328 |
|
|
329 |
|
if (!initialized_) initialize(); |
330 |
< |
|
330 |
> |
|
331 |
|
set<NonBondedInteraction*>::iterator it; |
332 |
|
|
333 |
< |
for (it = interactions_[ *(idat.atypes) ].begin(); |
334 |
< |
it != interactions_[ *(idat.atypes) ].end(); ++it){ |
333 |
> |
for (it = interactions_[ idat.atypes ].begin(); |
334 |
> |
it != interactions_[ idat.atypes ].end(); ++it){ |
335 |
|
if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) { |
336 |
|
dynamic_cast<ElectrostaticInteraction*>(*it)->calcSkipCorrection(idat); |
337 |
|
} |