ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/nonbonded/InteractionManager.cpp
Revision: 2005
Committed: Wed Jun 11 18:37:56 2014 UTC (10 years, 10 months ago) by gezelter
File size: 19963 byte(s)
Log Message:
Fixed a strange non-bonded reporting bug.

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

Properties

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