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 1479 by gezelter, Mon Jul 26 19:00:48 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;
60 +
61    
62    EAM* EAM::_instance = NULL;
63  
# Line 112 | Line 116 | namespace OpenMD {
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) {    
# Line 190 | Line 199 | namespace OpenMD {
199    void EAM::initialize() {
200  
201      // set up the mixing method:
202 <    ForceFieldOptions ffo = forceField_->getForceFieldOptions();
203 <    string EAMMixMeth = toUpperCopy(ffo.getEAMMixingMethod());
202 >    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
203 >    string EAMMixMeth = toUpperCopy(fopts.getEAMMixingMethod());
204  
205      if (EAMMixMeth == "JOHNSON")
206        mixMeth_ = eamJohnson;    
# 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 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 >    std::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);
# 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;
# 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,
352 <                        RealType rho_j_at_i) {
344 <
350 >  void EAM::calcDensity(AtomType* at1, AtomType* at2, const RealType rij,
351 >                        RealType &rho_i_at_j, RealType &rho_j_at_i) {
352 >    
353      if (!initialized_) initialize();
354 <
354 >    
355      EAMAtomData data1 = EAMMap[at1];
356      EAMAtomData data2 = EAMMap[at2];
357  
# Line 352 | Line 360 | namespace OpenMD {
360      return;
361    }
362  
363 <  void EAM::calcFunctional(AtomType* at1, RealType rho, RealType frho,
364 <                           RealType dfrhodrho) {
363 >  void EAM::calcFunctional(AtomType* at1, RealType rho, RealType &frho,
364 >                           RealType &dfrhodrho) {
365  
366      if (!initialized_) initialize();
367  
# Line 370 | Line 378 | namespace OpenMD {
378    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 rho1, RealType rho2, RealType dfrho1,
382 <                      RealType dfrho2, RealType fshift1, RealType fshift2) {
381 >                      RealType rho_i, RealType rho_j,
382 >                      RealType dfrhodrho_i, RealType dfrhodrho_j,
383 >                      RealType &fshift_i, RealType &fshift_j) {
384  
385      if (!initialized_) initialize();
386      
# Line 483 | Line 492 | namespace OpenMD {
492    }
493  
494  
495 <  void EAM::calc_eam_prepair_rho(int *atid1, int *atid2, RealType *d,
487 <                                 RealType *rij, RealType *r2,
495 >  void EAM::calc_eam_prepair_rho(int *atid1, int *atid2, RealType *rij,
496                                   RealType* rho_i_at_j, RealType* rho_j_at_i){
489 –    if (!initialized_) initialize();
497  
498 +    if (!initialized_) initialize();
499 +    
500      AtomType* atype1 = EAMlist[*atid1];
501      AtomType* atype2 = EAMlist[*atid2];
502      
503 <    Vector3d disp(d[0], d[1], d[2]);
503 >    calcDensity(atype1, atype2, *rij, *rho_i_at_j, *rho_j_at_i);
504  
496 –    calcDensity(atype1, atype2, disp, *rij, *r2, *rho_i_at_j, *rho_j_at_i);
497 –
505      return;    
506    }
507  
# Line 509 | Line 516 | namespace OpenMD {
516      
517      return;    
518    }
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,
# Line 518 | Line 533 | namespace OpenMD {
533  
534      if (!initialized_) initialize();
535      
536 <    AtomType* atype1 = EAMMap[*atid1];
537 <    AtomType* atype2 = EAMMap[*atid2];
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]);
# Line 535 | Line 550 | namespace OpenMD {
550    }
551    
552    void EAM::setCutoffEAM(RealType *thisRcut) {
553 <    eamRcut_ = thisRcut;
553 >    eamRcut_ = *thisRcut;
554    }
555   }
556  
# Line 545 | Line 560 | extern "C" {
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 <  
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) {
563 > #define fortranGetEAMcut FC_FUNC(geteamcut, GETEAMCUT)
564  
565 <    return OpenMD::EAM::Instance()->calc_eam_prepair_rho(*atid1, *atid2, *d,
566 <                                                         *rij, *r2,
567 <                                                         *rho_i_at_j,  
568 <                                                         *rho_j_at_i);
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 <  RealType fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho,
574 <                                 RealType *dfrhodrho) {  
575 <
576 <    return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(*atid1,
577 <                                                           *rho,
578 <                                                           *frho,
564 <                                                           *dfrhodrho);
565 <
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 fortranSetEAMCutoff(RealType *rcut) {
580 >  void fortranSetCutoffEAM(RealType *rcut) {
581      return OpenMD::EAM::Instance()->setCutoffEAM(rcut);
582    }
583 <  void fortranDoEAMPair(int *atid1, int *atid2, RealType *d, RealType *rij,
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);
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