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 1536 by gezelter, Wed Jan 5 14:49:05 2011 UTC vs.
Revision 1583 by gezelter, Thu Jun 16 22:00:08 2011 UTC

# Line 50 | Line 50 | namespace OpenMD {
50  
51   namespace OpenMD {
52  
53 <  MAW::MAW() : name_("MAW"), initialized_(false), forceField_(NULL),
54 <                   shiftedPot_(false), shiftedFrc_(false) {}
53 >  MAW::MAW() : name_("MAW"), initialized_(false), forceField_(NULL) {}
54    
55    void MAW::initialize() {    
56  
# Line 128 | Line 127 | namespace OpenMD {
127  
128      if (!initialized_) initialize();
129      
131 –    pair<AtomType*, AtomType*> key = make_pair(idat.atype1, idat.atype2);
130      map<pair<AtomType*, AtomType*>, MAWInteractionData>::iterator it;
131 <    it = MixingMap.find(key);
131 >    it = MixingMap.find( idat.atypes );
132      if (it != MixingMap.end()) {
133        MAWInteractionData mixer = (*it).second;
134        
# Line 145 | Line 143 | namespace OpenMD {
143        RealType ca1 = mixer.ca1;
144        RealType cb1 = mixer.cb1;
145  
146 <      bool j_is_Metal = idat.atype2->isMetal();
146 >      bool j_is_Metal = idat.atypes.second->isMetal();
147  
148        Vector3d r;
149        RotMat3x3d Atrans;
150        if (j_is_Metal) {
151          // rotate the inter-particle separation into the two different
152          // body-fixed coordinate systems:
153 <        r = idat.A1 * idat.d;
154 <        Atrans = idat.A1.transpose();
153 >        r = *(idat.A1) * *(idat.d);
154 >        Atrans = idat.A1->transpose();
155        } else {
156          // negative sign because this is the vector from j to i:      
157 <        r = -idat.A2 * idat.d;
158 <        Atrans = idat.A2.transpose();
157 >        r = -*(idat.A2) * *(idat.d);
158 >        Atrans = idat.A2->transpose();
159        }
160        
161        // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
162        
163 <      RealType expt     = -beta*(idat.rij - R_e);
163 >      RealType expt     = -beta*( *(idat.rij) - R_e);
164        RealType expfnc   = exp(expt);
165        RealType expfnc2  = expfnc*expfnc;
166        
# Line 173 | Line 171 | namespace OpenMD {
171        myPot  = D_e * (expfnc2  - 2.0 * expfnc);
172        myDeriv   = 2.0 * D_e * beta * (expfnc - expfnc2);
173        
174 <      if (MAW::shiftedPot_ || MAW::shiftedFrc_) {
175 <        exptC     = -beta*(idat.rcut - R_e);
174 >      if (idat.shiftedPot || idat.shiftedForce) {
175 >        exptC     = -beta*( *(idat.rcut)  - R_e);
176          expfncC   = exp(exptC);
177          expfnc2C  = expfncC*expfncC;
178        }
179        
180 <      if (MAW::shiftedPot_) {
180 >      if (idat.shiftedPot) {
181          myPotC = D_e * (expfnc2C - 2.0 * expfncC);
182          myDerivC = 0.0;
183 <      } else if (MAW::shiftedFrc_) {
183 >      } else if (idat.shiftedForce) {
184          myPotC = D_e * (expfnc2C - 2.0 * expfncC);
185          myDerivC  = 2.0 * D_e * beta * (expfnc2C - expfnc2C);
186 <        myPotC += myDerivC * (idat.rij - idat.rcut);
186 >        myPotC += myDerivC * ( *(idat.rij)  -  *(idat.rcut) );
187        } else {
188          myPotC = 0.0;
189          myDerivC = 0.0;
# Line 198 | Line 196 | namespace OpenMD {
196        RealType y2 = y * y;
197        RealType z2 = z * z;
198  
199 <      RealType r3 = idat.r2 * idat.rij;
200 <      RealType r4 = idat.r2 * idat.r2;
199 >      RealType r3 = *(idat.r2) *  *(idat.rij) ;
200 >      RealType r4 = *(idat.r2) *  *(idat.r2);
201  
202        // angular modulation of morse part of potential to approximate
203        // the squares of the two HOMO lone pair orbitals in water:
# Line 217 | Line 215 | namespace OpenMD {
215        // Vmorse(r)*[a*x2/r2 + b*z/r + (1-a-b)]
216        
217        RealType Vmorse = (myPot - myPotC);
218 <      RealType Vang = ca1 * x2 / idat.r2 + cb1 * z / idat.rij + (0.8 - ca1 / 3.0);
219 <        
220 <      RealType pot_temp = idat.vdwMult * Vmorse * Vang;
221 <      idat.vpair[0] += pot_temp;
222 <      idat.pot[0] += idat.sw * pot_temp;
218 >      RealType Vang = ca1 * x2 / *(idat.r2) +
219 >        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))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
224            
225 <      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);
225 >      Vector3d dVmorsedr = (myDeriv - myDerivC) * Vector3d(x, y, z) /  *(idat.rij) ;
226        
227 +      Vector3d dVangdr = Vector3d(-2.0 * ca1 * x2 * x / r4 + 2.0 * ca1 * x / *(idat.r2) - cb1 * x * z / r3,
228 +                                  -2.0 * ca1 * x2 * y / r4 - cb1 * z * y / r3,
229 +                                  -2.0 * ca1 * x2 * z / r4 + cb1 /  *(idat.rij)           - cb1 * z2  / r3);
230 +      
231        // chain rule to put these back on x, y, z
232  
233        Vector3d dvdr = Vang * dVmorsedr + Vmorse * dVangdr;
# Line 236 | Line 235 | namespace OpenMD {
235        // Torques for Vang using method of Price:
236        // S. L. Price, A. J. Stone, and M. Alderton, Mol. Phys. 52, 987 (1984).
237        
238 <      Vector3d dVangdu = Vector3d(cb1 * y / idat.rij,
239 <                                  2.0 * ca1 * x * z / idat.r2 - cb1 * x / idat.rij,
240 <                                  -2.0 * ca1 * y * x / idat.r2);
238 >      Vector3d dVangdu = Vector3d(cb1 * y /  *(idat.rij) ,
239 >                                  2.0 * ca1 * x * z / *(idat.r2) - cb1 * x /  *(idat.rij),
240 >                                  -2.0 * ca1 * y * x / *(idat.r2));
241  
242        // do the torques first since they are easy:
243        // remember that these are still in the body fixed axes    
244  
245 <      Vector3d trq = idat.vdwMult * Vmorse * dVangdu * idat.sw;
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;
250 >        *(idat.t1) += Atrans * trq;
251        } else {
252 <        idat.t2 += Atrans * trq;
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;
257 >      Vector3d ftmp = *(idat.vdwMult) * *(idat.sw) * dvdr;
258  
259        // rotate the terms back into the lab frame:
260        Vector3d flab;
# Line 265 | Line 264 | namespace OpenMD {
264          flab = - Atrans * ftmp;
265        }
266        
267 <      idat.f1 += flab;
267 >      *(idat.f1) += flab;
268      }
269      return;
270      
271    }
272      
273 <  RealType MAW::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
273 >  RealType MAW::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
274      if (!initialized_) initialize();  
276 –    pair<AtomType*, AtomType*> key = make_pair(at1, at2);
275      map<pair<AtomType*, AtomType*>, MAWInteractionData>::iterator it;
276 <    it = MixingMap.find(key);
276 >    it = MixingMap.find(atypes);
277      if (it == MixingMap.end())
278        return 0.0;
279      else  {

Diff Legend

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