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

Comparing:
branches/development/src/nonbonded/MAW.cpp (file contents), Revision 1545 by gezelter, Fri Apr 8 21:25:19 2011 UTC vs.
branches/devel_omp/src/nonbonded/MAW.cpp (file contents), Revision 1614 by mciznick, Tue Aug 23 20:55:51 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 124 | Line 123 | namespace OpenMD {
123      }    
124    }
125    
126 +  void MAW::initForce() {
127 +        if (!initialized_) initialize();
128 +  }
129 +
130    void MAW::calcForce(InteractionData &idat) {
131  
132      if (!initialized_) initialize();
133      
134      map<pair<AtomType*, AtomType*>, MAWInteractionData>::iterator it;
135 <    it = MixingMap.find(idat.atypes);
135 >    it = MixingMap.find( idat.atypes );
136      if (it != MixingMap.end()) {
137        MAWInteractionData mixer = (*it).second;
138        
# Line 151 | Line 154 | namespace OpenMD {
154        if (j_is_Metal) {
155          // rotate the inter-particle separation into the two different
156          // body-fixed coordinate systems:
157 <        r = idat.A1 * idat.d;
158 <        Atrans = idat.A1.transpose();
157 >        r = *(idat.A1) * *(idat.d);
158 >        Atrans = idat.A1->transpose();
159        } else {
160          // negative sign because this is the vector from j to i:      
161 <        r = -idat.A2 * idat.d;
162 <        Atrans = idat.A2.transpose();
161 >        r = -*(idat.A2) * *(idat.d);
162 >        Atrans = idat.A2->transpose();
163        }
164        
165        // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
166        
167 <      RealType expt     = -beta*(idat.rij - R_e);
167 >      RealType expt     = -beta*( *(idat.rij) - R_e);
168        RealType expfnc   = exp(expt);
169        RealType expfnc2  = expfnc*expfnc;
170        
# Line 172 | Line 175 | namespace OpenMD {
175        myPot  = D_e * (expfnc2  - 2.0 * expfnc);
176        myDeriv   = 2.0 * D_e * beta * (expfnc - expfnc2);
177        
178 <      if (MAW::shiftedPot_ || MAW::shiftedFrc_) {
179 <        exptC     = -beta*(idat.rcut - R_e);
178 >      if (idat.shiftedPot || idat.shiftedForce) {
179 >        exptC     = -beta*( *(idat.rcut)  - R_e);
180          expfncC   = exp(exptC);
181          expfnc2C  = expfncC*expfncC;
182        }
183        
184 <      if (MAW::shiftedPot_) {
184 >      if (idat.shiftedPot) {
185          myPotC = D_e * (expfnc2C - 2.0 * expfncC);
186          myDerivC = 0.0;
187 <      } else if (MAW::shiftedFrc_) {
187 >      } else if (idat.shiftedForce) {
188          myPotC = D_e * (expfnc2C - 2.0 * expfncC);
189          myDerivC  = 2.0 * D_e * beta * (expfnc2C - expfnc2C);
190 <        myPotC += myDerivC * (idat.rij - idat.rcut);
190 >        myPotC += myDerivC * ( *(idat.rij)  -  *(idat.rcut) );
191        } else {
192          myPotC = 0.0;
193          myDerivC = 0.0;
# Line 197 | Line 200 | namespace OpenMD {
200        RealType y2 = y * y;
201        RealType z2 = z * z;
202  
203 <      RealType r3 = idat.r2 * idat.rij;
204 <      RealType r4 = idat.r2 * idat.r2;
203 >      RealType r3 = *(idat.r2) *  *(idat.rij) ;
204 >      RealType r4 = *(idat.r2) *  *(idat.r2);
205  
206        // angular modulation of morse part of potential to approximate
207        // the squares of the two HOMO lone pair orbitals in water:
# Line 216 | Line 219 | namespace OpenMD {
219        // Vmorse(r)*[a*x2/r2 + b*z/r + (1-a-b)]
220        
221        RealType Vmorse = (myPot - myPotC);
222 <      RealType Vang = ca1 * x2 / idat.r2 + cb1 * z / idat.rij + (0.8 - ca1 / 3.0);
223 <        
224 <      RealType pot_temp = idat.vdwMult * Vmorse * Vang;
225 <      idat.vpair[0] += pot_temp;
226 <      idat.pot[0] += idat.sw * pot_temp;
222 >      RealType Vang = ca1 * x2 / *(idat.r2) +
223 >        cb1 * z /  *(idat.rij)  + (0.8 - ca1 / 3.0);
224 >      
225 >      RealType pot_temp = *(idat.vdwMult) * Vmorse * Vang;
226 >      *(idat.vpair) += pot_temp;
227 >      (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
228            
229 <      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);
229 >      Vector3d dVmorsedr = (myDeriv - myDerivC) * Vector3d(x, y, z) /  *(idat.rij) ;
230        
231 +      Vector3d dVangdr = Vector3d(-2.0 * ca1 * x2 * x / r4 + 2.0 * ca1 * x / *(idat.r2) - cb1 * x * z / r3,
232 +                                  -2.0 * ca1 * x2 * y / r4 - cb1 * z * y / r3,
233 +                                  -2.0 * ca1 * x2 * z / r4 + cb1 /  *(idat.rij)           - cb1 * z2  / r3);
234 +      
235        // chain rule to put these back on x, y, z
236  
237        Vector3d dvdr = Vang * dVmorsedr + Vmorse * dVangdr;
# Line 235 | Line 239 | namespace OpenMD {
239        // Torques for Vang using method of Price:
240        // S. L. Price, A. J. Stone, and M. Alderton, Mol. Phys. 52, 987 (1984).
241        
242 <      Vector3d dVangdu = Vector3d(cb1 * y / idat.rij,
243 <                                  2.0 * ca1 * x * z / idat.r2 - cb1 * x / idat.rij,
244 <                                  -2.0 * ca1 * y * x / idat.r2);
242 >      Vector3d dVangdu = Vector3d(cb1 * y /  *(idat.rij) ,
243 >                                  2.0 * ca1 * x * z / *(idat.r2) - cb1 * x /  *(idat.rij),
244 >                                  -2.0 * ca1 * y * x / *(idat.r2));
245  
246        // do the torques first since they are easy:
247        // remember that these are still in the body fixed axes    
248  
249 <      Vector3d trq = idat.vdwMult * Vmorse * dVangdu * idat.sw;
249 >      Vector3d trq = *(idat.vdwMult) * Vmorse * dVangdu * *(idat.sw);
250  
251        // go back to lab frame using transpose of rotation matrix:
252        
253        if (j_is_Metal) {
254 <        idat.t1 += Atrans * trq;
254 >        *(idat.t1) += Atrans * trq;
255        } else {
256 <        idat.t2 += Atrans * trq;
256 >        *(idat.t2) += Atrans * trq;
257        }
258  
259        // Now, on to the forces (still in body frame of water)
260  
261 <      Vector3d ftmp = idat.vdwMult * idat.sw * dvdr;
261 >      Vector3d ftmp = *(idat.vdwMult) * *(idat.sw) * dvdr;
262  
263        // rotate the terms back into the lab frame:
264        Vector3d flab;
# Line 264 | Line 268 | namespace OpenMD {
268          flab = - Atrans * ftmp;
269        }
270        
271 <      idat.f1 += flab;
271 >      *(idat.f1) += flab;
272      }
273      return;
274      

Diff Legend

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