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 1482 by gezelter, Tue Jul 27 14:55:15 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 +  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;
57 >  map<int, AtomType*> EAM::EAMlist;
58 >  map<AtomType*, EAMAtomData> EAM::EAMMap;
59 >  map<pair<AtomType*, AtomType*>, EAMInteractionData> EAM::MixingMap;
60 >
61    
62    EAM* EAM::_instance = NULL;
63  
# Line 107 | Line 111 | namespace OpenMD {
111      RealType dr = eamParam.dr;
112      vector<RealType> rvals;
113      
114 <    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
114 >    for (int i = 0; i < nr; i++) rvals.push_back(RealType(i) * dr);
115        
116      CubicSpline* cs = new CubicSpline();
117      cs->addPoints(rvals, eamParam.Z);
118      return cs;
119    }
120  
121 +  RealType EAM::getRcut(AtomType* atomType) {    
122 +    EAMParam eamParam = getEAMParam(atomType);
123 +    return eamParam.rcut;
124 +  }
125 +
126    CubicSpline* EAM::getRho(AtomType* atomType) {    
127      EAMParam eamParam = getEAMParam(atomType);
128      int nr = eamParam.nr;
129      RealType dr = eamParam.dr;
130      vector<RealType> rvals;
131      
132 <    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
132 >    for (int i = 0; i < nr; i++) rvals.push_back(RealType(i) * dr);
133        
134      CubicSpline* cs = new CubicSpline();
135      cs->addPoints(rvals, eamParam.rho);
# Line 135 | Line 144 | namespace OpenMD {
144      vector<RealType> scaledF;
145      
146      for (int i = 0; i < nrho; i++) {
147 <      rhovals.push_back(i * drho);
147 >      rhovals.push_back(RealType(i) * drho);
148        scaledF.push_back( eamParam.F[i] * 23.06054 );
149      }
150        
151      CubicSpline* cs = new CubicSpline();
152 <    cs->addPoints(rhovals, eamParam.F);
152 >    cs->addPoints(rhovals, scaledF);
153      return cs;
154    }
155    
# Line 152 | Line 161 | namespace OpenMD {
161  
162      // make the r grid:
163  
155 –    // set rcut to be the smaller of the two atomic rcuts
164  
165 <    RealType rcut = eamParam1.rcut < eamParam2.rcut ?
166 <      eamParam1.rcut : eamParam2.rcut;
165 >    // we need phi out to the largest value we'll encounter in the radial space;
166 >    
167 >    RealType rmax = 0.0;
168 >    rmax = max(rmax, eamParam1.rcut);
169 >    rmax = max(rmax, eamParam1.nr * eamParam1.dr);
170  
171 +    rmax = max(rmax, eamParam2.rcut);
172 +    rmax = max(rmax, eamParam2.nr * eamParam2.dr);
173 +
174      // use the smallest dr (finest grid) to build our grid:
175  
176 <    RealType dr = eamParam1.dr < eamParam2.dr ? eamParam1.dr : eamParam2.dr;
177 <    int nr = int(rcut/dr);
176 >    RealType dr = min(eamParam1.dr, eamParam2.dr);
177 >
178 >    int nr = int(rmax/dr + 0.5);
179 >
180      vector<RealType> rvals;
181 <    for (int i = 0; i < nr; i++) rvals.push_back(i*dr);
181 >    for (int i = 0; i < nr; i++) rvals.push_back(RealType(i*dr));
182  
183      // construct the pair potential:
184  
# Line 175 | Line 191 | namespace OpenMD {
191  
192      for (int i = 1; i < rvals.size(); i++ ) {
193        r = rvals[i];
178 –      zi = z1->getValueAt(r);
179 –      zj = z2->getValueAt(r);
194  
195 +      // only use z(r) if we're inside this atoms cutoff radius, otherwise, we'll use zero for the charge.
196 +      // This effectively means that our phi grid goes out beyond the cutoff of the pair potential
197 +
198 +      zi = r <= eamParam1.rcut ? z1->getValueAt(r) : 0.0;
199 +      zj = r <= eamParam2.rcut ? z2->getValueAt(r) : 0.0;
200 +
201        phi = 331.999296 * (zi * zj) / r;
202 +
203        phivals.push_back(phi);
204      }
205        
# Line 190 | Line 211 | namespace OpenMD {
211    void EAM::initialize() {
212  
213      // set up the mixing method:
214 <    ForceFieldOptions ffo = forceField_->getForceFieldOptions();
215 <    string EAMMixMeth = toUpperCopy(ffo.getEAMMixingMethod());
216 <
214 >    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
215 >    string EAMMixMeth = fopts.getEAMMixingMethod();
216 >    toUpper(EAMMixMeth);
217 >  
218      if (EAMMixMeth == "JOHNSON")
219        mixMeth_ = eamJohnson;    
220      else if (EAMMixMeth == "DAW")
# Line 222 | Line 244 | namespace OpenMD {
244        
245        if (nbt->isEAM()) {
246          
247 <        std::pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
247 >        pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
248          
249          GenericData* data = nbt->getPropertyByName("EAM");
250          if (data == NULL) {
# Line 247 | Line 269 | namespace OpenMD {
269            simError();          
270          }
271          
272 <        EAMMix eamParam = eamData->getData();
272 >        EAMMixingParam eamParam = eamData->getData();
273  
274 <        vector<RealType> phiAB = eamParam.phiAB;
274 >        vector<RealType> phiAB = eamParam.phi;
275          RealType dr = eamParam.dr;
276          int nr = eamParam.nr;
277  
# Line 264 | Line 286 | namespace OpenMD {
286    void EAM::addType(AtomType* atomType){
287  
288      EAMAtomData eamAtomData;
289 <
289 >    
290      eamAtomData.rho = getRho(atomType);
291      eamAtomData.F = getF(atomType);
292      eamAtomData.Z = getZ(atomType);
# Line 273 | Line 295 | namespace OpenMD {
295      // add it to the map:
296      AtomTypeProperties atp = atomType->getATP();    
297  
298 <    std::pair<std::map<int,AtomType*>::iterator,bool> ret;    
299 <    ret = EAMlist.insert( std::pair<int, AtomType*>(atp.ident, atomType) );
298 >    pair<map<int,AtomType*>::iterator,bool> ret;    
299 >    ret = EAMlist.insert( pair<int, AtomType*>(atp.ident, atomType) );
300      if (ret.second == false) {
301        sprintf( painCave.errMsg,
302                 "EAM already had a previous entry with ident %d\n",
# Line 288 | Line 310 | namespace OpenMD {
310      
311      // Now, iterate over all known types and add to the mixing map:
312      
313 <    std::map<int, AtomType*>::iterator it;
313 >    map<AtomType*, EAMAtomData>::iterator it;
314      for( it = EAMMap.begin(); it != EAMMap.end(); ++it) {
315        
316 <      AtomType* atype2 = (*it).second;
316 >      AtomType* atype2 = (*it).first;
317  
318        EAMInteractionData mixer;
319        mixer.phi = getPhi(atomType, atype2);
320        mixer.explicitlySet = false;
321  
322 <      std::pair<AtomType*, AtomType*> key1, key2;
323 <      key1 = std::make_pair(atomType, atype2);
324 <      key2 = std::make_pair(atype2, atomType);
322 >      pair<AtomType*, AtomType*> key1, key2;
323 >      key1 = make_pair(atomType, atype2);
324 >      key2 = make_pair(atype2, atomType);
325        
326        MixingMap[key1] = mixer;
327        if (key2 != key1) {
# Line 319 | Line 341 | namespace OpenMD {
341  
342      EAMInteractionData mixer;
343      CubicSpline* cs = new CubicSpline();
344 <    vector<RealType> rvals;
344 >    vector<RealType> rVals;
345  
346 <    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
346 >    for (int i = 0; i < nr; i++) rVals.push_back(i * dr);
347  
348      cs->addPoints(rVals, phiVals);
349      mixer.phi = cs;
350      mixer.explicitlySet = true;
351  
352 <    std::pair<AtomType*, AtomType*> key1, key2;
353 <    key1 = std::make_pair(atype1, atype2);
354 <    key2 = std::make_pair(atype2, atype1);
352 >    pair<AtomType*, AtomType*> key1, key2;
353 >    key1 = make_pair(atype1, atype2);
354 >    key2 = make_pair(atype2, atype1);
355      
356      MixingMap[key1] = mixer;
357      if (key2 != key1) {
# Line 338 | Line 360 | namespace OpenMD {
360      return;
361    }
362  
363 <  void EAM::calcDensity(AtomType* at1, AtomType* at2, Vector3d d,
364 <                        RealType rij, RealType r2, RealType rho_i_at_j,
365 <                        RealType rho_j_at_i) {
344 <
363 >  void EAM::calcDensity(AtomType* at1, AtomType* at2, const RealType rij,
364 >                        RealType &rho_i_at_j, RealType &rho_j_at_i) {
365 >    
366      if (!initialized_) initialize();
367 <
367 >    
368      EAMAtomData data1 = EAMMap[at1];
369      EAMAtomData data2 = EAMMap[at2];
370  
# Line 352 | Line 373 | namespace OpenMD {
373      return;
374    }
375  
376 <  void EAM::calcFunctional(AtomType* at1, RealType rho, RealType frho,
377 <                           RealType dfrhodrho) {
376 >  void EAM::calcFunctional(AtomType* at1, RealType rho, RealType &frho,
377 >                           RealType &dfrhodrho) {
378  
379      if (!initialized_) initialize();
380  
# Line 370 | Line 391 | namespace OpenMD {
391    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 rho1, RealType rho2, RealType dfrho1,
395 <                      RealType dfrho2, RealType fshift1, RealType fshift2) {
394 >                      RealType rho_i, RealType rho_j,
395 >                      RealType dfrhodrho_i, RealType dfrhodrho_j,
396 >                      RealType &fshift_i, RealType &fshift_j) {
397  
398      if (!initialized_) initialize();
399 <    
399 >
400      pair<RealType, RealType> res;
401      
402      if (rij < eamRcut_) {
# Line 433 | Line 455 | namespace OpenMD {
455          break;
456  
457        case eamDaw:
436 –                
458          res = MixingMap[make_pair(at1,at2)].phi->getValueAndDerivativeAt(rij);
459          phab = res.first;
460          dvpdr = res.second;
# Line 483 | Line 504 | namespace OpenMD {
504    }
505  
506  
507 <  void EAM::calc_eam_prepair_rho(int *atid1, int *atid2, RealType *d,
487 <                                 RealType *rij, RealType *r2,
507 >  void EAM::calc_eam_prepair_rho(int *atid1, int *atid2, RealType *rij,
508                                   RealType* rho_i_at_j, RealType* rho_j_at_i){
489 –    if (!initialized_) initialize();
509  
510 +    if (!initialized_) initialize();
511 +    
512      AtomType* atype1 = EAMlist[*atid1];
513      AtomType* atype2 = EAMlist[*atid2];
514      
515 <    Vector3d disp(d[0], d[1], d[2]);
515 >    calcDensity(atype1, atype2, *rij, *rho_i_at_j, *rho_j_at_i);
516  
496 –    calcDensity(atype1, atype2, disp, *rij, *r2, *rho_i_at_j, *rho_j_at_i);
497 –
517      return;    
518    }
519  
# Line 509 | Line 528 | namespace OpenMD {
528      
529      return;    
530    }
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,
# Line 518 | Line 545 | namespace OpenMD {
545  
546      if (!initialized_) initialize();
547      
548 <    AtomType* atype1 = EAMMap[*atid1];
549 <    AtomType* atype2 = EAMMap[*atid2];
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]);
# Line 535 | Line 562 | namespace OpenMD {
562    }
563    
564    void EAM::setCutoffEAM(RealType *thisRcut) {
565 <    eamRcut_ = thisRcut;
565 >    eamRcut_ = *thisRcut;
566    }
567   }
568  
# Line 545 | Line 572 | extern "C" {
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 <  
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) {
575 > #define fortranGetEAMcut FC_FUNC(geteamcut, GETEAMCUT)
576  
577 <    return OpenMD::EAM::Instance()->calc_eam_prepair_rho(*atid1, *atid2, *d,
578 <                                                         *rij, *r2,
579 <                                                         *rho_i_at_j,  
580 <                                                         *rho_j_at_i);
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 <  RealType fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho,
586 <                                 RealType *dfrhodrho) {  
587 <
588 <    return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(*atid1,
589 <                                                           *rho,
590 <                                                           *frho,
564 <                                                           *dfrhodrho);
565 <
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 fortranSetEAMCutoff(RealType *rcut) {
592 >  void fortranSetCutoffEAM(RealType *rcut) {
593      return OpenMD::EAM::Instance()->setCutoffEAM(rcut);
594    }
595 <  void fortranDoEAMPair(int *atid1, int *atid2, RealType *d, RealType *rij,
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);
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