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_); |
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: |
85 |
– |
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: |
311 |
|
} |
312 |
|
|
313 |
|
void InteractionManager::setCutoffRadius(RealType rcut) { |
314 |
+ |
|
315 |
|
electrostatic_->setCutoffRadius(rcut); |
316 |
|
eam_->setCutoffRadius(rcut); |
317 |
|
} |
318 |
|
|
278 |
– |
void InteractionManager::setSwitchingRadius(RealType rswitch) { |
279 |
– |
electrostatic_->setSwitchingRadius(rswitch); |
280 |
– |
} |
281 |
– |
|
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) |
322 |
< |
(*it)->calcForce(idat); |
323 |
< |
|
324 |
< |
return; |
325 |
< |
} |
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 |
< |
|
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); |
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); |