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 1473 by gezelter, Tue Jul 20 15:43:00 2010 UTC vs.
Revision 1502 by gezelter, Sat Oct 2 19:53:32 2010 UTC

# Line 46 | Line 46
46   #include "nonbonded/LJ.hpp"
47   #include "utils/simError.h"
48  
49 –
49   namespace OpenMD {
50  
51 <  bool LJ::initialized_ = false;
52 <  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;
51 >  LJ::LJ() : name_("LJ"), initialized_(false), shiftedPot_(false),
52 >             shiftedFrc_(false), forceField_(NULL) {}
53  
61 –  LJ* LJ::Instance() {
62 –    if (!_instance) {
63 –      _instance = new LJ();
64 –    }
65 –    return _instance;
66 –  }
67 –
54    LJParam LJ::getLJParam(AtomType* atomType) {
55      
56      // Do sanity checking on the AtomType we were passed before
# Line 106 | Line 92 | namespace OpenMD {
92      return ljParam.sigma;
93    }
94  
109 –  RealType LJ::getSigma(int atid) {
110 –    if (!initialized_) initialize();
111 –    std::map<int, AtomType*> :: const_iterator it;
112 –    it = LJMap.find(atid);
113 –    if (it == LJMap.end()) {
114 –      sprintf( painCave.errMsg,
115 –               "LJ::getSigma could not find atid %d in LJMap\n",
116 –               (atid));
117 –      painCave.severity = OPENMD_ERROR;
118 –      painCave.isFatal = 1;
119 –      simError();          
120 –    }
121 –    AtomType* atype = it->second;
122 –
123 –    return getSigma(atype);    
124 –  }
125 –
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 <    std::string DistanceMix = fopts.getDistanceMixingRule();
100 >    string DistanceMix = fopts.getDistanceMixingRule();
101      toUpper(DistanceMix);
102  
103      if (DistanceMix == "GEOMETRIC")
# Line 140 | Line 109 | namespace OpenMD {
109    RealType LJ::getEpsilon(AtomType* atomType) {    
110      LJParam ljParam = getLJParam(atomType);
111      return ljParam.epsilon;
143 –  }
144 –
145 –  RealType LJ::getEpsilon(int atid) {    
146 –    if (!initialized_) initialize();
147 –    std::map<int, AtomType*> :: const_iterator it;
148 –    it = LJMap.find(atid);
149 –    if (it == LJMap.end()) {
150 –      sprintf( painCave.errMsg,
151 –               "LJ::getEpsilon could not find atid %d in LJMap\n",
152 –               (atid));
153 –      painCave.severity = OPENMD_ERROR;
154 –      painCave.isFatal = 1;
155 –      simError();          
156 –    }
157 –    AtomType* atype = it->second;
158 –
159 –    return getEpsilon(atype);    
112    }
113  
114    RealType LJ::getEpsilon(AtomType* atomType1, AtomType* atomType2) {    
# Line 166 | Line 118 | namespace OpenMD {
118    }
119  
120    void LJ::initialize() {    
121 <    ForceField::AtomTypeContainer atomTypes = forceField_->getAtomTypes();
121 >    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
122      ForceField::AtomTypeContainer::MapTypeIterator i;
123      AtomType* at;
124  
125 <    for (at = atomTypes.beginType(i); at != NULL;
126 <         at = atomTypes.nextType(i)) {
125 >    for (at = atomTypes->beginType(i); at != NULL;
126 >         at = atomTypes->nextType(i)) {
127        
128        if (at->isLennardJones())
129          addType(at);
130      }
131  
132 <    ForceField::NonBondedInteractionTypeContainer nbiTypes = forceField_->getNonBondedInteractionTypes();
132 >    ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
133      ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
134      NonBondedInteractionType* nbt;
135  
136 <    for (nbt = nbiTypes.beginType(j); nbt != NULL;
137 <         nbt = nbiTypes.nextType(j)) {
136 >    for (nbt = nbiTypes->beginType(j); nbt != NULL;
137 >         nbt = nbiTypes->nextType(j)) {
138        
139        if (nbt->isLennardJones()) {
140          
141 <        std::pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
141 >        pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
142          
143          GenericData* data = nbt->getPropertyByName("LennardJones");
144          if (data == NULL) {
# Line 227 | Line 179 | namespace OpenMD {
179    void LJ::addType(AtomType* atomType){
180      RealType sigma1 = getSigma(atomType);
181      RealType epsilon1 = getEpsilon(atomType);
182 <
182 >    
183      // add it to the map:
184      AtomTypeProperties atp = atomType->getATP();    
185  
186 <    std::pair<std::map<int,AtomType*>::iterator,bool> ret;    
187 <    ret = LJMap.insert( std::pair<int, AtomType*>(atp.ident, atomType) );
186 >    pair<map<int,AtomType*>::iterator,bool> ret;    
187 >    ret = LJMap.insert( pair<int, AtomType*>(atp.ident, atomType) );
188      if (ret.second == false) {
189        sprintf( painCave.errMsg,
190                 "LJ already had a previous entry with ident %d\n",
# Line 288 | Line 240 | namespace OpenMD {
240      }    
241    }
242  
243 <  void LJ::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
244 <                     RealType rij, RealType r2, RealType rcut, RealType sw,
293 <                     RealType vdwMult, RealType &vpair, RealType &pot,
294 <                     Vector3d &f1) {
295 <
243 >  void LJ::calcForce(InteractionData idat) {
244 >    
245      if (!initialized_) initialize();
246      
247      RealType ros;
# Line 302 | Line 251 | namespace OpenMD {
251      RealType myDeriv = 0.0;
252      RealType myDerivC = 0.0;
253  
254 <    std::pair<AtomType*, AtomType*> key = std::make_pair(at1, at2);
254 >    std::pair<AtomType*, AtomType*> key = std::make_pair(idat.atype1,
255 >                                                         idat.atype2);
256      LJInteractionData mixer = MixingMap[key];
257  
258      RealType sigmai = mixer.sigmai;
259      RealType epsilon = mixer.epsilon;
260      
261  
262 <    ros = rij * sigmai;
262 >    ros = idat.rij * sigmai;
263  
264      getLJfunc(ros, myPot, myDeriv);
265  
266      if (shiftedPot_) {
267 <      rcos = rcut * sigmai;
267 >      rcos = idat.rcut * sigmai;
268        getLJfunc(rcos, myPotC, myDerivC);
269        myDerivC = 0.0;
270      } else if (LJ::shiftedFrc_) {
271 <      rcos = rcut * sigmai;
271 >      rcos = idat.rcut * sigmai;
272        getLJfunc(rcos, myPotC, myDerivC);
273 <      myPotC = myPotC + myDerivC * (rij - rcut) * sigmai;
273 >      myPotC = myPotC + myDerivC * (idat.rij - idat.rcut) * sigmai;
274      } else {
275        myPotC = 0.0;
276        myDerivC = 0.0;
277      }
278  
279 <    RealType pot_temp = vdwMult * epsilon * (myPot - myPotC);
280 <    vpair += pot_temp;
279 >    RealType pot_temp = idat.vdwMult * epsilon * (myPot - myPotC);
280 >    idat.vpair += pot_temp;
281  
282 <    RealType dudr = sw * vdwMult * epsilon * (myDeriv - myDerivC)*sigmai;
282 >    RealType dudr = idat.sw * idat.vdwMult * epsilon * (myDeriv -
283 >                                                        myDerivC)*sigmai;
284      
285 <    pot += sw * pot_temp;
286 <    f1 = d * dudr / rij;
285 >    idat.pot += idat.sw * pot_temp;
286 >    idat.f1 = idat.d * dudr / idat.rij;
287  
288      return;
289  
339 –
340 –
290    }
342 –
343 –  void LJ::do_lj_pair(int *atid1, int *atid2, RealType *d, RealType *rij,
344 –                      RealType *r2, RealType *rcut, RealType *sw,
345 –                      RealType *vdwMult,
346 –                      RealType *vpair, RealType *pot, RealType *f1) {
347 –
348 –    if (!initialized_) initialize();
349 –    
350 –    AtomType* atype1 = LJMap[*atid1];
351 –    AtomType* atype2 = LJMap[*atid2];
352 –    
353 –    Vector3d disp(d[0], d[1], d[2]);
354 –    Vector3d frc(f1[0], f1[1], f1[2]);
355 –    
356 –    calcForce(atype1, atype2, disp, *rij, *r2, *rcut, *sw, *vdwMult, *vpair,
357 –              *pot, frc);
358 –      
359 –    f1[0] = frc.x();
360 –    f1[1] = frc.y();
361 –    f1[2] = frc.z();
362 –
363 –    return;    
364 –  }
291    
292    void LJ::getLJfunc(RealType r, RealType &pot, RealType &deriv) {
293  
# Line 379 | Line 305 | namespace OpenMD {
305    }
306    
307  
382 –  void LJ::setLJDefaultCutoff(RealType *thisRcut, int *sP, int *sF) {
383 –    shiftedPot_ = (bool)(*sP);
384 –    shiftedFrc_ = (bool)(*sF);
385 –  }
308   }
387 –
388 – extern "C" {
389 –  
390 – #define fortranGetSigma FC_FUNC(getsigma, GETSIGMA)
391 – #define fortranGetEpsilon FC_FUNC(getepsilon, GETEPSILON)
392 – #define fortranSetLJCutoff FC_FUNC(setljdefaultcutoff, SETLJDEFAULTCUTOFF)
393 – #define fortranDoLJPair FC_FUNC(do_lj_pair, DO_LJ_PAIR)
394 –  
395 –  RealType fortranGetSigma(int* atid) {
396 –    return OpenMD::LJ::Instance()->getSigma(*atid);
397 –  }
398 –  RealType fortranGetEpsilon(int* atid) {  
399 –    return OpenMD::LJ::Instance()->getEpsilon(*atid);
400 –  }
401 –  void fortranSetLJCutoff(RealType *rcut, int *shiftedPot, int *shiftedFrc) {
402 –    return OpenMD::LJ::Instance()->setLJDefaultCutoff(rcut, shiftedPot,
403 –                                                      shiftedFrc);
404 –  }
405 –  void fortranDoLJPair(int *atid1, int *atid2, RealType *d, RealType *rij,
406 –                       RealType *r2, RealType *rcut, RealType *sw,
407 –                       RealType *vdwMult, RealType* vpair, RealType* pot,
408 –                       RealType *f1){
409 –    
410 –    return OpenMD::LJ::Instance()->do_lj_pair(atid1, atid2, d, rij, r2, rcut,
411 –                                              sw, vdwMult, vpair, pot, f1);
412 –  }
413 – }

Diff Legend

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