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 1549 by gezelter, Wed Apr 27 18:38:15 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).          
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 50 | Line 51 | namespace OpenMD {
51  
52   namespace OpenMD {
53  
54 <  MAW::MAW() : name_("MAW"), initialized_(false), forceField_(NULL),
54 <                   shiftedPot_(false), shiftedFrc_(false) {}
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;
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)) {
71        
72        if (nbt->isMAW()) {
73 <        pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
74 <        
75 <        GenericData* data = nbt->getPropertyByName("MAW");
76 <        if (data == NULL) {
77 <          sprintf( painCave.errMsg, "MAW::initialize could not find\n"
78 <                   "\tMAW parameters for %s - %s interaction.\n",
79 <                   atypes.first->getName().c_str(),
80 <                   atypes.second->getName().c_str());
81 <          painCave.severity = OPENMD_ERROR;
75 <          painCave.isFatal = 1;
76 <          simError();
73 >        keys = nbiTypes->getKeys(j);
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 <        
84 <        MAWData* mawData = dynamic_cast<MAWData*>(data);
85 <        if (mawData == NULL) {
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) {
93            sprintf( painCave.errMsg,
94 <                   "MAW::initialize could not convert GenericData to\n"
95 <                   "\tMAWData for %s - %s interaction.\n",
96 <                   atypes.first->getName().c_str(),
97 <                   atypes.second->getName().c_str());
94 >                   "MAW::initialize could not convert NonBondedInteractionType\n"
95 >                   "\tto MAWInteractionType for %s - %s interaction.\n",
96 >                   at1->getName().c_str(),
97 >                   at2->getName().c_str());
98            painCave.severity = OPENMD_ERROR;
99            painCave.isFatal = 1;
100            simError();          
101          }
102          
103 <        MAWParam mawParam = mawData->getData();
104 <
105 <        RealType De = mawParam.De;
106 <        RealType beta = mawParam.beta;
107 <        RealType Re = mawParam.Re;
96 <        RealType ca1 = mawParam.ca1;
97 <        RealType cb1 = mawParam.cb1;
103 >        RealType De = mit->getD();
104 >        RealType beta = mit->getBeta();
105 >        RealType Re = mit->getR();
106 >        RealType ca1 = mit->getCA1();
107 >        RealType cb1 = mit->getCB1();
108          
109 <        addExplicitInteraction(atypes.first, atypes.second,
109 >        addExplicitInteraction(at1, at2,
110                                 De, beta, Re, ca1, cb1);
111        }
112      }  
# Line 113 | 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();
130 –    
131 –    map<pair<AtomType*, AtomType*>, MAWInteractionData>::iterator it;
132 –    it = MixingMap.find(idat.atypes);
133 –    if (it != MixingMap.end()) {
134 –      MAWInteractionData mixer = (*it).second;
135 –      
136 –      RealType myPot = 0.0;
137 –      RealType myPotC = 0.0;
138 –      RealType myDeriv = 0.0;
139 –      RealType myDerivC = 0.0;
140 –      
141 –      RealType D_e = mixer.De;
142 –      RealType R_e = mixer.Re;
143 –      RealType beta = mixer.beta;
144 –      RealType ca1 = mixer.ca1;
145 –      RealType cb1 = mixer.cb1;
146  
147 <      bool j_is_Metal = idat.atypes.second->isMetal();
147 >    MAWInteractionData &mixer = MixingMap[MAWtids[idat.atid1]][MAWtids[idat.atid2]];
148  
149 <      Vector3d r;
150 <      RotMat3x3d Atrans;
151 <      if (j_is_Metal) {
152 <        // rotate the inter-particle separation into the two different
153 <        // body-fixed coordinate systems:
154 <        r = idat.A1 * idat.d;
155 <        Atrans = idat.A1.transpose();
156 <      } else {
157 <        // negative sign because this is the vector from j to i:      
158 <        r = -idat.A2 * idat.d;
159 <        Atrans = idat.A2.transpose();
160 <      }
161 <      
162 <      // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
163 <      
164 <      RealType expt     = -beta*(idat.rij - R_e);
165 <      RealType expfnc   = exp(expt);
166 <      RealType expfnc2  = expfnc*expfnc;
167 <      
168 <      RealType exptC = 0.0;
169 <      RealType expfncC = 0.0;
170 <      RealType expfnc2C = 0.0;
171 <
172 <      myPot  = D_e * (expfnc2  - 2.0 * expfnc);
173 <      myDeriv   = 2.0 * D_e * beta * (expfnc - expfnc2);
174 <      
175 <      if (MAW::shiftedPot_ || MAW::shiftedFrc_) {
176 <        exptC     = -beta*(idat.rcut - R_e);
177 <        expfncC   = exp(exptC);
178 <        expfnc2C  = expfncC*expfncC;
179 <      }
180 <      
181 <      if (MAW::shiftedPot_) {
182 <        myPotC = D_e * (expfnc2C - 2.0 * expfncC);
183 <        myDerivC = 0.0;
184 <      } else if (MAW::shiftedFrc_) {
185 <        myPotC = D_e * (expfnc2C - 2.0 * expfncC);
186 <        myDerivC  = 2.0 * D_e * beta * (expfnc2C - expfnc2C);
187 <        myPotC += myDerivC * (idat.rij - idat.rcut);
188 <      } else {
189 <        myPotC = 0.0;
190 <        myDerivC = 0.0;
191 <      }
192 <
193 <      RealType x = r.x();
194 <      RealType y = r.y();
195 <      RealType z = r.z();
196 <      RealType x2 = x * x;
197 <      RealType y2 = y * y;
198 <      RealType z2 = z * z;
199 <
200 <      RealType r3 = idat.r2 * idat.rij;
201 <      RealType r4 = idat.r2 * idat.r2;
202 <
203 <      // angular modulation of morse part of potential to approximate
204 <      // the squares of the two HOMO lone pair orbitals in water:
205 <      //********************** old form*************************
206 <      // s = 1 / (4 pi)
207 <      // ta1 = (s - pz)^2 = (1 - sqrt(3)*cos(theta))^2 / (4 pi)
208 <      // b1 = px^2 = 3 * (sin(theta)*cos(phi))^2 / (4 pi)  
209 <      //********************** old form*************************
210 <      // we'll leave out the 4 pi for now
211 <      
212 <      // new functional form just using the p orbitals.
213 <      // Vmorse(r)*[a*p_x + b p_z + (1-a-b)]
214 <      // which is
215 <      // Vmorse(r)*[a sin^2(theta) cos^2(phi) + b cos(theta) + (1-a-b)]
216 <      // Vmorse(r)*[a*x2/r2 + b*z/r + (1-a-b)]
217 <      
218 <      RealType Vmorse = (myPot - myPotC);
219 <      RealType Vang = ca1 * x2 / idat.r2 + cb1 * z / idat.rij + (0.8 - ca1 / 3.0);
220 <        
221 <      RealType pot_temp = idat.vdwMult * Vmorse * Vang;
222 <      idat.vpair += pot_temp;
223 <      idat.pot[0] += idat.sw * pot_temp;
224 <          
225 <      Vector3d dVmorsedr = (myDeriv - myDerivC) * Vector3d(x, y, z) / idat.rij;
149 >    RealType myPot = 0.0;
150 >    RealType myPotC = 0.0;
151 >    RealType myDeriv = 0.0;
152 >    RealType myDerivC = 0.0;
153      
154 <      Vector3d dVangdr = Vector3d(-2.0 * ca1 * x2 * x / r4 + 2.0 * ca1 * x / idat.r2 - cb1 * x * z / r3,
155 <                                  -2.0 * ca1 * x2 * y / r4                           - cb1 * z * y / r3,
156 <                                  -2.0 * ca1 * x2 * z / r4 + cb1 / idat.rij          - cb1 * z2  / r3);
157 <      
158 <      // chain rule to put these back on x, y, z
159 <
160 <      Vector3d dvdr = Vang * dVmorsedr + Vmorse * dVangdr;
161 <
162 <      // Torques for Vang using method of Price:
163 <      // S. L. Price, A. J. Stone, and M. Alderton, Mol. Phys. 52, 987 (1984).
164 <      
165 <      Vector3d dVangdu = Vector3d(cb1 * y / idat.rij,
166 <                                  2.0 * ca1 * x * z / idat.r2 - cb1 * x / idat.rij,
167 <                                  -2.0 * ca1 * y * x / idat.r2);
168 <
169 <      // do the torques first since they are easy:
170 <      // remember that these are still in the body fixed axes    
244 <
245 <      Vector3d trq = idat.vdwMult * Vmorse * dVangdu * idat.sw;
246 <
247 <      // go back to lab frame using transpose of rotation matrix:
248 <      
249 <      if (j_is_Metal) {
250 <        idat.t1 += Atrans * trq;
251 <      } else {
252 <        idat.t2 += Atrans * trq;
253 <      }
254 <
255 <      // Now, on to the forces (still in body frame of water)
256 <
257 <      Vector3d ftmp = idat.vdwMult * idat.sw * dvdr;
258 <
259 <      // rotate the terms back into the lab frame:
260 <      Vector3d flab;
261 <      if (j_is_Metal) {
262 <        flab = Atrans * ftmp;
263 <      } else {
264 <        flab = - Atrans * ftmp;
265 <      }
266 <      
267 <      idat.f1 += flab;
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      }
269 –    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