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

Diff Legend

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