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 1665 by gezelter, Tue Nov 22 20:38:56 2011 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 306 | Line 307 | namespace OpenMD {
307      
308      if (!initialized_) initialize();
309      
310 <    SCInteractionData mixer = MixingMap[ *(idat.atypes) ];
310 >    SCInteractionData mixer = MixingMap[ idat.atypes ];
311  
312      RealType rcij = mixer.rCut;
313  
314      if ( *(idat.rij)  < rcij) {
315 <      *(idat.rho_i_at_j) = mixer.phi->getValueAt( *(idat.rij) );
316 <      *(idat.rho_j_at_i) = *(idat.rho_i_at_j);
317 <    } else {
318 <      *(idat.rho_i_at_j) = 0.0;
318 <      *(idat.rho_j_at_i) = 0.0;
319 <    }
315 >      RealType rho = mixer.phi->getValueAt( *(idat.rij) );
316 >      *(idat.rho1) += rho;
317 >      *(idat.rho2) += rho;
318 >    }
319      
320      return;
321    }
# Line 326 | Line 325 | namespace OpenMD {
325      if (!initialized_) initialize();
326  
327      SCAtomData data1 = SCMap[sdat.atype];
328 <    
329 <    *(sdat.frho) = - data1.c * data1.epsilon * sqrt( *(sdat.rho) );
328 >  
329 >    RealType u = - data1.c * data1.epsilon * sqrt( *(sdat.rho) );
330 >    *(sdat.frho) = u;
331      *(sdat.dfrhodrho) = 0.5 * *(sdat.frho) / *(sdat.rho);
332 +
333 +    (*(sdat.pot))[METALLIC_FAMILY] += u;
334 +    *(sdat.particlePot) += u;
335      
336      return;
337    }
# Line 338 | Line 341 | namespace OpenMD {
341      
342      if (!initialized_) initialize();
343      
344 <    SCAtomData data1 = SCMap[idat.atypes->first];
345 <    SCAtomData data2 = SCMap[idat.atypes->second];
344 >    SCAtomData data1 = SCMap[idat.atypes.first];
345 >    SCAtomData data2 = SCMap[idat.atypes.second];
346  
347 <    SCInteractionData mixer = MixingMap[*(idat.atypes)];
347 >    SCInteractionData mixer = MixingMap[idat.atypes];
348  
349      RealType rcij = mixer.rCut;
350  
# Line 365 | Line 368 | namespace OpenMD {
368        
369        *(idat.f1) += *(idat.d) * dudr / *(idat.rij) ;
370          
371 <      // particle_pot is the difference between the full potential
372 <      // and the full potential without the presence of a particular
371 >      // particlePot is the difference between the full potential and
372 >      // the full potential without the presence of a particular
373        // particle (atom1).
374        //
375 <      // This reduces the density at other particle locations, so
376 <      // we need to recompute the density at atom2 assuming atom1
377 <      // didn't contribute.  This then requires recomputing the
378 <      // density functional for atom2 as well.
379 <      //
380 <      // Most of the particle_pot heavy lifting comes from the
381 <      // pair interaction, and will be handled by vpair.
375 >      // This reduces the density at other particle locations, so we
376 >      // need to recompute the density at atom2 assuming atom1 didn't
377 >      // contribute.  This then requires recomputing the density
378 >      // functional for atom2 as well.
379 >          
380 >      *(idat.particlePot1) -= data2.c * data2.epsilon *
381 >        sqrt( *(idat.rho2) - rhtmp) + *(idat.frho2);
382 >
383 >      *(idat.particlePot2) -= data1.c * data1.epsilon *
384 >        sqrt( *(idat.rho1) - rhtmp) + *(idat.frho1);
385        
386 <      *(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;
386 >      (*(idat.pot))[METALLIC_FAMILY] += pot_temp;
387      }
388        
389      return;    

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines