ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/nonbonded/InteractionManager.cpp
Revision: 1927
Committed: Wed Aug 14 20:19:19 2013 UTC (11 years, 8 months ago) by gezelter
File size: 20326 byte(s)
Log Message:
FlucRho/EAM initial commit

File Contents

# Content
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 sHash_.resize(nTypes);
94 iHash_.resize(nTypes);
95 interactions_.resize(nTypes);
96 ForceField::AtomTypeContainer::MapTypeIterator i1, i2;
97 AtomType* atype1;
98 AtomType* atype2;
99 int atid1, atid2;
100
101 // We only need to worry about the types that are actually in the simulation:
102
103 set<AtomType*> atypes = info_->getSimulatedAtomTypes();
104
105 lj_->setSimulatedAtomTypes(atypes);
106 gb_->setSimulatedAtomTypes(atypes);
107 sticky_->setSimulatedAtomTypes(atypes);
108 eam_->setSimulatedAtomTypes(atypes);
109 sc_->setSimulatedAtomTypes(atypes);
110 morse_->setSimulatedAtomTypes(atypes);
111 electrostatic_->setSimInfo(info_);
112 electrostatic_->setSimulatedAtomTypes(atypes);
113 maw_->setSimulatedAtomTypes(atypes);
114 repulsivePower_->setSimulatedAtomTypes(atypes);
115
116 set<AtomType*>::iterator at;
117
118 for (at = atypes.begin(); at != atypes.end(); ++at) {
119
120 //for (atype1 = atomTypes->beginType(i1); atype1 != NULL;
121 // atype1 = atomTypes->nextType(i1)) {
122
123 atype1 = *at;
124 atid1 = atype1->getIdent();
125 iHash_[atid1].resize(nTypes);
126 interactions_[atid1].resize(nTypes);
127
128 // add it to the map:
129 pair<map<int,AtomType*>::iterator,bool> ret;
130 ret = typeMap_.insert( pair<int, AtomType*>(atid1, atype1) );
131 if (ret.second == false) {
132 sprintf( painCave.errMsg,
133 "InteractionManager already had a previous entry with ident %d\n",
134 atype1->getIdent());
135 painCave.severity = OPENMD_INFO;
136 painCave.isFatal = 0;
137 simError();
138 }
139 }
140
141 if (atype1->isLennardJones()) {
142 sHash_[atid1] |= LJ_INTERACTION;
143 }
144 if (atype1->isElectrostatic()) {
145 sHash_[atid1] |= ELECTROSTATIC_INTERACTION;
146 }
147 if (atype1->isSticky()) {
148 sHash_[atid1] |= STICKY_INTERACTION;
149 }
150 if (atype1->isStickyPower()) {
151 sHash_[atid1] |= STICKY_INTERACTION;
152 }
153 if (atype1->isEAM()) {
154 sHash_[atid1] |= EAM_INTERACTION;
155 }
156 if (atype1->isSC()) {
157 sHash_[atid1] |= SC_INTERACTION;
158 }
159 if (atype1->isGayBerne()) {
160 sHash_[atid1] |= GB_INTERACTION;
161 }
162
163 // Now, iterate over all known types and add to the interaction map:
164
165 map<int, AtomType*>::iterator it1, it2;
166 for (it1 = typeMap_.begin(); it1 != typeMap_.end(); ++it1) {
167 atype1 = (*it1).second;
168 atid1 = atype1->getIdent();
169
170 for( it2 = typeMap_.begin(); it2 != typeMap_.end(); ++it2) {
171 atype2 = (*it2).second;
172 atid2 = atype2->getIdent();
173
174 iHash_[atid1][atid2] = 0;
175
176 if (atype1->isLennardJones() && atype2->isLennardJones()) {
177 interactions_[atid1][atid2].insert(lj_);
178 iHash_[atid1][atid2] |= LJ_INTERACTION;
179 }
180 if (atype1->isElectrostatic() && atype2->isElectrostatic() ) {
181 // Pairs of fluctuating density EAM atoms have their
182 // interactions handled via the EAM routines. All other
183 // interactions with these atoms are handled via normal
184 // electrostatic channels:
185 if (!(atype1->isEAM() && atype2->isEAM())) {
186 interactions_[atid1][atid2].insert(electrostatic_);
187 iHash_[atid1][atid2] |= ELECTROSTATIC_INTERACTION;
188 }
189 }
190 if (atype1->isSticky() && atype2->isSticky() ) {
191 interactions_[atid1][atid2].insert(sticky_);
192 iHash_[atid1][atid2] |= STICKY_INTERACTION;
193 }
194 if (atype1->isStickyPower() && atype2->isStickyPower() ) {
195 interactions_[atid1][atid2].insert(sticky_);
196 iHash_[atid1][atid2] |= STICKY_INTERACTION;
197 }
198 if (atype1->isEAM() && atype2->isEAM() ) {
199 interactions_[atid1][atid2].insert(eam_);
200 iHash_[atid1][atid2] |= EAM_INTERACTION;
201 }
202 if (atype1->isSC() && atype2->isSC() ) {
203 interactions_[atid1][atid2].insert(sc_);
204 iHash_[atid1][atid2] |= SC_INTERACTION;
205 }
206 if (atype1->isGayBerne() && atype2->isGayBerne() ) {
207 interactions_[atid1][atid2].insert(gb_);
208 iHash_[atid1][atid2] |= GB_INTERACTION;
209 }
210 if ((atype1->isGayBerne() && atype2->isLennardJones())
211 || (atype1->isLennardJones() && atype2->isGayBerne())) {
212 interactions_[atid1][atid2].insert(gb_);
213 iHash_[atid1][atid2] |= GB_INTERACTION;
214 }
215
216 // look for an explicitly-set non-bonded interaction type using the
217 // two atom types.
218 NonBondedInteractionType* nbiType = forceField_->getNonBondedInteractionType(atype1->getName(), atype2->getName());
219
220 if (nbiType != NULL) {
221
222 bool vdwExplicit = false;
223 bool metExplicit = false;
224 // bool hbExplicit = false;
225
226 if (nbiType->isLennardJones()) {
227 // We found an explicit Lennard-Jones interaction.
228 // override all other vdw entries for this pair of atom types:
229 set<NonBondedInteraction*>::iterator it;
230 for (it = interactions_[atid1][atid2].begin();
231 it != interactions_[atid1][atid2].end(); ++it) {
232 InteractionFamily ifam = (*it)->getFamily();
233 if (ifam == VANDERWAALS_FAMILY) {
234 interactions_[atid1][atid2].erase(*it);
235 iHash_[atid1][atid2] ^= (*it)->getHash();
236 }
237 }
238 interactions_[atid1][atid2].insert(lj_);
239 iHash_[atid1][atid2] |= LJ_INTERACTION;
240 vdwExplicit = true;
241 }
242
243 if (nbiType->isMorse()) {
244 if (vdwExplicit) {
245 sprintf( painCave.errMsg,
246 "InteractionManager::initialize found more than one "
247 "explicit \n"
248 "\tvan der Waals interaction for atom types %s - %s\n",
249 atype1->getName().c_str(), atype2->getName().c_str());
250 painCave.severity = OPENMD_ERROR;
251 painCave.isFatal = 1;
252 simError();
253 }
254 // We found an explicit Morse interaction.
255 // override all other vdw entries for this pair of atom types:
256 set<NonBondedInteraction*>::iterator it;
257 for (it = interactions_[atid1][atid2].begin();
258 it != interactions_[atid1][atid2].end(); ++it) {
259 InteractionFamily ifam = (*it)->getFamily();
260 if (ifam == VANDERWAALS_FAMILY) {
261 interactions_[atid1][atid2].erase(*it);
262 iHash_[atid1][atid2] ^= (*it)->getHash();
263 }
264 }
265 interactions_[atid1][atid2].insert(morse_);
266 iHash_[atid1][atid2] |= MORSE_INTERACTION;
267 vdwExplicit = true;
268 }
269
270 if (nbiType->isRepulsivePower()) {
271 if (vdwExplicit) {
272 sprintf( painCave.errMsg,
273 "InteractionManager::initialize found more than one "
274 "explicit \n"
275 "\tvan der Waals interaction for atom types %s - %s\n",
276 atype1->getName().c_str(), atype2->getName().c_str());
277 painCave.severity = OPENMD_ERROR;
278 painCave.isFatal = 1;
279 simError();
280 }
281 // We found an explicit RepulsivePower interaction.
282 // override all other vdw entries for this pair of atom types:
283 set<NonBondedInteraction*>::iterator it;
284 for (it = interactions_[atid1][atid2].begin();
285 it != interactions_[atid1][atid2].end(); ++it) {
286 InteractionFamily ifam = (*it)->getFamily();
287 if (ifam == VANDERWAALS_FAMILY) {
288 interactions_[atid1][atid2].erase(*it);
289 iHash_[atid1][atid2] ^= (*it)->getHash();
290 }
291 }
292 interactions_[atid1][atid2].insert(repulsivePower_);
293 iHash_[atid1][atid2] |= REPULSIVEPOWER_INTERACTION;
294 vdwExplicit = true;
295 }
296
297
298 if (nbiType->isEAM()) {
299 // We found an explicit EAM interaction.
300 // override all other metallic entries for this pair of atom types:
301 set<NonBondedInteraction*>::iterator it;
302 for (it = interactions_[atid1][atid2].begin();
303 it != interactions_[atid1][atid2].end(); ++it) {
304 InteractionFamily ifam = (*it)->getFamily();
305 if (ifam == METALLIC_FAMILY) {
306 interactions_[atid1][atid2].erase(*it);
307 iHash_[atid1][atid2] ^= (*it)->getHash();
308 }
309 }
310 interactions_[atid1][atid2].insert(eam_);
311 iHash_[atid1][atid2] |= EAM_INTERACTION;
312 metExplicit = true;
313 }
314
315 if (nbiType->isSC()) {
316 if (metExplicit) {
317 sprintf( painCave.errMsg,
318 "InteractionManager::initialize found more than one "
319 "explicit\n"
320 "\tmetallic interaction for atom types %s - %s\n",
321 atype1->getName().c_str(), atype2->getName().c_str());
322 painCave.severity = OPENMD_ERROR;
323 painCave.isFatal = 1;
324 simError();
325 }
326 // We found an explicit Sutton-Chen interaction.
327 // override all other metallic entries for this pair of atom types:
328 set<NonBondedInteraction*>::iterator it;
329 for (it = interactions_[atid1][atid2].begin();
330 it != interactions_[atid1][atid2].end(); ++it) {
331 InteractionFamily ifam = (*it)->getFamily();
332 if (ifam == METALLIC_FAMILY) {
333 interactions_[atid1][atid2].erase(*it);
334 iHash_[atid1][atid2] ^= (*it)->getHash();
335 }
336 }
337 interactions_[atid1][atid2].insert(sc_);
338 iHash_[atid1][atid2] |= SC_INTERACTION;
339 metExplicit = true;
340 }
341
342 if (nbiType->isMAW()) {
343 if (vdwExplicit) {
344 sprintf( painCave.errMsg,
345 "InteractionManager::initialize found more than one "
346 "explicit\n"
347 "\tvan der Waals interaction for atom types %s - %s\n",
348 atype1->getName().c_str(), atype2->getName().c_str());
349 painCave.severity = OPENMD_ERROR;
350 painCave.isFatal = 1;
351 simError();
352 }
353 // We found an explicit MAW interaction.
354 // override all other vdw entries for this pair of atom types:
355 set<NonBondedInteraction*>::iterator it;
356 for (it = interactions_[atid1][atid2].begin();
357 it != interactions_[atid1][atid2].end(); ++it) {
358 InteractionFamily ifam = (*it)->getFamily();
359 if (ifam == VANDERWAALS_FAMILY) {
360 interactions_[atid1][atid2].erase(*it);
361 iHash_[atid1][atid2] ^= (*it)->getHash();
362 }
363 }
364 interactions_[atid1][atid2].insert(maw_);
365 iHash_[atid1][atid2] |= MAW_INTERACTION;
366 vdwExplicit = true;
367 }
368 }
369 }
370 }
371
372
373 // Make sure every pair of atom types in this simulation has a
374 // non-bonded interaction. If not, just inform the user.
375
376 set<AtomType*> simTypes = info_->getSimulatedAtomTypes();
377 set<AtomType*>::iterator it, jt;
378
379 for (it = simTypes.begin(); it != simTypes.end(); ++it) {
380 atype1 = (*it);
381 atid1 = atype1->getIdent();
382 for (jt = it; jt != simTypes.end(); ++jt) {
383 atype2 = (*jt);
384 atid1 = atype1->getIdent();
385
386 if (interactions_[atid1][atid2].size() == 0) {
387 sprintf( painCave.errMsg,
388 "InteractionManager could not find a matching non-bonded\n"
389 "\tinteraction for atom types %s - %s\n"
390 "\tProceeding without this interaction.\n",
391 atype1->getName().c_str(), atype2->getName().c_str());
392 painCave.severity = OPENMD_INFO;
393 painCave.isFatal = 0;
394 simError();
395 }
396 }
397 }
398
399 initialized_ = true;
400 }
401
402 void InteractionManager::setCutoffRadius(RealType rcut) {
403
404 electrostatic_->setCutoffRadius(rcut);
405 eam_->setCutoffRadius(rcut);
406 }
407
408 void InteractionManager::doPrePair(InteractionData &idat){
409
410 if (!initialized_) initialize();
411
412 // excluded interaction, so just return
413 if (idat.excluded) return;
414
415 int& iHash = iHash_[idat.atid1][idat.atid2];
416
417 if ((iHash & EAM_INTERACTION) != 0) eam_->calcDensity(idat);
418 if ((iHash & SC_INTERACTION) != 0) sc_->calcDensity(idat);
419
420 // set<NonBondedInteraction*>::iterator it;
421
422 // for (it = interactions_[ idat.atypes ].begin();
423 // it != interactions_[ idat.atypes ].end(); ++it){
424 // if ((*it)->getFamily() == METALLIC_FAMILY) {
425 // dynamic_cast<MetallicInteraction*>(*it)->calcDensity(idat);
426 // }
427 // }
428
429 return;
430 }
431
432 void InteractionManager::doPreForce(SelfData &sdat){
433
434 if (!initialized_) initialize();
435
436 int& sHash = sHash_[sdat.atid];
437
438 if ((sHash & EAM_INTERACTION) != 0) eam_->calcFunctional(sdat);
439 if ((sHash & SC_INTERACTION) != 0) sc_->calcFunctional(sdat);
440
441 // set<NonBondedInteraction*>::iterator it;
442
443 // for (it = interactions_[atid1][atid2].begin(); it != interactions_[atid1][atid2].end(); ++it){
444 // if ((*it)->getFamily() == METALLIC_FAMILY) {
445 // dynamic_cast<MetallicInteraction*>(*it)->calcFunctional(sdat);
446 // }
447 // }
448
449 return;
450 }
451
452 void InteractionManager::doPair(InteractionData &idat){
453
454 if (!initialized_) initialize();
455
456 int& iHash = iHash_[idat.atid1][idat.atid2];
457
458 if ((iHash & ELECTROSTATIC_INTERACTION) != 0) electrostatic_->calcForce(idat);
459
460 // electrostatics still has to worry about indirect
461 // contributions from excluded pairs of atoms, but nothing else does:
462
463 if (idat.excluded) return;
464
465 if ((iHash & LJ_INTERACTION) != 0) lj_->calcForce(idat);
466 if ((iHash & GB_INTERACTION) != 0) gb_->calcForce(idat);
467 if ((iHash & STICKY_INTERACTION) != 0) sticky_->calcForce(idat);
468 if ((iHash & MORSE_INTERACTION) != 0) morse_->calcForce(idat);
469 if ((iHash & REPULSIVEPOWER_INTERACTION) != 0) repulsivePower_->calcForce(idat);
470 if ((iHash & EAM_INTERACTION) != 0) eam_->calcForce(idat);
471 if ((iHash & SC_INTERACTION) != 0) sc_->calcForce(idat);
472 if ((iHash & MAW_INTERACTION) != 0) maw_->calcForce(idat);
473
474 // set<NonBondedInteraction*>::iterator it;
475
476 // for (it = interactions_[ idat.atypes ].begin();
477 // it != interactions_[ idat.atypes ].end(); ++it) {
478
479 // // electrostatics still has to worry about indirect
480 // // contributions from excluded pairs of atoms:
481
482 // if (!idat.excluded || (*it)->getFamily() == ELECTROSTATIC_FAMILY) {
483 // (*it)->calcForce(idat);
484 // }
485 // }
486
487 return;
488 }
489
490 void InteractionManager::doSelfCorrection(SelfData &sdat){
491
492 if (!initialized_) initialize();
493
494 int& sHash = sHash_[sdat.atid];
495
496 if ((sHash & ELECTROSTATIC_INTERACTION) != 0) electrostatic_->calcSelfCorrection(sdat);
497
498 // set<NonBondedInteraction*>::iterator it;
499
500 // for (it = interactions_[atid1][atid2].begin(); it != interactions_[atid1][atid2].end(); ++it){
501 // if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
502 // dynamic_cast<ElectrostaticInteraction*>(*it)->calcSelfCorrection(sdat);
503 // }
504 // }
505
506 return;
507 }
508
509 void InteractionManager::doReciprocalSpaceSum(RealType &pot){
510 if (!initialized_) initialize();
511 electrostatic_->ReciprocalSpaceSum(pot);
512 }
513
514 RealType InteractionManager::getSuggestedCutoffRadius(int *atid) {
515 if (!initialized_) initialize();
516
517 AtomType* atype = typeMap_[*atid];
518
519 set<NonBondedInteraction*>::iterator it;
520 RealType cutoff = 0.0;
521
522 for (it = interactions_[*atid][*atid].begin();
523 it != interactions_[*atid][*atid].end();
524 ++it)
525 cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(make_pair(atype, atype)));
526 return cutoff;
527 }
528
529 RealType InteractionManager::getSuggestedCutoffRadius(AtomType* atype) {
530 if (!initialized_) initialize();
531
532 int atid = atype->getIdent();
533
534 set<NonBondedInteraction*>::iterator it;
535 RealType cutoff = 0.0;
536
537 for (it = interactions_[atid][atid].begin();
538 it != interactions_[atid][atid].end(); ++it)
539 cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(make_pair(atype, atype)));
540 return cutoff;
541 }
542 } //end namespace OpenMD

Properties

Name Value
svn:eol-style native
svn:executable *