ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/nonbonded/InteractionManager.cpp
Revision: 1576
Committed: Wed Jun 8 16:05:07 2011 UTC (13 years, 10 months ago) by gezelter
File size: 13838 byte(s)
Log Message:
migrating cutoff information from InteractionManager to ForceManager

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:
246
247 set<AtomType*> simTypes = info_->getSimulatedAtomTypes();
248 set<AtomType*>::iterator it, jt;
249 for (it = simTypes.begin(); it != simTypes.end(); ++it) {
250 atype1 = (*it);
251 for (jt = simTypes.begin(); jt != simTypes.end(); ++jt) {
252 atype2 = (*jt);
253 key = make_pair(atype1, atype2);
254
255 if (interactions_[key].size() == 0) {
256 sprintf( painCave.errMsg,
257 "InteractionManager unable to find an appropriate non-bonded\n"
258 "\tinteraction for atom types %s - %s\n",
259 atype1->getName().c_str(), atype2->getName().c_str());
260 painCave.severity = OPENMD_INFO;
261 painCave.isFatal = 1;
262 simError();
263 }
264 }
265 }
266
267 initialized_ = true;
268 }
269
270 void InteractionManager::doPrePair(InteractionData idat){
271
272 if (!initialized_) initialize();
273
274 set<NonBondedInteraction*>::iterator it;
275
276 for (it = interactions_[ idat.atypes ].begin();
277 it != interactions_[ idat.atypes ].end(); ++it){
278 if ((*it)->getFamily() == METALLIC_FAMILY) {
279 dynamic_cast<MetallicInteraction*>(*it)->calcDensity(idat);
280 }
281 }
282
283 return;
284 }
285
286 void InteractionManager::doPreForce(SelfData sdat){
287
288 if (!initialized_) initialize();
289
290 pair<AtomType*, AtomType*> key = make_pair(sdat.atype, sdat.atype);
291 set<NonBondedInteraction*>::iterator it;
292
293 for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
294 if ((*it)->getFamily() == METALLIC_FAMILY) {
295 dynamic_cast<MetallicInteraction*>(*it)->calcFunctional(sdat);
296 }
297 }
298
299 return;
300 }
301
302 void InteractionManager::doPair(InteractionData idat){
303
304 if (!initialized_) initialize();
305
306 set<NonBondedInteraction*>::iterator it;
307
308 for (it = interactions_[ idat.atypes ].begin();
309 it != interactions_[ idat.atypes ].end(); ++it)
310 (*it)->calcForce(idat);
311
312 return;
313 }
314
315 void InteractionManager::doSkipCorrection(InteractionData idat){
316
317 if (!initialized_) initialize();
318
319 set<NonBondedInteraction*>::iterator it;
320
321 for (it = interactions_[ idat.atypes ].begin();
322 it != interactions_[ idat.atypes ].end(); ++it){
323 if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
324 dynamic_cast<ElectrostaticInteraction*>(*it)->calcSkipCorrection(idat);
325 }
326 }
327
328 return;
329 }
330
331 void InteractionManager::doSelfCorrection(SelfData sdat){
332
333 if (!initialized_) initialize();
334
335 pair<AtomType*, AtomType*> key = make_pair(sdat.atype, sdat.atype);
336 set<NonBondedInteraction*>::iterator it;
337
338 for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
339 if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
340 dynamic_cast<ElectrostaticInteraction*>(*it)->calcSelfCorrection(sdat);
341 }
342 }
343
344 return;
345 }
346
347 RealType InteractionManager::getSuggestedCutoffRadius(int *atid) {
348 if (!initialized_) initialize();
349
350 AtomType* atype = typeMap_[*atid];
351
352 pair<AtomType*, AtomType*> key = make_pair(atype, atype);
353 set<NonBondedInteraction*>::iterator it;
354 RealType cutoff = 0.0;
355
356 for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
357 cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key));
358 return cutoff;
359 }
360
361 RealType InteractionManager::getSuggestedCutoffRadius(AtomType* atype) {
362 if (!initialized_) initialize();
363
364 pair<AtomType*, AtomType*> key = make_pair(atype, atype);
365 set<NonBondedInteraction*>::iterator it;
366 RealType cutoff = 0.0;
367
368 for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
369 cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key));
370 return cutoff;
371 }
372 } //end namespace OpenMD

Properties

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