ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/nonbonded/InteractionManager.cpp
Revision: 1929
Committed: Mon Aug 19 13:12:00 2013 UTC (11 years, 8 months ago) by gezelter
File size: 20022 byte(s)
Log Message:
Backing out fluc-rho and putting back the Electrostatic fluctuating
charge with coulomb integrals for atoms within a region.

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

Properties

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