ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/nonbonded/LJ.cpp
(Generate patch)

Comparing branches/development/src/nonbonded/LJ.cpp (file contents):
Revision 1472 by gezelter, Mon Jul 19 20:42:55 2010 UTC vs.
Revision 1683 by jmichalk, Wed Feb 29 20:33:01 2012 UTC

# Line 36 | Line 36
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).                        
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include <stdio.h>
# Line 45 | Line 46
46   #include <cmath>
47   #include "nonbonded/LJ.hpp"
48   #include "utils/simError.h"
49 + #include "types/LennardJonesInteractionType.hpp"
50  
49 –
51   namespace OpenMD {
52  
53 <  bool LJ::initialized_ = false;
53 <  bool LJ::shiftedPot_ = false;
54 <  bool LJ::shiftedFrc_ = false;
55 <  ForceField* LJ::forceField_ = NULL;
56 <  std::map<int, AtomType*> LJ::LJMap;
57 <  std::map<std::pair<AtomType*, AtomType*>, LJInteractionData> LJ::MixingMap;
58 <  
59 <  LJ* LJ::_instance = NULL;
53 >  LJ::LJ() : name_("LJ"), initialized_(false), forceField_(NULL) {}
54  
61 –  LJ* LJ::Instance() {
62 –    if (!_instance) {
63 –      _instance = new LJ();
64 –    }
65 –    return _instance;
66 –  }
67 –
55    LJParam LJ::getLJParam(AtomType* atomType) {
56      
57      // Do sanity checking on the AtomType we were passed before
# Line 106 | Line 93 | namespace OpenMD {
93      return ljParam.sigma;
94    }
95  
109 –  RealType LJ::getSigma(int atid) {
110 –    std::map<int, AtomType*> :: const_iterator it;
111 –    it = LJMap.find(atid);
112 –    if (it == LJMap.end()) {
113 –      sprintf( painCave.errMsg,
114 –               "LJ::getSigma could not find atid %d in LJMap\n",
115 –               (atid));
116 –      painCave.severity = OPENMD_ERROR;
117 –      painCave.isFatal = 1;
118 –      simError();          
119 –    }
120 –    AtomType* atype = it->second;
121 –
122 –    return getSigma(atype);    
123 –  }
124 –
96    RealType LJ::getSigma(AtomType* atomType1, AtomType* atomType2) {    
97      RealType sigma1 = getSigma(atomType1);
98      RealType sigma2 = getSigma(atomType2);
99      
100      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
101 <    std::string DistanceMix = fopts.getDistanceMixingRule();
101 >    string DistanceMix = fopts.getDistanceMixingRule();
102      toUpper(DistanceMix);
103  
104      if (DistanceMix == "GEOMETRIC")
# Line 141 | Line 112 | namespace OpenMD {
112      return ljParam.epsilon;
113    }
114  
144 –  RealType LJ::getEpsilon(int atid) {    
145 –    std::map<int, AtomType*> :: const_iterator it;
146 –    it = LJMap.find(atid);
147 –    if (it == LJMap.end()) {
148 –      sprintf( painCave.errMsg,
149 –               "LJ::getEpsilon could not find atid %d in LJMap\n",
150 –               (atid));
151 –      painCave.severity = OPENMD_ERROR;
152 –      painCave.isFatal = 1;
153 –      simError();          
154 –    }
155 –    AtomType* atype = it->second;
156 –
157 –    return getEpsilon(atype);    
158 –  }
159 –
115    RealType LJ::getEpsilon(AtomType* atomType1, AtomType* atomType2) {    
116      RealType epsilon1 = getEpsilon(atomType1);
117      RealType epsilon2 = getEpsilon(atomType2);
# Line 164 | Line 119 | namespace OpenMD {
119    }
120  
121    void LJ::initialize() {    
122 <    ForceField::AtomTypeContainer atomTypes = forceField_->getAtomTypes();
122 >    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
123      ForceField::AtomTypeContainer::MapTypeIterator i;
124      AtomType* at;
125  
126 <    for (at = atomTypes.beginType(i); at != NULL;
127 <         at = atomTypes.nextType(i)) {
128 <      
129 <      if (at->isLennardJones())
130 <        addType(at);
126 >    for (at = atomTypes->beginType(i); at != NULL;
127 >         at = atomTypes->nextType(i)) {
128 >      if (at->isLennardJones()){
129 >         addType(at);
130 >      }
131      }
132 <
178 <    ForceField::NonBondedInteractionTypeContainer nbiTypes = forceField_->getNonBondedInteractionTypes();
132 >    ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
133      ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
134      NonBondedInteractionType* nbt;
135 +    ForceField::NonBondedInteractionTypeContainer::KeyType keys;
136  
137 <    for (nbt = nbiTypes.beginType(j); nbt != NULL;
138 <         nbt = nbiTypes.nextType(j)) {
139 <      
137 >
138 >    for (nbt = nbiTypes->beginType(j); nbt != NULL;
139 >         nbt = nbiTypes->nextType(j)) {
140 >
141        if (nbt->isLennardJones()) {
142 <        
143 <        std::pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
144 <        
145 <        GenericData* data = nbt->getPropertyByName("LennardJones");
146 <        if (data == NULL) {
147 <          sprintf( painCave.errMsg, "LJ::rebuildMixingMap could not find\n"
148 <               "\tLennard-Jones parameters for %s - %s interaction.\n",
193 <                   atypes.first->getName().c_str(),
194 <                   atypes.second->getName().c_str());
195 <          painCave.severity = OPENMD_ERROR;
196 <          painCave.isFatal = 1;
197 <          simError();
198 <        }
199 <    
200 <        LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
201 <        if (ljData == NULL) {
142 >        keys = nbiTypes->getKeys(j);
143 >        AtomType* at1 = forceField_->getAtomType(keys[0]);
144 >        AtomType* at2 = forceField_->getAtomType(keys[1]);
145 >
146 >        LennardJonesInteractionType* ljit = dynamic_cast<LennardJonesInteractionType*>(nbt);
147 >
148 >        if (ljit == NULL) {
149            sprintf( painCave.errMsg,
150 <                   "LJ::rebuildMixingMap could not convert GenericData to\n"
151 <                   "\tLJParam for %s - %s interaction.\n",
152 <                   atypes.first->getName().c_str(),
153 <                   atypes.second->getName().c_str());
150 >                   "LJ::initialize could not convert NonBondedInteractionType\n"
151 >                   "\tto LennardJonesInteractionType for %s - %s interaction.\n",
152 >                   at1->getName().c_str(),
153 >                   at2->getName().c_str());
154            painCave.severity = OPENMD_ERROR;
155            painCave.isFatal = 1;
156            simError();          
157          }
211 –        
212 –        LJParam ljParam = ljData->getData();
158  
159 <        RealType sigma = ljParam.sigma;
160 <        RealType epsilon = ljParam.epsilon;
161 <
217 <        addExplicitInteraction(atypes.first, atypes.second, sigma, epsilon);
159 >        RealType sigma = ljit->getSigma();
160 >        RealType epsilon = ljit->getEpsilon();
161 >        addExplicitInteraction(at1, at2, sigma, epsilon);
162        }
163      }  
164      initialized_ = true;
# Line 225 | Line 169 | namespace OpenMD {
169    void LJ::addType(AtomType* atomType){
170      RealType sigma1 = getSigma(atomType);
171      RealType epsilon1 = getEpsilon(atomType);
172 <
172 >    
173      // add it to the map:
174      AtomTypeProperties atp = atomType->getATP();    
175 <    std::pair<std::map<int,AtomType*>::iterator,bool> ret;
176 <    ret = LJMap.insert( std::pair<int, AtomType*>(atp.ident, atomType) );
175 >
176 >    pair<map<int,AtomType*>::iterator,bool> ret;    
177 >    ret = LJMap.insert( pair<int, AtomType*>(atp.ident, atomType) );
178      if (ret.second == false) {
179        sprintf( painCave.errMsg,
180                 "LJ already had a previous entry with ident %d\n",
# Line 265 | Line 210 | namespace OpenMD {
210    
211    void LJ::addExplicitInteraction(AtomType* atype1, AtomType* atype2, RealType sigma, RealType epsilon){
212      
268 –    // in case these weren't already in the map
269 –    addType(atype1);
270 –    addType(atype2);
271 –
213      LJInteractionData mixer;
214      mixer.sigma = sigma;
215      mixer.epsilon = epsilon;
# Line 285 | Line 226 | namespace OpenMD {
226      }    
227    }
228  
229 <  void LJ::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
289 <                     RealType rij, RealType r2, RealType rcut, RealType sw,
290 <                     RealType vdwMult, RealType &vpair, RealType &pot,
291 <                     Vector3d &f1) {
292 <
229 >  void LJ::calcForce(InteractionData &idat) {
230      if (!initialized_) initialize();
231 +    map<pair<AtomType*, AtomType*>, LJInteractionData>::iterator it;
232 +    it = MixingMap.find( idat.atypes );
233      
234 <    RealType ros;
235 <    RealType rcos;
236 <    RealType myPot = 0.0;
237 <    RealType myPotC = 0.0;
238 <    RealType myDeriv = 0.0;
239 <    RealType myDerivC = 0.0;
234 >    if (it != MixingMap.end())  {
235 >      
236 >      LJInteractionData mixer = (*it).second;
237 >      
238 >      RealType sigmai = mixer.sigmai;
239 >      RealType epsilon = mixer.epsilon;
240 >      
241 >      RealType ros;
242 >      RealType rcos;
243 >      RealType myPot = 0.0;
244 >      RealType myPotC = 0.0;
245 >      RealType myDeriv = 0.0;
246 >      RealType myDerivC = 0.0;
247 >    
248 >      ros = *(idat.rij) * sigmai;    
249 >      
250 >      getLJfunc(ros, myPot, myDeriv);
251 >      
252 >      if (idat.shiftedPot) {
253 >        rcos = *(idat.rcut) * sigmai;
254 >        getLJfunc(rcos, myPotC, myDerivC);
255 >        myDerivC = 0.0;
256 >      } else if (idat.shiftedForce) {
257 >        rcos = *(idat.rcut) * sigmai;
258 >        getLJfunc(rcos, myPotC, myDerivC);
259 >        myPotC = myPotC + myDerivC * (*(idat.rij) - *(idat.rcut)) * sigmai;
260 >      } else {
261 >        myPotC = 0.0;
262 >        myDerivC = 0.0;        
263 >      }
264  
265 <    std::pair<AtomType*, AtomType*> key = std::make_pair(at1, at2);
266 <    LJInteractionData mixer = MixingMap[key];
265 >      RealType pot_temp = *(idat.vdwMult) * epsilon * (myPot - myPotC);
266 >      *(idat.vpair) += pot_temp;
267 >      
268 >      RealType dudr = *(idat.sw) * *(idat.vdwMult) * epsilon * (myDeriv -
269 >                                                                myDerivC)*sigmai;      
270  
271 <    RealType sigmai = mixer.sigmai;
272 <    RealType epsilon = mixer.epsilon;
307 <    
308 <    ros = rij * sigmai;
309 <
310 <    getLJfunc(ros, myPot, myDeriv);
311 <
312 <    if (shiftedPot_) {
313 <      rcos = rcut * sigmai;
314 <      getLJfunc(rcos, myPotC, myDerivC);
315 <      myDerivC = 0.0;
316 <    } else if (LJ::shiftedFrc_) {
317 <      rcos = rcut * sigmai;
318 <      getLJfunc(rcos, myPotC, myDerivC);
319 <      myPotC = myPotC + myDerivC * (rij - rcut) * sigmai;
320 <    } else {
321 <      myPotC = 0.0;
322 <      myDerivC = 0.0;
271 >      (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
272 >      *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
273      }
324 –
325 –    RealType pot_temp = vdwMult * epsilon * (myPot - myPotC);
326 –    vpair += pot_temp;
327 –
328 –    RealType dudr = sw * vdwMult * epsilon * (myDeriv - myDerivC)*sigmai;
329 –    
330 –    pot += sw * pot_temp;
331 –    f1 = d * dudr / rij;
332 –
274      return;
275    }
335 –
336 –  void LJ::do_lj_pair(int *atid1, int *atid2, RealType *d, RealType *rij,
337 –                      RealType *r2, RealType *rcut, RealType *sw,
338 –                      RealType *vdwMult,
339 –                      RealType *vpair, RealType *pot, RealType *f1) {
340 –
341 –    if (!initialized_) initialize();
342 –    
343 –    AtomType* atype1 = LJMap[*atid1];
344 –    AtomType* atype2 = LJMap[*atid2];
345 –    
346 –    Vector3d disp(d[0], d[1], d[2]);
347 –    Vector3d frc(f1[0], f1[1], f1[2]);
348 –    
349 –    calcForce(atype1, atype2, disp, *rij, *r2, *rcut, *sw, *vdwMult, *vpair,
350 –              *pot, frc);
351 –      
352 –    f1[0] = frc.x();
353 –    f1[1] = frc.y();
354 –    f1[2] = frc.z();
355 –    
356 –    return;    
357 –  }
276    
277    void LJ::getLJfunc(RealType r, RealType &pot, RealType &deriv) {
278  
# Line 371 | Line 289 | namespace OpenMD {
289      return;
290    }
291    
292 <
293 <  void LJ::setLJDefaultCutoff(RealType *thisRcut, int *sP, int *sF) {
294 <    shiftedPot_ = (bool)(*sP);
295 <    shiftedFrc_ = (bool)(*sF);
292 >  RealType LJ::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
293 >    if (!initialized_) initialize();  
294 >    map<pair<AtomType*, AtomType*>, LJInteractionData>::iterator it;
295 >    it = MixingMap.find(atypes);
296 >    if (it == MixingMap.end())
297 >      return 0.0;
298 >    else  {
299 >      LJInteractionData mixer = (*it).second;
300 >      return 2.5 * mixer.sigma;
301 >    }
302    }
379 – }
303  
381 – extern "C" {
382 –  
383 – #define fortranGetSigma FC_FUNC(getsigma, GETSIGMA)
384 – #define fortranGetEpsilon FC_FUNC(getepsilon, GETEPSILON)
385 – #define fortranSetLJCutoff FC_FUNC(setljdefaultcutoff, SETLJDEFAULTCUTOFF)
386 – #define fortranDoLJPair FC_FUNC(do_lj_pair, DO_LJ_PAIR)
387 –  
388 –  RealType fortranGetSigma(int* atid) {
389 –    return OpenMD::LJ::Instance()->getSigma(*atid);
390 –  }
391 –  RealType fortranGetEpsilon(int* atid) {  
392 –    return OpenMD::LJ::Instance()->getEpsilon(*atid);
393 –  }
394 –  void fortranSetLJCutoff(RealType *rcut, int *shiftedPot, int *shiftedFrc) {
395 –    return OpenMD::LJ::Instance()->setLJDefaultCutoff(rcut, shiftedPot,
396 –                                                      shiftedFrc);
397 –  }
398 –  void fortranDoLJPair(int *atid1, int *atid2, RealType *d, RealType *rij,
399 –                       RealType *r2, RealType *rcut, RealType *sw,
400 –                       RealType *vdwMult, RealType* vpair, RealType* pot,
401 –                       RealType *f1){
402 –    
403 –    return OpenMD::LJ::Instance()->do_lj_pair(atid1, atid2, d, rij, r2, rcut,
404 –                                              sw, vdwMult, vpair, pot, f1);
405 –  }
304   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines