ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/nonbonded/LJ.cpp
Revision: 1665
Committed: Tue Nov 22 20:38:56 2011 UTC (13 years, 6 months ago) by gezelter
File size: 10632 byte(s)
Log Message:
updated copyright notices

File Contents

# User Rev Content
1 gezelter 1471 /*
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 gezelter 1665 * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40     * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41 gezelter 1471 */
42    
43     #include <stdio.h>
44     #include <string.h>
45    
46     #include <cmath>
47     #include "nonbonded/LJ.hpp"
48     #include "utils/simError.h"
49    
50     namespace OpenMD {
51    
52 gezelter 1583 LJ::LJ() : name_("LJ"), initialized_(false), forceField_(NULL) {}
53 gezelter 1471
54     LJParam LJ::getLJParam(AtomType* atomType) {
55    
56     // Do sanity checking on the AtomType we were passed before
57     // building any data structures:
58     if (!atomType->isLennardJones()) {
59     sprintf( painCave.errMsg,
60     "LJ::getLJParam was passed an atomType (%s) that does not\n"
61     "\tappear to be a Lennard-Jones atom.\n",
62     atomType->getName().c_str());
63     painCave.severity = OPENMD_ERROR;
64     painCave.isFatal = 1;
65     simError();
66     }
67    
68     GenericData* data = atomType->getPropertyByName("LennardJones");
69     if (data == NULL) {
70     sprintf( painCave.errMsg, "LJ::getLJParam could not find Lennard-Jones\n"
71     "\tparameters for atomType %s.\n", atomType->getName().c_str());
72     painCave.severity = OPENMD_ERROR;
73     painCave.isFatal = 1;
74     simError();
75     }
76    
77     LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
78     if (ljData == NULL) {
79     sprintf( painCave.errMsg,
80     "LJ::getLJParam could not convert GenericData to LJParam for\n"
81     "\tatom type %s\n", atomType->getName().c_str());
82     painCave.severity = OPENMD_ERROR;
83     painCave.isFatal = 1;
84     simError();
85     }
86    
87     return ljData->getData();
88     }
89    
90     RealType LJ::getSigma(AtomType* atomType) {
91     LJParam ljParam = getLJParam(atomType);
92     return ljParam.sigma;
93     }
94    
95     RealType LJ::getSigma(AtomType* atomType1, AtomType* atomType2) {
96     RealType sigma1 = getSigma(atomType1);
97     RealType sigma2 = getSigma(atomType2);
98    
99     ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
100 gezelter 1502 string DistanceMix = fopts.getDistanceMixingRule();
101 gezelter 1471 toUpper(DistanceMix);
102    
103     if (DistanceMix == "GEOMETRIC")
104     return sqrt(sigma1 * sigma2);
105     else
106     return 0.5 * (sigma1 + sigma2);
107     }
108    
109     RealType LJ::getEpsilon(AtomType* atomType) {
110     LJParam ljParam = getLJParam(atomType);
111     return ljParam.epsilon;
112     }
113    
114     RealType LJ::getEpsilon(AtomType* atomType1, AtomType* atomType2) {
115     RealType epsilon1 = getEpsilon(atomType1);
116     RealType epsilon2 = getEpsilon(atomType2);
117     return sqrt(epsilon1 * epsilon2);
118     }
119    
120     void LJ::initialize() {
121 gezelter 1476 ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
122 gezelter 1471 ForceField::AtomTypeContainer::MapTypeIterator i;
123     AtomType* at;
124    
125 gezelter 1476 for (at = atomTypes->beginType(i); at != NULL;
126     at = atomTypes->nextType(i)) {
127 gezelter 1471
128     if (at->isLennardJones())
129     addType(at);
130     }
131    
132 gezelter 1476 ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
133 gezelter 1471 ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
134     NonBondedInteractionType* nbt;
135    
136 gezelter 1476 for (nbt = nbiTypes->beginType(j); nbt != NULL;
137     nbt = nbiTypes->nextType(j)) {
138 gezelter 1471
139     if (nbt->isLennardJones()) {
140    
141 gezelter 1502 pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
142 gezelter 1471
143     GenericData* data = nbt->getPropertyByName("LennardJones");
144     if (data == NULL) {
145     sprintf( painCave.errMsg, "LJ::rebuildMixingMap could not find\n"
146     "\tLennard-Jones parameters for %s - %s interaction.\n",
147     atypes.first->getName().c_str(),
148     atypes.second->getName().c_str());
149     painCave.severity = OPENMD_ERROR;
150     painCave.isFatal = 1;
151     simError();
152     }
153    
154     LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
155     if (ljData == NULL) {
156     sprintf( painCave.errMsg,
157     "LJ::rebuildMixingMap could not convert GenericData to\n"
158     "\tLJParam for %s - %s interaction.\n",
159     atypes.first->getName().c_str(),
160     atypes.second->getName().c_str());
161     painCave.severity = OPENMD_ERROR;
162     painCave.isFatal = 1;
163     simError();
164     }
165    
166     LJParam ljParam = ljData->getData();
167    
168     RealType sigma = ljParam.sigma;
169     RealType epsilon = ljParam.epsilon;
170    
171     addExplicitInteraction(atypes.first, atypes.second, sigma, epsilon);
172     }
173     }
174     initialized_ = true;
175     }
176    
177    
178    
179     void LJ::addType(AtomType* atomType){
180     RealType sigma1 = getSigma(atomType);
181     RealType epsilon1 = getEpsilon(atomType);
182 gezelter 1502
183 gezelter 1471 // add it to the map:
184     AtomTypeProperties atp = atomType->getATP();
185 gezelter 1473
186 gezelter 1502 pair<map<int,AtomType*>::iterator,bool> ret;
187     ret = LJMap.insert( pair<int, AtomType*>(atp.ident, atomType) );
188 gezelter 1471 if (ret.second == false) {
189     sprintf( painCave.errMsg,
190     "LJ already had a previous entry with ident %d\n",
191     atp.ident);
192     painCave.severity = OPENMD_INFO;
193     painCave.isFatal = 0;
194     simError();
195     }
196    
197     // Now, iterate over all known types and add to the mixing map:
198    
199     std::map<int, AtomType*>::iterator it;
200     for( it = LJMap.begin(); it != LJMap.end(); ++it) {
201    
202     AtomType* atype2 = (*it).second;
203    
204     LJInteractionData mixer;
205     mixer.sigma = getSigma(atomType, atype2);
206     mixer.epsilon = getEpsilon(atomType, atype2);
207     mixer.sigmai = 1.0 / mixer.sigma;
208     mixer.explicitlySet = false;
209    
210     std::pair<AtomType*, AtomType*> key1, key2;
211     key1 = std::make_pair(atomType, atype2);
212     key2 = std::make_pair(atype2, atomType);
213    
214     MixingMap[key1] = mixer;
215     if (key2 != key1) {
216     MixingMap[key2] = mixer;
217     }
218     }
219     }
220    
221     void LJ::addExplicitInteraction(AtomType* atype1, AtomType* atype2, RealType sigma, RealType epsilon){
222    
223     // in case these weren't already in the map
224     addType(atype1);
225     addType(atype2);
226    
227     LJInteractionData mixer;
228     mixer.sigma = sigma;
229     mixer.epsilon = epsilon;
230     mixer.sigmai = 1.0 / mixer.sigma;
231     mixer.explicitlySet = true;
232    
233     std::pair<AtomType*, AtomType*> key1, key2;
234     key1 = std::make_pair(atype1, atype2);
235     key2 = std::make_pair(atype2, atype1);
236    
237     MixingMap[key1] = mixer;
238     if (key2 != key1) {
239     MixingMap[key2] = mixer;
240     }
241     }
242    
243 gezelter 1536 void LJ::calcForce(InteractionData &idat) {
244 gezelter 1502
245 gezelter 1471 if (!initialized_) initialize();
246 gezelter 1587
247 gezelter 1505 map<pair<AtomType*, AtomType*>, LJInteractionData>::iterator it;
248 gezelter 1571 it = MixingMap.find( idat.atypes );
249 gezelter 1471
250 gezelter 1505 if (it != MixingMap.end()) {
251    
252     LJInteractionData mixer = (*it).second;
253    
254     RealType sigmai = mixer.sigmai;
255     RealType epsilon = mixer.epsilon;
256    
257     RealType ros;
258     RealType rcos;
259     RealType myPot = 0.0;
260     RealType myPotC = 0.0;
261     RealType myDeriv = 0.0;
262     RealType myDerivC = 0.0;
263    
264 gezelter 1582 ros = *(idat.rij) * sigmai;
265 gezelter 1505
266     getLJfunc(ros, myPot, myDeriv);
267    
268 gezelter 1583 if (idat.shiftedPot) {
269 gezelter 1554 rcos = *(idat.rcut) * sigmai;
270 gezelter 1505 getLJfunc(rcos, myPotC, myDerivC);
271     myDerivC = 0.0;
272 gezelter 1583 } else if (idat.shiftedForce) {
273 gezelter 1554 rcos = *(idat.rcut) * sigmai;
274 gezelter 1505 getLJfunc(rcos, myPotC, myDerivC);
275 gezelter 1554 myPotC = myPotC + myDerivC * (*(idat.rij) - *(idat.rcut)) * sigmai;
276 gezelter 1505 } else {
277     myPotC = 0.0;
278     myDerivC = 0.0;
279     }
280 gezelter 1473
281 gezelter 1554 RealType pot_temp = *(idat.vdwMult) * epsilon * (myPot - myPotC);
282     *(idat.vpair) += pot_temp;
283 gezelter 1505
284 gezelter 1554 RealType dudr = *(idat.sw) * *(idat.vdwMult) * epsilon * (myDeriv -
285 gezelter 1582 myDerivC)*sigmai;
286 gezelter 1583
287 gezelter 1582 (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
288 gezelter 1629 *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
289 gezelter 1471 }
290     return;
291     }
292    
293 gezelter 1472 void LJ::getLJfunc(RealType r, RealType &pot, RealType &deriv) {
294    
295 gezelter 1471 RealType ri = 1.0 / r;
296     RealType ri2 = ri * ri;
297     RealType ri6 = ri2 * ri2 * ri2;
298     RealType ri7 = ri6 * ri;
299     RealType ri12 = ri6 * ri6;
300     RealType ri13 = ri12 * ri;
301 gezelter 1472
302 gezelter 1471 pot = 4.0 * (ri12 - ri6);
303     deriv = 24.0 * (ri7 - 2.0 * ri13);
304 gezelter 1472
305 gezelter 1471 return;
306     }
307 gezelter 1472
308 gezelter 1545 RealType LJ::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
309 gezelter 1505 if (!initialized_) initialize();
310     map<pair<AtomType*, AtomType*>, LJInteractionData>::iterator it;
311 gezelter 1545 it = MixingMap.find(atypes);
312 gezelter 1505 if (it == MixingMap.end())
313     return 0.0;
314     else {
315     LJInteractionData mixer = (*it).second;
316     return 2.5 * mixer.sigma;
317     }
318     }
319 gezelter 1471
320     }

Properties

Name Value
svn:eol-style native