| 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(); |
| 69 |
|
eam_->setForceField(forceField_); |
| 70 |
|
sc_->setForceField(forceField_); |
| 71 |
|
morse_->setForceField(forceField_); |
| 72 |
+ |
electrostatic_->setSimInfo(info_); |
| 73 |
|
electrostatic_->setForceField(forceField_); |
| 74 |
|
maw_->setForceField(forceField_); |
| 75 |
+ |
repulsivePower_->setForceField(forceField_); |
| 76 |
|
|
| 77 |
|
ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes(); |
| 78 |
|
ForceField::AtomTypeContainer::MapTypeIterator i1, i2; |
| 179 |
|
if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it); |
| 180 |
|
} |
| 181 |
|
interactions_[key].insert(morse_); |
| 182 |
+ |
vdwExplicit = true; |
| 183 |
+ |
} |
| 184 |
+ |
|
| 185 |
+ |
if (nbiType->isRepulsivePower()) { |
| 186 |
+ |
if (vdwExplicit) { |
| 187 |
+ |
sprintf( painCave.errMsg, |
| 188 |
+ |
"InteractionManager::initialize found more than one " |
| 189 |
+ |
"explicit \n" |
| 190 |
+ |
"\tvan der Waals interaction for atom types %s - %s\n", |
| 191 |
+ |
atype1->getName().c_str(), atype2->getName().c_str()); |
| 192 |
+ |
painCave.severity = OPENMD_ERROR; |
| 193 |
+ |
painCave.isFatal = 1; |
| 194 |
+ |
simError(); |
| 195 |
+ |
} |
| 196 |
+ |
// We found an explicit RepulsivePower interaction. |
| 197 |
+ |
// override all other vdw entries for this pair of atom types: |
| 198 |
+ |
set<NonBondedInteraction*>::iterator it; |
| 199 |
+ |
for (it = interactions_[key].begin(); |
| 200 |
+ |
it != interactions_[key].end(); ++it) { |
| 201 |
+ |
InteractionFamily ifam = (*it)->getFamily(); |
| 202 |
+ |
if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it); |
| 203 |
+ |
} |
| 204 |
+ |
interactions_[key].insert(repulsivePower_); |
| 205 |
|
vdwExplicit = true; |
| 206 |
|
} |
| 207 |
|
|
| 208 |
+ |
|
| 209 |
|
if (nbiType->isEAM()) { |
| 210 |
|
// We found an explicit EAM interaction. |
| 211 |
|
// override all other metallic entries for this pair of atom types: |
| 269 |
|
} |
| 270 |
|
|
| 271 |
|
|
| 272 |
< |
// make sure every pair of atom types in this simulation has a |
| 273 |
< |
// non-bonded interaction: |
| 272 |
> |
// Make sure every pair of atom types in this simulation has a |
| 273 |
> |
// non-bonded interaction. If not, just inform the user. |
| 274 |
|
|
| 275 |
|
set<AtomType*> simTypes = info_->getSimulatedAtomTypes(); |
| 276 |
|
set<AtomType*>::iterator it, jt; |
| 277 |
+ |
|
| 278 |
|
for (it = simTypes.begin(); it != simTypes.end(); ++it) { |
| 279 |
|
atype1 = (*it); |
| 280 |
< |
for (jt = simTypes.begin(); jt != simTypes.end(); ++jt) { |
| 280 |
> |
for (jt = it; jt != simTypes.end(); ++jt) { |
| 281 |
|
atype2 = (*jt); |
| 282 |
|
key = make_pair(atype1, atype2); |
| 283 |
|
|
| 284 |
|
if (interactions_[key].size() == 0) { |
| 285 |
|
sprintf( painCave.errMsg, |
| 286 |
< |
"InteractionManager unable to find an appropriate non-bonded\n" |
| 287 |
< |
"\tinteraction for atom types %s - %s\n", |
| 286 |
> |
"InteractionManager could not find a matching non-bonded\n" |
| 287 |
> |
"\tinteraction for atom types %s - %s\n" |
| 288 |
> |
"\tProceeding without this interaction.\n", |
| 289 |
|
atype1->getName().c_str(), atype2->getName().c_str()); |
| 290 |
|
painCave.severity = OPENMD_INFO; |
| 291 |
< |
painCave.isFatal = 1; |
| 291 |
> |
painCave.isFatal = 0; |
| 292 |
|
simError(); |
| 293 |
|
} |
| 294 |
|
} |
| 296 |
|
|
| 297 |
|
initialized_ = true; |
| 298 |
|
} |
| 299 |
+ |
|
| 300 |
+ |
void InteractionManager::setCutoffRadius(RealType rcut) { |
| 301 |
+ |
|
| 302 |
+ |
electrostatic_->setCutoffRadius(rcut); |
| 303 |
+ |
eam_->setCutoffRadius(rcut); |
| 304 |
+ |
} |
| 305 |
+ |
|
| 306 |
+ |
void InteractionManager::setSwitchingRadius(RealType rswitch) { |
| 307 |
+ |
electrostatic_->setSwitchingRadius(rswitch); |
| 308 |
+ |
} |
| 309 |
|
|
| 310 |
|
void InteractionManager::doPrePair(InteractionData idat){ |
| 311 |
|
|
| 312 |
|
if (!initialized_) initialize(); |
| 313 |
|
|
| 314 |
+ |
// excluded interaction, so just return |
| 315 |
+ |
if (idat.excluded) return; |
| 316 |
+ |
|
| 317 |
|
set<NonBondedInteraction*>::iterator it; |
| 318 |
|
|
| 319 |
|
for (it = interactions_[ idat.atypes ].begin(); |
| 345 |
|
void InteractionManager::doPair(InteractionData idat){ |
| 346 |
|
|
| 347 |
|
if (!initialized_) initialize(); |
| 348 |
< |
|
| 348 |
> |
|
| 349 |
|
set<NonBondedInteraction*>::iterator it; |
| 350 |
|
|
| 351 |
|
for (it = interactions_[ idat.atypes ].begin(); |
| 352 |
< |
it != interactions_[ idat.atypes ].end(); ++it) |
| 310 |
< |
(*it)->calcForce(idat); |
| 311 |
< |
|
| 312 |
< |
return; |
| 313 |
< |
} |
| 352 |
> |
it != interactions_[ idat.atypes ].end(); ++it) { |
| 353 |
|
|
| 354 |
< |
void InteractionManager::doSkipCorrection(InteractionData idat){ |
| 354 |
> |
// electrostatics still has to worry about indirect |
| 355 |
> |
// contributions from excluded pairs of atoms: |
| 356 |
|
|
| 357 |
< |
if (!initialized_) initialize(); |
| 358 |
< |
|
| 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); |
| 357 |
> |
if (!idat.excluded || (*it)->getFamily() == ELECTROSTATIC_FAMILY) { |
| 358 |
> |
(*it)->calcForce(idat); |
| 359 |
|
} |
| 360 |
|
} |
| 361 |
|
|
