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 1481 by gezelter, Mon Jul 26 21:55:18 2010 UTC vs.
Revision 1586 by gezelter, Tue Jun 21 06:34:35 2011 UTC

# Line 50 | Line 50 | namespace OpenMD {
50  
51   namespace OpenMD {
52  
53 <  bool EAM::initialized_ = false;
54 <  RealType EAM::eamRcut_ = 0.0;
55 <  EAMMixingMethod EAM::mixMeth_ = eamJohnson;
56 <  ForceField* EAM::forceField_ = NULL;
57 <  map<int, AtomType*> EAM::EAMlist;
58 <  map<AtomType*, EAMAtomData> EAM::EAMMap;
59 <  map<pair<AtomType*, AtomType*>, EAMInteractionData> EAM::MixingMap;
60 <
53 >  EAM::EAM() : name_("EAM"), initialized_(false), forceField_(NULL),
54 >               mixMeth_(eamJohnson), eamRcut_(0.0), haveCutoffRadius_(false) {}
55    
62 –  EAM* EAM::_instance = NULL;
63 –
64 –  EAM* EAM::Instance() {
65 –    if (!_instance) {
66 –      _instance = new EAM();
67 –    }
68 –    return _instance;
69 –  }
70 –  
56    EAMParam EAM::getEAMParam(AtomType* atomType) {
57      
58      // Do sanity checking on the AtomType we were passed before
# Line 111 | Line 96 | namespace OpenMD {
96      RealType dr = eamParam.dr;
97      vector<RealType> rvals;
98      
99 <    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
99 >    for (int i = 0; i < nr; i++) rvals.push_back(RealType(i) * dr);
100        
101      CubicSpline* cs = new CubicSpline();
102      cs->addPoints(rvals, eamParam.Z);
# Line 129 | Line 114 | namespace OpenMD {
114      RealType dr = eamParam.dr;
115      vector<RealType> rvals;
116      
117 <    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
117 >    for (int i = 0; i < nr; i++) rvals.push_back(RealType(i) * dr);
118        
119      CubicSpline* cs = new CubicSpline();
120      cs->addPoints(rvals, eamParam.rho);
# Line 144 | Line 129 | namespace OpenMD {
129      vector<RealType> scaledF;
130      
131      for (int i = 0; i < nrho; i++) {
132 <      rhovals.push_back(i * drho);
132 >      rhovals.push_back(RealType(i) * drho);
133        scaledF.push_back( eamParam.F[i] * 23.06054 );
134      }
135        
136      CubicSpline* cs = new CubicSpline();
137 <    cs->addPoints(rhovals, eamParam.F);
137 >    cs->addPoints(rhovals, scaledF);
138      return cs;
139    }
140    
# Line 192 | Line 177 | namespace OpenMD {
177      for (int i = 1; i < rvals.size(); i++ ) {
178        r = rvals[i];
179  
180 <      // only use z(r) if we're inside this atoms cutoff radius, otherwise, we'll use zero for the charge.
181 <      // This effectively means that our phi grid goes out beyond the cutoff of the pair potential
180 >      // only use z(r) if we're inside this atom's cutoff radius,
181 >      // otherwise, we'll use zero for the charge.  This effectively
182 >      // means that our phi grid goes out beyond the cutoff of the
183 >      // pair potential
184  
185        zi = r <= eamParam1.rcut ? z1->getValueAt(r) : 0.0;
186        zj = r <= eamParam2.rcut ? z2->getValueAt(r) : 0.0;
# Line 208 | 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 360 | Line 352 | namespace OpenMD {
352      return;
353    }
354  
355 <  void EAM::calcDensity(AtomType* at1, AtomType* at2, const RealType rij,
364 <                        RealType &rho_i_at_j, RealType &rho_j_at_i) {
355 >  void EAM::calcDensity(InteractionData &idat) {
356      
357      if (!initialized_) initialize();
358      
359 <    EAMAtomData data1 = EAMMap[at1];
360 <    EAMAtomData data2 = EAMMap[at2];
361 <
362 <    if (rij < data1.rcut) rho_i_at_j = data1.rho->getValueAt(rij);
363 <    if (rij < data2.rcut) rho_j_at_i = data2.rho->getValueAt(rij);
364 <    return;
365 <  }
366 <
367 <  void EAM::calcFunctional(AtomType* at1, RealType rho, RealType &frho,
368 <                           RealType &dfrhodrho) {
369 <
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 >    
378      if (!initialized_) initialize();
379  
380 <    EAMAtomData data1 = EAMMap[at1];
380 >    EAMAtomData data1 = EAMMap[ sdat.atype ];
381          
382 <    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt(rho);
382 >    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt( *(sdat.rho) );
383  
384 <    frho = result.first;
385 <    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(AtomType* at1, AtomType* at2, Vector3d d,
392 <                      RealType rij, RealType r2, RealType sw,
393 <                      RealType &vpair, RealType &pot, Vector3d &f1,
394 <                      RealType rho_i, RealType rho_j,
395 <                      RealType dfrhodrho_i, RealType dfrhodrho_j,
396 <                      RealType &fshift_i, RealType &fshift_j) {
394 >  void EAM::calcForce(InteractionData &idat) {
395  
396      if (!initialized_) initialize();
397  
400 –    pair<RealType, RealType> res;
401 –    
402 –    if (rij < eamRcut_) {
398  
404 –      EAMAtomData data1 = EAMMap[at1];
405 –      EAMAtomData data2 = EAMMap[at2];
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 (rij < rci) {
435 <        res = data1.rho->getValueAndDerivativeAt(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(at1, at1)].phi->getValueAndDerivativeAt(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 <
448 <      if (rij < rcj) {
449 <        res = data2.rho->getValueAndDerivativeAt(rij);
450 <        rhb = res.first;
451 <        drhb = res.second;
452 <
453 <        res = MixingMap[make_pair(at2, at2)].phi->getValueAndDerivativeAt(rij);
433 <        phb = res.first;
434 <        dphb = res.second;
435 <      }
436 <
437 <      phab = 0.0;
438 <      dvpdr = 0.0;
439 <
440 <      switch(mixMeth_) {
441 <      case eamJohnson:
442 <      
443 <        if (rij < rci) {
444 <          phab = phab + 0.5 * (rhb / rha) * pha;
445 <          dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
446 <                               pha*((drhb/rha) - (rhb*drha/rha/rha)));
447 <        }
448 <
449 <        if (rij < rcj) {
450 <          phab = phab + 0.5 * (rha / rhb) * phb;
451 <          dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
452 <                                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
453 <        }
454 <
455 <        break;
456 <
457 <      case eamDaw:
458 <        res = MixingMap[make_pair(at1,at2)].phi->getValueAndDerivativeAt(rij);
459 <        phab = res.first;
460 <        dvpdr = res.second;
461 <
462 <        break;
463 <      case eamUnknown:
464 <      default:
465 <
466 <        sprintf(painCave.errMsg,
467 <                "EAM::calcForce hit a mixing method it doesn't know about!\n"
468 <                );
469 <        painCave.severity = OPENMD_ERROR;
470 <        painCave.isFatal = 1;
471 <        simError();        
472 <          
447 >      
448 >      
449 >      
450 >      if ( *(idat.rij) < rcj) {
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        }
455        
456 <      drhoidr = drha;
457 <      drhojdr = drhb;
458 <
459 <      dudr = drhojdr*dfrhodrho_i + drhoidr*dfrhodrho_j + dvpdr;
460 <
461 <      f1 = d * dudr / 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.
493 <    
494 <      fshift_i = data1.F->getValueAt( rho_i - rhb );
495 <      fshift_j = data1.F->getValueAt( rho_j - rha );
496 <
497 <      pot += phab;
498 <
499 <      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      }
501 –
502 –    return;
475      
476 <  }
477 <
506 <
507 <  void EAM::calc_eam_prepair_rho(int *atid1, int *atid2, RealType *rij,
508 <                                 RealType* rho_i_at_j, RealType* rho_j_at_i){
509 <
510 <    if (!initialized_) initialize();
476 >    drhoidr = drha;
477 >    drhojdr = drhb;
478      
479 <    AtomType* atype1 = EAMlist[*atid1];
513 <    AtomType* atype2 = EAMlist[*atid2];
479 >    dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
480      
481 <    calcDensity(atype1, atype2, *rij, *rho_i_at_j, *rho_j_at_i);
482 <
483 <    return;    
484 <  }
485 <
486 <  void EAM::calc_eam_preforce_Frho(int *atid1, RealType *rho, RealType *frho,
487 <                                   RealType *dfrhodrho) {
488 <
489 <    if (!initialized_) initialize();
490 <
525 <    AtomType* atype1 = EAMlist[*atid1];  
526 <
527 <    calcFunctional(atype1, *rho, *frho, *dfrhodrho);
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 <    return;    
493 <  }
531 <  RealType EAM::getEAMcut(int *atid1) {
532 <
533 <    if (!initialized_) initialize();
492 >    *(idat.particlePot1) += data2.F->getValueAt( *(idat.rho2) - rha )
493 >      - *(idat.frho2);
494      
495 <    AtomType* atype1 = EAMlist[*atid1];  
496 <      
537 <    return getRcut(atype1);
538 <  }
539 <
540 <  void EAM::do_eam_pair(int *atid1, int *atid2, RealType *d, RealType *rij,
541 <                        RealType *r2, RealType *sw, RealType *vpair,
542 <                        RealType *pot, RealType *f1, RealType *rho1,
543 <                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
544 <                        RealType *fshift1, RealType *fshift2) {
545 <
546 <    if (!initialized_) initialize();
495 >    *(idat.particlePot2) += data1.F->getValueAt( *(idat.rho1) - rhb)
496 >      - *(idat.frho1);
497      
498 <    AtomType* atype1 = EAMlist[*atid1];
549 <    AtomType* atype2 = EAMlist[*atid2];
498 >    (*(idat.pot))[METALLIC_FAMILY] += phab;
499      
500 <    Vector3d disp(d[0], d[1], d[2]);
501 <    Vector3d frc(f1[0], f1[1], f1[2]);
500 >    *(idat.vpair) += phab;
501 >  
502 >    return;
503      
554 –    calcForce(atype1, atype2, disp, *rij, *r2, *sw, *vpair,  *pot, frc,
555 –              *rho1, *rho2, *dfrho1, *dfrho2, *fshift1, *fshift2);
556 –      
557 –    f1[0] = frc.x();
558 –    f1[1] = frc.y();
559 –    f1[2] = frc.z();
560 –
561 –    return;    
504    }
563 –  
564 –  void EAM::setCutoffEAM(RealType *thisRcut) {
565 –    eamRcut_ = *thisRcut;
566 –  }
567 – }
505  
506 < extern "C" {
507 <  
571 < #define fortranCalcDensity FC_FUNC(calc_eam_prepair_rho, CALC_EAM_PREPAIR_RHO)
572 < #define fortranCalcFunctional FC_FUNC(calc_eam_preforce_frho, CALC_EAM_PREFORCE_FRHO)
573 < #define fortranCalcForce FC_FUNC(do_eam_pair, DO_EAM_PAIR)
574 < #define fortranSetCutoffEAM FC_FUNC(setcutoffeam, SETCUTOFFEAM)
575 < #define fortranGetEAMcut FC_FUNC(geteamcut, GETEAMCUT)
506 >  RealType EAM::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
507 >    if (!initialized_) initialize();  
508  
509 <  
578 <  void fortranCalcDensity(int *atid1, int *atid2, RealType *rij,
579 <                          RealType *rho_i_at_j, RealType *rho_j_at_i) {
580 <    
581 <    return OpenMD::EAM::Instance()->calc_eam_prepair_rho(atid1, atid2, rij,
582 <                                                         rho_i_at_j,  
583 <                                                         rho_j_at_i);
584 <  }
585 <  void fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho,
586 <                             RealType *dfrhodrho) {  
587 <    
588 <    return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(atid1, rho, frho,
589 <                                                           dfrhodrho);
590 <    
591 <  }
592 <  void fortranSetCutoffEAM(RealType *rcut) {
593 <    return OpenMD::EAM::Instance()->setCutoffEAM(rcut);
594 <  }
595 <  void fortranCalcForce(int *atid1, int *atid2, RealType *d, RealType *rij,
596 <                        RealType *r2, RealType *sw, RealType *vpair,
597 <                        RealType *pot, RealType *f1, RealType *rho1,
598 <                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
599 <                        RealType *fshift1, RealType *fshift2){
600 <    
601 <    return OpenMD::EAM::Instance()->do_eam_pair(atid1, atid2, d, rij,
602 <                                                r2, sw, vpair,
603 <                                                pot, f1, rho1,
604 <                                                rho2, dfrho1, dfrho2,
605 <                                                fshift1,  fshift2);
606 <  }
607 <  RealType fortranGetEAMcut(int* atid) {
608 <    return OpenMD::EAM::Instance()->getEAMcut(atid);
609 <  }
509 >    RealType cut = 0.0;
510  
511 +    map<AtomType*, EAMAtomData>::iterator it;
512 +
513 +    it = EAMMap.find(atypes.first);
514 +    if (it != EAMMap.end()) {
515 +      EAMAtomData data1 = (*it).second;
516 +      cut = data1.rcut;
517 +    }
518 +
519 +    it = EAMMap.find(atypes.second);
520 +    if (it != EAMMap.end()) {
521 +      EAMAtomData data2 = (*it).second;
522 +      if (data2.rcut > cut)
523 +        cut = data2.rcut;
524 +    }
525 +
526 +    return cut;
527 +  }
528   }
529 +

Diff Legend

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