ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/nonbonded/InteractionManager.cpp
Revision: 1583
Committed: Thu Jun 16 22:00:08 2011 UTC (13 years, 10 months ago) by gezelter
File size: 13914 byte(s)
Log Message:
Bug squashing

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, 24107 (2008).
39 * [4] Vardeman & Gezelter, in progress (2009).
40 */
41
42 #include "nonbonded/InteractionManager.hpp"
43
44 namespace OpenMD {
45
46 InteractionManager::InteractionManager() {
47
48 initialized_ = false;
49
50 lj_ = new LJ();
51 gb_ = new GB();
52 sticky_ = new Sticky();
53 morse_ = new Morse();
54 eam_ = new EAM();
55 sc_ = new SC();
56 electrostatic_ = new Electrostatic();
57 maw_ = new MAW();
58 }
59
60 void InteractionManager::initialize() {
61
62 ForceField* forceField_ = info_->getForceField();
63
64 lj_->setForceField(forceField_);
65 gb_->setForceField(forceField_);
66 sticky_->setForceField(forceField_);
67 eam_->setForceField(forceField_);
68 sc_->setForceField(forceField_);
69 morse_->setForceField(forceField_);
70 electrostatic_->setForceField(forceField_);
71 maw_->setForceField(forceField_);
72
73 ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
74 ForceField::AtomTypeContainer::MapTypeIterator i1, i2;
75 AtomType* atype1;
76 AtomType* atype2;
77 pair<AtomType*, AtomType*> key;
78 pair<set<NonBondedInteraction*>::iterator, bool> ret;
79
80 for (atype1 = atomTypes->beginType(i1); atype1 != NULL;
81 atype1 = atomTypes->nextType(i1)) {
82
83 // add it to the map:
84 AtomTypeProperties atp = atype1->getATP();
85
86 pair<map<int,AtomType*>::iterator,bool> ret;
87 ret = typeMap_.insert( pair<int, AtomType*>(atp.ident, atype1) );
88 if (ret.second == false) {
89 sprintf( painCave.errMsg,
90 "InteractionManager already had a previous entry with ident %d\n",
91 atp.ident);
92 painCave.severity = OPENMD_INFO;
93 painCave.isFatal = 0;
94 simError();
95 }
96 }
97
98 // Now, iterate over all known types and add to the interaction map:
99
100 map<int, AtomType*>::iterator it1, it2;
101 for (it1 = typeMap_.begin(); it1 != typeMap_.end(); ++it1) {
102 atype1 = (*it1).second;
103
104 for( it2 = typeMap_.begin(); it2 != typeMap_.end(); ++it2) {
105 atype2 = (*it2).second;
106
107 bool vdwExplicit = false;
108 bool metExplicit = false;
109 bool hbExplicit = false;
110
111 key = make_pair(atype1, atype2);
112
113 if (atype1->isLennardJones() && atype2->isLennardJones()) {
114 interactions_[key].insert(lj_);
115 }
116 if (atype1->isElectrostatic() && atype2->isElectrostatic() ) {
117 interactions_[key].insert(electrostatic_);
118 }
119 if (atype1->isSticky() && atype2->isSticky() ) {
120 interactions_[key].insert(sticky_);
121 }
122 if (atype1->isStickyPower() && atype2->isStickyPower() ) {
123 interactions_[key].insert(sticky_);
124 }
125 if (atype1->isEAM() && atype2->isEAM() ) {
126 interactions_[key].insert(eam_);
127 }
128 if (atype1->isSC() && atype2->isSC() ) {
129 interactions_[key].insert(sc_);
130 }
131 if (atype1->isGayBerne() && atype2->isGayBerne() ) {
132 interactions_[key].insert(gb_);
133 }
134 if ((atype1->isGayBerne() && atype2->isLennardJones())
135 || (atype1->isLennardJones() && atype2->isGayBerne())) {
136 interactions_[key].insert(gb_);
137 }
138
139 // look for an explicitly-set non-bonded interaction type using the
140 // two atom types.
141 NonBondedInteractionType* nbiType = forceField_->getNonBondedInteractionType(atype1->getName(), atype2->getName());
142
143 if (nbiType != NULL) {
144
145 if (nbiType->isLennardJones()) {
146 // We found an explicit Lennard-Jones interaction.
147 // override all other vdw entries for this pair of atom types:
148 set<NonBondedInteraction*>::iterator it;
149 for (it = interactions_[key].begin();
150 it != interactions_[key].end(); ++it) {
151 InteractionFamily ifam = (*it)->getFamily();
152 if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
153 }
154 interactions_[key].insert(lj_);
155 vdwExplicit = true;
156 }
157
158 if (nbiType->isMorse()) {
159 if (vdwExplicit) {
160 sprintf( painCave.errMsg,
161 "InteractionManager::initialize found more than one "
162 "explicit \n"
163 "\tvan der Waals interaction for atom types %s - %s\n",
164 atype1->getName().c_str(), atype2->getName().c_str());
165 painCave.severity = OPENMD_ERROR;
166 painCave.isFatal = 1;
167 simError();
168 }
169 // We found an explicit Morse interaction.
170 // override all other vdw entries for this pair of atom types:
171 set<NonBondedInteraction*>::iterator it;
172 for (it = interactions_[key].begin();
173 it != interactions_[key].end(); ++it) {
174 InteractionFamily ifam = (*it)->getFamily();
175 if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
176 }
177 interactions_[key].insert(morse_);
178 vdwExplicit = true;
179 }
180
181 if (nbiType->isEAM()) {
182 // We found an explicit EAM interaction.
183 // override all other metallic entries for this pair of atom types:
184 set<NonBondedInteraction*>::iterator it;
185 for (it = interactions_[key].begin();
186 it != interactions_[key].end(); ++it) {
187 InteractionFamily ifam = (*it)->getFamily();
188 if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it);
189 }
190 interactions_[key].insert(eam_);
191 metExplicit = true;
192 }
193
194 if (nbiType->isSC()) {
195 if (metExplicit) {
196 sprintf( painCave.errMsg,
197 "InteractionManager::initialize found more than one "
198 "explicit\n"
199 "\tmetallic interaction for atom types %s - %s\n",
200 atype1->getName().c_str(), atype2->getName().c_str());
201 painCave.severity = OPENMD_ERROR;
202 painCave.isFatal = 1;
203 simError();
204 }
205 // We found an explicit Sutton-Chen interaction.
206 // override all other metallic entries for this pair of atom types:
207 set<NonBondedInteraction*>::iterator it;
208 for (it = interactions_[key].begin();
209 it != interactions_[key].end(); ++it) {
210 InteractionFamily ifam = (*it)->getFamily();
211 if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it);
212 }
213 interactions_[key].insert(sc_);
214 metExplicit = true;
215 }
216
217 if (nbiType->isMAW()) {
218 if (vdwExplicit) {
219 sprintf( painCave.errMsg,
220 "InteractionManager::initialize found more than one "
221 "explicit\n"
222 "\tvan der Waals interaction for atom types %s - %s\n",
223 atype1->getName().c_str(), atype2->getName().c_str());
224 painCave.severity = OPENMD_ERROR;
225 painCave.isFatal = 1;
226 simError();
227 }
228 // We found an explicit MAW interaction.
229 // override all other vdw entries for this pair of atom types:
230 set<NonBondedInteraction*>::iterator it;
231 for (it = interactions_[key].begin();
232 it != interactions_[key].end(); ++it) {
233 InteractionFamily ifam = (*it)->getFamily();
234 if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
235 }
236 interactions_[key].insert(maw_);
237 vdwExplicit = true;
238 }
239 }
240 }
241 }
242
243
244 // Make sure every pair of atom types in this simulation has a
245 // non-bonded interaction. If not, just inform the user.
246
247 set<AtomType*> simTypes = info_->getSimulatedAtomTypes();
248 set<AtomType*>::iterator it, jt;
249
250 for (it = simTypes.begin(); it != simTypes.end(); ++it) {
251 atype1 = (*it);
252 for (jt = it; jt != simTypes.end(); ++jt) {
253 atype2 = (*jt);
254 key = make_pair(atype1, atype2);
255
256 if (interactions_[key].size() == 0) {
257 sprintf( painCave.errMsg,
258 "InteractionManager could not find a matching non-bonded\n"
259 "\tinteraction for atom types %s - %s\n"
260 "\tProceeding without this interaction.\n",
261 atype1->getName().c_str(), atype2->getName().c_str());
262 painCave.severity = OPENMD_INFO;
263 painCave.isFatal = 0;
264 simError();
265 }
266 }
267 }
268
269 initialized_ = true;
270 }
271
272 void InteractionManager::doPrePair(InteractionData idat){
273
274 if (!initialized_) initialize();
275
276 set<NonBondedInteraction*>::iterator it;
277
278 for (it = interactions_[ idat.atypes ].begin();
279 it != interactions_[ idat.atypes ].end(); ++it){
280 if ((*it)->getFamily() == METALLIC_FAMILY) {
281 dynamic_cast<MetallicInteraction*>(*it)->calcDensity(idat);
282 }
283 }
284
285 return;
286 }
287
288 void InteractionManager::doPreForce(SelfData sdat){
289
290 if (!initialized_) initialize();
291
292 pair<AtomType*, AtomType*> key = make_pair(sdat.atype, sdat.atype);
293 set<NonBondedInteraction*>::iterator it;
294
295 for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
296 if ((*it)->getFamily() == METALLIC_FAMILY) {
297 dynamic_cast<MetallicInteraction*>(*it)->calcFunctional(sdat);
298 }
299 }
300
301 return;
302 }
303
304 void InteractionManager::doPair(InteractionData idat){
305
306 if (!initialized_) initialize();
307
308 set<NonBondedInteraction*>::iterator it;
309
310 for (it = interactions_[ idat.atypes ].begin();
311 it != interactions_[ idat.atypes ].end(); ++it)
312 (*it)->calcForce(idat);
313
314 return;
315 }
316
317 void InteractionManager::doSkipCorrection(InteractionData idat){
318
319 if (!initialized_) initialize();
320
321 set<NonBondedInteraction*>::iterator it;
322
323 for (it = interactions_[ idat.atypes ].begin();
324 it != interactions_[ idat.atypes ].end(); ++it){
325 if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
326 dynamic_cast<ElectrostaticInteraction*>(*it)->calcSkipCorrection(idat);
327 }
328 }
329
330 return;
331 }
332
333 void InteractionManager::doSelfCorrection(SelfData sdat){
334
335 if (!initialized_) initialize();
336
337 pair<AtomType*, AtomType*> key = make_pair(sdat.atype, sdat.atype);
338 set<NonBondedInteraction*>::iterator it;
339
340 for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
341 if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
342 dynamic_cast<ElectrostaticInteraction*>(*it)->calcSelfCorrection(sdat);
343 }
344 }
345
346 return;
347 }
348
349 RealType InteractionManager::getSuggestedCutoffRadius(int *atid) {
350 if (!initialized_) initialize();
351
352 AtomType* atype = typeMap_[*atid];
353
354 pair<AtomType*, AtomType*> key = make_pair(atype, atype);
355 set<NonBondedInteraction*>::iterator it;
356 RealType cutoff = 0.0;
357
358 for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
359 cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key));
360 return cutoff;
361 }
362
363 RealType InteractionManager::getSuggestedCutoffRadius(AtomType* atype) {
364 if (!initialized_) initialize();
365
366 pair<AtomType*, AtomType*> key = make_pair(atype, atype);
367 set<NonBondedInteraction*>::iterator it;
368 RealType cutoff = 0.0;
369
370 for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
371 cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key));
372 return cutoff;
373 }
374 } //end namespace OpenMD

Properties

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