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

Comparing:
branches/development/src/nonbonded/Sticky.cpp (file contents), Revision 1549 by gezelter, Wed Apr 27 18:38:15 2011 UTC vs.
branches/devel_omp/src/nonbonded/Sticky.cpp (file contents), Revision 1614 by mciznick, Tue Aug 23 20:55:51 2011 UTC

# Line 177 | Line 177 | namespace OpenMD {
177      }
178    }
179  
180 +  void Sticky::initForce() {
181 +        if (!initialized_) initialize();
182 +  }
183 +
184    /**
185     * This function does the sticky portion of the SSD potential
186     * [Chandra and Ichiye, Journal of Chemical Physics 111, 2701
# Line 206 | Line 210 | namespace OpenMD {
210        RealType rbig = mixer.rbig;
211        bool isPower = mixer.isPower;
212        
213 <      if (idat.rij <= rbig) {
213 >      if ( *(idat.rij) <= rbig) {
214          
215 <        RealType r3 = idat.r2 * idat.rij;
216 <        RealType r5 = r3 * idat.r2;
215 >        RealType r3 = *(idat.r2) * *(idat.rij);
216 >        RealType r5 = r3 * *(idat.r2);
217          
218 <        RotMat3x3d A1trans = idat.A1.transpose();
219 <        RotMat3x3d A2trans = idat.A2.transpose();
218 >        RotMat3x3d A1trans = idat.A1->transpose();
219 >        RotMat3x3d A2trans = idat.A2->transpose();
220          
221          // rotate the inter-particle separation into the two different
222          // body-fixed coordinate systems:
223          
224 <        Vector3d ri = idat.A1 * idat.d;
224 >        Vector3d ri = *(idat.A1) * *(idat.d);
225          
226          // negative sign because this is the vector from j to i:
227          
228 <        Vector3d rj = - idat.A2 * idat.d;
228 >        Vector3d rj = - *(idat.A2) * *(idat.d);
229          
230          RealType xi = ri.x();
231          RealType yi = ri.y();
# Line 246 | Line 250 | namespace OpenMD {
250          RealType sp = 0.0;
251          RealType dspdr = 0.0;
252          
253 <        if (idat.rij < ru) {
254 <          if (idat.rij < rl) {
253 >        if ( *(idat.rij) < ru) {
254 >          if ( *(idat.rij) < rl) {
255              s = 1.0;
256              dsdr = 0.0;
257            } else {          
258              // we are in the switching region
259              
260 <            pair<RealType, RealType> res = mixer.s->getValueAndDerivativeAt(idat.rij);
260 >            pair<RealType, RealType> res = mixer.s->getValueAndDerivativeAt(*(idat.rij));
261              s = res.first;
262              dsdr = res.second;
263            }
264          }
265          
266 <        if (idat.rij < rup) {
267 <          if (idat.rij < rlp) {
266 >        if (*(idat.rij) < rup) {
267 >          if ( *(idat.rij) < rlp) {
268              sp = 1.0;
269              dspdr = 0.0;
270            } else {
271              // we are in the switching region
272              
273 <            pair<RealType, RealType> res =mixer.sp->getValueAndDerivativeAt(idat.rij);
273 >            pair<RealType, RealType> res =mixer.sp->getValueAndDerivativeAt( *(idat.rij));
274              sp = res.first;
275              dspdr = res.second;
276            }
# Line 277 | Line 281 | namespace OpenMD {
281          RealType w = wi+wj;
282          
283          
284 <        RealType zif = zi/idat.rij - 0.6;
285 <        RealType zis = zi/idat.rij + 0.8;
284 >        RealType zif = zi/ *(idat.rij)  - 0.6;
285 >        RealType zis = zi/ *(idat.rij)  + 0.8;
286          
287 <        RealType zjf = zj/idat.rij - 0.6;
288 <        RealType zjs = zj/idat.rij + 0.8;
287 >        RealType zjf = zj/ *(idat.rij)  - 0.6;
288 >        RealType zjs = zj/ *(idat.rij)  + 0.8;
289          
290          RealType wip = zif*zif*zis*zis - w0;
291          RealType wjp = zjf*zjf*zjs*zjs - w0;
# Line 300 | Line 304 | namespace OpenMD {
304  
305          Vector3d dwip(-2.0*xi*zi*uglyi/r3,
306                        -2.0*yi*zi*uglyi/r3,
307 <                      2.0*(1.0/idat.rij - zi2/r3)*uglyi);
307 >                      2.0*(1.0/ *(idat.rij)  - zi2/r3)*uglyi);
308          
309          Vector3d dwjp(-2.0*xj*zj*uglyj/r3,
310                        -2.0*yj*zj*uglyj/r3,
311 <                      2.0*(1.0/idat.rij - zj2/r3)*uglyj);
311 >                      2.0*(1.0/ *(idat.rij)  - zj2/r3)*uglyj);
312          
313          Vector3d dwidu(4.0*(yi*zi2 + 0.5*yi*(xi2-yi2))/r3,
314                         4.0*(xi*zi2 - 0.5*xi*(xi2-yi2))/r3,
# Line 314 | Line 318 | namespace OpenMD {
318                         4.0*(xj*zj2 - 0.5*xj*(xj2-yj2))/r3,
319                         - 8.0*xj*yj*zj/r3);
320          
321 <        Vector3d dwipdu(2.0*yi*uglyi/idat.rij,
322 <                        -2.0*xi*uglyi/idat.rij,
321 >        Vector3d dwipdu(2.0*yi*uglyi/ *(idat.rij) ,
322 >                        -2.0*xi*uglyi/ *(idat.rij) ,
323                          0.0);
324          
325 <        Vector3d dwjpdu(2.0*yj*uglyj/idat.rij,
326 <                        -2.0*xj*uglyj/idat.rij,
325 >        Vector3d dwjpdu(2.0*yj*uglyj/ *(idat.rij) ,
326 >                        -2.0*xj*uglyj/ *(idat.rij) ,
327                          0.0);
328          
329          if (isPower) {
# Line 342 | Line 346 | namespace OpenMD {
346            dspdr = 0.0;
347          }
348          
349 <        idat.vpair += 0.5*(v0*s*w + v0p*sp*wp);
350 <        idat.pot[2] += 0.5*(v0*s*w + v0p*sp*wp)*idat.sw;
349 >        *(idat.vpair) += 0.5*(v0*s*w + v0p*sp*wp);
350 >        (*(idat.pot))[HYDROGENBONDING_FAMILY] += 0.5*(v0*s*w + v0p*sp*wp)* *(idat.sw) ;
351          
352          // do the torques first since they are easy:
353          // remember that these are still in the body-fixed axes
354          
355 <        Vector3d ti = 0.5*idat.sw*(v0*s*dwidu + v0p*sp*dwipdu);
356 <        Vector3d tj = 0.5*idat.sw*(v0*s*dwjdu + v0p*sp*dwjpdu);
355 >        Vector3d ti = 0.5* *(idat.sw) *(v0*s*dwidu + v0p*sp*dwipdu);
356 >        Vector3d tj = 0.5* *(idat.sw) *(v0*s*dwjdu + v0p*sp*dwjpdu);
357          
358          // go back to lab frame using transpose of rotation matrix:
359          
360 <        idat.t1 += A1trans * ti;
361 <        idat.t2 += A2trans * tj;
360 >        *(idat.t1) += A1trans * ti;
361 >        *(idat.t2) += A2trans * tj;
362          
363          // Now, on to the forces:
364          
365          // first rotate the i terms back into the lab frame:
366          
367 <        Vector3d radcomi = (v0 * s * dwi + v0p * sp * dwip) * idat.sw;
368 <        Vector3d radcomj = (v0 * s * dwj + v0p * sp * dwjp) * idat.sw;
367 >        Vector3d radcomi = (v0 * s * dwi + v0p * sp * dwip) *  *(idat.sw);
368 >        Vector3d radcomj = (v0 * s * dwj + v0p * sp * dwjp) *  *(idat.sw);
369          
370          Vector3d fii = A1trans * radcomi;
371          Vector3d fjj = A2trans * radcomj;
372          
373          // now assemble these with the radial-only terms:
374          
375 <        idat.f1 += 0.5 * ((v0*dsdr*w + v0p*dspdr*wp) * idat.d /
376 <                          idat.rij + fii - fjj);
375 >        *(idat.f1) += 0.5 * ((v0*dsdr*w + v0p*dspdr*wp) * *(idat.d) /
376 >                             *(idat.rij)  + fii - fjj);
377          
378        }
379      }

Diff Legend

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