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

# User Rev Content
1 gezelter 1502 /*
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 gezelter 1535
46 gezelter 1576 InteractionManager::InteractionManager() {
47 gezelter 1535
48 gezelter 1576 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 gezelter 1502 }
59    
60     void InteractionManager::initialize() {
61    
62 gezelter 1535 ForceField* forceField_ = info_->getForceField();
63    
64 gezelter 1502 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 gezelter 1532 maw_->setForceField(forceField_);
72 gezelter 1502
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 gezelter 1535
143     if (nbiType != NULL) {
144 gezelter 1502
145 gezelter 1535 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 gezelter 1502 }
157 gezelter 1535
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 gezelter 1502 }
180 gezelter 1535
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 gezelter 1502 }
193 gezelter 1535
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 gezelter 1502 }
216 gezelter 1535
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 gezelter 1502 }
240     }
241     }
242    
243 gezelter 1535
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 gezelter 1502 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 gezelter 1535
267     initialized_ = true;
268     }
269 gezelter 1502
270 gezelter 1545 void InteractionManager::doPrePair(InteractionData idat){
271 gezelter 1502
272     if (!initialized_) initialize();
273 gezelter 1545
274 gezelter 1504 set<NonBondedInteraction*>::iterator it;
275    
276 gezelter 1571 for (it = interactions_[ idat.atypes ].begin();
277     it != interactions_[ idat.atypes ].end(); ++it){
278 gezelter 1504 if ((*it)->getFamily() == METALLIC_FAMILY) {
279 gezelter 1545 dynamic_cast<MetallicInteraction*>(*it)->calcDensity(idat);
280 gezelter 1504 }
281     }
282 gezelter 1502
283     return;
284     }
285 gezelter 1504
286 gezelter 1545 void InteractionManager::doPreForce(SelfData sdat){
287 gezelter 1502
288 gezelter 1504 if (!initialized_) initialize();
289 gezelter 1545
290     pair<AtomType*, AtomType*> key = make_pair(sdat.atype, sdat.atype);
291 gezelter 1504 set<NonBondedInteraction*>::iterator it;
292 gezelter 1502
293 gezelter 1504 for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
294     if ((*it)->getFamily() == METALLIC_FAMILY) {
295 gezelter 1545 dynamic_cast<MetallicInteraction*>(*it)->calcFunctional(sdat);
296 gezelter 1504 }
297     }
298    
299 gezelter 1502 return;
300     }
301    
302 gezelter 1545 void InteractionManager::doPair(InteractionData idat){
303 gezelter 1504
304 gezelter 1502 if (!initialized_) initialize();
305 gezelter 1545
306 gezelter 1502 set<NonBondedInteraction*>::iterator it;
307    
308 gezelter 1571 for (it = interactions_[ idat.atypes ].begin();
309     it != interactions_[ idat.atypes ].end(); ++it)
310 gezelter 1502 (*it)->calcForce(idat);
311    
312     return;
313     }
314    
315 gezelter 1545 void InteractionManager::doSkipCorrection(InteractionData idat){
316 gezelter 1502
317 gezelter 1545 if (!initialized_) initialize();
318 gezelter 1504
319     set<NonBondedInteraction*>::iterator it;
320    
321 gezelter 1571 for (it = interactions_[ idat.atypes ].begin();
322     it != interactions_[ idat.atypes ].end(); ++it){
323 gezelter 1504 if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
324 gezelter 1545 dynamic_cast<ElectrostaticInteraction*>(*it)->calcSkipCorrection(idat);
325 gezelter 1504 }
326     }
327 gezelter 1502
328 gezelter 1504 return;
329 gezelter 1502 }
330    
331 gezelter 1545 void InteractionManager::doSelfCorrection(SelfData sdat){
332 gezelter 1502
333     if (!initialized_) initialize();
334 gezelter 1504
335 gezelter 1545 pair<AtomType*, AtomType*> key = make_pair(sdat.atype, sdat.atype);
336 gezelter 1504 set<NonBondedInteraction*>::iterator it;
337 gezelter 1502
338 gezelter 1504 for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
339     if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
340 gezelter 1545 dynamic_cast<ElectrostaticInteraction*>(*it)->calcSelfCorrection(sdat);
341 gezelter 1504 }
342     }
343 gezelter 1536
344 gezelter 1504 return;
345 gezelter 1502 }
346    
347 gezelter 1505 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 gezelter 1545 cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key));
358 gezelter 1505 return cutoff;
359     }
360    
361 gezelter 1528 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 gezelter 1545 cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key));
370 gezelter 1528 return cutoff;
371     }
372 gezelter 1502 } //end namespace OpenMD

Properties

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