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 1850 by gezelter, Wed Feb 20 15:39:39 2013 UTC

# Line 35 | Line 35
35   *                                                                      
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).                        
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
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/Morse.hpp"
48   #include "utils/simError.h"
49 < #include "types/NonBondedInteractionType.hpp"
49 > #include "types/MorseInteractionType.hpp"
50  
51   using namespace std;
52  
53   namespace OpenMD {
54  
55 <  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_;
55 >  Morse::Morse() : name_("Morse"), initialized_(false), forceField_(NULL) {}
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  
73 –    stringToEnumMap_["shiftedMorse"] = shiftedMorse;
74 –    stringToEnumMap_["repulsiveMorse"] = repulsiveMorse;
75 –
59      ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
60      ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
61      NonBondedInteractionType* nbt;
62 +    ForceField::NonBondedInteractionTypeContainer::KeyType keys;
63  
64      for (nbt = nbiTypes->beginType(j); nbt != NULL;
65           nbt = nbiTypes->nextType(j)) {
66        
67        if (nbt->isMorse()) {
68 <        
69 <        pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
70 <        
87 <        GenericData* data = nbt->getPropertyByName("Morse");
88 <        if (data == NULL) {
89 <          sprintf( painCave.errMsg, "Morse::initialize could not find\n"
90 <                   "\tMorse parameters for %s - %s interaction.\n",
91 <                   atypes.first->getName().c_str(),
92 <                   atypes.second->getName().c_str());
93 <          painCave.severity = OPENMD_ERROR;
94 <          painCave.isFatal = 1;
95 <          simError();
96 <        }
97 <        
98 <        MorseData* morseData = dynamic_cast<MorseData*>(data);
99 <        if (morseData == NULL) {
100 <          sprintf( painCave.errMsg,
101 <                   "Morse::initialize could not convert GenericData to\n"
102 <                   "\tMorseData for %s - %s interaction.\n",
103 <                   atypes.first->getName().c_str(),
104 <                   atypes.second->getName().c_str());
105 <          painCave.severity = OPENMD_ERROR;
106 <          painCave.isFatal = 1;
107 <          simError();          
108 <        }
109 <        
110 <        MorseParam morseParam = morseData->getData();
68 >        keys = nbiTypes->getKeys(j);
69 >        AtomType* at1 = forceField_->getAtomType(keys[0]);
70 >        AtomType* at2 = forceField_->getAtomType(keys[1]);
71  
72 <        RealType De = morseParam.De;
113 <        RealType Re = morseParam.Re;
114 <        RealType beta = morseParam.beta;
115 <        string interactionType = morseParam.interactionType;
72 >        MorseInteractionType* mit = dynamic_cast<MorseInteractionType*>(nbt);
73  
74 <        toUpper(interactionType);
118 <        map<string, MorseInteractionType>::iterator i;
119 <        i = stringToEnumMap_.find(interactionType);
120 <        if (i != stringToEnumMap_.end()) {
121 <          addExplicitInteraction(atypes.first, atypes.second,
122 <                                 De, Re, beta, i->second );
123 <        } else {
74 >        if (mit == NULL) {
75            sprintf( painCave.errMsg,
76 <                   "Morse::initialize found unknown Morse interaction type\n"
77 <                   "\t(%s) for %s - %s interaction.\n",
78 <                   morseParam.interactionType.c_str(),
79 <                   atypes.first->getName().c_str(),
129 <                   atypes.second->getName().c_str());
76 >                   "Morse::initialize could not convert NonBondedInteractionType\n"
77 >                   "\tto MorseInteractionType for %s - %s interaction.\n",
78 >                   at1->getName().c_str(),
79 >                   at2->getName().c_str());
80            painCave.severity = OPENMD_ERROR;
81            painCave.isFatal = 1;
82            simError();          
83          }
84 +        
85 +        RealType De = mit->getD();
86 +        RealType Re = mit->getR();
87 +        RealType beta = mit->getBeta();
88 +  
89 +        MorseType variant = mit->getInteractionType();
90 +        addExplicitInteraction(at1, at2, De, Re, beta, variant );
91        }
92      }  
93      initialized_ = true;
# Line 138 | Line 95 | namespace OpenMD {
95        
96    void Morse::addExplicitInteraction(AtomType* atype1, AtomType* atype2,
97                                       RealType De, RealType Re, RealType beta,
98 <                                     MorseInteractionType mit) {
98 >                                     MorseType mt) {
99  
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 –
100      MorseInteractionData mixer;
101      mixer.De = De;
102      mixer.Re = Re;
103      mixer.beta = beta;
104 <    mixer.interactionType = mit;
104 >    mixer.variant = mt;
105  
106      pair<AtomType*, AtomType*> key1, key2;
107      key1 = make_pair(atype1, atype2);
# Line 162 | Line 113 | namespace OpenMD {
113      }    
114    }
115    
116 <  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) {
116 >  void Morse::calcForce(InteractionData &idat) {
117  
118      if (!initialized_) initialize();
119      
120 <    RealType myPot = 0.0;
121 <    RealType myPotC = 0.0;
122 <    RealType myDeriv = 0.0;
123 <    RealType myDerivC = 0.0;
124 <
125 <    pair<AtomType*, AtomType*> key = make_pair(at1, at2);
126 <    MorseInteractionData mixer = MixingMap[key];
127 <
128 <    RealType De = mixer.De;
129 <    RealType Re = mixer.Re;
130 <    RealType beta = mixer.beta;
131 <    MorseInteractionType interactionType = mixer.interactionType;
132 <  
133 <    // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
134 <    
135 <    RealType expt     = -beta*(r - Re);
136 <    RealType expfnc   = exp(expt);
137 <    RealType expfnc2  = expfnc*expfnc;
138 <
139 <    RealType exptC = 0.0;
140 <    RealType expfncC = 0.0;
141 <    RealType expfnc2C = 0.0;
142 <
143 <    if (Morse::shiftedPot_ || Morse::shiftedFrc_) {
144 <      exptC     = -beta*(rcut - Re);
145 <      expfncC   = exp(exptC);
146 <      expfnc2C  = expfncC*expfncC;
147 <    }
148 <
201 <
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;
120 >    map<pair<AtomType*, AtomType*>, MorseInteractionData>::iterator it;
121 >    it = MixingMap.find( idat.atypes );
122 >    if (it != MixingMap.end()) {
123 >      MorseInteractionData mixer = (*it).second;
124 >      
125 >      RealType myPot = 0.0;
126 >      RealType myPotC = 0.0;
127 >      RealType myDeriv = 0.0;
128 >      RealType myDerivC = 0.0;
129 >      
130 >      RealType De = mixer.De;
131 >      RealType Re = mixer.Re;
132 >      RealType beta = mixer.beta;
133 >      MorseType variant = mixer.variant;
134 >      
135 >      // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
136 >      
137 >      RealType expt     = -beta*( *(idat.rij) - Re);
138 >      RealType expfnc   = exp(expt);
139 >      RealType expfnc2  = expfnc*expfnc;
140 >      
141 >      RealType exptC = 0.0;
142 >      RealType expfncC = 0.0;
143 >      RealType expfnc2C = 0.0;
144 >      
145 >      if (idat.shiftedPot || idat.shiftedForce) {
146 >        exptC     = -beta*( *(idat.rcut) - Re);
147 >        expfncC   = exp(exptC);
148 >        expfnc2C  = expfncC*expfncC;
149        }
150 <
151 <      break;
152 <    }
153 <    case repulsiveMorse : {
154 <
155 <      myPot  = De * expfnc2;
156 <      myDeriv  = -2.0 * De * beta * expfnc2;
157 <
158 <      if (Morse::shiftedPot_) {
159 <        myPotC = De * expfnc2C;
160 <        myDerivC = 0.0;
161 <      } else if (Morse::shiftedFrc_) {
162 <        myPotC = De * expfnc2C;
163 <        myDerivC = -2.0 * De * beta * expfnc2C;
164 <        myPotC += myDerivC * (r - rcut);
165 <      } else {
166 <        myPotC = 0.0;
167 <        myDerivC = 0.0;
150 >      
151 >      
152 >      switch(variant) {
153 >      case mtShifted : {
154 >        
155 >        myPot  = De * (expfnc2  - 2.0 * expfnc);
156 >        myDeriv   = 2.0 * De * beta * (expfnc - expfnc2);
157 >        
158 >        if (idat.shiftedPot) {
159 >          myPotC = De * (expfnc2C - 2.0 * expfncC);
160 >          myDerivC = 0.0;
161 >        } else if (idat.shiftedForce) {
162 >          myPotC = De * (expfnc2C - 2.0 * expfncC);
163 >          myDerivC  = 2.0 * De * beta * (expfnc2C - expfnc2C);
164 >          myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut) );
165 >        } else {
166 >          myPotC = 0.0;
167 >          myDerivC = 0.0;
168 >        }
169 >        
170 >        break;
171        }
172 <
173 <      break;
172 >      case mtRepulsive : {
173 >        
174 >        myPot  = De * expfnc2;
175 >        myDeriv  = -2.0 * De * beta * expfnc2;
176 >        
177 >        if (idat.shiftedPot) {
178 >          myPotC = De * expfnc2C;
179 >          myDerivC = 0.0;
180 >        } else if (idat.shiftedForce) {
181 >          myPotC = De * expfnc2C;
182 >          myDerivC = -2.0 * De * beta * expfnc2C;
183 >          myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut));
184 >        } else {
185 >          myPotC = 0.0;
186 >          myDerivC = 0.0;
187 >        }
188 >        
189 >        break;
190 >      }
191 >      case mtUnknown: {
192 >        // don't know what to do so don't do anything
193 >        break;
194 >      }
195 >      }
196 >      
197 >      RealType pot_temp = *(idat.vdwMult) * (myPot - myPotC);
198 >      *(idat.vpair) += pot_temp;
199 >      
200 >      RealType dudr = *(idat.sw) * *(idat.vdwMult) * (myDeriv - myDerivC);
201 >      
202 >      (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
203 >      *(idat.f1) = *(idat.d) * dudr / *(idat.rij);
204      }
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 –
205      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();
206      
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;    
207    }
208      
209 <  void Morse::setMorseDefaultCutoff(RealType *thisRcut, int *sP, int *sF) {
210 <    shiftedPot_ = (bool)(*sP);
211 <    shiftedFrc_ = (bool)(*sF);
212 <  }
213 < }
209 >  RealType Morse::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
210 >    if (!initialized_) initialize();  
211 >    map<pair<AtomType*, AtomType*>, MorseInteractionData>::iterator it;
212 >    it = MixingMap.find(atypes);
213 >    if (it == MixingMap.end())
214 >      return 0.0;
215 >    else  {
216 >      MorseInteractionData mixer = (*it).second;
217  
218 < extern "C" {
219 <  
220 < #define fortranSetMorseCutoff FC_FUNC(setmorsedefaultcutoff, SETMORSEDEFAULTCUTOFF)
221 < #define fortranDoMorsePair FC_FUNC(do_morse_pair, DO_MORSE_PAIR)
222 <  
223 <  void fortranSetMorseCutoff(RealType *rcut, int *shiftedPot, int *shiftedFrc) {
224 <    return OpenMD::Morse::Instance()->setMorseDefaultCutoff(rcut, shiftedPot,
225 <                                                            shiftedFrc);
218 >      RealType Re = mixer.Re;
219 >      RealType beta = mixer.beta;
220 >      // This value of the r corresponds to an energy about 1.48% of
221 >      // the energy at the bottom of the Morse well.  For comparison, the
222 >      // Lennard-Jones function is about 1.63% of it's minimum value at
223 >      // a distance of 2.5 sigma.
224 >      return (4.9 + beta * Re) / beta;
225 >    }
226    }
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 –  }
227   }
228 +

Diff Legend

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