| 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). |
| 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(); |
| 60 |
|
} |
| 61 |
|
|
| 62 |
|
void InteractionManager::initialize() { |
| 63 |
< |
|
| 63 |
> |
|
| 64 |
> |
if (initialized_) return; |
| 65 |
> |
|
| 66 |
|
ForceField* forceField_ = info_->getForceField(); |
| 67 |
|
|
| 68 |
|
lj_->setForceField(forceField_); |
| 74 |
|
electrostatic_->setSimInfo(info_); |
| 75 |
|
electrostatic_->setForceField(forceField_); |
| 76 |
|
maw_->setForceField(forceField_); |
| 77 |
+ |
repulsivePower_->setForceField(forceField_); |
| 78 |
|
|
| 79 |
|
ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes(); |
| 80 |
|
ForceField::AtomTypeContainer::MapTypeIterator i1, i2; |
| 87 |
|
atype1 = atomTypes->nextType(i1)) { |
| 88 |
|
|
| 89 |
|
// add it to the map: |
| 85 |
– |
AtomTypeProperties atp = atype1->getATP(); |
| 90 |
|
|
| 91 |
|
pair<map<int,AtomType*>::iterator,bool> ret; |
| 92 |
< |
ret = typeMap_.insert( pair<int, AtomType*>(atp.ident, atype1) ); |
| 92 |
> |
ret = typeMap_.insert( pair<int, AtomType*>(atype1->getIdent(), atype1) ); |
| 93 |
|
if (ret.second == false) { |
| 94 |
|
sprintf( painCave.errMsg, |
| 95 |
|
"InteractionManager already had a previous entry with ident %d\n", |
| 96 |
< |
atp.ident); |
| 96 |
> |
atype1->getIdent()); |
| 97 |
|
painCave.severity = OPENMD_INFO; |
| 98 |
|
painCave.isFatal = 0; |
| 99 |
|
simError(); |
| 182 |
|
interactions_[key].insert(morse_); |
| 183 |
|
vdwExplicit = true; |
| 184 |
|
} |
| 185 |
+ |
|
| 186 |
+ |
if (nbiType->isRepulsivePower()) { |
| 187 |
+ |
if (vdwExplicit) { |
| 188 |
+ |
sprintf( painCave.errMsg, |
| 189 |
+ |
"InteractionManager::initialize found more than one " |
| 190 |
+ |
"explicit \n" |
| 191 |
+ |
"\tvan der Waals interaction for atom types %s - %s\n", |
| 192 |
+ |
atype1->getName().c_str(), atype2->getName().c_str()); |
| 193 |
+ |
painCave.severity = OPENMD_ERROR; |
| 194 |
+ |
painCave.isFatal = 1; |
| 195 |
+ |
simError(); |
| 196 |
+ |
} |
| 197 |
+ |
// We found an explicit RepulsivePower interaction. |
| 198 |
+ |
// override all other vdw entries for this pair of atom types: |
| 199 |
+ |
set<NonBondedInteraction*>::iterator it; |
| 200 |
+ |
for (it = interactions_[key].begin(); |
| 201 |
+ |
it != interactions_[key].end(); ++it) { |
| 202 |
+ |
InteractionFamily ifam = (*it)->getFamily(); |
| 203 |
+ |
if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it); |
| 204 |
+ |
} |
| 205 |
+ |
interactions_[key].insert(repulsivePower_); |
| 206 |
+ |
vdwExplicit = true; |
| 207 |
+ |
} |
| 208 |
|
|
| 209 |
+ |
|
| 210 |
|
if (nbiType->isEAM()) { |
| 211 |
|
// We found an explicit EAM interaction. |
| 212 |
|
// override all other metallic entries for this pair of atom types: |
| 299 |
|
} |
| 300 |
|
|
| 301 |
|
void InteractionManager::setCutoffRadius(RealType rcut) { |
| 302 |
+ |
|
| 303 |
|
electrostatic_->setCutoffRadius(rcut); |
| 304 |
|
eam_->setCutoffRadius(rcut); |
| 305 |
|
} |
| 312 |
|
|
| 313 |
|
if (!initialized_) initialize(); |
| 314 |
|
|
| 315 |
+ |
// excluded interaction, so just return |
| 316 |
+ |
if (idat.excluded) return; |
| 317 |
+ |
|
| 318 |
|
set<NonBondedInteraction*>::iterator it; |
| 319 |
|
|
| 320 |
|
for (it = interactions_[ idat.atypes ].begin(); |
| 346 |
|
void InteractionManager::doPair(InteractionData idat){ |
| 347 |
|
|
| 348 |
|
if (!initialized_) initialize(); |
| 349 |
< |
|
| 349 |
> |
|
| 350 |
|
set<NonBondedInteraction*>::iterator it; |
| 351 |
|
|
| 352 |
|
for (it = interactions_[ idat.atypes ].begin(); |
| 353 |
< |
it != interactions_[ idat.atypes ].end(); ++it) |
| 322 |
< |
(*it)->calcForce(idat); |
| 323 |
< |
|
| 324 |
< |
return; |
| 325 |
< |
} |
| 353 |
> |
it != interactions_[ idat.atypes ].end(); ++it) { |
| 354 |
|
|
| 355 |
< |
void InteractionManager::doSkipCorrection(InteractionData idat){ |
| 355 |
> |
// electrostatics still has to worry about indirect |
| 356 |
> |
// contributions from excluded pairs of atoms: |
| 357 |
|
|
| 358 |
< |
if (!initialized_) initialize(); |
| 359 |
< |
|
| 331 |
< |
set<NonBondedInteraction*>::iterator it; |
| 332 |
< |
|
| 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); |
| 358 |
> |
if (!idat.excluded || (*it)->getFamily() == ELECTROSTATIC_FAMILY) { |
| 359 |
> |
(*it)->calcForce(idat); |
| 360 |
|
} |
| 361 |
|
} |
| 362 |
|
|
| 381 |
|
|
| 382 |
|
RealType InteractionManager::getSuggestedCutoffRadius(int *atid) { |
| 383 |
|
if (!initialized_) initialize(); |
| 384 |
< |
|
| 384 |
> |
|
| 385 |
|
AtomType* atype = typeMap_[*atid]; |
| 386 |
|
|
| 387 |
|
pair<AtomType*, AtomType*> key = make_pair(atype, atype); |
