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 1664 by gezelter, Tue Nov 22 14:37:41 2011 UTC

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

Diff Legend

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