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 1478 by gezelter, Fri Jul 23 20:45:40 2010 UTC vs.
Revision 1505 by gezelter, Sun Oct 3 22:18:59 2010 UTC

# Line 45 | Line 45
45   #include <cmath>
46   #include "nonbonded/EAM.hpp"
47   #include "utils/simError.h"
48 + #include "types/NonBondedInteractionType.hpp"
49  
50  
51   namespace OpenMD {
52  
53 <  bool EAM::initialized_ = false;
54 <  ForceField* EAM::forceField_ = NULL;
54 <  std::map<int, AtomType*> EAM::EAMlist;
55 <  std::map<AtomType*, EAMAtomData> EAM::EAMMap;
56 <  std::map<std::pair<AtomType*, AtomType*>, EAMInteractionData> EAM::MixingMap;
53 >  EAM::EAM() : name_("EAM"), initialized_(false), forceField_(NULL),
54 >               mixMeth_(eamJohnson), eamRcut_(0.0) {}
55    
58 –  EAM* EAM::_instance = NULL;
59 –
60 –  EAM* EAM::Instance() {
61 –    if (!_instance) {
62 –      _instance = new EAM();
63 –    }
64 –    return _instance;
65 –  }
66 –  
56    EAMParam EAM::getEAMParam(AtomType* atomType) {
57      
58      // Do sanity checking on the AtomType we were passed before
# Line 107 | 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);
103      return cs;
104    }
105  
106 +  RealType EAM::getRcut(AtomType* atomType) {    
107 +    EAMParam eamParam = getEAMParam(atomType);
108 +    return eamParam.rcut;
109 +  }
110 +
111    CubicSpline* EAM::getRho(AtomType* atomType) {    
112      EAMParam eamParam = getEAMParam(atomType);
113      int nr = eamParam.nr;
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 135 | 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 152 | Line 146 | namespace OpenMD {
146  
147      // make the r grid:
148  
155 –    // 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 <    // use the smallest dr (finest grid) to build our grid:
156 >    rmax = max(rmax, eamParam2.rcut);
157 >    rmax = max(rmax, eamParam2.nr * eamParam2.dr);
158  
159 <    RealType dr = eamParam1.dr < eamParam2.dr ? eamParam1.dr : eamParam2.dr;
160 <    int nr = int(rcut/dr);
159 >    // use the smallest dr (finest grid) to build our grid:
160 >
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 175 | Line 176 | namespace OpenMD {
176  
177      for (int i = 1; i < rvals.size(); i++ ) {
178        r = rvals[i];
178 –      zi = z1->getValueAt(r);
179 –      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 190 | Line 198 | namespace OpenMD {
198    void EAM::initialize() {
199  
200      // set up the mixing method:
201 <    ForceFieldOptions ffo = forceField_->getForceFieldOptions();
202 <    string EAMMixMeth = toUpperCopy(ffo.getEAMMixingMethod());
203 <
201 >    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
202 >    string EAMMixMeth = fopts.getEAMMixingMethod();
203 >    toUpper(EAMMixMeth);
204 >  
205      if (EAMMixMeth == "JOHNSON")
206        mixMeth_ = eamJohnson;    
207      else if (EAMMixMeth == "DAW")
# Line 222 | 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 247 | Line 256 | namespace OpenMD {
256            simError();          
257          }
258          
259 <        EAMMix eamParam = eamData->getData();
259 >        EAMMixingParam eamParam = eamData->getData();
260  
261 <        vector<RealType> phiAB = eamParam.phiAB;
261 >        vector<RealType> phiAB = eamParam.phi;
262          RealType dr = eamParam.dr;
263          int nr = eamParam.nr;
264  
# Line 264 | Line 273 | namespace OpenMD {
273    void EAM::addType(AtomType* atomType){
274  
275      EAMAtomData eamAtomData;
276 <
276 >    
277      eamAtomData.rho = getRho(atomType);
278      eamAtomData.F = getF(atomType);
279      eamAtomData.Z = getZ(atomType);
# Line 273 | 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 288 | Line 297 | namespace OpenMD {
297      
298      // Now, iterate over all known types and add to the mixing map:
299      
300 <    std::map<int, AtomType*>::iterator it;
300 >    map<AtomType*, EAMAtomData>::iterator it;
301      for( it = EAMMap.begin(); it != EAMMap.end(); ++it) {
302        
303 <      AtomType* atype2 = (*it).second;
303 >      AtomType* atype2 = (*it).first;
304  
305        EAMInteractionData mixer;
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 319 | Line 328 | namespace OpenMD {
328  
329      EAMInteractionData mixer;
330      CubicSpline* cs = new CubicSpline();
331 <    vector<RealType> rvals;
331 >    vector<RealType> rVals;
332  
333 <    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
333 >    for (int i = 0; i < nr; i++) rVals.push_back(i * dr);
334  
335      cs->addPoints(rVals, phiVals);
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 338 | Line 347 | namespace OpenMD {
347      return;
348    }
349  
350 <  void EAM::calcDensity(AtomType* at1, AtomType* at2, Vector3d d,
351 <                        RealType rij, RealType r2, RealType rho_i_at_j,
343 <                        RealType rho_j_at_i) {
344 <
350 >  void EAM::calcDensity(DensityData ddat) {
351 >    
352      if (!initialized_) initialize();
353 +    
354 +    EAMAtomData data1 = EAMMap[ddat.atype1];
355 +    EAMAtomData data2 = EAMMap[ddat.atype2];
356  
357 <    EAMAtomData data1 = EAMMap[at1];
358 <    EAMAtomData data2 = EAMMap[at2];
357 >    if (ddat.rij < data1.rcut)
358 >      ddat.rho_i_at_j = data1.rho->getValueAt(ddat.rij);
359  
360 <    if (rij < data1.rcut) rho_i_at_j = data1.rho->getValueAt(rij);
361 <    if (rij < data2.rcut) rho_j_at_i = data2.rho->getValueAt(rij);
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,
356 <                           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,
371 <                      RealType rij, RealType r2, RealType sw,
372 <                      RealType &vpair, RealType &pot, Vector3d &f1,
373 <                      RealType rho1, RealType rho2, RealType dfrho1,
374 <                      RealType dfrho2, RealType fshift1, RealType fshift2) {
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 392 | 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 418 | 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 433 | Line 439 | namespace OpenMD {
439          break;
440  
441        case eamDaw:
442 <                
437 <        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 454 | 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 470 | 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 *d,
487 <                                 RealType *rij, RealType *r2,
488 <                                 RealType* rho_i_at_j, RealType* rho_j_at_i){
489 <    if (!initialized_) initialize();
493 >    RealType cut = 0.0;
494  
495 <    AtomType* atype1 = EAMlist[*atid1];
492 <    AtomType* atype2 = EAMlist[*atid2];
493 <    
494 <    Vector3d disp(d[0], d[1], d[2]);
495 >    map<AtomType*, EAMAtomData>::iterator it;
496  
497 <    calcDensity(atype1, atype2, disp, *rij, *r2, *rho_i_at_j, *rho_j_at_i);
497 >    it = EAMMap.find(at1);
498 >    if (it != EAMMap.end()) {
499 >      EAMAtomData data1 = (*it).second;
500 >      cut = data1.rcut;
501 >    }
502  
503 <    return;    
504 <  }
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 <  void EAM::calc_eam_preforce_Frho(int *atid1, RealType *rho, RealType *frho,
502 <                                   RealType *dfrhodrho) {
503 <
504 <    if (!initialized_) initialize();
505 <
506 <    AtomType* atype1 = EAMlist[*atid1];  
507 <
508 <    calcFunctional(atype1, *rho, *frho, *dfrhodrho);
509 <    
510 <    return;    
510 >    return cut;
511    }
512 –
513 –  void EAM::do_eam_pair(int *atid1, int *atid2, RealType *d, RealType *rij,
514 –                        RealType *r2, RealType *sw, RealType *vpair,
515 –                        RealType *pot, RealType *f1, RealType *rho1,
516 –                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
517 –                        RealType *fshift1, RealType *fshift2) {
518 –
519 –    if (!initialized_) initialize();
520 –    
521 –    AtomType* atype1 = EAMMap[*atid1];
522 –    AtomType* atype2 = EAMMap[*atid2];
523 –    
524 –    Vector3d disp(d[0], d[1], d[2]);
525 –    Vector3d frc(f1[0], f1[1], f1[2]);
526 –    
527 –    calcForce(atype1, atype2, disp, *rij, *r2, *sw, *vpair,  *pot, frc,
528 –              *rho1, *rho2, *dfrho1, *dfrho2, *fshift1, *fshift2);
529 –      
530 –    f1[0] = frc.x();
531 –    f1[1] = frc.y();
532 –    f1[2] = frc.z();
533 –
534 –    return;    
535 –  }
536 –  
537 –  void EAM::setCutoffEAM(RealType *thisRcut) {
538 –    eamRcut_ = thisRcut;
539 –  }
512   }
513  
542 – extern "C" {
543 –  
544 – #define fortranCalcDensity FC_FUNC(calc_eam_prepair_rho, CALC_EAM_PREPAIR_RHO)
545 – #define fortranCalcFunctional FC_FUNC(calc_eam_preforce_frho, CALC_EAM_PREFORCE_FRHO)
546 – #define fortranCalcForce FC_FUNC(do_eam_pair, DO_EAM_PAIR)
547 – #define fortranSetCutoffEAM FC_FUNC(setcutoffeam, SETCUTOFFEAM)
548 –  
549 –  RealType fortranCalcDensity(int *atid1, int *atid2, RealType *d,
550 –                              RealType *rij, RealType *r2,
551 –                              RealType *rho_i_at_j, RealType *rho_j_at_i) {
552 –
553 –    return OpenMD::EAM::Instance()->calc_eam_prepair_rho(*atid1, *atid2, *d,
554 –                                                         *rij, *r2,
555 –                                                         *rho_i_at_j,  
556 –                                                         *rho_j_at_i);
557 –  }
558 –  RealType fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho,
559 –                                 RealType *dfrhodrho) {  
560 –
561 –    return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(*atid1,
562 –                                                           *rho,
563 –                                                           *frho,
564 –                                                           *dfrhodrho);
565 –
566 –  }
567 –  void fortranSetEAMCutoff(RealType *rcut) {
568 –    return OpenMD::EAM::Instance()->setCutoffEAM(rcut);
569 –  }
570 –  void fortranDoEAMPair(int *atid1, int *atid2, RealType *d, RealType *rij,
571 –                        RealType *r2, RealType *sw, RealType *vpair,
572 –                        RealType *pot, RealType *f1, RealType *rho1,
573 –                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
574 –                        RealType *fshift1, RealType *fshift2){
575 –    
576 –    return OpenMD::EAM::Instance()->do_eam_pair(*atid1, *atid2, *d, *rij,
577 –                                                *r2, *sw,  *vpair,
578 –                                                *pot, *f1,  *rho1,
579 –                                                *rho2,  *dfrho1,  *dfrho2,
580 –                                                *fshift1,  *fshift2);
581 –  }
582 – }

Diff Legend

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