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 1586 by gezelter, Tue Jun 21 06:34:35 2011 UTC

# Line 51 | Line 51 | namespace OpenMD {
51   namespace OpenMD {
52  
53    EAM::EAM() : name_("EAM"), initialized_(false), forceField_(NULL),
54 <               mixMeth_(eamJohnson), eamRcut_(0.0) {}
54 >               mixMeth_(eamJohnson), eamRcut_(0.0), haveCutoffRadius_(false) {}
55    
56    EAMParam EAM::getEAMParam(AtomType* atomType) {
57      
# Line 195 | Line 195 | namespace OpenMD {
195      return cs;
196    }
197  
198 +  void EAM::setCutoffRadius( RealType rCut ) {
199 +    eamRcut_ = rCut;
200 +    haveCutoffRadius_ = true;
201 +  }
202 +
203    void EAM::initialize() {
204  
205      // set up the mixing method:
# Line 351 | Line 356 | namespace OpenMD {
356      
357      if (!initialized_) initialize();
358      
359 <    EAMAtomData data1 = EAMMap[idat.atypes->first];
360 <    EAMAtomData data2 = EAMMap[idat.atypes->second];
361 <
362 <    if ( *(idat.rij) < data1.rcut)
363 <      *(idat.rho_i_at_j) = data1.rho->getValueAt( *(idat.rij));
364 <
365 <    if ( *(idat.rij) < data2.rcut)
366 <      *(idat.rho_j_at_i) = data2.rho->getValueAt( *(idat.rij));
367 <
368 <    return;
369 <  }
370 <
359 >    EAMAtomData data1 = EAMMap[idat.atypes.first];
360 >    EAMAtomData data2 = EAMMap[idat.atypes.second];
361 >    
362 >    if (haveCutoffRadius_)
363 >      if ( *(idat.rij) > eamRcut_) return;
364 >    
365 >    if ( *(idat.rij) < data1.rcut) {
366 >      *(idat.rho1) += data1.rho->getValueAt( *(idat.rij));
367 >      
368 >      
369 >      if ( *(idat.rij) < data2.rcut)
370 >        *(idat.rho2) += data2.rho->getValueAt( *(idat.rij));
371 >      
372 >      return;
373 >    }
374 >  }
375 >  
376    void EAM::calcFunctional(SelfData &sdat) {
377 <
377 >    
378      if (!initialized_) initialize();
379  
380      EAMAtomData data1 = EAMMap[ sdat.atype ];
# Line 373 | Line 383 | namespace OpenMD {
383  
384      *(sdat.frho) = result.first;
385      *(sdat.dfrhodrho) = result.second;
386 +
387 +    (*(sdat.pot))[METALLIC_FAMILY] += result.first;
388 +    *(sdat.particlePot) += result.first;
389 +
390      return;
391    }
392  
# Line 381 | Line 395 | namespace OpenMD {
395  
396      if (!initialized_) initialize();
397  
384 –    pair<RealType, RealType> res;
385 –    
386 –    if ( *(idat.rij) < eamRcut_) {
398  
388 –      EAMAtomData data1 = EAMMap[idat.atypes->first];
389 –      EAMAtomData data2 = EAMMap[idat.atypes->second];
399  
400 <      // get type-specific cutoff radii
401 <
402 <      RealType rci = data1.rcut;
403 <      RealType rcj = data2.rcut;
400 >    if (haveCutoffRadius_)
401 >      if ( *(idat.rij) > eamRcut_) return;
402 >  
403 >    pair<RealType, RealType> res;
404 >    
405 >    
406 >    EAMAtomData data1 = EAMMap[idat.atypes.first];
407 >    EAMAtomData data2 = EAMMap[idat.atypes.second];
408 >    
409 >    // get type-specific cutoff radii
410 >    
411 >    RealType rci = data1.rcut;
412 >    RealType rcj = data2.rcut;
413 >    
414 >    RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0);
415 >    RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0);
416 >    RealType phab(0.0), dvpdr(0.0);
417 >    RealType drhoidr, drhojdr, dudr;
418 >    
419 >    if ( *(idat.rij) < rci) {
420 >      res = data1.rho->getValueAndDerivativeAt( *(idat.rij));
421 >      rha = res.first;
422 >      drha = res.second;
423        
424 <      RealType rha, drha, rhb, drhb;
425 <      RealType pha, dpha, phb, dphb;
426 <      RealType phab, dvpdr;
427 <      RealType drhoidr, drhojdr, dudr;
424 >      res = MixingMap[make_pair(idat.atypes.first, idat.atypes.first)].phi->getValueAndDerivativeAt( *(idat.rij) );
425 >      pha = res.first;
426 >      dpha = res.second;
427 >    }
428 >    
429 >    if ( *(idat.rij) < rcj) {
430 >      res = data2.rho->getValueAndDerivativeAt( *(idat.rij) );
431 >      rhb = res.first;
432 >      drhb = res.second;
433        
434 <      if ( *(idat.rij) < rci) {
435 <        res = data1.rho->getValueAndDerivativeAt( *(idat.rij));
436 <        rha = res.first;
437 <        drha = res.second;
434 >      res = MixingMap[make_pair(idat.atypes.second, idat.atypes.second)].phi->getValueAndDerivativeAt( *(idat.rij) );
435 >      phb = res.first;
436 >      dphb = res.second;
437 >    }
438  
439 <        res = MixingMap[make_pair(idat.atypes->first, idat.atypes->first)].phi->getValueAndDerivativeAt( *(idat.rij) );
440 <        pha = res.first;
441 <        dpha = res.second;
439 >    switch(mixMeth_) {
440 >    case eamJohnson:
441 >      
442 >      if ( *(idat.rij) < rci) {
443 >        phab = phab + 0.5 * (rhb / rha) * pha;
444 >        dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
445 >                             pha*((drhb/rha) - (rhb*drha/rha/rha)));
446        }
447 <
447 >      
448 >      
449 >      
450        if ( *(idat.rij) < rcj) {
451 <        res = data2.rho->getValueAndDerivativeAt( *(idat.rij) );
452 <        rhb = res.first;
453 <        drhb = res.second;
415 <
416 <        res = MixingMap[make_pair(idat.atypes->second, idat.atypes->second)].phi->getValueAndDerivativeAt( *(idat.rij) );
417 <        phb = res.first;
418 <        dphb = res.second;
451 >        phab = phab + 0.5 * (rha / rhb) * phb;
452 >        dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
453 >                               phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
454        }
420 –
421 –      phab = 0.0;
422 –      dvpdr = 0.0;
423 –
424 –      switch(mixMeth_) {
425 –      case eamJohnson:
426 –      
427 –        if ( *(idat.rij) < rci) {
428 –          phab = phab + 0.5 * (rhb / rha) * pha;
429 –          dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
430 –                               pha*((drhb/rha) - (rhb*drha/rha/rha)));
431 –        }
432 –
433 –        if ( *(idat.rij) < rcj) {
434 –          phab = phab + 0.5 * (rha / rhb) * phb;
435 –          dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
436 –                                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
437 –        }
438 –
439 –        break;
440 –
441 –      case eamDaw:
442 –        res = MixingMap[*(idat.atypes)].phi->getValueAndDerivativeAt( *(idat.rij));
443 –        phab = res.first;
444 –        dvpdr = res.second;
445 –
446 –        break;
447 –      case eamUnknown:
448 –      default:
449 –
450 –        sprintf(painCave.errMsg,
451 –                "EAM::calcForce hit a mixing method it doesn't know about!\n"
452 –                );
453 –        painCave.severity = OPENMD_ERROR;
454 –        painCave.isFatal = 1;
455 –        simError();        
456 –          
457 –      }
455        
456 <      drhoidr = drha;
457 <      drhojdr = drhb;
458 <
459 <      dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
460 <
461 <      *(idat.f1) = *(idat.d) * dudr / *(idat.rij);
462 <        
463 <      // particle_pot is the difference between the full potential
464 <      // and the full potential without the presence of a particular
465 <      // particle (atom1).
466 <      //
467 <      // This reduces the density at other particle locations, so
468 <      // we need to recompute the density at atom2 assuming atom1
469 <      // didn't contribute.  This then requires recomputing the
470 <      // density functional for atom2 as well.
471 <      //
472 <      // Most of the particle_pot heavy lifting comes from the
473 <      // 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 );
480 <
481 <      idat.pot[METALLIC_FAMILY] += phab;
482 <
483 <      *(idat.vpair) += phab;
456 >      break;
457 >      
458 >    case eamDaw:
459 >      res = MixingMap[idat.atypes].phi->getValueAndDerivativeAt( *(idat.rij));
460 >      phab = res.first;
461 >      dvpdr = res.second;
462 >      
463 >      break;
464 >    case eamUnknown:
465 >    default:
466 >      
467 >      sprintf(painCave.errMsg,
468 >              "EAM::calcForce hit a mixing method it doesn't know about!\n"
469 >              );
470 >      painCave.severity = OPENMD_ERROR;
471 >      painCave.isFatal = 1;
472 >      simError();        
473 >      
474      }
475 <
475 >    
476 >    drhoidr = drha;
477 >    drhojdr = drhb;
478 >    
479 >    dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
480 >    
481 >    *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
482 >        
483 >    // particlePot is the difference between the full potential and
484 >    // the full potential without the presence of a particular
485 >    // particle (atom1).
486 >    //
487 >    // This reduces the density at other particle locations, so we
488 >    // need to recompute the density at atom2 assuming atom1 didn't
489 >    // contribute.  This then requires recomputing the density
490 >    // functional for atom2 as well.
491 >    
492 >    *(idat.particlePot1) += data2.F->getValueAt( *(idat.rho2) - rha )
493 >      - *(idat.frho2);
494 >    
495 >    *(idat.particlePot2) += data1.F->getValueAt( *(idat.rho1) - rhb)
496 >      - *(idat.frho1);
497 >    
498 >    (*(idat.pot))[METALLIC_FAMILY] += phab;
499 >    
500 >    *(idat.vpair) += phab;
501 >  
502      return;
503      
504    }

Diff Legend

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