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 1895 by gezelter, Mon Jul 1 21:09:37 2013 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    
# Line 117 | 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;
127 <      RealType myDeriv = 0.0;
128 <      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 <      MorseType variant = mixer.variant;
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;
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 <      
191 <      switch(variant) {
192 <      case mtShifted : {
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 (idat.shiftedPot) {
198 <          myPotC = De * (expfnc2C - 2.0 * expfncC);
199 <          myDerivC = 0.0;
200 <        } else if (idat.shiftedForce) {
201 <          myPotC = De * (expfnc2C - 2.0 * expfncC);
202 <          myDerivC  = 2.0 * De * beta * (expfncC - expfnc2C);
203 <          myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut) );
204 <        } else {
205 <          myPotC = 0.0;
206 <          myDerivC = 0.0;
168 <        }
169 <        
170 <        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        }
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 –      }
208        
209 <      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);
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 >    RealType pot_temp = *(idat.vdwMult) * (myPot - myPotC);
218 >    *(idat.vpair) += pot_temp;
219      
220 +    RealType dudr = *(idat.sw) * *(idat.vdwMult) * (myDeriv - myDerivC);
221 +    
222 +    (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
223 +    *(idat.f1) = *(idat.d) * dudr / *(idat.rij);
224 +    
225 +    return;    
226 +  }
227 +
228    RealType Morse::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
229      if (!initialized_) initialize();  
230 <    map<pair<AtomType*, AtomType*>, MorseInteractionData>::iterator it;
231 <    it = MixingMap.find(atypes);
232 <    if (it == MixingMap.end())
233 <      return 0.0;
234 <    else  {
235 <      MorseInteractionData mixer = (*it).second;
236 <
230 >    
231 >    int atid1 = atypes.first->getIdent();
232 >    int atid2 = atypes.second->getIdent();
233 >    int mtid1 = Mtids[atid1];
234 >    int mtid2 = Mtids[atid2];
235 >    
236 >    if ( mtid1 == -1 || mtid2 == -1) return 0.0;
237 >    else {      
238 >      MorseInteractionData mixer = MixingMap[mtid1][mtid2];
239        RealType Re = mixer.Re;
240        RealType beta = mixer.beta;
241        // This value of the r corresponds to an energy about 1.48% of

Diff Legend

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