40 |
|
*/ |
41 |
|
|
42 |
|
#include "nonbonded/InteractionManager.hpp" |
43 |
– |
#include "UseTheForce/doForces_interface.h" |
43 |
|
|
44 |
|
namespace OpenMD { |
45 |
|
|
46 |
< |
InteractionManager* InteractionManager::_instance = NULL; |
48 |
< |
SimInfo* InteractionManager::info_ = NULL; |
49 |
< |
bool InteractionManager::initialized_ = false; |
46 |
> |
InteractionManager::InteractionManager() { |
47 |
|
|
48 |
< |
RealType InteractionManager::rCut_ = 0.0; |
49 |
< |
RealType InteractionManager::rSwitch_ = 0.0; |
50 |
< |
RealType InteractionManager::skinThickness_ = 0.0; |
51 |
< |
RealType InteractionManager::listRadius_ = 0.0; |
52 |
< |
CutoffMethod InteractionManager::cutoffMethod_ = SHIFTED_FORCE; |
53 |
< |
SwitchingFunctionType InteractionManager::sft_ = cubic; |
54 |
< |
RealType InteractionManager::vdwScale_[4] = {1.0, 0.0, 0.0, 0.0}; |
55 |
< |
RealType InteractionManager::electrostaticScale_[4] = {1.0, 0.0, 0.0, 0.0}; |
56 |
< |
|
57 |
< |
map<int, AtomType*> InteractionManager::typeMap_; |
61 |
< |
map<pair<AtomType*, AtomType*>, set<NonBondedInteraction*> > InteractionManager::interactions_; |
62 |
< |
|
63 |
< |
LJ* InteractionManager::lj_ = new LJ(); |
64 |
< |
GB* InteractionManager::gb_ = new GB(); |
65 |
< |
Sticky* InteractionManager::sticky_ = new Sticky(); |
66 |
< |
Morse* InteractionManager::morse_ = new Morse(); |
67 |
< |
EAM* InteractionManager::eam_ = new EAM(); |
68 |
< |
SC* InteractionManager::sc_ = new SC(); |
69 |
< |
Electrostatic* InteractionManager::electrostatic_ = new Electrostatic(); |
70 |
< |
MAW* InteractionManager::maw_ = new MAW(); |
71 |
< |
SwitchingFunction* InteractionManager::switcher_ = new SwitchingFunction(); |
72 |
< |
|
73 |
< |
InteractionManager* InteractionManager::Instance() { |
74 |
< |
if (!_instance) { |
75 |
< |
_instance = new InteractionManager(); |
76 |
< |
} |
77 |
< |
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() { |
70 |
|
electrostatic_->setForceField(forceField_); |
71 |
|
maw_->setForceField(forceField_); |
72 |
|
|
93 |
– |
ForceFieldOptions& fopts = forceField_->getForceFieldOptions(); |
94 |
– |
|
95 |
– |
// Force fields can set options on how to scale van der Waals and electrostatic |
96 |
– |
// interactions for atoms connected via bonds, bends and torsions |
97 |
– |
// in this case the topological distance between atoms is: |
98 |
– |
// 0 = topologically unconnected |
99 |
– |
// 1 = bonded together |
100 |
– |
// 2 = connected via a bend |
101 |
– |
// 3 = connected via a torsion |
102 |
– |
|
103 |
– |
vdwScale_[0] = 1.0; |
104 |
– |
vdwScale_[1] = fopts.getvdw12scale(); |
105 |
– |
vdwScale_[2] = fopts.getvdw13scale(); |
106 |
– |
vdwScale_[3] = fopts.getvdw14scale(); |
107 |
– |
|
108 |
– |
electrostaticScale_[0] = 1.0; |
109 |
– |
electrostaticScale_[1] = fopts.getelectrostatic12scale(); |
110 |
– |
electrostaticScale_[2] = fopts.getelectrostatic13scale(); |
111 |
– |
electrostaticScale_[3] = fopts.getelectrostatic14scale(); |
112 |
– |
|
73 |
|
ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes(); |
74 |
|
ForceField::AtomTypeContainer::MapTypeIterator i1, i2; |
75 |
|
AtomType* atype1; |
263 |
|
} |
264 |
|
} |
265 |
|
} |
306 |
– |
|
307 |
– |
setupCutoffs(); |
308 |
– |
setupSwitching(); |
309 |
– |
setupNeighborlists(); |
266 |
|
|
311 |
– |
//int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; |
312 |
– |
//int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; |
313 |
– |
//notifyFortranCutoffs(&rCut_, &rSwitch_, &ljsp, &ljsf); |
314 |
– |
//notifyFortranSkinThickness(&skinThickness_); |
315 |
– |
|
267 |
|
initialized_ = true; |
268 |
|
} |
269 |
|
|
319 |
– |
/** |
320 |
– |
* setupCutoffs |
321 |
– |
* |
322 |
– |
* Sets the values of cutoffRadius and cutoffMethod |
323 |
– |
* |
324 |
– |
* cutoffRadius : realType |
325 |
– |
* If the cutoffRadius was explicitly set, use that value. |
326 |
– |
* If the cutoffRadius was not explicitly set: |
327 |
– |
* Are there electrostatic atoms? Use 12.0 Angstroms. |
328 |
– |
* No electrostatic atoms? Poll the atom types present in the |
329 |
– |
* simulation for suggested cutoff values (e.g. 2.5 * sigma). |
330 |
– |
* Use the maximum suggested value that was found. |
331 |
– |
* |
332 |
– |
* cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL) |
333 |
– |
* If cutoffMethod was explicitly set, use that choice. |
334 |
– |
* If cutoffMethod was not explicitly set, use SHIFTED_FORCE |
335 |
– |
*/ |
336 |
– |
void InteractionManager::setupCutoffs() { |
337 |
– |
|
338 |
– |
Globals* simParams_ = info_->getSimParams(); |
339 |
– |
|
340 |
– |
if (simParams_->haveCutoffRadius()) { |
341 |
– |
rCut_ = simParams_->getCutoffRadius(); |
342 |
– |
} else { |
343 |
– |
if (info_->usesElectrostaticAtoms()) { |
344 |
– |
sprintf(painCave.errMsg, |
345 |
– |
"InteractionManager::setupCutoffs: No value was set for the cutoffRadius.\n" |
346 |
– |
"\tOpenMD will use a default value of 12.0 angstroms" |
347 |
– |
"\tfor the cutoffRadius.\n"); |
348 |
– |
painCave.isFatal = 0; |
349 |
– |
painCave.severity = OPENMD_INFO; |
350 |
– |
simError(); |
351 |
– |
rCut_ = 12.0; |
352 |
– |
} else { |
353 |
– |
RealType thisCut; |
354 |
– |
set<AtomType*>::iterator i; |
355 |
– |
set<AtomType*> atomTypes; |
356 |
– |
atomTypes = info_->getSimulatedAtomTypes(); |
357 |
– |
for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
358 |
– |
thisCut = getSuggestedCutoffRadius((*i)); |
359 |
– |
rCut_ = max(thisCut, rCut_); |
360 |
– |
} |
361 |
– |
sprintf(painCave.errMsg, |
362 |
– |
"InteractionManager::setupCutoffs: No value was set for the cutoffRadius.\n" |
363 |
– |
"\tOpenMD will use %lf angstroms.\n", |
364 |
– |
rCut_); |
365 |
– |
painCave.isFatal = 0; |
366 |
– |
painCave.severity = OPENMD_INFO; |
367 |
– |
simError(); |
368 |
– |
} |
369 |
– |
} |
370 |
– |
|
371 |
– |
map<string, CutoffMethod> stringToCutoffMethod; |
372 |
– |
stringToCutoffMethod["HARD"] = HARD; |
373 |
– |
stringToCutoffMethod["SWITCHED"] = SWITCHED; |
374 |
– |
stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL; |
375 |
– |
stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE; |
376 |
– |
|
377 |
– |
if (simParams_->haveCutoffMethod()) { |
378 |
– |
string cutMeth = toUpperCopy(simParams_->getCutoffMethod()); |
379 |
– |
map<string, CutoffMethod>::iterator i; |
380 |
– |
i = stringToCutoffMethod.find(cutMeth); |
381 |
– |
if (i == stringToCutoffMethod.end()) { |
382 |
– |
sprintf(painCave.errMsg, |
383 |
– |
"InteractionManager::setupCutoffs: Could not find chosen cutoffMethod %s\n" |
384 |
– |
"\tShould be one of: " |
385 |
– |
"HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n", |
386 |
– |
cutMeth.c_str()); |
387 |
– |
painCave.isFatal = 1; |
388 |
– |
painCave.severity = OPENMD_ERROR; |
389 |
– |
simError(); |
390 |
– |
} else { |
391 |
– |
cutoffMethod_ = i->second; |
392 |
– |
} |
393 |
– |
} else { |
394 |
– |
sprintf(painCave.errMsg, |
395 |
– |
"InteractionManager::setupCutoffs: No value was set for the cutoffMethod.\n" |
396 |
– |
"\tOpenMD will use SHIFTED_FORCE.\n"); |
397 |
– |
painCave.isFatal = 0; |
398 |
– |
painCave.severity = OPENMD_INFO; |
399 |
– |
simError(); |
400 |
– |
cutoffMethod_ = SHIFTED_FORCE; |
401 |
– |
} |
402 |
– |
} |
403 |
– |
|
404 |
– |
|
405 |
– |
/** |
406 |
– |
* setupSwitching |
407 |
– |
* |
408 |
– |
* Sets the values of switchingRadius and |
409 |
– |
* If the switchingRadius was explicitly set, use that value (but check it) |
410 |
– |
* If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_ |
411 |
– |
*/ |
412 |
– |
void InteractionManager::setupSwitching() { |
413 |
– |
Globals* simParams_ = info_->getSimParams(); |
414 |
– |
|
415 |
– |
if (simParams_->haveSwitchingRadius()) { |
416 |
– |
rSwitch_ = simParams_->getSwitchingRadius(); |
417 |
– |
if (rSwitch_ > rCut_) { |
418 |
– |
sprintf(painCave.errMsg, |
419 |
– |
"InteractionManager::setupSwitching: switchingRadius (%f) is larger than cutoffRadius(%f)\n", |
420 |
– |
rSwitch_, rCut_); |
421 |
– |
painCave.isFatal = 1; |
422 |
– |
painCave.severity = OPENMD_ERROR; |
423 |
– |
simError(); |
424 |
– |
} |
425 |
– |
} else { |
426 |
– |
rSwitch_ = 0.85 * rCut_; |
427 |
– |
sprintf(painCave.errMsg, |
428 |
– |
"InteractionManager::setupSwitching: No value was set for the switchingRadius.\n" |
429 |
– |
"\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" |
430 |
– |
"\tswitchingRadius = %f. for this simulation\n", rSwitch_); |
431 |
– |
painCave.isFatal = 0; |
432 |
– |
painCave.severity = OPENMD_WARNING; |
433 |
– |
simError(); |
434 |
– |
} |
435 |
– |
|
436 |
– |
if (simParams_->haveSwitchingFunctionType()) { |
437 |
– |
string funcType = simParams_->getSwitchingFunctionType(); |
438 |
– |
toUpper(funcType); |
439 |
– |
if (funcType == "CUBIC") { |
440 |
– |
sft_ = cubic; |
441 |
– |
} else { |
442 |
– |
if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
443 |
– |
sft_ = fifth_order_poly; |
444 |
– |
} else { |
445 |
– |
// throw error |
446 |
– |
sprintf( painCave.errMsg, |
447 |
– |
"InteractionManager::setupSwitching : Unknown switchingFunctionType. (Input file specified %s .)\n" |
448 |
– |
"\tswitchingFunctionType must be one of: " |
449 |
– |
"\"cubic\" or \"fifth_order_polynomial\".", |
450 |
– |
funcType.c_str() ); |
451 |
– |
painCave.isFatal = 1; |
452 |
– |
painCave.severity = OPENMD_ERROR; |
453 |
– |
simError(); |
454 |
– |
} |
455 |
– |
} |
456 |
– |
} |
457 |
– |
|
458 |
– |
switcher_->setSwitchType(sft_); |
459 |
– |
switcher_->setSwitch(rSwitch_, rCut_); |
460 |
– |
} |
461 |
– |
|
462 |
– |
/** |
463 |
– |
* setupNeighborlists |
464 |
– |
* |
465 |
– |
* If the skinThickness was explicitly set, use that value (but check it) |
466 |
– |
* If the skinThickness was not explicitly set: use 1.0 angstroms |
467 |
– |
*/ |
468 |
– |
void InteractionManager::setupNeighborlists() { |
469 |
– |
|
470 |
– |
Globals* simParams_ = info_->getSimParams(); |
471 |
– |
|
472 |
– |
if (simParams_->haveSkinThickness()) { |
473 |
– |
skinThickness_ = simParams_->getSkinThickness(); |
474 |
– |
} else { |
475 |
– |
skinThickness_ = 1.0; |
476 |
– |
sprintf(painCave.errMsg, |
477 |
– |
"InteractionManager::setupNeighborlists: No value was set for the skinThickness.\n" |
478 |
– |
"\tOpenMD will use a default value of %f Angstroms\n" |
479 |
– |
"\tfor this simulation\n", skinThickness_); |
480 |
– |
painCave.severity = OPENMD_INFO; |
481 |
– |
painCave.isFatal = 0; |
482 |
– |
simError(); |
483 |
– |
} |
484 |
– |
|
485 |
– |
listRadius_ = rCut_ + skinThickness_; |
486 |
– |
} |
487 |
– |
|
488 |
– |
|
270 |
|
void InteractionManager::doPrePair(InteractionData idat){ |
271 |
|
|
272 |
|
if (!initialized_) initialize(); |
273 |
|
|
274 |
|
set<NonBondedInteraction*>::iterator it; |
275 |
|
|
276 |
< |
for (it = interactions_[idat.atypes].begin(); |
277 |
< |
it != interactions_[idat.atypes].end(); ++it){ |
276 |
> |
for (it = interactions_[ idat.atypes ].begin(); |
277 |
> |
it != interactions_[ idat.atypes ].end(); ++it){ |
278 |
|
if ((*it)->getFamily() == METALLIC_FAMILY) { |
279 |
|
dynamic_cast<MetallicInteraction*>(*it)->calcDensity(idat); |
280 |
|
} |
305 |
|
|
306 |
|
set<NonBondedInteraction*>::iterator it; |
307 |
|
|
308 |
< |
for (it = interactions_[idat.atypes].begin(); |
309 |
< |
it != interactions_[idat.atypes].end(); ++it) |
308 |
> |
for (it = interactions_[ idat.atypes ].begin(); |
309 |
> |
it != interactions_[ idat.atypes ].end(); ++it) |
310 |
|
(*it)->calcForce(idat); |
311 |
|
|
312 |
|
return; |
318 |
|
|
319 |
|
set<NonBondedInteraction*>::iterator it; |
320 |
|
|
321 |
< |
for (it = interactions_[idat.atypes].begin(); |
322 |
< |
it != interactions_[idat.atypes].end(); ++it){ |
321 |
> |
for (it = interactions_[ idat.atypes ].begin(); |
322 |
> |
it != interactions_[ idat.atypes ].end(); ++it){ |
323 |
|
if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) { |
324 |
|
dynamic_cast<ElectrostaticInteraction*>(*it)->calcSkipCorrection(idat); |
325 |
|
} |