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

Comparing branches/development/src/nonbonded/Morse.cpp (file contents):
Revision 1502 by gezelter, Sat Oct 2 19:53:32 2010 UTC vs.
Revision 1583 by gezelter, Thu Jun 16 22:00:08 2011 UTC

# Line 51 | Line 51 | namespace OpenMD {
51  
52   namespace OpenMD {
53  
54 <  Morse::Morse() : name_("Morse"), initialized_(false), forceField_(NULL),
55 <                   shiftedPot_(false), shiftedFrc_(false) {}
54 >  Morse::Morse() : name_("Morse"), initialized_(false), forceField_(NULL) {}
55    
56    void Morse::initialize() {    
57  
# Line 142 | Line 141 | namespace OpenMD {
141      }    
142    }
143    
144 <  void Morse::calcForce(InteractionData idat) {
144 >  void Morse::calcForce(InteractionData &idat) {
145  
146      if (!initialized_) initialize();
147      
148 <    RealType myPot = 0.0;
149 <    RealType myPotC = 0.0;
150 <    RealType myDeriv = 0.0;
151 <    RealType myDerivC = 0.0;
152 <
153 <    pair<AtomType*, AtomType*> key = make_pair(idat.atype1, idat.atype2);
154 <    MorseInteractionData mixer = MixingMap[key];
155 <
156 <    RealType De = mixer.De;
157 <    RealType Re = mixer.Re;
158 <    RealType beta = mixer.beta;
159 <    MorseInteractionType interactionType = mixer.interactionType;
160 <  
161 <    // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
162 <    
163 <    RealType expt     = -beta*(idat.rij - Re);
164 <    RealType expfnc   = exp(expt);
165 <    RealType expfnc2  = expfnc*expfnc;
166 <
167 <    RealType exptC = 0.0;
168 <    RealType expfncC = 0.0;
169 <    RealType expfnc2C = 0.0;
170 <
171 <    if (Morse::shiftedPot_ || Morse::shiftedFrc_) {
172 <      exptC     = -beta*(idat.rcut - Re);
173 <      expfncC   = exp(exptC);
174 <      expfnc2C  = expfncC*expfncC;
175 <    }
176 <
178 <
179 <    switch(interactionType) {
180 <    case shiftedMorse : {
181 <
182 <      myPot  = De * (expfnc2  - 2.0 * expfnc);
183 <      myDeriv   = 2.0 * De * beta * (expfnc - expfnc2);
184 <
185 <      if (Morse::shiftedPot_) {
186 <        myPotC = De * (expfnc2C - 2.0 * expfncC);
187 <        myDerivC = 0.0;
188 <      } else if (Morse::shiftedFrc_) {
189 <        myPotC = De * (expfnc2C - 2.0 * expfncC);
190 <        myDerivC  = 2.0 * De * beta * (expfnc2C - expfnc2C);
191 <        myPotC += myDerivC * (idat.rij - idat.rcut);
192 <      } else {
193 <        myPotC = 0.0;
194 <        myDerivC = 0.0;
148 >    map<pair<AtomType*, AtomType*>, MorseInteractionData>::iterator it;
149 >    it = MixingMap.find( idat.atypes );
150 >    if (it != MixingMap.end()) {
151 >      MorseInteractionData mixer = (*it).second;
152 >      
153 >      RealType myPot = 0.0;
154 >      RealType myPotC = 0.0;
155 >      RealType myDeriv = 0.0;
156 >      RealType myDerivC = 0.0;
157 >      
158 >      RealType De = mixer.De;
159 >      RealType Re = mixer.Re;
160 >      RealType beta = mixer.beta;
161 >      MorseInteractionType interactionType = mixer.interactionType;
162 >      
163 >      // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
164 >      
165 >      RealType expt     = -beta*( *(idat.rij) - Re);
166 >      RealType expfnc   = exp(expt);
167 >      RealType expfnc2  = expfnc*expfnc;
168 >      
169 >      RealType exptC = 0.0;
170 >      RealType expfncC = 0.0;
171 >      RealType expfnc2C = 0.0;
172 >      
173 >      if (idat.shiftedPot || idat.shiftedForce) {
174 >        exptC     = -beta*( *(idat.rcut) - Re);
175 >        expfncC   = exp(exptC);
176 >        expfnc2C  = expfncC*expfncC;
177        }
178 <
179 <      break;
180 <    }
181 <    case repulsiveMorse : {
182 <
183 <      myPot  = De * expfnc2;
184 <      myDeriv  = -2.0 * De * beta * expfnc2;
185 <
186 <      if (Morse::shiftedPot_) {
187 <        myPotC = De * expfnc2C;
188 <        myDerivC = 0.0;
189 <      } else if (Morse::shiftedFrc_) {
190 <        myPotC = De * expfnc2C;
191 <        myDerivC = -2.0 * De * beta * expfnc2C;
192 <        myPotC += myDerivC * (idat.rij - idat.rcut);
193 <      } else {
194 <        myPotC = 0.0;
195 <        myDerivC = 0.0;
178 >      
179 >      
180 >      switch(interactionType) {
181 >      case shiftedMorse : {
182 >        
183 >        myPot  = De * (expfnc2  - 2.0 * expfnc);
184 >        myDeriv   = 2.0 * De * beta * (expfnc - expfnc2);
185 >        
186 >        if (idat.shiftedPot) {
187 >          myPotC = De * (expfnc2C - 2.0 * expfncC);
188 >          myDerivC = 0.0;
189 >        } else if (idat.shiftedForce) {
190 >          myPotC = De * (expfnc2C - 2.0 * expfncC);
191 >          myDerivC  = 2.0 * De * beta * (expfnc2C - expfnc2C);
192 >          myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut) );
193 >        } else {
194 >          myPotC = 0.0;
195 >          myDerivC = 0.0;
196 >        }
197 >        
198 >        break;
199        }
200 <
201 <      break;
200 >      case repulsiveMorse : {
201 >        
202 >        myPot  = De * expfnc2;
203 >        myDeriv  = -2.0 * De * beta * expfnc2;
204 >        
205 >        if (idat.shiftedPot) {
206 >          myPotC = De * expfnc2C;
207 >          myDerivC = 0.0;
208 >        } else if (idat.shiftedForce) {
209 >          myPotC = De * expfnc2C;
210 >          myDerivC = -2.0 * De * beta * expfnc2C;
211 >          myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut));
212 >        } else {
213 >          myPotC = 0.0;
214 >          myDerivC = 0.0;
215 >        }
216 >        
217 >        break;
218 >      }
219 >      }
220 >      
221 >      RealType pot_temp = *(idat.vdwMult) * (myPot - myPotC);
222 >      *(idat.vpair) += pot_temp;
223 >      
224 >      RealType dudr = *(idat.sw) * *(idat.vdwMult) * (myDeriv - myDerivC);
225 >      
226 >      (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
227 >      *(idat.f1) = *(idat.d) * dudr / *(idat.rij);
228      }
218 –    }
219 –
220 –    RealType pot_temp = idat.vdwMult * (myPot - myPotC);
221 –    idat.vpair += pot_temp;
222 –
223 –    RealType dudr = idat.sw * idat.vdwMult * (myDeriv - myDerivC);
224 –    
225 –    idat.pot += idat.sw * pot_temp;
226 –    idat.f1 = idat.d * dudr / idat.rij;
227 –
229      return;
230 <
230 >    
231    }
232      
233 +  RealType Morse::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
234 +    if (!initialized_) initialize();  
235 +    map<pair<AtomType*, AtomType*>, MorseInteractionData>::iterator it;
236 +    it = MixingMap.find(atypes);
237 +    if (it == MixingMap.end())
238 +      return 0.0;
239 +    else  {
240 +      MorseInteractionData mixer = (*it).second;
241 +
242 +      RealType Re = mixer.Re;
243 +      RealType beta = mixer.beta;
244 +      // This value of the r corresponds to an energy about 1.48% of
245 +      // the energy at the bottom of the Morse well.  For comparison, the
246 +      // Lennard-Jones function is about 1.63% of it's minimum value at
247 +      // a distance of 2.5 sigma.
248 +      return (4.9 + beta * Re) / beta;
249 +    }
250 +  }
251   }
252  

Diff Legend

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