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

Comparing branches/development/src/nonbonded/SC.cpp (file contents):
Revision 1549 by gezelter, Wed Apr 27 18:38:15 2011 UTC vs.
Revision 1586 by gezelter, Tue Jun 21 06:34:35 2011 UTC

# Line 306 | Line 306 | namespace OpenMD {
306      
307      if (!initialized_) initialize();
308      
309 <    SCInteractionData mixer = MixingMap[idat.atypes];
309 >    SCInteractionData mixer = MixingMap[ idat.atypes ];
310  
311      RealType rcij = mixer.rCut;
312  
313 <    if (idat.rij < rcij) {
314 <      idat.rho_i_at_j = mixer.phi->getValueAt(idat.rij);
315 <      idat.rho_j_at_i = idat.rho_i_at_j;
316 <    } else {
317 <      idat.rho_i_at_j = 0.0;
318 <      idat.rho_j_at_i = 0.0;
319 <    }
320 <
313 >    if ( *(idat.rij)  < rcij) {
314 >      RealType rho = mixer.phi->getValueAt( *(idat.rij) );
315 >      *(idat.rho1) += rho;
316 >      *(idat.rho2) += rho;
317 >    }
318 >    
319      return;
320    }
321  
# Line 326 | Line 324 | namespace OpenMD {
324      if (!initialized_) initialize();
325  
326      SCAtomData data1 = SCMap[sdat.atype];
327 <    
328 <    sdat.frho = - data1.c * data1.epsilon * sqrt(sdat.rho);
329 <    sdat.dfrhodrho = 0.5 * sdat.frho / sdat.rho;
327 >  
328 >    RealType u = - data1.c * data1.epsilon * sqrt( *(sdat.rho) );
329 >    *(sdat.frho) = u;
330 >    *(sdat.dfrhodrho) = 0.5 * *(sdat.frho) / *(sdat.rho);
331 >
332 >    (*(sdat.pot))[METALLIC_FAMILY] += u;
333 >    *(sdat.particlePot) += u;
334      
335      return;
336    }
# Line 345 | Line 347 | namespace OpenMD {
347  
348      RealType rcij = mixer.rCut;
349  
350 <    if (idat.rij < rcij) {
350 >    if ( *(idat.rij)  < rcij) {
351        RealType vcij = mixer.vCut;
352        
353        pair<RealType, RealType> res;
354        
355 <      res = mixer.phi->getValueAndDerivativeAt(idat.rij);
355 >      res = mixer.phi->getValueAndDerivativeAt( *(idat.rij) );
356        RealType rhtmp = res.first;
357        RealType drhodr = res.second;
358        
359 <      res = mixer.V->getValueAndDerivativeAt(idat.rij);
359 >      res = mixer.V->getValueAndDerivativeAt( *(idat.rij) );
360        RealType vptmp = res.first;
361        RealType dvpdr = res.second;
362        
363        RealType pot_temp = vptmp - vcij;
364 <      idat.vpair += pot_temp;
363 <      
364 <      RealType dudr = drhodr * (idat.dfrho1 + idat.dfrho2) + dvpdr;
364 >      *(idat.vpair) += pot_temp;
365        
366 <      idat.f1 += idat.d * dudr / idat.rij;
366 >      RealType dudr = drhodr * ( *(idat.dfrho1) + *(idat.dfrho2) ) + dvpdr;
367 >      
368 >      *(idat.f1) += *(idat.d) * dudr / *(idat.rij) ;
369          
370 <      // particle_pot is the difference between the full potential
371 <      // and the full potential without the presence of a particular
370 >      // particlePot is the difference between the full potential and
371 >      // the full potential without the presence of a particular
372        // particle (atom1).
373        //
374 <      // This reduces the density at other particle locations, so
375 <      // we need to recompute the density at atom2 assuming atom1
376 <      // didn't contribute.  This then requires recomputing the
377 <      // density functional for atom2 as well.
378 <      //
379 <      // Most of the particle_pot heavy lifting comes from the
380 <      // pair interaction, and will be handled by vpair.
374 >      // This reduces the density at other particle locations, so we
375 >      // need to recompute the density at atom2 assuming atom1 didn't
376 >      // contribute.  This then requires recomputing the density
377 >      // functional for atom2 as well.
378 >          
379 >      *(idat.particlePot1) -= data2.c * data2.epsilon *
380 >        sqrt( *(idat.rho2) - rhtmp) + *(idat.frho2);
381 >
382 >      *(idat.particlePot2) -= data1.c * data1.epsilon *
383 >        sqrt( *(idat.rho1) - rhtmp) + *(idat.frho1);
384        
385 <      idat.fshift1 = - data1.c * data1.epsilon * sqrt(idat.rho1 - rhtmp);
381 <      idat.fshift2 = - data2.c * data2.epsilon * sqrt(idat.rho2 - rhtmp);
382 <      
383 <      idat.pot[3] += pot_temp;
385 >      (*(idat.pot))[METALLIC_FAMILY] += pot_temp;
386      }
387        
388      return;    

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines