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 1471 by gezelter, Mon Jul 19 18:59:59 2010 UTC vs.
Revision 1710 by gezelter, Fri May 18 21:44:02 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/LennardJonesAdapter.hpp"
50 + #include "types/LennardJonesInteractionType.hpp"
51  
49
52   namespace OpenMD {
53  
54 <  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;
54 >  LJ::LJ() : name_("LJ"), initialized_(false), forceField_(NULL) {}
55  
56 <  LJ* LJ::Instance() {
62 <    if (!_instance) {
63 <      _instance = new LJ();
64 <    }
65 <    return _instance;
66 <  }
56 >  RealType LJ::getSigma(AtomType* atomType1, AtomType* atomType2) {
57  
58 <  LJParam LJ::getLJParam(AtomType* atomType) {
58 >    LennardJonesAdapter lja1 = LennardJonesAdapter(atomType1);
59 >    LennardJonesAdapter lja2 = LennardJonesAdapter(atomType2);
60 >    RealType sigma1 = lja1.getSigma();
61 >    RealType sigma2 = lja2.getSigma();
62      
70    // Do sanity checking on the AtomType we were passed before
71    // building any data structures:
72    if (!atomType->isLennardJones()) {
73      sprintf( painCave.errMsg,
74               "LJ::getLJParam was passed an atomType (%s) that does not\n"
75               "\tappear to be a Lennard-Jones atom.\n",
76               atomType->getName().c_str());
77      painCave.severity = OPENMD_ERROR;
78      painCave.isFatal = 1;
79      simError();
80    }
81    
82    GenericData* data = atomType->getPropertyByName("LennardJones");
83    if (data == NULL) {
84      sprintf( painCave.errMsg, "LJ::getLJParam could not find Lennard-Jones\n"
85               "\tparameters for atomType %s.\n", atomType->getName().c_str());
86      painCave.severity = OPENMD_ERROR;
87      painCave.isFatal = 1;
88      simError();
89    }
90    
91    LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
92    if (ljData == NULL) {
93      sprintf( painCave.errMsg,
94               "LJ::getLJParam could not convert GenericData to LJParam for\n"
95               "\tatom type %s\n", atomType->getName().c_str());
96      painCave.severity = OPENMD_ERROR;
97      painCave.isFatal = 1;
98      simError();          
99    }
100    
101    return ljData->getData();
102  }
103
104  RealType LJ::getSigma(AtomType* atomType) {    
105    LJParam ljParam = getLJParam(atomType);
106    return ljParam.sigma;
107  }
108
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
125  RealType LJ::getSigma(AtomType* atomType1, AtomType* atomType2) {    
126    RealType sigma1 = getSigma(atomType1);
127    RealType sigma2 = getSigma(atomType2);
128    
63      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
64 <    std::string DistanceMix = fopts.getDistanceMixingRule();
64 >    string DistanceMix = fopts.getDistanceMixingRule();
65      toUpper(DistanceMix);
66  
67      if (DistanceMix == "GEOMETRIC")
# Line 136 | Line 70 | namespace OpenMD {
70        return 0.5 * (sigma1 + sigma2);
71    }
72  
73 <  RealType LJ::getEpsilon(AtomType* atomType) {    
74 <    LJParam ljParam = getLJParam(atomType);
75 <    return ljParam.epsilon;
76 <  }
77 <
78 <  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 <
160 <  RealType LJ::getEpsilon(AtomType* atomType1, AtomType* atomType2) {    
161 <    RealType epsilon1 = getEpsilon(atomType1);
162 <    RealType epsilon2 = getEpsilon(atomType2);
73 >  RealType LJ::getEpsilon(AtomType* atomType1, AtomType* atomType2) {  
74 >    LennardJonesAdapter lja1 = LennardJonesAdapter(atomType1);
75 >    LennardJonesAdapter lja2 = LennardJonesAdapter(atomType2);
76 >  
77 >    RealType epsilon1 = lja1.getEpsilon();
78 >    RealType epsilon2 = lja2.getEpsilon();
79      return sqrt(epsilon1 * epsilon2);
80    }
81  
82    void LJ::initialize() {    
83 <    ForceField::AtomTypeContainer atomTypes = forceField_->getAtomTypes();
83 >    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
84      ForceField::AtomTypeContainer::MapTypeIterator i;
85      AtomType* at;
86  
87 <    for (at = atomTypes.beginType(i); at != NULL;
88 <         at = atomTypes.nextType(i)) {
89 <      
90 <      if (at->isLennardJones())
91 <        addType(at);
87 >    for (at = atomTypes->beginType(i); at != NULL;
88 >         at = atomTypes->nextType(i)) {
89 >      LennardJonesAdapter lja = LennardJonesAdapter(at);
90 >      if (lja.isLennardJones()){
91 >         addType(at);
92 >      }
93      }
94 <
178 <    ForceField::NonBondedInteractionTypeContainer nbiTypes = forceField_->getNonBondedInteractionTypes();
94 >    ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
95      ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
96      NonBondedInteractionType* nbt;
97 +    ForceField::NonBondedInteractionTypeContainer::KeyType keys;
98  
99 <    for (nbt = nbiTypes.beginType(j); nbt != NULL;
100 <         nbt = nbiTypes.nextType(j)) {
101 <      
99 >
100 >    for (nbt = nbiTypes->beginType(j); nbt != NULL;
101 >         nbt = nbiTypes->nextType(j)) {
102 >
103        if (nbt->isLennardJones()) {
104 <        
105 <        std::pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
106 <        
107 <        GenericData* data = nbt->getPropertyByName("LennardJones");
108 <        if (data == NULL) {
109 <          sprintf( painCave.errMsg, "LJ::rebuildMixingMap could not find\n"
110 <               "\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) {
104 >        keys = nbiTypes->getKeys(j);
105 >        AtomType* at1 = forceField_->getAtomType(keys[0]);
106 >        AtomType* at2 = forceField_->getAtomType(keys[1]);
107 >
108 >        LennardJonesInteractionType* ljit = dynamic_cast<LennardJonesInteractionType*>(nbt);
109 >
110 >        if (ljit == NULL) {
111            sprintf( painCave.errMsg,
112 <                   "LJ::rebuildMixingMap could not convert GenericData to\n"
113 <                   "\tLJParam for %s - %s interaction.\n",
114 <                   atypes.first->getName().c_str(),
115 <                   atypes.second->getName().c_str());
112 >                   "LJ::initialize could not convert NonBondedInteractionType\n"
113 >                   "\tto LennardJonesInteractionType for %s - %s interaction.\n",
114 >                   at1->getName().c_str(),
115 >                   at2->getName().c_str());
116            painCave.severity = OPENMD_ERROR;
117            painCave.isFatal = 1;
118            simError();          
119          }
211        
212        LJParam ljParam = ljData->getData();
120  
121 <        RealType sigma = ljParam.sigma;
122 <        RealType epsilon = ljParam.epsilon;
123 <
217 <        addExplicitInteraction(atypes.first, atypes.second, sigma, epsilon);
121 >        RealType sigma = ljit->getSigma();
122 >        RealType epsilon = ljit->getEpsilon();
123 >        addExplicitInteraction(at1, at2, sigma, epsilon);
124        }
125      }  
126      initialized_ = true;
# Line 223 | Line 129 | namespace OpenMD {
129  
130  
131    void LJ::addType(AtomType* atomType){
132 <    RealType sigma1 = getSigma(atomType);
227 <    RealType epsilon1 = getEpsilon(atomType);
132 >    LennardJonesAdapter lja1 = LennardJonesAdapter(atomType);
133  
134 +    RealType sigma1 = lja1.getSigma();
135 +    RealType epsilon1 = lja1.getEpsilon();
136 +    
137      // add it to the map:
138 <    AtomTypeProperties atp = atomType->getATP();    
139 <    std::pair<std::map<int,AtomType*>::iterator,bool> ret;
140 <    ret = LJMap.insert( std::pair<int, AtomType*>(atp.ident, atomType) );
138 >
139 >    pair<map<int,AtomType*>::iterator,bool> ret;    
140 >    ret = LJMap.insert( pair<int, AtomType*>(atomType->getIdent(), atomType) );
141      if (ret.second == false) {
142        sprintf( painCave.errMsg,
143                 "LJ already had a previous entry with ident %d\n",
144 <               atp.ident);
144 >               atomType->getIdent());
145        painCave.severity = OPENMD_INFO;
146        painCave.isFatal = 0;
147        simError();        
# Line 265 | Line 173 | namespace OpenMD {
173    
174    void LJ::addExplicitInteraction(AtomType* atype1, AtomType* atype2, RealType sigma, RealType epsilon){
175      
268    // in case these weren't already in the map
269    addType(atype1);
270    addType(atype2);
271
176      LJInteractionData mixer;
177      mixer.sigma = sigma;
178      mixer.epsilon = epsilon;
# Line 285 | Line 189 | namespace OpenMD {
189      }    
190    }
191  
192 <  void LJ::calcForce(AtomType* at1, AtomType* at2, Vector3d d, RealType rij,
289 <                     RealType r2, RealType rcut, RealType sw, RealType vdwMult,
290 <                     RealType vpair, RealType pot, Vector3d f1) {
291 <
192 >  void LJ::calcForce(InteractionData &idat) {
193      if (!initialized_) initialize();
194 +    map<pair<AtomType*, AtomType*>, LJInteractionData>::iterator it;
195 +    it = MixingMap.find( idat.atypes );
196      
197 <    RealType ros;
198 <    RealType rcos;
199 <    RealType myPot = 0.0;
200 <    RealType myPotC = 0.0;
201 <    RealType myDeriv = 0.0;
202 <    RealType myDerivC = 0.0;
197 >    if (it != MixingMap.end())  {
198 >      
199 >      LJInteractionData mixer = (*it).second;
200 >      
201 >      RealType sigmai = mixer.sigmai;
202 >      RealType epsilon = mixer.epsilon;
203 >      
204 >      RealType ros;
205 >      RealType rcos;
206 >      RealType myPot = 0.0;
207 >      RealType myPotC = 0.0;
208 >      RealType myDeriv = 0.0;
209 >      RealType myDerivC = 0.0;
210 >    
211 >      ros = *(idat.rij) * sigmai;    
212 >      
213 >      getLJfunc(ros, myPot, myDeriv);
214 >      
215 >      if (idat.shiftedPot) {
216 >        rcos = *(idat.rcut) * sigmai;
217 >        getLJfunc(rcos, myPotC, myDerivC);
218 >        myDerivC = 0.0;
219 >      } else if (idat.shiftedForce) {
220 >        rcos = *(idat.rcut) * sigmai;
221 >        getLJfunc(rcos, myPotC, myDerivC);
222 >        myPotC = myPotC + myDerivC * (*(idat.rij) - *(idat.rcut)) * sigmai;
223 >      } else {
224 >        myPotC = 0.0;
225 >        myDerivC = 0.0;        
226 >      }
227  
228 <    std::pair<AtomType*, AtomType*> key = std::make_pair(at1, at2);
229 <    LJInteractionData mixer = MixingMap[key];
228 >      RealType pot_temp = *(idat.vdwMult) * epsilon * (myPot - myPotC);
229 >      *(idat.vpair) += pot_temp;
230 >      
231 >      RealType dudr = *(idat.sw) * *(idat.vdwMult) * epsilon * (myDeriv -
232 >                                                                myDerivC)*sigmai;      
233  
234 <    RealType sigmai = mixer.sigmai;
235 <    RealType epsilon = mixer.epsilon;
306 <    
307 <    ros = rij * sigmai;
308 <
309 <    getLJfunc(ros, myPot, myDeriv);
310 <
311 <    if (shiftedPot_) {
312 <      rcos = rcut * sigmai;
313 <      getLJfunc(rcos, myPotC, myDerivC);
314 <      myDerivC = 0.0;
315 <    } else if (LJ::shiftedFrc_) {
316 <      rcos = rcut * sigmai;
317 <      getLJfunc(rcos, myPotC, myDerivC);
318 <      myPotC = myPotC + myDerivC * (rij - rcut) * sigmai;
319 <    } else {
320 <      myPotC = 0.0;
321 <      myDerivC = 0.0;
234 >      (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
235 >      *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
236      }
323
324    RealType pot_temp = vdwMult * epsilon * (myPot - myPotC);
325    vpair += pot_temp;
326
327    RealType dudr = sw * vdwMult * epsilon * (myDeriv - myDerivC)*sigmai;
328    
329    pot += sw * pot_temp;
330    f1 = d * dudr / rij;
331
237      return;
238    }
334
335  void LJ::do_lj_pair(int *atid1, int *atid2, RealType *d, RealType *rij,
336                      RealType *r2, RealType *rcut, RealType *sw,
337                      RealType *vdwMult,
338                      RealType *vpair, RealType *pot, RealType *f1) {
339
340    if (!initialized_) initialize();
341    
342    AtomType* atype1 = LJMap[*atid1];
343    AtomType* atype2 = LJMap[*atid2];
344    
345    Vector3d disp(d[0], d[1], d[2]);
346    Vector3d frc(f1[0], f1[1], f1[2]);
347    
348    calcForce(atype1, atype2, disp, *rij, *r2, *rcut, *sw, *vdwMult, *vpair,
349              *pot, frc);
350    return;    
351  }
239    
240 <  void LJ::getLJfunc(const RealType r, RealType pot, RealType deriv) {
240 >  void LJ::getLJfunc(RealType r, RealType &pot, RealType &deriv) {
241 >
242      RealType ri = 1.0 / r;
243      RealType ri2 = ri * ri;
244      RealType ri6 = ri2 * ri2 * ri2;
245      RealType ri7 = ri6 * ri;
246      RealType ri12 = ri6 * ri6;
247      RealType ri13 = ri12 * ri;
248 <
248 >    
249      pot = 4.0 * (ri12 - ri6);
250      deriv = 24.0 * (ri7 - 2.0 * ri13);
251 +
252      return;
253    }
365
366
367  void LJ::setLJDefaultCutoff(RealType *thisRcut, int *sP, int *sF) {
368    shiftedPot_ = (bool)(*sP);
369    shiftedFrc_ = (bool)(*sF);
370  }
371 }
372
373 extern "C" {
254    
255 < #define fortranGetSigma FC_FUNC(getsigma, GETSIGMA)
256 < #define fortranGetEpsilon FC_FUNC(getepsilon, GETEPSILON)
257 < #define fortranSetLJCutoff FC_FUNC(setljdefaultcutoff, SETLJDEFAULTCUTOFF)
258 < #define fortranDoLJPair FC_FUNC(do_lj_pair, DO_LJ_PAIR)
259 <  
260 <  RealType fortranGetSigma(int* atid) {
261 <    return OpenMD::LJ::Instance()->getSigma(*atid);
255 >  RealType LJ::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
256 >    if (!initialized_) initialize();  
257 >    map<pair<AtomType*, AtomType*>, LJInteractionData>::iterator it;
258 >    it = MixingMap.find(atypes);
259 >    if (it == MixingMap.end())
260 >      return 0.0;
261 >    else  {
262 >      LJInteractionData mixer = (*it).second;
263 >      return 2.5 * mixer.sigma;
264 >    }
265    }
266 <  RealType fortranGetEpsilon(int* atid) {  
384 <    return OpenMD::LJ::Instance()->getEpsilon(*atid);
385 <  }
386 <  void fortranSetLJCutoff(RealType *rcut, int *shiftedPot, int *shiftedFrc) {
387 <    return OpenMD::LJ::Instance()->setLJDefaultCutoff(rcut, shiftedPot,
388 <                                                      shiftedFrc);
389 <  }
390 <  void fortranDoLJPair(int *atid1, int *atid2, RealType *d, RealType *rij,
391 <                       RealType *r2, RealType *rcut, RealType *sw,
392 <                       RealType *vdwMult, RealType* vpair, RealType* pot,
393 <                       RealType *f1){
394 <    
395 <    return OpenMD::LJ::Instance()->do_lj_pair(atid1, atid2, d, rij, r2, rcut,
396 <                                              sw, vdwMult, vpair, pot, f1);
397 <  }
266 >
267   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines