| 35 |
|
* |
| 36 |
|
* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). |
| 37 |
|
* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). |
| 38 |
< |
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). |
| 39 |
< |
* [4] Vardeman & Gezelter, in progress (2009). |
| 38 |
> |
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). |
| 39 |
> |
* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
| 40 |
> |
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
| 41 |
|
*/ |
| 42 |
|
|
| 43 |
|
#include "nonbonded/InteractionManager.hpp" |
| 52 |
|
gb_ = new GB(); |
| 53 |
|
sticky_ = new Sticky(); |
| 54 |
|
morse_ = new Morse(); |
| 55 |
+ |
repulsivePower_ = new RepulsivePower(); |
| 56 |
|
eam_ = new EAM(); |
| 57 |
|
sc_ = new SC(); |
| 58 |
|
electrostatic_ = new Electrostatic(); |
| 59 |
|
maw_ = new MAW(); |
| 60 |
|
} |
| 61 |
|
|
| 62 |
+ |
InteractionManager::~InteractionManager() { |
| 63 |
+ |
delete lj_; |
| 64 |
+ |
delete gb_; |
| 65 |
+ |
delete sticky_; |
| 66 |
+ |
delete morse_; |
| 67 |
+ |
delete repulsivePower_; |
| 68 |
+ |
delete eam_; |
| 69 |
+ |
delete sc_; |
| 70 |
+ |
delete electrostatic_; |
| 71 |
+ |
delete maw_; |
| 72 |
+ |
} |
| 73 |
+ |
|
| 74 |
|
void InteractionManager::initialize() { |
| 75 |
< |
|
| 75 |
> |
|
| 76 |
> |
if (initialized_) return; |
| 77 |
> |
|
| 78 |
|
ForceField* forceField_ = info_->getForceField(); |
| 79 |
|
|
| 80 |
|
lj_->setForceField(forceField_); |
| 83 |
|
eam_->setForceField(forceField_); |
| 84 |
|
sc_->setForceField(forceField_); |
| 85 |
|
morse_->setForceField(forceField_); |
| 86 |
+ |
electrostatic_->setSimInfo(info_); |
| 87 |
|
electrostatic_->setForceField(forceField_); |
| 88 |
|
maw_->setForceField(forceField_); |
| 89 |
+ |
repulsivePower_->setForceField(forceField_); |
| 90 |
|
|
| 91 |
|
ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes(); |
| 92 |
|
ForceField::AtomTypeContainer::MapTypeIterator i1, i2; |
| 99 |
|
atype1 = atomTypes->nextType(i1)) { |
| 100 |
|
|
| 101 |
|
// add it to the map: |
| 84 |
– |
AtomTypeProperties atp = atype1->getATP(); |
| 102 |
|
|
| 103 |
|
pair<map<int,AtomType*>::iterator,bool> ret; |
| 104 |
< |
ret = typeMap_.insert( pair<int, AtomType*>(atp.ident, atype1) ); |
| 104 |
> |
ret = typeMap_.insert( pair<int, AtomType*>(atype1->getIdent(), atype1) ); |
| 105 |
|
if (ret.second == false) { |
| 106 |
|
sprintf( painCave.errMsg, |
| 107 |
|
"InteractionManager already had a previous entry with ident %d\n", |
| 108 |
< |
atp.ident); |
| 108 |
> |
atype1->getIdent()); |
| 109 |
|
painCave.severity = OPENMD_INFO; |
| 110 |
|
painCave.isFatal = 0; |
| 111 |
|
simError(); |
| 123 |
|
|
| 124 |
|
bool vdwExplicit = false; |
| 125 |
|
bool metExplicit = false; |
| 126 |
< |
bool hbExplicit = false; |
| 126 |
> |
// bool hbExplicit = false; |
| 127 |
|
|
| 128 |
|
key = make_pair(atype1, atype2); |
| 129 |
|
|
| 194 |
|
interactions_[key].insert(morse_); |
| 195 |
|
vdwExplicit = true; |
| 196 |
|
} |
| 197 |
+ |
|
| 198 |
+ |
if (nbiType->isRepulsivePower()) { |
| 199 |
+ |
if (vdwExplicit) { |
| 200 |
+ |
sprintf( painCave.errMsg, |
| 201 |
+ |
"InteractionManager::initialize found more than one " |
| 202 |
+ |
"explicit \n" |
| 203 |
+ |
"\tvan der Waals interaction for atom types %s - %s\n", |
| 204 |
+ |
atype1->getName().c_str(), atype2->getName().c_str()); |
| 205 |
+ |
painCave.severity = OPENMD_ERROR; |
| 206 |
+ |
painCave.isFatal = 1; |
| 207 |
+ |
simError(); |
| 208 |
+ |
} |
| 209 |
+ |
// We found an explicit RepulsivePower interaction. |
| 210 |
+ |
// override all other vdw entries for this pair of atom types: |
| 211 |
+ |
set<NonBondedInteraction*>::iterator it; |
| 212 |
+ |
for (it = interactions_[key].begin(); |
| 213 |
+ |
it != interactions_[key].end(); ++it) { |
| 214 |
+ |
InteractionFamily ifam = (*it)->getFamily(); |
| 215 |
+ |
if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it); |
| 216 |
+ |
} |
| 217 |
+ |
interactions_[key].insert(repulsivePower_); |
| 218 |
+ |
vdwExplicit = true; |
| 219 |
+ |
} |
| 220 |
|
|
| 221 |
+ |
|
| 222 |
|
if (nbiType->isEAM()) { |
| 223 |
|
// We found an explicit EAM interaction. |
| 224 |
|
// override all other metallic entries for this pair of atom types: |
| 282 |
|
} |
| 283 |
|
|
| 284 |
|
|
| 285 |
< |
// make sure every pair of atom types in this simulation has a |
| 286 |
< |
// non-bonded interaction: |
| 285 |
> |
// Make sure every pair of atom types in this simulation has a |
| 286 |
> |
// non-bonded interaction. If not, just inform the user. |
| 287 |
|
|
| 288 |
|
set<AtomType*> simTypes = info_->getSimulatedAtomTypes(); |
| 289 |
|
set<AtomType*>::iterator it, jt; |
| 290 |
+ |
|
| 291 |
|
for (it = simTypes.begin(); it != simTypes.end(); ++it) { |
| 292 |
|
atype1 = (*it); |
| 293 |
< |
for (jt = simTypes.begin(); jt != simTypes.end(); ++jt) { |
| 293 |
> |
for (jt = it; jt != simTypes.end(); ++jt) { |
| 294 |
|
atype2 = (*jt); |
| 295 |
|
key = make_pair(atype1, atype2); |
| 296 |
|
|
| 297 |
|
if (interactions_[key].size() == 0) { |
| 298 |
|
sprintf( painCave.errMsg, |
| 299 |
< |
"InteractionManager unable to find an appropriate non-bonded\n" |
| 300 |
< |
"\tinteraction for atom types %s - %s\n", |
| 299 |
> |
"InteractionManager could not find a matching non-bonded\n" |
| 300 |
> |
"\tinteraction for atom types %s - %s\n" |
| 301 |
> |
"\tProceeding without this interaction.\n", |
| 302 |
|
atype1->getName().c_str(), atype2->getName().c_str()); |
| 303 |
|
painCave.severity = OPENMD_INFO; |
| 304 |
< |
painCave.isFatal = 1; |
| 304 |
> |
painCave.isFatal = 0; |
| 305 |
|
simError(); |
| 306 |
|
} |
| 307 |
|
} |
| 309 |
|
|
| 310 |
|
initialized_ = true; |
| 311 |
|
} |
| 312 |
< |
|
| 312 |
> |
|
| 313 |
> |
void InteractionManager::setCutoffRadius(RealType rcut) { |
| 314 |
> |
|
| 315 |
> |
electrostatic_->setCutoffRadius(rcut); |
| 316 |
> |
eam_->setCutoffRadius(rcut); |
| 317 |
> |
} |
| 318 |
> |
|
| 319 |
|
void InteractionManager::doPrePair(InteractionData idat){ |
| 320 |
|
|
| 321 |
|
if (!initialized_) initialize(); |
| 322 |
|
|
| 323 |
+ |
// excluded interaction, so just return |
| 324 |
+ |
if (idat.excluded) return; |
| 325 |
+ |
|
| 326 |
|
set<NonBondedInteraction*>::iterator it; |
| 327 |
|
|
| 328 |
|
for (it = interactions_[ idat.atypes ].begin(); |
| 354 |
|
void InteractionManager::doPair(InteractionData idat){ |
| 355 |
|
|
| 356 |
|
if (!initialized_) initialize(); |
| 357 |
< |
|
| 357 |
> |
|
| 358 |
|
set<NonBondedInteraction*>::iterator it; |
| 359 |
|
|
| 360 |
|
for (it = interactions_[ idat.atypes ].begin(); |
| 361 |
< |
it != interactions_[ idat.atypes ].end(); ++it) |
| 310 |
< |
(*it)->calcForce(idat); |
| 311 |
< |
|
| 312 |
< |
return; |
| 313 |
< |
} |
| 361 |
> |
it != interactions_[ idat.atypes ].end(); ++it) { |
| 362 |
|
|
| 363 |
< |
void InteractionManager::doSkipCorrection(InteractionData idat){ |
| 363 |
> |
// electrostatics still has to worry about indirect |
| 364 |
> |
// contributions from excluded pairs of atoms: |
| 365 |
|
|
| 366 |
< |
if (!initialized_) initialize(); |
| 367 |
< |
|
| 319 |
< |
set<NonBondedInteraction*>::iterator it; |
| 320 |
< |
|
| 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); |
| 366 |
> |
if (!idat.excluded || (*it)->getFamily() == ELECTROSTATIC_FAMILY) { |
| 367 |
> |
(*it)->calcForce(idat); |
| 368 |
|
} |
| 369 |
|
} |
| 370 |
|
|
| 389 |
|
|
| 390 |
|
RealType InteractionManager::getSuggestedCutoffRadius(int *atid) { |
| 391 |
|
if (!initialized_) initialize(); |
| 392 |
< |
|
| 392 |
> |
|
| 393 |
|
AtomType* atype = typeMap_[*atid]; |
| 394 |
|
|
| 395 |
|
pair<AtomType*, AtomType*> key = make_pair(atype, atype); |
