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

Comparing trunk/src/nonbonded/Morse.cpp (file contents):
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 2033 by gezelter, Sat Nov 1 14:12:16 2014 UTC

# Line 56 | Line 56 | namespace OpenMD {
56    
57    void Morse::initialize() {    
58  
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;
61 –    NonBondedInteractionType* nbt;
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)) {
# Line 69 | Line 75 | namespace OpenMD {
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 +        int atid2 = at2->getIdent();
85 +        if (Mtids[atid2] == -1) {
86 +          mtid2 = Mtypes.size();
87 +          Mtypes.insert(atid2);
88 +          Mtids[atid2] = mtid2;
89 +        }
90 +    
91          MorseInteractionType* mit = dynamic_cast<MorseInteractionType*>(nbt);
92  
93          if (mit == NULL) {
# Line 103 | Line 122 | namespace OpenMD {
122      mixer.beta = beta;
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    
139    void Morse::calcForce(InteractionData &idat) {
140  
141      if (!initialized_) initialize();
142 +
143 +  
144 +    MorseInteractionData &mixer = MixingMap[Mtids[idat.atid1]][Mtids[idat.atid2]];
145 +
146 +    RealType myPot = 0.0;
147 +    RealType myPotC = 0.0;
148 +    RealType myDeriv = 0.0;
149 +    RealType myDerivC = 0.0;
150      
151 <    map<pair<AtomType*, AtomType*>, MorseInteractionData>::iterator it;
152 <    it = MixingMap.find( idat.atypes );
153 <    if (it != MixingMap.end()) {
154 <      MorseInteractionData mixer = (*it).second;
151 >    RealType De = mixer.De;
152 >    RealType Re = mixer.Re;
153 >    RealType beta = mixer.beta;
154 >    MorseType variant = mixer.variant;
155 >    
156 >    // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
157 >    
158 >    RealType expt     = -beta*( *(idat.rij) - Re);
159 >    RealType expfnc   = exp(expt);
160 >    RealType expfnc2  = expfnc*expfnc;
161 >    
162 >    RealType exptC = 0.0;
163 >    RealType expfncC = 0.0;
164 >    RealType expfnc2C = 0.0;
165 >    
166 >    if (idat.shiftedPot || idat.shiftedForce) {
167 >      exptC     = -beta*( *(idat.rcut) - Re);
168 >      expfncC   = exp(exptC);
169 >      expfnc2C  = expfncC*expfncC;
170 >    }
171 >    
172 >    
173 >    switch(variant) {
174 >    case mtShifted : {
175 >      myPot  = De * (expfnc2  - 2.0 * expfnc);
176 >      myDeriv   = 2.0 * De * beta * (expfnc - expfnc2);
177        
178 <      RealType myPot = 0.0;
179 <      RealType myPotC = 0.0;
180 <      RealType myDeriv = 0.0;
181 <      RealType myDerivC = 0.0;
182 <      
183 <      RealType De = mixer.De;
184 <      RealType Re = mixer.Re;
185 <      RealType beta = mixer.beta;
186 <      MorseType variant = mixer.variant;
187 <      
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;
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 +      break;
191 +    }
192 +    case mtRepulsive : {
193 +      myPot  = De * expfnc2;
194 +      myDeriv  = -2.0 * De * beta * expfnc2;
195        
196 <      switch(variant) {
197 <      case mtShifted : {
198 <        
199 <        myPot  = De * (expfnc2  - 2.0 * expfnc);
200 <        myDeriv   = 2.0 * De * beta * (expfnc - expfnc2);
201 <        
202 <        if (idat.shiftedPot) {
203 <          myPotC = De * (expfnc2C - 2.0 * expfncC);
204 <          myDerivC = 0.0;
205 <        } else if (idat.shiftedForce) {
162 <          myPotC = De * (expfnc2C - 2.0 * expfncC);
163 <          myDerivC  = 2.0 * De * beta * (expfncC - expfnc2C);
164 <          myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut) );
165 <        } else {
166 <          myPotC = 0.0;
167 <          myDerivC = 0.0;
168 <        }
169 <        
170 <        break;
196 >      if (idat.shiftedPot) {
197 >        myPotC = De * expfnc2C;
198 >        myDerivC = 0.0;
199 >      } else if (idat.shiftedForce) {
200 >        myPotC = De * expfnc2C;
201 >        myDerivC = -2.0 * De * beta * expfnc2C;
202 >        myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut));
203 >      } else {
204 >        myPotC = 0.0;
205 >        myDerivC = 0.0;
206        }
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 –      }
207        
208 <      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);
208 >      break;
209      }
210 <    return;
210 >    case mtUnknown: {
211 >      // don't know what to do so don't do anything
212 >      break;
213 >    }
214 >    }
215      
216 <  }
216 >    RealType pot_temp = *(idat.vdwMult) * (myPot - myPotC);
217 >    *(idat.vpair) += pot_temp;
218      
219 +    RealType dudr = *(idat.sw) * *(idat.vdwMult) * (myDeriv - myDerivC);
220 +    
221 +    (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
222 +    *(idat.f1) = *(idat.d) * dudr / *(idat.rij);
223 +
224 +    return;    
225 +  }
226 +
227    RealType Morse::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
228      if (!initialized_) initialize();  
229 <    map<pair<AtomType*, AtomType*>, MorseInteractionData>::iterator it;
230 <    it = MixingMap.find(atypes);
231 <    if (it == MixingMap.end())
232 <      return 0.0;
233 <    else  {
234 <      MorseInteractionData mixer = (*it).second;
235 <
229 >    
230 >    int atid1 = atypes.first->getIdent();
231 >    int atid2 = atypes.second->getIdent();
232 >    int mtid1 = Mtids[atid1];
233 >    int mtid2 = Mtids[atid2];
234 >    
235 >    if ( mtid1 == -1 || mtid2 == -1) return 0.0;
236 >    else {      
237 >      MorseInteractionData mixer = MixingMap[mtid1][mtid2];
238        RealType Re = mixer.Re;
239        RealType beta = mixer.beta;
240        // This value of the r corresponds to an energy about 1.48% of

Diff Legend

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