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

Comparing branches/development/src/nonbonded/Sticky.cpp (file contents):
Revision 1505 by gezelter, Sun Oct 3 22:18:59 2010 UTC vs.
Revision 1668 by gezelter, Fri Jan 6 19:03:05 2012 UTC

# Line 36 | Line 36
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).                        
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 186 | Line 187 | namespace OpenMD {
187     * idat structure.
188     */
189    
190 <  void Sticky::calcForce(InteractionData idat) {
190 >  void Sticky::calcForce(InteractionData &idat) {
191    
192      if (!initialized_) initialize();
193      
193 –    pair<AtomType*, AtomType*> key = make_pair(idat.atype1, idat.atype2);
194      map<pair<AtomType*, AtomType*>, StickyInteractionData>::iterator it;
195 <    it = MixingMap.find(key);
195 >    it = MixingMap.find(idat.atypes);
196      if (it != MixingMap.end()) {
197  
198        StickyInteractionData mixer = (*it).second;
# Line 207 | Line 207 | namespace OpenMD {
207        RealType rbig = mixer.rbig;
208        bool isPower = mixer.isPower;
209        
210 <      if (idat.rij <= rbig) {
210 >      if ( *(idat.rij) <= rbig) {
211          
212 <        RealType r3 = idat.r2 * idat.rij;
213 <        RealType r5 = r3 * idat.r2;
212 >        RealType r3 = *(idat.r2) * *(idat.rij);
213 >        RealType r5 = r3 * *(idat.r2);
214          
215 <        RotMat3x3d A1trans = idat.A1.transpose();
216 <        RotMat3x3d A2trans = idat.A2.transpose();
215 >        RotMat3x3d A1trans = idat.A1->transpose();
216 >        RotMat3x3d A2trans = idat.A2->transpose();
217          
218          // rotate the inter-particle separation into the two different
219          // body-fixed coordinate systems:
220          
221 <        Vector3d ri = idat.A1 * idat.d;
221 >        Vector3d ri = *(idat.A1) * *(idat.d);
222          
223          // negative sign because this is the vector from j to i:
224          
225 <        Vector3d rj = - idat.A2 * idat.d;
225 >        Vector3d rj = - *(idat.A2) * *(idat.d);
226          
227          RealType xi = ri.x();
228          RealType yi = ri.y();
# Line 247 | Line 247 | namespace OpenMD {
247          RealType sp = 0.0;
248          RealType dspdr = 0.0;
249          
250 <        if (idat.rij < ru) {
251 <          if (idat.rij < rl) {
250 >        if ( *(idat.rij) < ru) {
251 >          if ( *(idat.rij) < rl) {
252              s = 1.0;
253              dsdr = 0.0;
254            } else {          
255              // we are in the switching region
256              
257 <            pair<RealType, RealType> res = mixer.s->getValueAndDerivativeAt(idat.rij);
257 >            pair<RealType, RealType> res = mixer.s->getValueAndDerivativeAt(*(idat.rij));
258              s = res.first;
259              dsdr = res.second;
260            }
261          }
262          
263 <        if (idat.rij < rup) {
264 <          if (idat.rij < rlp) {
263 >        if (*(idat.rij) < rup) {
264 >          if ( *(idat.rij) < rlp) {
265              sp = 1.0;
266              dspdr = 0.0;
267            } else {
268              // we are in the switching region
269              
270 <            pair<RealType, RealType> res =mixer.sp->getValueAndDerivativeAt(idat.rij);
270 >            pair<RealType, RealType> res =mixer.sp->getValueAndDerivativeAt( *(idat.rij));
271              sp = res.first;
272              dspdr = res.second;
273            }
# Line 278 | Line 278 | namespace OpenMD {
278          RealType w = wi+wj;
279          
280          
281 <        RealType zif = zi/idat.rij - 0.6;
282 <        RealType zis = zi/idat.rij + 0.8;
281 >        RealType zif = zi/ *(idat.rij)  - 0.6;
282 >        RealType zis = zi/ *(idat.rij)  + 0.8;
283          
284 <        RealType zjf = zj/idat.rij - 0.6;
285 <        RealType zjs = zj/idat.rij + 0.8;
284 >        RealType zjf = zj/ *(idat.rij)  - 0.6;
285 >        RealType zjs = zj/ *(idat.rij)  + 0.8;
286          
287          RealType wip = zif*zif*zis*zis - w0;
288          RealType wjp = zjf*zjf*zjs*zjs - w0;
# Line 301 | Line 301 | namespace OpenMD {
301  
302          Vector3d dwip(-2.0*xi*zi*uglyi/r3,
303                        -2.0*yi*zi*uglyi/r3,
304 <                      2.0*(1.0/idat.rij - zi2/r3)*uglyi);
304 >                      2.0*(1.0/ *(idat.rij)  - zi2/r3)*uglyi);
305          
306          Vector3d dwjp(-2.0*xj*zj*uglyj/r3,
307                        -2.0*yj*zj*uglyj/r3,
308 <                      2.0*(1.0/idat.rij - zj2/r3)*uglyj);
308 >                      2.0*(1.0/ *(idat.rij)  - zj2/r3)*uglyj);
309          
310          Vector3d dwidu(4.0*(yi*zi2 + 0.5*yi*(xi2-yi2))/r3,
311                         4.0*(xi*zi2 - 0.5*xi*(xi2-yi2))/r3,
# Line 315 | Line 315 | namespace OpenMD {
315                         4.0*(xj*zj2 - 0.5*xj*(xj2-yj2))/r3,
316                         - 8.0*xj*yj*zj/r3);
317          
318 <        Vector3d dwipdu(2.0*yi*uglyi/idat.rij,
319 <                        -2.0*xi*uglyi/idat.rij,
318 >        Vector3d dwipdu(2.0*yi*uglyi/ *(idat.rij) ,
319 >                        -2.0*xi*uglyi/ *(idat.rij) ,
320                          0.0);
321          
322 <        Vector3d dwjpdu(2.0*yj*uglyj/idat.rij,
323 <                        -2.0*xj*uglyj/idat.rij,
322 >        Vector3d dwjpdu(2.0*yj*uglyj/ *(idat.rij) ,
323 >                        -2.0*xj*uglyj/ *(idat.rij) ,
324                          0.0);
325          
326          if (isPower) {
# Line 331 | Line 331 | namespace OpenMD {
331            // sticky power has no w' function:
332            w = frac1 * wi * wi2 + frac2*wi + frac1*wj*wj2 + frac2*wj + v0p;
333            wp = 0.0;
334 <          dwi = frac1*3.0*wi2*dwi + frac2*dwi;
335 <          dwj = frac1*3.0*wj2*dwi + frac2*dwi;
334 >          dwi = frac1*RealType(3.0)*wi2*dwi + frac2*dwi;
335 >          dwj = frac1*RealType(3.0)*wj2*dwi + frac2*dwi;
336            dwip = V3Zero;
337            dwjp = V3Zero;
338 <          dwidu = frac1*3.0*wi2*dwidu + frac2*dwidu;
339 <          dwidu = frac1*3.0*wj2*dwjdu + frac2*dwjdu;
338 >          dwidu = frac1*RealType(3.0)*wi2*dwidu + frac2*dwidu;
339 >          dwidu = frac1*RealType(3.0)*wj2*dwjdu + frac2*dwjdu;
340            dwipdu = V3Zero;
341            dwjpdu = V3Zero;
342            sp = 0.0;
343            dspdr = 0.0;
344          }
345          
346 <        idat.vpair += 0.5*(v0*s*w + v0p*sp*wp);
347 <        idat.pot += 0.5*(v0*s*w + v0p*sp*wp)*idat.sw;
346 >        *(idat.vpair) += RealType(0.5)*(v0*s*w + v0p*sp*wp);
347 >        (*(idat.pot))[HYDROGENBONDING_FAMILY] += RealType(0.5)*(v0*s*w + v0p*sp*wp)* *(idat.sw) ;
348          
349          // do the torques first since they are easy:
350          // remember that these are still in the body-fixed axes
351          
352 <        Vector3d ti = 0.5*idat.sw*(v0*s*dwidu + v0p*sp*dwipdu);
353 <        Vector3d tj = 0.5*idat.sw*(v0*s*dwjdu + v0p*sp*dwjpdu);
352 >        Vector3d ti = RealType(0.5)* *(idat.sw) *(v0*s*dwidu + v0p*sp*dwipdu);
353 >        Vector3d tj = RealType(0.5)* *(idat.sw) *(v0*s*dwjdu + v0p*sp*dwjpdu);
354          
355          // go back to lab frame using transpose of rotation matrix:
356          
357 <        idat.t1 += A1trans * ti;
358 <        idat.t2 += A2trans * tj;
357 >        *(idat.t1) += A1trans * ti;
358 >        *(idat.t2) += A2trans * tj;
359          
360          // Now, on to the forces:
361          
362          // first rotate the i terms back into the lab frame:
363          
364 <        Vector3d radcomi = (v0 * s * dwi + v0p * sp * dwip) * idat.sw;
365 <        Vector3d radcomj = (v0 * s * dwj + v0p * sp * dwjp) * idat.sw;
364 >        Vector3d radcomi = (v0 * s * dwi + v0p * sp * dwip) *  *(idat.sw);
365 >        Vector3d radcomj = (v0 * s * dwj + v0p * sp * dwjp) *  *(idat.sw);
366          
367          Vector3d fii = A1trans * radcomi;
368          Vector3d fjj = A2trans * radcomj;
369          
370          // now assemble these with the radial-only terms:
371          
372 <        idat.f1 += 0.5 * ((v0*dsdr*w + v0p*dspdr*wp) * idat.d /
373 <                          idat.rij + fii - fjj);
372 >        *(idat.f1) += RealType(0.5) * ((v0*dsdr*w + v0p*dspdr*wp) * *(idat.d) /
373 >                                       *(idat.rij)  + fii - fjj);
374          
375        }
376      }
# Line 378 | Line 378 | namespace OpenMD {
378      return;      
379    }
380  
381 <  RealType Sticky::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
381 >  RealType Sticky::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
382      if (!initialized_) initialize();  
383 –    pair<AtomType*, AtomType*> key = make_pair(at1, at2);
383      map<pair<AtomType*, AtomType*>, StickyInteractionData>::iterator it;
384 <    it = MixingMap.find(key);
384 >    it = MixingMap.find(atypes);
385      if (it == MixingMap.end())
386        return 0.0;
387      else  {

Diff Legend

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