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

Diff Legend

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