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 1479 by gezelter, Mon Jul 26 19:00:48 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 <  std::map<int, AtomType*> EAM::EAMlist;
58 <  std::map<AtomType*, EAMAtomData> EAM::EAMMap;
59 <  std::map<std::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 161 | Line 146 | namespace OpenMD {
146  
147      // make the r grid:
148  
164 –    // set rcut to be the smaller of the two atomic rcuts
149  
150 <    RealType rcut = eamParam1.rcut < eamParam2.rcut ?
151 <      eamParam1.rcut : eamParam2.rcut;
150 >    // we need phi out to the largest value we'll encounter in the radial space;
151 >    
152 >    RealType rmax = 0.0;
153 >    rmax = max(rmax, eamParam1.rcut);
154 >    rmax = max(rmax, eamParam1.nr * eamParam1.dr);
155  
156 +    rmax = max(rmax, eamParam2.rcut);
157 +    rmax = max(rmax, eamParam2.nr * eamParam2.dr);
158 +
159      // use the smallest dr (finest grid) to build our grid:
160  
161 <    RealType dr = eamParam1.dr < eamParam2.dr ? eamParam1.dr : eamParam2.dr;
162 <    int nr = int(rcut/dr);
161 >    RealType dr = min(eamParam1.dr, eamParam2.dr);
162 >
163 >    int nr = int(rmax/dr + 0.5);
164 >
165      vector<RealType> rvals;
166 <    for (int i = 0; i < nr; i++) rvals.push_back(i*dr);
166 >    for (int i = 0; i < nr; i++) rvals.push_back(RealType(i*dr));
167  
168      // construct the pair potential:
169  
# Line 184 | Line 176 | namespace OpenMD {
176  
177      for (int i = 1; i < rvals.size(); i++ ) {
178        r = rvals[i];
187 –      zi = z1->getValueAt(r);
188 –      zj = z2->getValueAt(r);
179  
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;
187 +
188        phi = 331.999296 * (zi * zj) / r;
189 +
190        phivals.push_back(phi);
191      }
192        
# Line 200 | Line 199 | namespace OpenMD {
199  
200      // set up the mixing method:
201      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
202 <    string EAMMixMeth = toUpperCopy(fopts.getEAMMixingMethod());
203 <
202 >    string EAMMixMeth = fopts.getEAMMixingMethod();
203 >    toUpper(EAMMixMeth);
204 >  
205      if (EAMMixMeth == "JOHNSON")
206        mixMeth_ = eamJohnson;    
207      else if (EAMMixMeth == "DAW")
# Line 231 | Line 231 | namespace OpenMD {
231        
232        if (nbt->isEAM()) {
233          
234 <        std::pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
234 >        pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
235          
236          GenericData* data = nbt->getPropertyByName("EAM");
237          if (data == NULL) {
# Line 282 | Line 282 | namespace OpenMD {
282      // add it to the map:
283      AtomTypeProperties atp = atomType->getATP();    
284  
285 <    std::pair<std::map<int,AtomType*>::iterator,bool> ret;    
286 <    ret = EAMlist.insert( std::pair<int, AtomType*>(atp.ident, atomType) );
285 >    pair<map<int,AtomType*>::iterator,bool> ret;    
286 >    ret = EAMlist.insert( pair<int, AtomType*>(atp.ident, atomType) );
287      if (ret.second == false) {
288        sprintf( painCave.errMsg,
289                 "EAM already had a previous entry with ident %d\n",
# Line 297 | Line 297 | namespace OpenMD {
297      
298      // Now, iterate over all known types and add to the mixing map:
299      
300 <    std::map<AtomType*, EAMAtomData>::iterator it;
300 >    map<AtomType*, EAMAtomData>::iterator it;
301      for( it = EAMMap.begin(); it != EAMMap.end(); ++it) {
302        
303        AtomType* atype2 = (*it).first;
# Line 306 | Line 306 | namespace OpenMD {
306        mixer.phi = getPhi(atomType, atype2);
307        mixer.explicitlySet = false;
308  
309 <      std::pair<AtomType*, AtomType*> key1, key2;
310 <      key1 = std::make_pair(atomType, atype2);
311 <      key2 = std::make_pair(atype2, atomType);
309 >      pair<AtomType*, AtomType*> key1, key2;
310 >      key1 = make_pair(atomType, atype2);
311 >      key2 = make_pair(atype2, atomType);
312        
313        MixingMap[key1] = mixer;
314        if (key2 != key1) {
# Line 336 | Line 336 | namespace OpenMD {
336      mixer.phi = cs;
337      mixer.explicitlySet = true;
338  
339 <    std::pair<AtomType*, AtomType*> key1, key2;
340 <    key1 = std::make_pair(atype1, atype2);
341 <    key2 = std::make_pair(atype2, atype1);
339 >    pair<AtomType*, AtomType*> key1, key2;
340 >    key1 = make_pair(atype1, atype2);
341 >    key2 = make_pair(atype2, atype1);
342      
343      MixingMap[key1] = mixer;
344      if (key2 != key1) {
# Line 347 | Line 347 | namespace OpenMD {
347      return;
348    }
349  
350 <  void EAM::calcDensity(AtomType* at1, AtomType* at2, const RealType rij,
351 <                        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,
364 <                           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,
379 <                      RealType rij, RealType r2, RealType sw,
380 <                      RealType &vpair, RealType &pot, Vector3d &f1,
381 <                      RealType rho_i, RealType rho_j,
382 <                      RealType dfrhodrho_i, RealType dfrhodrho_j,
383 <                      RealType &fshift_i, RealType &fshift_j) {
380 >  void EAM::calcForce(InteractionData idat) {
381  
382      if (!initialized_) initialize();
383 <    
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 401 | 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 427 | 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 442 | Line 439 | namespace OpenMD {
439          break;
440  
441        case eamDaw:
442 <                
446 <        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 463 | 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 479 | 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,
496 <                                 RealType* rho_i_at_j, RealType* rho_j_at_i){
493 >    RealType cut = 0.0;
494  
495 <    if (!initialized_) initialize();
499 <    
500 <    AtomType* atype1 = EAMlist[*atid1];
501 <    AtomType* atype2 = EAMlist[*atid2];
502 <    
503 <    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();
512 <
513 <    AtomType* atype1 = EAMlist[*atid1];  
514 <
515 <    calcFunctional(atype1, *rho, *frho, *dfrhodrho);
516 <    
517 <    return;    
510 >    return cut;
511    }
519 –  RealType EAM::getEAMcut(int *atid1) {
520 –
521 –    if (!initialized_) initialize();
522 –    
523 –    AtomType* atype1 = EAMlist[*atid1];  
524 –      
525 –    return getRcut(atype1);
526 –  }
527 –
528 –  void EAM::do_eam_pair(int *atid1, int *atid2, RealType *d, RealType *rij,
529 –                        RealType *r2, RealType *sw, RealType *vpair,
530 –                        RealType *pot, RealType *f1, RealType *rho1,
531 –                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
532 –                        RealType *fshift1, RealType *fshift2) {
533 –
534 –    if (!initialized_) initialize();
535 –    
536 –    AtomType* atype1 = EAMlist[*atid1];
537 –    AtomType* atype2 = EAMlist[*atid2];
538 –    
539 –    Vector3d disp(d[0], d[1], d[2]);
540 –    Vector3d frc(f1[0], f1[1], f1[2]);
541 –    
542 –    calcForce(atype1, atype2, disp, *rij, *r2, *sw, *vpair,  *pot, frc,
543 –              *rho1, *rho2, *dfrho1, *dfrho2, *fshift1, *fshift2);
544 –      
545 –    f1[0] = frc.x();
546 –    f1[1] = frc.y();
547 –    f1[2] = frc.z();
548 –
549 –    return;    
550 –  }
551 –  
552 –  void EAM::setCutoffEAM(RealType *thisRcut) {
553 –    eamRcut_ = *thisRcut;
554 –  }
512   }
513  
557 – extern "C" {
558 –  
559 – #define fortranCalcDensity FC_FUNC(calc_eam_prepair_rho, CALC_EAM_PREPAIR_RHO)
560 – #define fortranCalcFunctional FC_FUNC(calc_eam_preforce_frho, CALC_EAM_PREFORCE_FRHO)
561 – #define fortranCalcForce FC_FUNC(do_eam_pair, DO_EAM_PAIR)
562 – #define fortranSetCutoffEAM FC_FUNC(setcutoffeam, SETCUTOFFEAM)
563 – #define fortranGetEAMcut FC_FUNC(geteamcut, GETEAMCUT)
564 –
565 –  
566 –  void fortranCalcDensity(int *atid1, int *atid2, RealType *rij,
567 –                          RealType *rho_i_at_j, RealType *rho_j_at_i) {
568 –    
569 –    return OpenMD::EAM::Instance()->calc_eam_prepair_rho(atid1, atid2, rij,
570 –                                                         rho_i_at_j,  
571 –                                                         rho_j_at_i);
572 –  }
573 –  void fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho,
574 –                             RealType *dfrhodrho) {  
575 –    
576 –    return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(atid1, rho, frho,
577 –                                                           dfrhodrho);
578 –    
579 –  }
580 –  void fortranSetCutoffEAM(RealType *rcut) {
581 –    return OpenMD::EAM::Instance()->setCutoffEAM(rcut);
582 –  }
583 –  void fortranCalcForce(int *atid1, int *atid2, RealType *d, RealType *rij,
584 –                        RealType *r2, RealType *sw, RealType *vpair,
585 –                        RealType *pot, RealType *f1, RealType *rho1,
586 –                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
587 –                        RealType *fshift1, RealType *fshift2){
588 –    
589 –    return OpenMD::EAM::Instance()->do_eam_pair(atid1, atid2, d, rij,
590 –                                                r2, sw, vpair,
591 –                                                pot, f1, rho1,
592 –                                                rho2, dfrho1, dfrho2,
593 –                                                fshift1,  fshift2);
594 –  }
595 –  RealType fortranGetEAMcut(int* atid) {
596 –    return OpenMD::EAM::Instance()->getEAMcut(atid);
597 –  }
598 –
599 – }

Diff Legend

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