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

Comparing:
branches/development/src/nonbonded/Morse.cpp (file contents), Revision 1571 by gezelter, Fri May 27 16:45:44 2011 UTC vs.
trunk/src/nonbonded/Morse.cpp (file contents), Revision 2017 by gezelter, Tue Sep 2 18:31:44 2014 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 <  Morse::Morse() : name_("Morse"), initialized_(false), forceField_(NULL),
55 <                   shiftedPot_(false), shiftedFrc_(false) {}
55 >  Morse::Morse() : name_("Morse"), initialized_(false), forceField_(NULL) {}
56    
57    void Morse::initialize() {    
58  
59 <    stringToEnumMap_["shiftedMorse"] = shiftedMorse;
60 <    stringToEnumMap_["repulsiveMorse"] = repulsiveMorse;
59 >    Mtypes.clear();
60 >    Mtids.clear();
61 >    MixingMap.clear();
62 >    Mtids.resize( forceField_->getNAtomType(), -1);
63  
64      ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
65      ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
66 +    ForceField::NonBondedInteractionTypeContainer::KeyType keys;
67      NonBondedInteractionType* nbt;
68 +    int mtid1, mtid2;
69  
70      for (nbt = nbiTypes->beginType(j); nbt != NULL;
71           nbt = nbiTypes->nextType(j)) {
72        
73        if (nbt->isMorse()) {
74 <        
75 <        pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
76 <        
77 <        GenericData* data = nbt->getPropertyByName("Morse");
78 <        if (data == NULL) {
79 <          sprintf( painCave.errMsg, "Morse::initialize could not find\n"
80 <                   "\tMorse parameters for %s - %s interaction.\n",
81 <                   atypes.first->getName().c_str(),
82 <                   atypes.second->getName().c_str());
79 <          painCave.severity = OPENMD_ERROR;
80 <          painCave.isFatal = 1;
81 <          simError();
74 >        keys = nbiTypes->getKeys(j);
75 >        AtomType* at1 = forceField_->getAtomType(keys[0]);
76 >        AtomType* at2 = forceField_->getAtomType(keys[1]);
77 >
78 >        int atid1 = at1->getIdent();
79 >        if (Mtids[atid1] == -1) {
80 >          mtid1 = Mtypes.size();
81 >          Mtypes.insert(atid1);
82 >          Mtids[atid1] = mtid1;        
83          }
84 <        
85 <        MorseData* morseData = dynamic_cast<MorseData*>(data);
86 <        if (morseData == NULL) {
87 <          sprintf( painCave.errMsg,
88 <                   "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();          
84 >        int atid2 = at2->getIdent();
85 >        if (Mtids[atid2] == -1) {
86 >          mtid2 = Mtypes.size();
87 >          Mtypes.insert(atid2);
88 >          Mtids[atid2] = mtid2;
89          }
90 <        
91 <        MorseParam morseParam = morseData->getData();
90 >    
91 >        MorseInteractionType* mit = dynamic_cast<MorseInteractionType*>(nbt);
92  
93 <        RealType De = morseParam.De;
99 <        RealType Re = morseParam.Re;
100 <        RealType beta = morseParam.beta;
101 <        string interactionType = morseParam.interactionType;
102 <
103 <        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 {
93 >        if (mit == NULL) {
94            sprintf( painCave.errMsg,
95 <                   "Morse::initialize found unknown Morse interaction type\n"
96 <                   "\t(%s) for %s - %s interaction.\n",
97 <                   morseParam.interactionType.c_str(),
98 <                   atypes.first->getName().c_str(),
115 <                   atypes.second->getName().c_str());
95 >                   "Morse::initialize could not convert NonBondedInteractionType\n"
96 >                   "\tto MorseInteractionType for %s - %s interaction.\n",
97 >                   at1->getName().c_str(),
98 >                   at2->getName().c_str());
99            painCave.severity = OPENMD_ERROR;
100            painCave.isFatal = 1;
101            simError();          
102          }
103 +        
104 +        RealType De = mit->getD();
105 +        RealType Re = mit->getR();
106 +        RealType beta = mit->getBeta();
107 +  
108 +        MorseType variant = mit->getInteractionType();
109 +        addExplicitInteraction(at1, at2, De, Re, beta, variant );
110        }
111      }  
112      initialized_ = true;
# Line 124 | Line 114 | namespace OpenMD {
114        
115    void Morse::addExplicitInteraction(AtomType* atype1, AtomType* atype2,
116                                       RealType De, RealType Re, RealType beta,
117 <                                     MorseInteractionType mit) {
117 >                                     MorseType mt) {
118  
119      MorseInteractionData mixer;
120      mixer.De = De;
121      mixer.Re = Re;
122      mixer.beta = beta;
123 <    mixer.interactionType = mit;
123 >    mixer.variant = mt;
124  
125 <    pair<AtomType*, AtomType*> key1, key2;
126 <    key1 = make_pair(atype1, atype2);
127 <    key2 = make_pair(atype2, atype1);
125 >    int mtid1 = Mtids[atype1->getIdent()];
126 >    int mtid2 = Mtids[atype2->getIdent()];
127 >    int nM = Mtypes.size();
128 >
129 >    MixingMap.resize(nM);
130 >    MixingMap[mtid1].resize(nM);
131      
132 <    MixingMap[key1] = mixer;
133 <    if (key2 != key1) {
134 <      MixingMap[key2] = mixer;
132 >    MixingMap[mtid1][mtid2] = mixer;
133 >    if (mtid2 != mtid1) {
134 >      MixingMap[mtid2].resize(nM);
135 >      MixingMap[mtid2][mtid1] = mixer;
136      }    
137    }
138    
# Line 146 | Line 140 | namespace OpenMD {
140  
141      if (!initialized_) initialize();
142      
143 <    map<pair<AtomType*, AtomType*>, MorseInteractionData>::iterator it;
144 <    it = MixingMap.find( idat.atypes );
145 <    if (it != MixingMap.end()) {
146 <      MorseInteractionData mixer = (*it).second;
143 >    MorseInteractionData &mixer = MixingMap[Mtids[idat.atid1]][Mtids[idat.atid2]];
144 >
145 >    RealType myPot = 0.0;
146 >    RealType myPotC = 0.0;
147 >    RealType myDeriv = 0.0;
148 >    RealType myDerivC = 0.0;
149 >    
150 >    RealType De = mixer.De;
151 >    RealType Re = mixer.Re;
152 >    RealType beta = mixer.beta;
153 >    MorseType variant = mixer.variant;
154 >    
155 >    // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
156 >    
157 >    RealType expt     = -beta*( *(idat.rij) - Re);
158 >    RealType expfnc   = exp(expt);
159 >    RealType expfnc2  = expfnc*expfnc;
160 >    
161 >    RealType exptC = 0.0;
162 >    RealType expfncC = 0.0;
163 >    RealType expfnc2C = 0.0;
164 >    
165 >    if (idat.shiftedPot || idat.shiftedForce) {
166 >      exptC     = -beta*( *(idat.rcut) - Re);
167 >      expfncC   = exp(exptC);
168 >      expfnc2C  = expfncC*expfncC;
169 >    }
170 >    
171 >    
172 >    switch(variant) {
173 >    case mtShifted : {
174        
175 <      RealType myPot = 0.0;
176 <      RealType myPotC = 0.0;
156 <      RealType myDeriv = 0.0;
157 <      RealType myDerivC = 0.0;
175 >      myPot  = De * (expfnc2  - 2.0 * expfnc);
176 >      myDeriv   = 2.0 * De * beta * (expfnc - expfnc2);
177        
178 <      RealType De = mixer.De;
179 <      RealType Re = mixer.Re;
180 <      RealType beta = mixer.beta;
181 <      MorseInteractionType interactionType = mixer.interactionType;
182 <      
183 <      // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
184 <      
185 <      RealType expt     = -beta*( *(idat.rij) - Re);
186 <      RealType expfnc   = exp(expt);
187 <      RealType expfnc2  = expfnc*expfnc;
169 <      
170 <      RealType exptC = 0.0;
171 <      RealType expfncC = 0.0;
172 <      RealType expfnc2C = 0.0;
173 <      
174 <      if (Morse::shiftedPot_ || Morse::shiftedFrc_) {
175 <        exptC     = -beta*( *(idat.rcut) - Re);
176 <        expfncC   = exp(exptC);
177 <        expfnc2C  = expfncC*expfncC;
178 >      if (idat.shiftedPot) {
179 >        myPotC = De * (expfnc2C - 2.0 * expfncC);
180 >        myDerivC = 0.0;
181 >      } else if (idat.shiftedForce) {
182 >        myPotC = De * (expfnc2C - 2.0 * expfncC);
183 >        myDerivC  = 2.0 * De * beta * (expfncC - expfnc2C);
184 >        myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut) );
185 >      } else {
186 >        myPotC = 0.0;
187 >        myDerivC = 0.0;
188        }
189        
190 <      
191 <      switch(interactionType) {
192 <      case shiftedMorse : {
190 >      break;
191 >    }
192 >    case mtRepulsive : {
193          
194 <        myPot  = De * (expfnc2  - 2.0 * expfnc);
195 <        myDeriv   = 2.0 * De * beta * (expfnc - expfnc2);
196 <        
197 <        if (Morse::shiftedPot_) {
198 <          myPotC = De * (expfnc2C - 2.0 * expfncC);
199 <          myDerivC = 0.0;
200 <        } else if (Morse::shiftedFrc_) {
201 <          myPotC = De * (expfnc2C - 2.0 * expfncC);
202 <          myDerivC  = 2.0 * De * beta * (expfnc2C - expfnc2C);
203 <          myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut) );
204 <        } else {
205 <          myPotC = 0.0;
206 <          myDerivC = 0.0;
197 <        }
198 <        
199 <        break;
194 >      myPot  = De * expfnc2;
195 >      myDeriv  = -2.0 * De * beta * expfnc2;
196 >      
197 >      if (idat.shiftedPot) {
198 >        myPotC = De * expfnc2C;
199 >        myDerivC = 0.0;
200 >      } else if (idat.shiftedForce) {
201 >        myPotC = De * expfnc2C;
202 >        myDerivC = -2.0 * De * beta * expfnc2C;
203 >        myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut));
204 >      } else {
205 >        myPotC = 0.0;
206 >        myDerivC = 0.0;
207        }
201 –      case repulsiveMorse : {
202 –        
203 –        myPot  = De * expfnc2;
204 –        myDeriv  = -2.0 * De * beta * expfnc2;
205 –        
206 –        if (Morse::shiftedPot_) {
207 –          myPotC = De * expfnc2C;
208 –          myDerivC = 0.0;
209 –        } else if (Morse::shiftedFrc_) {
210 –          myPotC = De * expfnc2C;
211 –          myDerivC = -2.0 * De * beta * expfnc2C;
212 –          myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut));
213 –        } else {
214 –          myPotC = 0.0;
215 –          myDerivC = 0.0;
216 –        }
217 –        
218 –        break;
219 –      }
220 –      }
208        
209 <      RealType pot_temp = *(idat.vdwMult) * (myPot - myPotC);
223 <      *(idat.vpair) += pot_temp;
224 <      
225 <      RealType dudr = *(idat.sw) * *(idat.vdwMult) * (myDeriv - myDerivC);
226 <      
227 <      idat.pot[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
228 <      *(idat.f1) = *(idat.d) * dudr / *(idat.rij);
209 >      break;
210      }
211 <    return;
211 >    case mtUnknown: {
212 >      // don't know what to do so don't do anything
213 >      break;
214 >    }
215 >    }
216      
217 <  }
217 >
218 >    RealType pot_temp = *(idat.vdwMult) * (myPot - myPotC);
219 >    *(idat.vpair) += pot_temp;
220      
221 +    RealType dudr = *(idat.sw) * *(idat.vdwMult) * (myDeriv - myDerivC);
222 +
223 +    
224 +    (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
225 +    *(idat.f1) = *(idat.d) * dudr / *(idat.rij);
226 +    
227 +    return;    
228 +  }
229 +
230    RealType Morse::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
231      if (!initialized_) initialize();  
232 <    map<pair<AtomType*, AtomType*>, MorseInteractionData>::iterator it;
233 <    it = MixingMap.find(atypes);
234 <    if (it == MixingMap.end())
235 <      return 0.0;
236 <    else  {
237 <      MorseInteractionData mixer = (*it).second;
238 <
232 >    
233 >    int atid1 = atypes.first->getIdent();
234 >    int atid2 = atypes.second->getIdent();
235 >    int mtid1 = Mtids[atid1];
236 >    int mtid2 = Mtids[atid2];
237 >    
238 >    if ( mtid1 == -1 || mtid2 == -1) return 0.0;
239 >    else {      
240 >      MorseInteractionData mixer = MixingMap[mtid1][mtid2];
241        RealType Re = mixer.Re;
242        RealType beta = mixer.beta;
243        // This value of the r corresponds to an energy about 1.48% of

Diff Legend

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