1 |
/* |
2 |
* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. |
3 |
* |
4 |
* The University of Notre Dame grants you ("Licensee") a |
5 |
* non-exclusive, royalty free, license to use, modify and |
6 |
* redistribute this software in source and binary code form, provided |
7 |
* that the following conditions are met: |
8 |
* |
9 |
* 1. Redistributions of source code must retain the above copyright |
10 |
* notice, this list of conditions and the following disclaimer. |
11 |
* |
12 |
* 2. Redistributions in binary form must reproduce the above copyright |
13 |
* notice, this list of conditions and the following disclaimer in the |
14 |
* documentation and/or other materials provided with the |
15 |
* distribution. |
16 |
* |
17 |
* This software is provided "AS IS," without a warranty of any |
18 |
* kind. All express or implied conditions, representations and |
19 |
* warranties, including any implied warranty of merchantability, |
20 |
* fitness for a particular purpose or non-infringement, are hereby |
21 |
* excluded. The University of Notre Dame and its licensors shall not |
22 |
* be liable for any damages suffered by licensee as a result of |
23 |
* using, modifying or distributing the software or its |
24 |
* derivatives. In no event will the University of Notre Dame or its |
25 |
* licensors be liable for any lost revenue, profit or data, or for |
26 |
* direct, indirect, special, consequential, incidental or punitive |
27 |
* damages, however caused and regardless of the theory of liability, |
28 |
* arising out of the use of or inability to use software, even if the |
29 |
* University of Notre Dame has been advised of the possibility of |
30 |
* such damages. |
31 |
* |
32 |
* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your |
33 |
* research, please cite the appropriate papers when you publish your |
34 |
* work. Good starting points are: |
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, 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" |
44 |
|
45 |
namespace OpenMD { |
46 |
|
47 |
InteractionManager::InteractionManager() { |
48 |
|
49 |
initialized_ = false; |
50 |
|
51 |
lj_ = new LJ(); |
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 |
|
76 |
if (initialized_) return; |
77 |
|
78 |
ForceField* forceField_ = info_->getForceField(); |
79 |
|
80 |
lj_->setForceField(forceField_); |
81 |
gb_->setForceField(forceField_); |
82 |
sticky_->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 |
int nTypes = atomTypes->size(); |
93 |
iHash_.resize(nTypes); |
94 |
interactions_.resize(nTypes); |
95 |
ForceField::AtomTypeContainer::MapTypeIterator i1, i2; |
96 |
AtomType* atype1; |
97 |
AtomType* atype2; |
98 |
int atid1, atid2; |
99 |
|
100 |
// We only need to worry about the types that are actually in the simulation: |
101 |
|
102 |
set<AtomType*> atypes = info_->getSimulatedAtomTypes(); |
103 |
|
104 |
lj_->setSimulatedAtomTypes(atypes); |
105 |
gb_->setSimulatedAtomTypes(atypes); |
106 |
sticky_->setSimulatedAtomTypes(atypes); |
107 |
eam_->setSimulatedAtomTypes(atypes); |
108 |
sc_->setSimulatedAtomTypes(atypes); |
109 |
morse_->setSimulatedAtomTypes(atypes); |
110 |
electrostatic_->setSimInfo(info_); |
111 |
electrostatic_->setSimulatedAtomTypes(atypes); |
112 |
maw_->setSimulatedAtomTypes(atypes); |
113 |
repulsivePower_->setSimulatedAtomTypes(atypes); |
114 |
|
115 |
set<AtomType*>::iterator at; |
116 |
|
117 |
for (at = atypes.begin(); at != atypes.end(); ++at) { |
118 |
|
119 |
//for (atype1 = atomTypes->beginType(i1); atype1 != NULL; |
120 |
// atype1 = atomTypes->nextType(i1)) { |
121 |
|
122 |
atype1 = *at; |
123 |
atid1 = atype1->getIdent(); |
124 |
iHash_[atid1].resize(nTypes); |
125 |
interactions_[atid1].resize(nTypes); |
126 |
|
127 |
// add it to the map: |
128 |
pair<map<int,AtomType*>::iterator,bool> ret; |
129 |
ret = typeMap_.insert( pair<int, AtomType*>(atid1, atype1) ); |
130 |
if (ret.second == false) { |
131 |
sprintf( painCave.errMsg, |
132 |
"InteractionManager already had a previous entry with ident %d\n", |
133 |
atype1->getIdent()); |
134 |
painCave.severity = OPENMD_INFO; |
135 |
painCave.isFatal = 0; |
136 |
simError(); |
137 |
} |
138 |
} |
139 |
|
140 |
// Now, iterate over all known types and add to the interaction map: |
141 |
|
142 |
map<int, AtomType*>::iterator it1, it2; |
143 |
for (it1 = typeMap_.begin(); it1 != typeMap_.end(); ++it1) { |
144 |
atype1 = (*it1).second; |
145 |
atid1 = atype1->getIdent(); |
146 |
|
147 |
for( it2 = typeMap_.begin(); it2 != typeMap_.end(); ++it2) { |
148 |
atype2 = (*it2).second; |
149 |
atid2 = atype2->getIdent(); |
150 |
|
151 |
iHash_[atid1][atid2] = 0; |
152 |
|
153 |
if (atype1->isLennardJones() && atype2->isLennardJones()) { |
154 |
interactions_[atid1][atid2].insert(lj_); |
155 |
iHash_[atid1][atid2] |= LJ_PAIR; |
156 |
} |
157 |
if (atype1->isElectrostatic() && atype2->isElectrostatic() ) { |
158 |
interactions_[atid1][atid2].insert(electrostatic_); |
159 |
iHash_[atid1][atid2] |= ELECTROSTATIC_PAIR; |
160 |
} |
161 |
if (atype1->isSticky() && atype2->isSticky() ) { |
162 |
interactions_[atid1][atid2].insert(sticky_); |
163 |
iHash_[atid1][atid2] |= STICKY_PAIR; |
164 |
} |
165 |
if (atype1->isStickyPower() && atype2->isStickyPower() ) { |
166 |
interactions_[atid1][atid2].insert(sticky_); |
167 |
iHash_[atid1][atid2] |= STICKY_PAIR; |
168 |
} |
169 |
if (atype1->isEAM() && atype2->isEAM() ) { |
170 |
interactions_[atid1][atid2].insert(eam_); |
171 |
iHash_[atid1][atid2] |= EAM_PAIR; |
172 |
} |
173 |
if (atype1->isSC() && atype2->isSC() ) { |
174 |
interactions_[atid1][atid2].insert(sc_); |
175 |
iHash_[atid1][atid2] |= SC_PAIR; |
176 |
} |
177 |
if (atype1->isGayBerne() && atype2->isGayBerne() ) { |
178 |
interactions_[atid1][atid2].insert(gb_); |
179 |
iHash_[atid1][atid2] |= GB_PAIR; |
180 |
} |
181 |
if ((atype1->isGayBerne() && atype2->isLennardJones()) |
182 |
|| (atype1->isLennardJones() && atype2->isGayBerne())) { |
183 |
interactions_[atid1][atid2].insert(gb_); |
184 |
iHash_[atid1][atid2] |= GB_PAIR; |
185 |
} |
186 |
|
187 |
// look for an explicitly-set non-bonded interaction type using the |
188 |
// two atom types. |
189 |
NonBondedInteractionType* nbiType = forceField_->getNonBondedInteractionType(atype1->getName(), atype2->getName()); |
190 |
|
191 |
if (nbiType != NULL) { |
192 |
|
193 |
bool vdwExplicit = false; |
194 |
bool metExplicit = false; |
195 |
// bool hbExplicit = false; |
196 |
|
197 |
if (nbiType->isLennardJones()) { |
198 |
// We found an explicit Lennard-Jones interaction. |
199 |
// override all other vdw entries for this pair of atom types: |
200 |
set<NonBondedInteraction*>::iterator it; |
201 |
for (it = interactions_[atid1][atid2].begin(); |
202 |
it != interactions_[atid1][atid2].end(); ++it) { |
203 |
InteractionFamily ifam = (*it)->getFamily(); |
204 |
if (ifam == VANDERWAALS_FAMILY) { |
205 |
interactions_[atid1][atid2].erase(*it); |
206 |
iHash_[atid1][atid2] ^= (*it)->getHash(); |
207 |
} |
208 |
} |
209 |
interactions_[atid1][atid2].insert(lj_); |
210 |
iHash_[atid1][atid2] |= LJ_PAIR; |
211 |
vdwExplicit = true; |
212 |
} |
213 |
|
214 |
if (nbiType->isMorse()) { |
215 |
if (vdwExplicit) { |
216 |
sprintf( painCave.errMsg, |
217 |
"InteractionManager::initialize found more than one " |
218 |
"explicit \n" |
219 |
"\tvan der Waals interaction for atom types %s - %s\n", |
220 |
atype1->getName().c_str(), atype2->getName().c_str()); |
221 |
painCave.severity = OPENMD_ERROR; |
222 |
painCave.isFatal = 1; |
223 |
simError(); |
224 |
} |
225 |
// We found an explicit Morse interaction. |
226 |
// override all other vdw entries for this pair of atom types: |
227 |
set<NonBondedInteraction*>::iterator it; |
228 |
for (it = interactions_[atid1][atid2].begin(); |
229 |
it != interactions_[atid1][atid2].end(); ++it) { |
230 |
InteractionFamily ifam = (*it)->getFamily(); |
231 |
if (ifam == VANDERWAALS_FAMILY) { |
232 |
interactions_[atid1][atid2].erase(*it); |
233 |
iHash_[atid1][atid2] ^= (*it)->getHash(); |
234 |
} |
235 |
} |
236 |
interactions_[atid1][atid2].insert(morse_); |
237 |
iHash_[atid1][atid2] |= MORSE_PAIR; |
238 |
vdwExplicit = true; |
239 |
} |
240 |
|
241 |
if (nbiType->isRepulsivePower()) { |
242 |
if (vdwExplicit) { |
243 |
sprintf( painCave.errMsg, |
244 |
"InteractionManager::initialize found more than one " |
245 |
"explicit \n" |
246 |
"\tvan der Waals interaction for atom types %s - %s\n", |
247 |
atype1->getName().c_str(), atype2->getName().c_str()); |
248 |
painCave.severity = OPENMD_ERROR; |
249 |
painCave.isFatal = 1; |
250 |
simError(); |
251 |
} |
252 |
// We found an explicit RepulsivePower interaction. |
253 |
// override all other vdw entries for this pair of atom types: |
254 |
set<NonBondedInteraction*>::iterator it; |
255 |
for (it = interactions_[atid1][atid2].begin(); |
256 |
it != interactions_[atid1][atid2].end(); ++it) { |
257 |
InteractionFamily ifam = (*it)->getFamily(); |
258 |
if (ifam == VANDERWAALS_FAMILY) { |
259 |
interactions_[atid1][atid2].erase(*it); |
260 |
iHash_[atid1][atid2] ^= (*it)->getHash(); |
261 |
} |
262 |
} |
263 |
interactions_[atid1][atid2].insert(repulsivePower_); |
264 |
iHash_[atid1][atid2] |= REPULSIVEPOWER_PAIR; |
265 |
vdwExplicit = true; |
266 |
} |
267 |
|
268 |
|
269 |
if (nbiType->isEAM()) { |
270 |
// We found an explicit EAM interaction. |
271 |
// override all other metallic entries for this pair of atom types: |
272 |
set<NonBondedInteraction*>::iterator it; |
273 |
for (it = interactions_[atid1][atid2].begin(); |
274 |
it != interactions_[atid1][atid2].end(); ++it) { |
275 |
InteractionFamily ifam = (*it)->getFamily(); |
276 |
if (ifam == METALLIC_FAMILY) { |
277 |
interactions_[atid1][atid2].erase(*it); |
278 |
iHash_[atid1][atid2] ^= (*it)->getHash(); |
279 |
} |
280 |
} |
281 |
interactions_[atid1][atid2].insert(eam_); |
282 |
iHash_[atid1][atid2] |= EAM_PAIR; |
283 |
metExplicit = true; |
284 |
} |
285 |
|
286 |
if (nbiType->isSC()) { |
287 |
if (metExplicit) { |
288 |
sprintf( painCave.errMsg, |
289 |
"InteractionManager::initialize found more than one " |
290 |
"explicit\n" |
291 |
"\tmetallic interaction for atom types %s - %s\n", |
292 |
atype1->getName().c_str(), atype2->getName().c_str()); |
293 |
painCave.severity = OPENMD_ERROR; |
294 |
painCave.isFatal = 1; |
295 |
simError(); |
296 |
} |
297 |
// We found an explicit Sutton-Chen interaction. |
298 |
// override all other metallic entries for this pair of atom types: |
299 |
set<NonBondedInteraction*>::iterator it; |
300 |
for (it = interactions_[atid1][atid2].begin(); |
301 |
it != interactions_[atid1][atid2].end(); ++it) { |
302 |
InteractionFamily ifam = (*it)->getFamily(); |
303 |
if (ifam == METALLIC_FAMILY) { |
304 |
interactions_[atid1][atid2].erase(*it); |
305 |
iHash_[atid1][atid2] ^= (*it)->getHash(); |
306 |
} |
307 |
} |
308 |
interactions_[atid1][atid2].insert(sc_); |
309 |
iHash_[atid1][atid2] |= SC_PAIR; |
310 |
metExplicit = true; |
311 |
} |
312 |
|
313 |
if (nbiType->isMAW()) { |
314 |
if (vdwExplicit) { |
315 |
sprintf( painCave.errMsg, |
316 |
"InteractionManager::initialize found more than one " |
317 |
"explicit\n" |
318 |
"\tvan der Waals interaction for atom types %s - %s\n", |
319 |
atype1->getName().c_str(), atype2->getName().c_str()); |
320 |
painCave.severity = OPENMD_ERROR; |
321 |
painCave.isFatal = 1; |
322 |
simError(); |
323 |
} |
324 |
// We found an explicit MAW interaction. |
325 |
// override all other vdw entries for this pair of atom types: |
326 |
set<NonBondedInteraction*>::iterator it; |
327 |
for (it = interactions_[atid1][atid2].begin(); |
328 |
it != interactions_[atid1][atid2].end(); ++it) { |
329 |
InteractionFamily ifam = (*it)->getFamily(); |
330 |
if (ifam == VANDERWAALS_FAMILY) { |
331 |
interactions_[atid1][atid2].erase(*it); |
332 |
iHash_[atid1][atid2] ^= (*it)->getHash(); |
333 |
} |
334 |
} |
335 |
interactions_[atid1][atid2].insert(maw_); |
336 |
iHash_[atid1][atid2] |= MAW_PAIR; |
337 |
vdwExplicit = true; |
338 |
} |
339 |
} |
340 |
} |
341 |
} |
342 |
|
343 |
|
344 |
// Make sure every pair of atom types in this simulation has a |
345 |
// non-bonded interaction. If not, just inform the user. |
346 |
|
347 |
set<AtomType*> simTypes = info_->getSimulatedAtomTypes(); |
348 |
set<AtomType*>::iterator it, jt; |
349 |
|
350 |
for (it = simTypes.begin(); it != simTypes.end(); ++it) { |
351 |
atype1 = (*it); |
352 |
atid1 = atype1->getIdent(); |
353 |
for (jt = it; jt != simTypes.end(); ++jt) { |
354 |
atype2 = (*jt); |
355 |
atid1 = atype1->getIdent(); |
356 |
|
357 |
if (interactions_[atid1][atid2].size() == 0) { |
358 |
sprintf( painCave.errMsg, |
359 |
"InteractionManager could not find a matching non-bonded\n" |
360 |
"\tinteraction for atom types %s - %s\n" |
361 |
"\tProceeding without this interaction.\n", |
362 |
atype1->getName().c_str(), atype2->getName().c_str()); |
363 |
painCave.severity = OPENMD_INFO; |
364 |
painCave.isFatal = 0; |
365 |
simError(); |
366 |
} |
367 |
} |
368 |
} |
369 |
|
370 |
initialized_ = true; |
371 |
} |
372 |
|
373 |
void InteractionManager::setCutoffRadius(RealType rcut) { |
374 |
|
375 |
electrostatic_->setCutoffRadius(rcut); |
376 |
eam_->setCutoffRadius(rcut); |
377 |
} |
378 |
|
379 |
void InteractionManager::doPrePair(InteractionData &idat){ |
380 |
|
381 |
if (!initialized_) initialize(); |
382 |
|
383 |
// excluded interaction, so just return |
384 |
if (idat.excluded) return; |
385 |
|
386 |
int& iHash = iHash_[idat.atid1][idat.atid2]; |
387 |
|
388 |
if ((iHash & EAM_PAIR) != 0) eam_->calcDensity(idat); |
389 |
if ((iHash & SC_PAIR) != 0) sc_->calcDensity(idat); |
390 |
|
391 |
// set<NonBondedInteraction*>::iterator it; |
392 |
|
393 |
// for (it = interactions_[ idat.atypes ].begin(); |
394 |
// it != interactions_[ idat.atypes ].end(); ++it){ |
395 |
// if ((*it)->getFamily() == METALLIC_FAMILY) { |
396 |
// dynamic_cast<MetallicInteraction*>(*it)->calcDensity(idat); |
397 |
// } |
398 |
// } |
399 |
|
400 |
return; |
401 |
} |
402 |
|
403 |
void InteractionManager::doPreForce(SelfData &sdat){ |
404 |
|
405 |
if (!initialized_) initialize(); |
406 |
|
407 |
int& iHash = iHash_[sdat.atid][sdat.atid]; |
408 |
|
409 |
if ((iHash & EAM_PAIR) != 0) eam_->calcFunctional(sdat); |
410 |
if ((iHash & SC_PAIR) != 0) sc_->calcFunctional(sdat); |
411 |
|
412 |
// set<NonBondedInteraction*>::iterator it; |
413 |
|
414 |
// for (it = interactions_[atid1][atid2].begin(); it != interactions_[atid1][atid2].end(); ++it){ |
415 |
// if ((*it)->getFamily() == METALLIC_FAMILY) { |
416 |
// dynamic_cast<MetallicInteraction*>(*it)->calcFunctional(sdat); |
417 |
// } |
418 |
// } |
419 |
|
420 |
return; |
421 |
} |
422 |
|
423 |
void InteractionManager::doPair(InteractionData &idat){ |
424 |
|
425 |
if (!initialized_) initialize(); |
426 |
|
427 |
int& iHash = iHash_[idat.atid1][idat.atid2]; |
428 |
|
429 |
if ((iHash & ELECTROSTATIC_PAIR) != 0) electrostatic_->calcForce(idat); |
430 |
|
431 |
// electrostatics still has to worry about indirect |
432 |
// contributions from excluded pairs of atoms, but nothing else does: |
433 |
|
434 |
if (idat.excluded) return; |
435 |
|
436 |
if ((iHash & LJ_PAIR) != 0) lj_->calcForce(idat); |
437 |
if ((iHash & GB_PAIR) != 0) gb_->calcForce(idat); |
438 |
if ((iHash & STICKY_PAIR) != 0) sticky_->calcForce(idat); |
439 |
if ((iHash & MORSE_PAIR) != 0) morse_->calcForce(idat); |
440 |
if ((iHash & REPULSIVEPOWER_PAIR) != 0) repulsivePower_->calcForce(idat); |
441 |
if ((iHash & EAM_PAIR) != 0) eam_->calcForce(idat); |
442 |
if ((iHash & SC_PAIR) != 0) sc_->calcForce(idat); |
443 |
if ((iHash & MAW_PAIR) != 0) maw_->calcForce(idat); |
444 |
|
445 |
// set<NonBondedInteraction*>::iterator it; |
446 |
|
447 |
// for (it = interactions_[ idat.atypes ].begin(); |
448 |
// it != interactions_[ idat.atypes ].end(); ++it) { |
449 |
|
450 |
// // electrostatics still has to worry about indirect |
451 |
// // contributions from excluded pairs of atoms: |
452 |
|
453 |
// if (!idat.excluded || (*it)->getFamily() == ELECTROSTATIC_FAMILY) { |
454 |
// (*it)->calcForce(idat); |
455 |
// } |
456 |
// } |
457 |
|
458 |
return; |
459 |
} |
460 |
|
461 |
void InteractionManager::doSelfCorrection(SelfData &sdat){ |
462 |
|
463 |
if (!initialized_) initialize(); |
464 |
|
465 |
int& iHash = iHash_[sdat.atid][sdat.atid]; |
466 |
|
467 |
if ((iHash & ELECTROSTATIC_PAIR) != 0) electrostatic_->calcSelfCorrection(sdat); |
468 |
|
469 |
// set<NonBondedInteraction*>::iterator it; |
470 |
|
471 |
// for (it = interactions_[atid1][atid2].begin(); it != interactions_[atid1][atid2].end(); ++it){ |
472 |
// if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) { |
473 |
// dynamic_cast<ElectrostaticInteraction*>(*it)->calcSelfCorrection(sdat); |
474 |
// } |
475 |
// } |
476 |
|
477 |
return; |
478 |
} |
479 |
|
480 |
void InteractionManager::doReciprocalSpaceSum(RealType &pot){ |
481 |
if (!initialized_) initialize(); |
482 |
electrostatic_->ReciprocalSpaceSum(pot); |
483 |
} |
484 |
|
485 |
RealType InteractionManager::getSuggestedCutoffRadius(int *atid) { |
486 |
if (!initialized_) initialize(); |
487 |
|
488 |
AtomType* atype = typeMap_[*atid]; |
489 |
|
490 |
set<NonBondedInteraction*>::iterator it; |
491 |
RealType cutoff = 0.0; |
492 |
|
493 |
for (it = interactions_[*atid][*atid].begin(); |
494 |
it != interactions_[*atid][*atid].end(); |
495 |
++it) |
496 |
cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(make_pair(atype, atype))); |
497 |
return cutoff; |
498 |
} |
499 |
|
500 |
RealType InteractionManager::getSuggestedCutoffRadius(AtomType* atype) { |
501 |
if (!initialized_) initialize(); |
502 |
|
503 |
int atid = atype->getIdent(); |
504 |
|
505 |
set<NonBondedInteraction*>::iterator it; |
506 |
RealType cutoff = 0.0; |
507 |
|
508 |
for (it = interactions_[atid][atid].begin(); |
509 |
it != interactions_[atid][atid].end(); ++it) |
510 |
cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(make_pair(atype, atype))); |
511 |
return cutoff; |
512 |
} |
513 |
} //end namespace OpenMD |