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

Comparing branches/development/src/nonbonded/EAM.cpp (file contents):
Revision 1554 by gezelter, Sat Apr 30 02:54:02 2011 UTC vs.
Revision 1582 by gezelter, Tue Jun 14 20:41:44 2011 UTC

# Line 351 | Line 351 | namespace OpenMD {
351      
352      if (!initialized_) initialize();
353      
354 <    EAMAtomData data1 = EAMMap[idat.atypes->first];
355 <    EAMAtomData data2 = EAMMap[idat.atypes->second];
354 >    EAMAtomData data1 = EAMMap[idat.atypes.first];
355 >    EAMAtomData data2 = EAMMap[idat.atypes.second];
356  
357      if ( *(idat.rij) < data1.rcut)
358 <      *(idat.rho_i_at_j) = data1.rho->getValueAt( *(idat.rij));
358 >      *(idat.rho1) += data1.rho->getValueAt( *(idat.rij));
359  
360 +
361      if ( *(idat.rij) < data2.rcut)
362 <      *(idat.rho_j_at_i) = data2.rho->getValueAt( *(idat.rij));
362 >      *(idat.rho2) += data2.rho->getValueAt( *(idat.rij));
363  
364      return;
365    }
# Line 373 | Line 374 | namespace OpenMD {
374  
375      *(sdat.frho) = result.first;
376      *(sdat.dfrhodrho) = result.second;
377 +
378 +    sdat.pot[METALLIC_FAMILY] += result.first;
379 +    *(sdat.particlePot) += result.first;
380 +
381      return;
382    }
383  
# Line 385 | Line 390 | namespace OpenMD {
390      
391      if ( *(idat.rij) < eamRcut_) {
392  
393 <      EAMAtomData data1 = EAMMap[idat.atypes->first];
394 <      EAMAtomData data2 = EAMMap[idat.atypes->second];
393 >      EAMAtomData data1 = EAMMap[idat.atypes.first];
394 >      EAMAtomData data2 = EAMMap[idat.atypes.second];
395  
396        // get type-specific cutoff radii
397  
398        RealType rci = data1.rcut;
399        RealType rcj = data2.rcut;
400        
401 <      RealType rha, drha, rhb, drhb;
402 <      RealType pha, dpha, phb, dphb;
403 <      RealType phab, dvpdr;
401 >      RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0);
402 >      RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0);
403 >      RealType phab(0.0), dvpdr(0.0);
404        RealType drhoidr, drhojdr, dudr;
405        
406        if ( *(idat.rij) < rci) {
# Line 403 | Line 408 | namespace OpenMD {
408          rha = res.first;
409          drha = res.second;
410  
411 <        res = MixingMap[make_pair(idat.atypes->first, idat.atypes->first)].phi->getValueAndDerivativeAt( *(idat.rij) );
411 >        res = MixingMap[make_pair(idat.atypes.first, idat.atypes.first)].phi->getValueAndDerivativeAt( *(idat.rij) );
412          pha = res.first;
413          dpha = res.second;
414        }
# Line 413 | Line 418 | namespace OpenMD {
418          rhb = res.first;
419          drhb = res.second;
420  
421 <        res = MixingMap[make_pair(idat.atypes->second, idat.atypes->second)].phi->getValueAndDerivativeAt( *(idat.rij) );
421 >        res = MixingMap[make_pair(idat.atypes.second, idat.atypes.second)].phi->getValueAndDerivativeAt( *(idat.rij) );
422          phb = res.first;
423          dphb = res.second;
424        }
425  
421 –      phab = 0.0;
422 –      dvpdr = 0.0;
423 –
426        switch(mixMeth_) {
427        case eamJohnson:
428        
# Line 429 | Line 431 | namespace OpenMD {
431            dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
432                                 pha*((drhb/rha) - (rhb*drha/rha/rha)));
433          }
434 +          
435 +          
436  
437          if ( *(idat.rij) < rcj) {
438            phab = phab + 0.5 * (rha / rhb) * phb;
# Line 439 | Line 443 | namespace OpenMD {
443          break;
444  
445        case eamDaw:
446 <        res = MixingMap[*(idat.atypes)].phi->getValueAndDerivativeAt( *(idat.rij));
446 >        res = MixingMap[idat.atypes].phi->getValueAndDerivativeAt( *(idat.rij));
447          phab = res.first;
448          dvpdr = res.second;
449  
# Line 463 | Line 467 | namespace OpenMD {
467  
468        *(idat.f1) = *(idat.d) * dudr / *(idat.rij);
469          
470 <      // particle_pot is the difference between the full potential
471 <      // and the full potential without the presence of a particular
470 >      // particlePot is the difference between the full potential and
471 >      // the full potential without the presence of a particular
472        // particle (atom1).
473        //
474 <      // This reduces the density at other particle locations, so
475 <      // we need to recompute the density at atom2 assuming atom1
476 <      // didn't contribute.  This then requires recomputing the
477 <      // density functional for atom2 as well.
474 <      //
475 <      // Most of the particle_pot heavy lifting comes from the
476 <      // pair interaction, and will be handled by vpair.
477 <    
478 <      *(idat.fshift1) = data1.F->getValueAt( *(idat.rho1) - rhb );
479 <      *(idat.fshift2) = data1.F->getValueAt( *(idat.rho2) - rha );
474 >      // This reduces the density at other particle locations, so we
475 >      // need to recompute the density at atom2 assuming atom1 didn't
476 >      // contribute.  This then requires recomputing the density
477 >      // functional for atom2 as well.
478  
479 <      idat.pot[METALLIC_FAMILY] += phab;
479 >      *(idat.particlePot1) += data2.F->getValueAt( *(idat.rho2) - rha )
480 >        - *(idat.frho2);
481  
482 +      *(idat.particlePot2) += data1.F->getValueAt( *(idat.rho1) - rhb)
483 +        - *(idat.frho1);
484 +
485 +      (*(idat.pot))[METALLIC_FAMILY] += phab;
486 +
487        *(idat.vpair) += phab;
488      }
489  

Diff Legend

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