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 1505 by gezelter, Sun Oct 3 22:18: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  
52   namespace OpenMD {
53  
54 <  LJ::LJ() : name_("LJ"), initialized_(false), shiftedPot_(false),
52 <             shiftedFrc_(false), forceField_(NULL) {}
54 >  LJ::LJ() : name_("LJ"), initialized_(false), forceField_(NULL) {}
55  
56 <  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 <  }
56 >  RealType LJ::getSigma(AtomType* atomType1, AtomType* atomType2) {
57  
58 <  RealType LJ::getSigma(AtomType* atomType) {    
59 <    LJParam ljParam = getLJParam(atomType);
60 <    return ljParam.sigma;
61 <  }
94 <
95 <  RealType LJ::getSigma(AtomType* atomType1, AtomType* atomType2) {    
96 <    RealType sigma1 = getSigma(atomType1);
97 <    RealType sigma2 = getSigma(atomType2);
58 >    LennardJonesAdapter lja1 = LennardJonesAdapter(atomType1);
59 >    LennardJonesAdapter lja2 = LennardJonesAdapter(atomType2);
60 >    RealType sigma1 = lja1.getSigma();
61 >    RealType sigma2 = lja2.getSigma();
62      
63      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
64      string DistanceMix = fopts.getDistanceMixingRule();
# Line 106 | 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(AtomType* atomType1, AtomType* atomType2) {    
115 <    RealType epsilon1 = getEpsilon(atomType1);
116 <    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  
# Line 124 | Line 86 | namespace OpenMD {
86  
87      for (at = atomTypes->beginType(i); at != NULL;
88           at = atomTypes->nextType(i)) {
89 <      
90 <      if (at->isLennardJones())
91 <        addType(at);
89 >      LennardJonesAdapter lja = LennardJonesAdapter(at);
90 >      if (lja.isLennardJones()){
91 >         addType(at);
92 >      }
93      }
131 –
94      ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
95      ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
96      NonBondedInteractionType* nbt;
97 +    ForceField::NonBondedInteractionTypeContainer::KeyType keys;
98  
99 +
100      for (nbt = nbiTypes->beginType(j); nbt != NULL;
101           nbt = nbiTypes->nextType(j)) {
102 <      
103 <      if (nbt->isLennardJones()) {
104 <        
105 <        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",
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) {
102 >
103 >      if (nbt->isLennardJones()) {
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          }
165 –        
166 –        LJParam ljParam = ljData->getData();
120  
121 <        RealType sigma = ljParam.sigma;
122 <        RealType epsilon = ljParam.epsilon;
123 <
171 <        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 177 | Line 129 | namespace OpenMD {
129  
130  
131    void LJ::addType(AtomType* atomType){
132 <    RealType sigma1 = getSigma(atomType);
133 <    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:
184 –    AtomTypeProperties atp = atomType->getATP();    
138  
139      pair<map<int,AtomType*>::iterator,bool> ret;    
140 <    ret = LJMap.insert( pair<int, AtomType*>(atp.ident, atomType) );
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 220 | Line 173 | namespace OpenMD {
173    
174    void LJ::addExplicitInteraction(AtomType* atype1, AtomType* atype2, RealType sigma, RealType epsilon){
175      
223 –    // in case these weren't already in the map
224 –    addType(atype1);
225 –    addType(atype2);
226 –
176      LJInteractionData mixer;
177      mixer.sigma = sigma;
178      mixer.epsilon = epsilon;
# Line 240 | Line 189 | namespace OpenMD {
189      }    
190    }
191  
192 <  void LJ::calcForce(InteractionData idat) {
244 <    
192 >  void LJ::calcForce(InteractionData &idat) {
193      if (!initialized_) initialize();
246 –    
247 –    pair<AtomType*, AtomType*> key = make_pair(idat.atype1, idat.atype2);
194      map<pair<AtomType*, AtomType*>, LJInteractionData>::iterator it;
195 <    it = MixingMap.find(key);
195 >    it = MixingMap.find( idat.atypes );
196      
197      if (it != MixingMap.end())  {
198        
# Line 262 | Line 208 | namespace OpenMD {
208        RealType myDeriv = 0.0;
209        RealType myDerivC = 0.0;
210      
211 <      ros = idat.rij * sigmai;
211 >      ros = *(idat.rij) * sigmai;    
212        
213        getLJfunc(ros, myPot, myDeriv);
214        
215 <      if (shiftedPot_) {
216 <        rcos = idat.rcut * sigmai;
215 >      if (idat.shiftedPot) {
216 >        rcos = *(idat.rcut) * sigmai;
217          getLJfunc(rcos, myPotC, myDerivC);
218          myDerivC = 0.0;
219 <      } else if (LJ::shiftedFrc_) {
220 <        rcos = idat.rcut * sigmai;
219 >      } else if (idat.shiftedForce) {
220 >        rcos = *(idat.rcut) * sigmai;
221          getLJfunc(rcos, myPotC, myDerivC);
222 <        myPotC = myPotC + myDerivC * (idat.rij - idat.rcut) * sigmai;
222 >        myPotC = myPotC + myDerivC * (*(idat.rij) - *(idat.rcut)) * sigmai;
223        } else {
224          myPotC = 0.0;
225          myDerivC = 0.0;        
226        }
227  
228 <      RealType pot_temp = idat.vdwMult * epsilon * (myPot - myPotC);
229 <      idat.vpair += pot_temp;
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 <      idat.pot += idat.sw * pot_temp;
235 <      idat.f1 = idat.d * dudr / idat.rij;
290 <      
231 >      RealType dudr = *(idat.sw) * *(idat.vdwMult) * epsilon * (myDeriv -
232 >                                                                myDerivC)*sigmai;      
233 >
234 >      (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
235 >      *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
236      }
237      return;
238    }
# Line 307 | Line 252 | namespace OpenMD {
252      return;
253    }
254    
255 <  RealType LJ::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
255 >  RealType LJ::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
256      if (!initialized_) initialize();  
312 –    pair<AtomType*, AtomType*> key = make_pair(at1, at2);
257      map<pair<AtomType*, AtomType*>, LJInteractionData>::iterator it;
258 <    it = MixingMap.find(key);
258 >    it = MixingMap.find(atypes);
259      if (it == MixingMap.end())
260        return 0.0;
261      else  {

Diff Legend

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