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 |