ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/nonbonded/InteractionManager.cpp
Revision: 1925
Committed: Wed Aug 7 15:24:16 2013 UTC (11 years, 8 months ago) by gezelter
File size: 19263 byte(s)
Log Message:
More ewald fixes, reporting reciprocal potential in stats.

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 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

Properties

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