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 1554 by gezelter, Sat Apr 30 02:54:02 2011 UTC

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

Diff Legend

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