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

Comparing:
branches/development/src/nonbonded/MAW.cpp (file contents), Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC vs.
trunk/src/nonbonded/MAW.cpp (file contents), Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 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).          
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   */
# Line 51 | Line 51 | namespace OpenMD {
51  
52   namespace OpenMD {
53  
54 <  MAW::MAW() : name_("MAW"), initialized_(false), forceField_(NULL) {}
54 >  MAW::MAW() : initialized_(false), forceField_(NULL), name_("MAW") {}
55    
56    void MAW::initialize() {    
57  
58 +    MAWtypes.clear();
59 +    MAWtids.clear();
60 +    MixingMap.clear();
61 +    MAWtids.resize( forceField_->getNAtomType(), -1);
62 +
63      ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
64      ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
60 –    NonBondedInteractionType* nbt;
65      ForceField::NonBondedInteractionTypeContainer::KeyType keys;
66 +    NonBondedInteractionType* nbt;
67 +    int mtid1, mtid2;
68  
69      for (nbt = nbiTypes->beginType(j); nbt != NULL;
70           nbt = nbiTypes->nextType(j)) {
# Line 68 | Line 74 | namespace OpenMD {
74          AtomType* at1 = forceField_->getAtomType(keys[0]);
75          AtomType* at2 = forceField_->getAtomType(keys[1]);
76  
77 +        int atid1 = at1->getIdent();
78 +        if (MAWtids[atid1] == -1) {
79 +          mtid1 = MAWtypes.size();
80 +          MAWtypes.insert(atid1);
81 +          MAWtids[atid1] = mtid1;        
82 +        }
83 +        int atid2 = at2->getIdent();
84 +        if (MAWtids[atid2] == -1) {
85 +          mtid2 = MAWtypes.size();
86 +          MAWtypes.insert(atid2);
87 +          MAWtids[atid2] = mtid2;
88 +        }
89 +
90          MAWInteractionType* mit = dynamic_cast<MAWInteractionType*>(nbt);
91  
92          if (mit == NULL) {
# Line 104 | Line 123 | namespace OpenMD {
123      mixer.Re = Re;
124      mixer.ca1 = ca1;
125      mixer.cb1 = cb1;
126 +    mixer.j_is_Metal = atype2->isMetal();
127  
128 <    pair<AtomType*, AtomType*> key1, key2;
129 <    key1 = make_pair(atype1, atype2);
130 <    key2 = make_pair(atype2, atype1);
128 >    int mtid1 = MAWtids[atype1->getIdent()];
129 >    int mtid2 = MAWtids[atype2->getIdent()];
130 >    int nM = MAWtypes.size();
131 >
132 >    MixingMap.resize(nM);
133 >    MixingMap[mtid1].resize(nM);
134      
135 <    MixingMap[key1] = mixer;
136 <    if (key2 != key1) {
137 <      MixingMap[key2] = mixer;
135 >    MixingMap[mtid1][mtid2] = mixer;
136 >    if (mtid2 != mtid1) {
137 >      MixingMap[mtid2].resize(nM);
138 >      mixer.j_is_Metal = atype1->isMetal();
139 >      MixingMap[mtid2][mtid1] = mixer;
140      }    
141    }
142    
143    void MAW::calcForce(InteractionData &idat) {
144  
145      if (!initialized_) initialize();
121 –    
122 –    map<pair<AtomType*, AtomType*>, MAWInteractionData>::iterator it;
123 –    it = MixingMap.find( idat.atypes );
124 –    if (it != MixingMap.end()) {
125 –      MAWInteractionData mixer = (*it).second;
126 –      
127 –      RealType myPot = 0.0;
128 –      RealType myPotC = 0.0;
129 –      RealType myDeriv = 0.0;
130 –      RealType myDerivC = 0.0;
131 –      
132 –      RealType D_e = mixer.De;
133 –      RealType R_e = mixer.Re;
134 –      RealType beta = mixer.beta;
135 –      RealType ca1 = mixer.ca1;
136 –      RealType cb1 = mixer.cb1;
146  
147 <      bool j_is_Metal = idat.atypes.second->isMetal();
139 <
140 <      Vector3d r;
141 <      RotMat3x3d Atrans;
142 <      if (j_is_Metal) {
143 <        // rotate the inter-particle separation into the two different
144 <        // body-fixed coordinate systems:
145 <        r = *(idat.A1) * *(idat.d);
146 <        Atrans = idat.A1->transpose();
147 <      } else {
148 <        // negative sign because this is the vector from j to i:      
149 <        r = -*(idat.A2) * *(idat.d);
150 <        Atrans = idat.A2->transpose();
151 <      }
152 <      
153 <      // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
154 <      
155 <      RealType expt     = -beta*( *(idat.rij) - R_e);
156 <      RealType expfnc   = exp(expt);
157 <      RealType expfnc2  = expfnc*expfnc;
158 <      
159 <      RealType exptC = 0.0;
160 <      RealType expfncC = 0.0;
161 <      RealType expfnc2C = 0.0;
162 <
163 <      myPot  = D_e * (expfnc2  - 2.0 * expfnc);
164 <      myDeriv   = 2.0 * D_e * beta * (expfnc - expfnc2);
165 <      
166 <      if (idat.shiftedPot || idat.shiftedForce) {
167 <        exptC     = -beta*( *(idat.rcut)  - R_e);
168 <        expfncC   = exp(exptC);
169 <        expfnc2C  = expfncC*expfncC;
170 <      }
171 <      
172 <      if (idat.shiftedPot) {
173 <        myPotC = D_e * (expfnc2C - 2.0 * expfncC);
174 <        myDerivC = 0.0;
175 <      } else if (idat.shiftedForce) {
176 <        myPotC = D_e * (expfnc2C - 2.0 * expfncC);
177 <        myDerivC  = 2.0 * D_e * beta * (expfnc2C - expfnc2C);
178 <        myPotC += myDerivC * ( *(idat.rij)  -  *(idat.rcut) );
179 <      } else {
180 <        myPotC = 0.0;
181 <        myDerivC = 0.0;
182 <      }
183 <
184 <      RealType x = r.x();
185 <      RealType y = r.y();
186 <      RealType z = r.z();
187 <      RealType x2 = x * x;
188 <      RealType y2 = y * y;
189 <      RealType z2 = z * z;
190 <
191 <      RealType r3 = *(idat.r2) *  *(idat.rij) ;
192 <      RealType r4 = *(idat.r2) *  *(idat.r2);
193 <
194 <      // angular modulation of morse part of potential to approximate
195 <      // the squares of the two HOMO lone pair orbitals in water:
196 <      //********************** old form*************************
197 <      // s = 1 / (4 pi)
198 <      // ta1 = (s - pz)^2 = (1 - sqrt(3)*cos(theta))^2 / (4 pi)
199 <      // b1 = px^2 = 3 * (sin(theta)*cos(phi))^2 / (4 pi)  
200 <      //********************** old form*************************
201 <      // we'll leave out the 4 pi for now
202 <      
203 <      // new functional form just using the p orbitals.
204 <      // Vmorse(r)*[a*p_x + b p_z + (1-a-b)]
205 <      // which is
206 <      // Vmorse(r)*[a sin^2(theta) cos^2(phi) + b cos(theta) + (1-a-b)]
207 <      // Vmorse(r)*[a*x2/r2 + b*z/r + (1-a-b)]
208 <      
209 <      RealType Vmorse = (myPot - myPotC);
210 <      RealType Vang = ca1 * x2 / *(idat.r2) +
211 <        cb1 * z /  *(idat.rij)  + (0.8 - ca1 / 3.0);
212 <      
213 <      RealType pot_temp = *(idat.vdwMult) * Vmorse * Vang;
214 <      *(idat.vpair) += pot_temp;
215 <      (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
216 <          
217 <      Vector3d dVmorsedr = (myDeriv - myDerivC) * Vector3d(x, y, z) /  *(idat.rij) ;
218 <      
219 <      Vector3d dVangdr = Vector3d(-2.0 * ca1 * x2 * x / r4 + 2.0 * ca1 * x / *(idat.r2) - cb1 * x * z / r3,
220 <                                  -2.0 * ca1 * x2 * y / r4 - cb1 * z * y / r3,
221 <                                  -2.0 * ca1 * x2 * z / r4 + cb1 /  *(idat.rij)           - cb1 * z2  / r3);
222 <      
223 <      // chain rule to put these back on x, y, z
224 <
225 <      Vector3d dvdr = Vang * dVmorsedr + Vmorse * dVangdr;
226 <
227 <      // Torques for Vang using method of Price:
228 <      // S. L. Price, A. J. Stone, and M. Alderton, Mol. Phys. 52, 987 (1984).
229 <      
230 <      Vector3d dVangdu = Vector3d(cb1 * y /  *(idat.rij) ,
231 <                                  2.0 * ca1 * x * z / *(idat.r2) - cb1 * x /  *(idat.rij),
232 <                                  -2.0 * ca1 * y * x / *(idat.r2));
233 <
234 <      // do the torques first since they are easy:
235 <      // remember that these are still in the body fixed axes    
236 <
237 <      Vector3d trq = *(idat.vdwMult) * Vmorse * dVangdu * *(idat.sw);
238 <
239 <      // go back to lab frame using transpose of rotation matrix:
240 <      
241 <      if (j_is_Metal) {
242 <        *(idat.t1) += Atrans * trq;
243 <      } else {
244 <        *(idat.t2) += Atrans * trq;
245 <      }
147 >    MAWInteractionData &mixer = MixingMap[MAWtids[idat.atid1]][MAWtids[idat.atid2]];
148  
149 <      // Now, on to the forces (still in body frame of water)
150 <
151 <      Vector3d ftmp = *(idat.vdwMult) * *(idat.sw) * dvdr;
152 <
153 <      // rotate the terms back into the lab frame:
154 <      Vector3d flab;
155 <      if (j_is_Metal) {
156 <        flab = Atrans * ftmp;
157 <      } else {
158 <        flab = - Atrans * ftmp;
159 <      }
160 <      
161 <      *(idat.f1) += flab;
149 >    RealType myPot = 0.0;
150 >    RealType myPotC = 0.0;
151 >    RealType myDeriv = 0.0;
152 >    RealType myDerivC = 0.0;
153 >    
154 >    RealType D_e = mixer.De;
155 >    RealType R_e = mixer.Re;
156 >    RealType beta = mixer.beta;
157 >    RealType ca1 = mixer.ca1;
158 >    RealType cb1 = mixer.cb1;  
159 >    
160 >    Vector3d r;
161 >    RotMat3x3d Atrans;
162 >    if (mixer.j_is_Metal) {
163 >      // rotate the inter-particle separation into the two different
164 >      // body-fixed coordinate systems:
165 >      r = *(idat.A1) * *(idat.d);
166 >      Atrans = idat.A1->transpose();
167 >    } else {
168 >      // negative sign because this is the vector from j to i:      
169 >      r = -*(idat.A2) * *(idat.d);
170 >      Atrans = idat.A2->transpose();
171      }
261 –    return;
172      
173 <  }
173 >    // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
174      
175 +    RealType expt     = -beta*( *(idat.rij) - R_e);
176 +    RealType expfnc   = exp(expt);
177 +    RealType expfnc2  = expfnc*expfnc;
178 +    
179 +    RealType exptC = 0.0;
180 +    RealType expfncC = 0.0;
181 +    RealType expfnc2C = 0.0;
182 +    
183 +    myPot  = D_e * (expfnc2  - 2.0 * expfnc);
184 +    myDeriv   = 2.0 * D_e * beta * (expfnc - expfnc2);
185 +    
186 +    if (idat.shiftedPot || idat.shiftedForce) {
187 +      exptC     = -beta*( *(idat.rcut)  - R_e);
188 +      expfncC   = exp(exptC);
189 +      expfnc2C  = expfncC*expfncC;
190 +    }
191 +    
192 +    if (idat.shiftedPot) {
193 +      myPotC = D_e * (expfnc2C - 2.0 * expfncC);
194 +      myDerivC = 0.0;
195 +    } else if (idat.shiftedForce) {
196 +      myPotC = D_e * (expfnc2C - 2.0 * expfncC);
197 +      myDerivC  = 2.0 * D_e * beta * (expfncC - expfnc2C);
198 +      myPotC += myDerivC * ( *(idat.rij)  -  *(idat.rcut) );
199 +    } else {
200 +      myPotC = 0.0;
201 +      myDerivC = 0.0;
202 +    }
203 +    
204 +    RealType x = r.x();
205 +    RealType y = r.y();
206 +    RealType z = r.z();
207 +    RealType x2 = x * x;
208 +    RealType z2 = z * z;
209 +    
210 +    RealType r3 = *(idat.r2) *  *(idat.rij) ;
211 +    RealType r4 = *(idat.r2) *  *(idat.r2);
212 +    
213 +    // angular modulation of morse part of potential to approximate
214 +    // the squares of the two HOMO lone pair orbitals in water:
215 +    //********************** old form*************************
216 +    // s = 1 / (4 pi)
217 +    // ta1 = (s - pz)^2 = (1 - sqrt(3)*cos(theta))^2 / (4 pi)
218 +    // b1 = px^2 = 3 * (sin(theta)*cos(phi))^2 / (4 pi)  
219 +    //********************** old form*************************
220 +    // we'll leave out the 4 pi for now
221 +    
222 +    // new functional form just using the p orbitals.
223 +    // Vmorse(r)*[a*p_x + b p_z + (1-a-b)]
224 +    // which is
225 +    // Vmorse(r)*[a sin^2(theta) cos^2(phi) + b cos(theta) + (1-a-b)]
226 +    // Vmorse(r)*[a*x2/r2 + b*z/r + (1-a-b)]
227 +    
228 +    RealType Vmorse = (myPot - myPotC);
229 +    RealType Vang = ca1 * x2 / *(idat.r2) +
230 +      cb1 * z /  *(idat.rij)  + (0.8 - ca1 / 3.0);
231 +    
232 +    RealType pot_temp = *(idat.vdwMult) * Vmorse * Vang;
233 +    *(idat.vpair) += pot_temp;
234 +    (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
235 +    
236 +    Vector3d dVmorsedr = (myDeriv - myDerivC) * Vector3d(x, y, z) /  *(idat.rij) ;
237 +    
238 +    Vector3d dVangdr = Vector3d(-2.0 * ca1 * x2 * x / r4 + 2.0 * ca1 * x / *(idat.r2) - cb1 * x * z / r3,
239 +                                -2.0 * ca1 * x2 * y / r4 - cb1 * z * y / r3,
240 +                                -2.0 * ca1 * x2 * z / r4 + cb1 /  *(idat.rij)           - cb1 * z2  / r3);
241 +    
242 +    // chain rule to put these back on x, y, z
243 +    
244 +    Vector3d dvdr = Vang * dVmorsedr + Vmorse * dVangdr;
245 +    
246 +    // Torques for Vang using method of Price:
247 +    // S. L. Price, A. J. Stone, and M. Alderton, Mol. Phys. 52, 987 (1984).
248 +    
249 +    Vector3d dVangdu = Vector3d(cb1 * y /  *(idat.rij) ,
250 +                                2.0 * ca1 * x * z / *(idat.r2) - cb1 * x /  *(idat.rij),
251 +                                -2.0 * ca1 * y * x / *(idat.r2));
252 +    
253 +    // do the torques first since they are easy:
254 +    // remember that these are still in the body fixed axes    
255 +    
256 +    Vector3d trq = *(idat.vdwMult) * Vmorse * dVangdu * *(idat.sw);
257 +    
258 +    // go back to lab frame using transpose of rotation matrix:
259 +    
260 +    if (mixer.j_is_Metal) {
261 +      *(idat.t1) += Atrans * trq;
262 +    } else {
263 +      *(idat.t2) += Atrans * trq;
264 +    }
265 +    
266 +    // Now, on to the forces (still in body frame of water)
267 +    
268 +    Vector3d ftmp = *(idat.vdwMult) * *(idat.sw) * dvdr;
269 +    
270 +    // rotate the terms back into the lab frame:
271 +    Vector3d flab;
272 +    if (mixer.j_is_Metal) {
273 +      flab = Atrans * ftmp;
274 +    } else {
275 +      flab = - Atrans * ftmp;
276 +    }
277 +    
278 +    *(idat.f1) += flab;
279 +    
280 +    return;    
281 +  }
282 +
283    RealType MAW::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
284      if (!initialized_) initialize();  
285 <    map<pair<AtomType*, AtomType*>, MAWInteractionData>::iterator it;
286 <    it = MixingMap.find(atypes);
287 <    if (it == MixingMap.end())
288 <      return 0.0;
289 <    else  {
290 <      MAWInteractionData mixer = (*it).second;
291 <
285 >    int atid1 = atypes.first->getIdent();
286 >    int atid2 = atypes.second->getIdent();
287 >    int mtid1 = MAWtids[atid1];
288 >    int mtid2 = MAWtids[atid2];
289 >    
290 >    if ( mtid1 == -1 || mtid2 == -1) return 0.0;
291 >    else {      
292 >      MAWInteractionData mixer = MixingMap[mtid1][mtid2];
293        RealType R_e = mixer.Re;
294        RealType beta = mixer.beta;
295        // This value of the r corresponds to an energy about 1.48% of

Diff Legend

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