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 1554 by gezelter, Sat Apr 30 02:54:02 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 <    }
314 >      RealType rho = mixer.phi->getValueAt( *(idat.rij) );
315 >      *(idat.rho1) += rho;
316 >      *(idat.rho2) += rho;
317 >    }
318      
319      return;
320    }
# 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) );
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 338 | Line 340 | namespace OpenMD {
340      
341      if (!initialized_) initialize();
342      
343 <    SCAtomData data1 = SCMap[idat.atypes->first];
344 <    SCAtomData data2 = SCMap[idat.atypes->second];
343 >    SCAtomData data1 = SCMap[idat.atypes.first];
344 >    SCAtomData data2 = SCMap[idat.atypes.second];
345  
346 <    SCInteractionData mixer = MixingMap[*(idat.atypes)];
346 >    SCInteractionData mixer = MixingMap[idat.atypes];
347  
348      RealType rcij = mixer.rCut;
349  
# Line 365 | Line 367 | namespace OpenMD {
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[METALLIC_FAMILY] += pot_temp;
385 >      (*(idat.pot))[METALLIC_FAMILY] += pot_temp;
386      }
387        
388      return;    

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines