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

Comparing branches/development/src/nonbonded/MAW.cpp (file contents):
Revision 1532 by gezelter, Wed Dec 29 19:59:21 2010 UTC vs.
Revision 1582 by gezelter, Tue Jun 14 20:41:44 2011 UTC

# Line 124 | Line 124 | namespace OpenMD {
124      }    
125    }
126    
127 <  void MAW::calcForce(InteractionData idat) {
127 >  void MAW::calcForce(InteractionData &idat) {
128  
129      if (!initialized_) initialize();
130      
131 –    pair<AtomType*, AtomType*> key = make_pair(idat.atype1, idat.atype2);
131      map<pair<AtomType*, AtomType*>, MAWInteractionData>::iterator it;
132 <    it = MixingMap.find(key);
132 >    it = MixingMap.find( idat.atypes );
133      if (it != MixingMap.end()) {
134        MAWInteractionData mixer = (*it).second;
135        
# Line 145 | Line 144 | namespace OpenMD {
144        RealType ca1 = mixer.ca1;
145        RealType cb1 = mixer.cb1;
146  
147 <      bool j_is_Metal = idat.atype2->isMetal();
147 >      bool j_is_Metal = idat.atypes.second->isMetal();
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();
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();
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);
164 >      RealType expt     = -beta*( *(idat.rij) - R_e);
165        RealType expfnc   = exp(expt);
166        RealType expfnc2  = expfnc*expfnc;
167        
# Line 174 | Line 173 | namespace OpenMD {
173        myDeriv   = 2.0 * D_e * beta * (expfnc - expfnc2);
174        
175        if (MAW::shiftedPot_ || MAW::shiftedFrc_) {
176 <        exptC     = -beta*(idat.rcut - R_e);
176 >        exptC     = -beta*( *(idat.rcut)  - R_e);
177          expfncC   = exp(exptC);
178          expfnc2C  = expfncC*expfncC;
179        }
# Line 185 | Line 184 | namespace OpenMD {
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);
187 >        myPotC += myDerivC * ( *(idat.rij)  -  *(idat.rcut) );
188        } else {
189          myPotC = 0.0;
190          myDerivC = 0.0;
# Line 198 | Line 197 | namespace OpenMD {
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;
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:
# Line 217 | Line 216 | namespace OpenMD {
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 += idat.sw * pot_temp;
219 >      RealType Vang = ca1 * x2 / *(idat.r2) +
220 >        cb1 * z /  *(idat.rij)  + (0.8 - ca1 / 3.0);
221 >      
222 >      RealType pot_temp = *(idat.vdwMult) * Vmorse * Vang;
223 >      *(idat.vpair) += pot_temp;
224 >      (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
225            
226 <      Vector3d dVmorsedr = (myDeriv - myDerivC) * Vector3d(x, y, z) / idat.rij;
227 <    
228 <      Vector3d dVangdr = Vector3d(-2.0 * ca1 * x2 * x / r4 + 2.0 * ca1 * x / idat.r2 - cb1 * x * z / r3,
229 <                                  -2.0 * ca1 * x2 * y / r4                           - cb1 * z * y / r3,
230 <                                  -2.0 * ca1 * x2 * z / r4 + cb1 / idat.rij          - cb1 * z2  / r3);
226 >      Vector3d dVmorsedr = (myDeriv - myDerivC) * Vector3d(x, y, z) /  *(idat.rij) ;
227        
228 +      Vector3d dVangdr = Vector3d(-2.0 * ca1 * x2 * x / r4 + 2.0 * ca1 * x / *(idat.r2) - cb1 * x * z / r3,
229 +                                  -2.0 * ca1 * x2 * y / r4 - cb1 * z * y / r3,
230 +                                  -2.0 * ca1 * x2 * z / r4 + cb1 /  *(idat.rij)           - cb1 * z2  / r3);
231 +      
232        // chain rule to put these back on x, y, z
233  
234        Vector3d dvdr = Vang * dVmorsedr + Vmorse * dVangdr;
# Line 236 | Line 236 | namespace OpenMD {
236        // Torques for Vang using method of Price:
237        // S. L. Price, A. J. Stone, and M. Alderton, Mol. Phys. 52, 987 (1984).
238        
239 <      Vector3d dVangdu = Vector3d(cb1 * y / idat.rij,
240 <                                  2.0 * ca1 * x * z / idat.r2 - cb1 * x / idat.rij,
241 <                                  -2.0 * ca1 * y * x / idat.r2);
239 >      Vector3d dVangdu = Vector3d(cb1 * y /  *(idat.rij) ,
240 >                                  2.0 * ca1 * x * z / *(idat.r2) - cb1 * x /  *(idat.rij),
241 >                                  -2.0 * ca1 * y * x / *(idat.r2));
242  
243        // do the torques first since they are easy:
244        // remember that these are still in the body fixed axes    
245  
246 <      Vector3d trq = idat.vdwMult * Vmorse * dVangdu * idat.sw;
246 >      Vector3d trq = *(idat.vdwMult) * Vmorse * dVangdu * *(idat.sw);
247  
248        // go back to lab frame using transpose of rotation matrix:
249        
250        if (j_is_Metal) {
251 <        idat.t1 += Atrans * trq;
251 >        *(idat.t1) += Atrans * trq;
252        } else {
253 <        idat.t2 += Atrans * trq;
253 >        *(idat.t2) += Atrans * trq;
254        }
255  
256        // Now, on to the forces (still in body frame of water)
257  
258 <      Vector3d ftmp = idat.vdwMult * idat.sw * dvdr;
258 >      Vector3d ftmp = *(idat.vdwMult) * *(idat.sw) * dvdr;
259  
260        // rotate the terms back into the lab frame:
261        Vector3d flab;
# Line 265 | Line 265 | namespace OpenMD {
265          flab = - Atrans * ftmp;
266        }
267        
268 <      idat.f1 += flab;
268 >      *(idat.f1) += flab;
269      }
270      return;
271      
272    }
273      
274 <  RealType MAW::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
274 >  RealType MAW::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
275      if (!initialized_) initialize();  
276 –    pair<AtomType*, AtomType*> key = make_pair(at1, at2);
276      map<pair<AtomType*, AtomType*>, MAWInteractionData>::iterator it;
277 <    it = MixingMap.find(key);
277 >    it = MixingMap.find(atypes);
278      if (it == MixingMap.end())
279        return 0.0;
280      else  {

Diff Legend

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