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 1501 by gezelter, Wed Sep 15 19:32:10 2010 UTC vs.
Revision 1505 by gezelter, Sun Oct 3 22:18:59 2010 UTC

# Line 51 | Line 51 | namespace OpenMD {
51  
52   namespace OpenMD {
53  
54 <  bool Morse::initialized_ = false;
55 <  bool Morse::shiftedPot_ = false;
56 <  bool Morse::shiftedFrc_ = false;
57 <  ForceField* Morse::forceField_ = NULL;
58 <  map<int, AtomType*> Morse::MorseMap;
59 <  map<pair<AtomType*, AtomType*>, MorseInteractionData> Morse::MixingMap;
60 <  map<string, MorseInteractionType> Morse::stringToEnumMap_;
54 >  Morse::Morse() : name_("Morse"), initialized_(false), forceField_(NULL),
55 >                   shiftedPot_(false), shiftedFrc_(false) {}
56    
62 –  Morse* Morse::_instance = NULL;
63 –
64 –  Morse* Morse::Instance() {
65 –    if (!_instance) {
66 –      _instance = new Morse();
67 –    }
68 –    return _instance;
69 –  }
70 –
57    void Morse::initialize() {    
58  
59      stringToEnumMap_["shiftedMorse"] = shiftedMorse;
# Line 140 | Line 126 | namespace OpenMD {
126                                       RealType De, RealType Re, RealType beta,
127                                       MorseInteractionType mit) {
128  
143 –    AtomTypeProperties atp = atype1->getATP();    
144 –    MorseMap.insert( pair<int, AtomType*>(atp.ident, atype1) );
145 –
146 –    atp = atype2->getATP();    
147 –    MorseMap.insert( pair<int, AtomType*>(atp.ident, atype2) );
148 –
129      MorseInteractionData mixer;
130      mixer.De = De;
131      mixer.Re = Re;
# Line 162 | Line 142 | namespace OpenMD {
142      }    
143    }
144    
145 <  void Morse::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
166 <                        RealType r, RealType r2, RealType rcut, RealType sw,
167 <                        RealType vdwMult, RealType &vpair, RealType &pot,
168 <                        Vector3d &f1) {
145 >  void Morse::calcForce(InteractionData idat) {
146  
147      if (!initialized_) initialize();
148      
149 <    RealType myPot = 0.0;
150 <    RealType myPotC = 0.0;
151 <    RealType myDeriv = 0.0;
152 <    RealType myDerivC = 0.0;
153 <
154 <    pair<AtomType*, AtomType*> key = make_pair(at1, at2);
155 <    MorseInteractionData mixer = MixingMap[key];
156 <
157 <    RealType De = mixer.De;
158 <    RealType Re = mixer.Re;
159 <    RealType beta = mixer.beta;
160 <    MorseInteractionType interactionType = mixer.interactionType;
161 <  
162 <    // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
163 <    
164 <    RealType expt     = -beta*(r - Re);
165 <    RealType expfnc   = exp(expt);
166 <    RealType expfnc2  = expfnc*expfnc;
167 <
168 <    RealType exptC = 0.0;
169 <    RealType expfncC = 0.0;
170 <    RealType expfnc2C = 0.0;
171 <
172 <    if (Morse::shiftedPot_ || Morse::shiftedFrc_) {
173 <      exptC     = -beta*(rcut - Re);
174 <      expfncC   = exp(exptC);
175 <      expfnc2C  = expfncC*expfncC;
176 <    }
177 <
178 <
202 <    switch(interactionType) {
203 <    case shiftedMorse : {
204 <
205 <      myPot  = De * (expfnc2  - 2.0 * expfnc);
206 <      myDeriv   = 2.0 * De * beta * (expfnc - expfnc2);
207 <
208 <      if (Morse::shiftedPot_) {
209 <        myPotC = De * (expfnc2C - 2.0 * expfncC);
210 <        myDerivC = 0.0;
211 <      } else if (Morse::shiftedFrc_) {
212 <        myPotC = De * (expfnc2C - 2.0 * expfncC);
213 <        myDerivC  = 2.0 * De * beta * (expfnc2C - expfnc2C);
214 <        myPotC += myDerivC * (r - rcut);
215 <      } else {
216 <        myPotC = 0.0;
217 <        myDerivC = 0.0;
149 >    pair<AtomType*, AtomType*> key = make_pair(idat.atype1, idat.atype2);
150 >    map<pair<AtomType*, AtomType*>, MorseInteractionData>::iterator it;
151 >    it = MixingMap.find(key);
152 >    if (it != MixingMap.end()) {
153 >      MorseInteractionData mixer = (*it).second;
154 >      
155 >      RealType myPot = 0.0;
156 >      RealType myPotC = 0.0;
157 >      RealType myDeriv = 0.0;
158 >      RealType myDerivC = 0.0;
159 >      
160 >      RealType De = mixer.De;
161 >      RealType Re = mixer.Re;
162 >      RealType beta = mixer.beta;
163 >      MorseInteractionType interactionType = mixer.interactionType;
164 >      
165 >      // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
166 >      
167 >      RealType expt     = -beta*(idat.rij - Re);
168 >      RealType expfnc   = exp(expt);
169 >      RealType expfnc2  = expfnc*expfnc;
170 >      
171 >      RealType exptC = 0.0;
172 >      RealType expfncC = 0.0;
173 >      RealType expfnc2C = 0.0;
174 >      
175 >      if (Morse::shiftedPot_ || Morse::shiftedFrc_) {
176 >        exptC     = -beta*(idat.rcut - Re);
177 >        expfncC   = exp(exptC);
178 >        expfnc2C  = expfncC*expfncC;
179        }
180 <
181 <      break;
182 <    }
183 <    case repulsiveMorse : {
184 <
185 <      myPot  = De * expfnc2;
186 <      myDeriv  = -2.0 * De * beta * expfnc2;
187 <
188 <      if (Morse::shiftedPot_) {
189 <        myPotC = De * expfnc2C;
190 <        myDerivC = 0.0;
191 <      } else if (Morse::shiftedFrc_) {
192 <        myPotC = De * expfnc2C;
193 <        myDerivC = -2.0 * De * beta * expfnc2C;
194 <        myPotC += myDerivC * (r - rcut);
195 <      } else {
196 <        myPotC = 0.0;
197 <        myDerivC = 0.0;
180 >      
181 >      
182 >      switch(interactionType) {
183 >      case shiftedMorse : {
184 >        
185 >        myPot  = De * (expfnc2  - 2.0 * expfnc);
186 >        myDeriv   = 2.0 * De * beta * (expfnc - expfnc2);
187 >        
188 >        if (Morse::shiftedPot_) {
189 >          myPotC = De * (expfnc2C - 2.0 * expfncC);
190 >          myDerivC = 0.0;
191 >        } else if (Morse::shiftedFrc_) {
192 >          myPotC = De * (expfnc2C - 2.0 * expfncC);
193 >          myDerivC  = 2.0 * De * beta * (expfnc2C - expfnc2C);
194 >          myPotC += myDerivC * (idat.rij - idat.rcut);
195 >        } else {
196 >          myPotC = 0.0;
197 >          myDerivC = 0.0;
198 >        }
199 >        
200 >        break;
201        }
202 <
203 <      break;
202 >      case repulsiveMorse : {
203 >        
204 >        myPot  = De * expfnc2;
205 >        myDeriv  = -2.0 * De * beta * expfnc2;
206 >        
207 >        if (Morse::shiftedPot_) {
208 >          myPotC = De * expfnc2C;
209 >          myDerivC = 0.0;
210 >        } else if (Morse::shiftedFrc_) {
211 >          myPotC = De * expfnc2C;
212 >          myDerivC = -2.0 * De * beta * expfnc2C;
213 >          myPotC += myDerivC * (idat.rij - idat.rcut);
214 >        } else {
215 >          myPotC = 0.0;
216 >          myDerivC = 0.0;
217 >        }
218 >        
219 >        break;
220 >      }
221 >      }
222 >      
223 >      RealType pot_temp = idat.vdwMult * (myPot - myPotC);
224 >      idat.vpair += pot_temp;
225 >      
226 >      RealType dudr = idat.sw * idat.vdwMult * (myDeriv - myDerivC);
227 >      
228 >      idat.pot += idat.sw * pot_temp;
229 >      idat.f1 = idat.d * dudr / idat.rij;
230      }
241 –    }
242 –
243 –    RealType pot_temp = vdwMult * (myPot - myPotC);
244 –    vpair += pot_temp;
245 –
246 –    RealType dudr = sw * vdwMult * (myDeriv - myDerivC);
247 –    
248 –    pot += sw * pot_temp;
249 –    f1 = d * dudr / r;
250 –
231      return;
252 –
253 –  }
254 –
255 –  void Morse::do_morse_pair(int *atid1, int *atid2, RealType *d,
256 –                            RealType *rij, RealType *r2, RealType *rcut,
257 –                            RealType *sw, RealType *vdwMult,
258 –                            RealType *vpair, RealType *pot, RealType *f1) {
259 –
260 –    if (!initialized_) initialize();
232      
262 –    AtomType* atype1 = MorseMap[*atid1];
263 –    AtomType* atype2 = MorseMap[*atid2];
264 –    
265 –    Vector3d disp(d[0], d[1], d[2]);
266 –    Vector3d frc(f1[0], f1[1], f1[2]);
267 –    
268 –    calcForce(atype1, atype2, disp, *rij, *r2, *rcut, *sw, *vdwMult, *vpair,
269 –              *pot, frc);
270 –      
271 –    f1[0] = frc.x();
272 –    f1[1] = frc.y();
273 –    f1[2] = frc.z();
274 –
275 –    return;    
233    }
234      
235 <  void Morse::setMorseDefaultCutoff(RealType *thisRcut, int *sP, int *sF) {
236 <    shiftedPot_ = (bool)(*sP);
237 <    shiftedFrc_ = (bool)(*sF);
238 <  }
239 < }
235 >  RealType Morse::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
236 >    if (!initialized_) initialize();  
237 >    pair<AtomType*, AtomType*> key = make_pair(at1, at2);
238 >    map<pair<AtomType*, AtomType*>, MorseInteractionData>::iterator it;
239 >    it = MixingMap.find(key);
240 >    if (it == MixingMap.end())
241 >      return 0.0;
242 >    else  {
243 >      MorseInteractionData mixer = (*it).second;
244  
245 < extern "C" {
246 <  
247 < #define fortranSetMorseCutoff FC_FUNC(setmorsedefaultcutoff, SETMORSEDEFAULTCUTOFF)
248 < #define fortranDoMorsePair FC_FUNC(do_morse_pair, DO_MORSE_PAIR)
249 <  
250 <  void fortranSetMorseCutoff(RealType *rcut, int *shiftedPot, int *shiftedFrc) {
251 <    return OpenMD::Morse::Instance()->setMorseDefaultCutoff(rcut, shiftedPot,
252 <                                                            shiftedFrc);
245 >      RealType Re = mixer.Re;
246 >      RealType beta = mixer.beta;
247 >      // This value of the r corresponds to an energy about 1.48% of
248 >      // the energy at the bottom of the Morse well.  For comparison, the
249 >      // Lennard-Jones function is about 1.63% of it's minimum value at
250 >      // a distance of 2.5 sigma.
251 >      return (4.9 + beta * Re) / beta;
252 >    }
253    }
293 –  void fortranDoMorsePair(int *atid1, int *atid2, RealType *d, RealType *rij,
294 –                          RealType *r2, RealType *rcut, RealType *sw,
295 –                          RealType *vdwMult, RealType* vpair, RealType* pot,
296 –                          RealType *f1){
297 –    
298 –    return OpenMD::Morse::Instance()->do_morse_pair(atid1, atid2, d, rij, r2,
299 –                                                    rcut, sw, vdwMult, vpair,
300 –                                                    pot, f1);
301 –  }
254   }
255 +

Diff Legend

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