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 1505 by gezelter, Sun Oct 3 22:18:59 2010 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) {}
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 360 | Line 347 | namespace OpenMD {
347      return;
348    }
349  
350 <  void EAM::calcDensity(AtomType* at1, AtomType* at2, const RealType rij,
364 <                        RealType &rho_i_at_j, RealType &rho_j_at_i) {
350 >  void EAM::calcDensity(DensityData ddat) {
351      
352      if (!initialized_) initialize();
353      
354 <    EAMAtomData data1 = EAMMap[at1];
355 <    EAMAtomData data2 = EAMMap[at2];
354 >    EAMAtomData data1 = EAMMap[ddat.atype1];
355 >    EAMAtomData data2 = EAMMap[ddat.atype2];
356  
357 <    if (rij < data1.rcut) rho_i_at_j = data1.rho->getValueAt(rij);
358 <    if (rij < data2.rcut) rho_j_at_i = data2.rho->getValueAt(rij);
357 >    if (ddat.rij < data1.rcut)
358 >      ddat.rho_i_at_j = data1.rho->getValueAt(ddat.rij);
359 >
360 >    if (ddat.rij < data2.rcut)
361 >      ddat.rho_j_at_i = data2.rho->getValueAt(ddat.rij);
362 >
363      return;
364    }
365  
366 <  void EAM::calcFunctional(AtomType* at1, RealType rho, RealType &frho,
377 <                           RealType &dfrhodrho) {
366 >  void EAM::calcFunctional(FunctionalData fdat) {
367  
368      if (!initialized_) initialize();
369  
370 <    EAMAtomData data1 = EAMMap[at1];
370 >    EAMAtomData data1 = EAMMap[fdat.atype];
371          
372 <    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt(rho);
372 >    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt(fdat.rho);
373  
374 <    frho = result.first;
375 <    dfrhodrho = result.second;
374 >    fdat.frho = result.first;
375 >    fdat.dfrhodrho = result.second;
376      return;
377    }
378  
379  
380 <  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) {
380 >  void EAM::calcForce(InteractionData idat) {
381  
382      if (!initialized_) initialize();
383  
384      pair<RealType, RealType> res;
385      
386 <    if (rij < eamRcut_) {
386 >    if (idat.rij < eamRcut_) {
387  
388 <      EAMAtomData data1 = EAMMap[at1];
389 <      EAMAtomData data2 = EAMMap[at2];
388 >      EAMAtomData data1 = EAMMap[idat.atype1];
389 >      EAMAtomData data2 = EAMMap[idat.atype2];
390  
391        // get type-specific cutoff radii
392  
# Line 414 | Line 398 | namespace OpenMD {
398        RealType phab, dvpdr;
399        RealType drhoidr, drhojdr, dudr;
400        
401 <      if (rij < rci) {
402 <        res = data1.rho->getValueAndDerivativeAt(rij);
401 >      if (idat.rij < rci) {
402 >        res = data1.rho->getValueAndDerivativeAt(idat.rij);
403          rha = res.first;
404          drha = res.second;
405  
406 <        res = MixingMap[make_pair(at1, at1)].phi->getValueAndDerivativeAt(rij);
406 >        res = MixingMap[make_pair(idat.atype1, idat.atype1)].phi->getValueAndDerivativeAt(idat.rij);
407          pha = res.first;
408          dpha = res.second;
409        }
410  
411 <      if (rij < rcj) {
412 <        res = data2.rho->getValueAndDerivativeAt(rij);
411 >      if (idat.rij < rcj) {
412 >        res = data2.rho->getValueAndDerivativeAt(idat.rij);
413          rhb = res.first;
414          drhb = res.second;
415  
416 <        res = MixingMap[make_pair(at2, at2)].phi->getValueAndDerivativeAt(rij);
416 >        res = MixingMap[make_pair(idat.atype2, idat.atype2)].phi->getValueAndDerivativeAt(idat.rij);
417          phb = res.first;
418          dphb = res.second;
419        }
# Line 440 | Line 424 | namespace OpenMD {
424        switch(mixMeth_) {
425        case eamJohnson:
426        
427 <        if (rij < rci) {
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 (rij < rcj) {
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)));
# Line 455 | Line 439 | namespace OpenMD {
439          break;
440  
441        case eamDaw:
442 <        res = MixingMap[make_pair(at1,at2)].phi->getValueAndDerivativeAt(rij);
442 >        res = MixingMap[make_pair(idat.atype1,idat.atype2)].phi->getValueAndDerivativeAt(idat.rij);
443          phab = res.first;
444          dvpdr = res.second;
445  
# Line 475 | Line 459 | namespace OpenMD {
459        drhoidr = drha;
460        drhojdr = drhb;
461  
462 <      dudr = drhojdr*dfrhodrho_i + drhoidr*dfrhodrho_j + dvpdr;
462 >      dudr = drhojdr*idat.dfrho1 + drhoidr*idat.dfrho2 + dvpdr;
463  
464 <      f1 = d * dudr / rij;
464 >      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
# Line 491 | Line 475 | namespace OpenMD {
475        // Most of the particle_pot heavy lifting comes from the
476        // pair interaction, and will be handled by vpair.
477      
478 <      fshift_i = data1.F->getValueAt( rho_i - rhb );
479 <      fshift_j = data1.F->getValueAt( rho_j - rha );
478 >      idat.fshift1 = data1.F->getValueAt( idat.rho1 - rhb );
479 >      idat.fshift2 = data1.F->getValueAt( idat.rho2 - rha );
480  
481 <      pot += phab;
481 >      idat.pot += phab;
482  
483 <      vpair += phab;
483 >      idat.vpair += phab;
484      }
485  
486      return;
487      
488    }
489  
490 +  RealType EAM::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
491 +    if (!initialized_) initialize();  
492  
493 <  void EAM::calc_eam_prepair_rho(int *atid1, int *atid2, RealType *rij,
508 <                                 RealType* rho_i_at_j, RealType* rho_j_at_i){
493 >    RealType cut = 0.0;
494  
495 <    if (!initialized_) initialize();
511 <    
512 <    AtomType* atype1 = EAMlist[*atid1];
513 <    AtomType* atype2 = EAMlist[*atid2];
514 <    
515 <    calcDensity(atype1, atype2, *rij, *rho_i_at_j, *rho_j_at_i);
495 >    map<AtomType*, EAMAtomData>::iterator it;
496  
497 <    return;    
498 <  }
497 >    it = EAMMap.find(at1);
498 >    if (it != EAMMap.end()) {
499 >      EAMAtomData data1 = (*it).second;
500 >      cut = data1.rcut;
501 >    }
502  
503 <  void EAM::calc_eam_preforce_Frho(int *atid1, RealType *rho, RealType *frho,
504 <                                   RealType *dfrhodrho) {
503 >    it = EAMMap.find(at2);
504 >    if (it != EAMMap.end()) {
505 >      EAMAtomData data2 = (*it).second;
506 >      if (data2.rcut > cut)
507 >        cut = data2.rcut;
508 >    }
509  
510 <    if (!initialized_) initialize();
524 <
525 <    AtomType* atype1 = EAMlist[*atid1];  
526 <
527 <    calcFunctional(atype1, *rho, *frho, *dfrhodrho);
528 <    
529 <    return;    
510 >    return cut;
511    }
531 –  RealType EAM::getEAMcut(int *atid1) {
532 –
533 –    if (!initialized_) initialize();
534 –    
535 –    AtomType* atype1 = EAMlist[*atid1];  
536 –      
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();
547 –    
548 –    AtomType* atype1 = EAMlist[*atid1];
549 –    AtomType* atype2 = EAMlist[*atid2];
550 –    
551 –    Vector3d disp(d[0], d[1], d[2]);
552 –    Vector3d frc(f1[0], f1[1], f1[2]);
553 –    
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;    
562 –  }
563 –  
564 –  void EAM::setCutoffEAM(RealType *thisRcut) {
565 –    eamRcut_ = *thisRcut;
566 –  }
512   }
513  
569 – extern "C" {
570 –  
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)
576 –
577 –  
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 –  }
610 –
611 – }

Diff Legend

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