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 1480 by gezelter, Mon Jul 26 19:50:53 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 114 | Line 118 | namespace OpenMD {
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;
# 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());
204 <
202 >    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
203 >    std::string EAMMixMeth = fopts.getEAMMixingMethod();
204 >    toUpper(EAMMixMeth);
205 >  
206      if (EAMMixMeth == "JOHNSON")
207        mixMeth_ = eamJohnson;    
208      else if (EAMMixMeth == "DAW")
# Line 247 | Line 257 | namespace OpenMD {
257            simError();          
258          }
259          
260 <        EAMMix eamParam = eamData->getData();
260 >        EAMMixingParam eamParam = eamData->getData();
261  
262 <        vector<RealType> phiAB = eamParam.phiAB;
262 >        vector<RealType> phiAB = eamParam.phi;
263          RealType dr = eamParam.dr;
264          int nr = eamParam.nr;
265  
# Line 264 | Line 274 | namespace OpenMD {
274    void EAM::addType(AtomType* atomType){
275  
276      EAMAtomData eamAtomData;
277 <
277 >    
278      eamAtomData.rho = getRho(atomType);
279      eamAtomData.F = getF(atomType);
280      eamAtomData.Z = getZ(atomType);
# Line 288 | Line 298 | namespace OpenMD {
298      
299      // Now, iterate over all known types and add to the mixing map:
300      
301 <    std::map<int, AtomType*>::iterator it;
301 >    std::map<AtomType*, EAMAtomData>::iterator it;
302      for( it = EAMMap.begin(); it != EAMMap.end(); ++it) {
303        
304 <      AtomType* atype2 = (*it).second;
304 >      AtomType* atype2 = (*it).first;
305  
306        EAMInteractionData mixer;
307        mixer.phi = getPhi(atomType, atype2);
# Line 319 | Line 329 | namespace OpenMD {
329  
330      EAMInteractionData mixer;
331      CubicSpline* cs = new CubicSpline();
332 <    vector<RealType> rvals;
332 >    vector<RealType> rVals;
333  
334 <    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
334 >    for (int i = 0; i < nr; i++) rVals.push_back(i * dr);
335  
336      cs->addPoints(rVals, phiVals);
337      mixer.phi = cs;
# Line 338 | Line 348 | namespace OpenMD {
348      return;
349    }
350  
351 <  void EAM::calcDensity(AtomType* at1, AtomType* at2, Vector3d d,
352 <                        RealType rij, RealType r2, RealType rho_i_at_j,
353 <                        RealType rho_j_at_i) {
344 <
351 >  void EAM::calcDensity(AtomType* at1, AtomType* at2, const RealType rij,
352 >                        RealType &rho_i_at_j, RealType &rho_j_at_i) {
353 >    
354      if (!initialized_) initialize();
355 <
355 >    
356      EAMAtomData data1 = EAMMap[at1];
357      EAMAtomData data2 = EAMMap[at2];
358  
# Line 352 | Line 361 | namespace OpenMD {
361      return;
362    }
363  
364 <  void EAM::calcFunctional(AtomType* at1, RealType rho, RealType frho,
365 <                           RealType dfrhodrho) {
364 >  void EAM::calcFunctional(AtomType* at1, RealType rho, RealType &frho,
365 >                           RealType &dfrhodrho) {
366  
367      if (!initialized_) initialize();
368  
# Line 370 | Line 379 | namespace OpenMD {
379    void EAM::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
380                        RealType rij, RealType r2, RealType sw,
381                        RealType &vpair, RealType &pot, Vector3d &f1,
382 <                      RealType rho1, RealType rho2, RealType dfrho1,
383 <                      RealType dfrho2, RealType fshift1, RealType fshift2) {
382 >                      RealType rho_i, RealType rho_j,
383 >                      RealType dfrhodrho_i, RealType dfrhodrho_j,
384 >                      RealType &fshift_i, RealType &fshift_j) {
385  
386      if (!initialized_) initialize();
387      
# Line 483 | Line 493 | namespace OpenMD {
493    }
494  
495  
496 <  void EAM::calc_eam_prepair_rho(int *atid1, int *atid2, RealType *d,
487 <                                 RealType *rij, RealType *r2,
496 >  void EAM::calc_eam_prepair_rho(int *atid1, int *atid2, RealType *rij,
497                                   RealType* rho_i_at_j, RealType* rho_j_at_i){
489 –    if (!initialized_) initialize();
498  
499 +    if (!initialized_) initialize();
500 +    
501      AtomType* atype1 = EAMlist[*atid1];
502      AtomType* atype2 = EAMlist[*atid2];
503      
504 <    Vector3d disp(d[0], d[1], d[2]);
504 >    calcDensity(atype1, atype2, *rij, *rho_i_at_j, *rho_j_at_i);
505  
496 –    calcDensity(atype1, atype2, disp, *rij, *r2, *rho_i_at_j, *rho_j_at_i);
497 –
506      return;    
507    }
508  
# Line 508 | Line 516 | namespace OpenMD {
516      calcFunctional(atype1, *rho, *frho, *dfrhodrho);
517      
518      return;    
519 +  }
520 +  RealType EAM::getEAMcut(int *atid1) {
521 +
522 +    if (!initialized_) initialize();
523 +    
524 +    AtomType* atype1 = EAMlist[*atid1];  
525 +      
526 +    return getRcut(atype1);
527    }
528  
529    void EAM::do_eam_pair(int *atid1, int *atid2, RealType *d, RealType *rij,
# Line 518 | Line 534 | namespace OpenMD {
534  
535      if (!initialized_) initialize();
536      
537 <    AtomType* atype1 = EAMMap[*atid1];
538 <    AtomType* atype2 = EAMMap[*atid2];
537 >    AtomType* atype1 = EAMlist[*atid1];
538 >    AtomType* atype2 = EAMlist[*atid2];
539      
540      Vector3d disp(d[0], d[1], d[2]);
541      Vector3d frc(f1[0], f1[1], f1[2]);
# Line 535 | Line 551 | namespace OpenMD {
551    }
552    
553    void EAM::setCutoffEAM(RealType *thisRcut) {
554 <    eamRcut_ = thisRcut;
554 >    eamRcut_ = *thisRcut;
555    }
556   }
557  
# Line 545 | Line 561 | extern "C" {
561   #define fortranCalcFunctional FC_FUNC(calc_eam_preforce_frho, CALC_EAM_PREFORCE_FRHO)
562   #define fortranCalcForce FC_FUNC(do_eam_pair, DO_EAM_PAIR)
563   #define fortranSetCutoffEAM FC_FUNC(setcutoffeam, SETCUTOFFEAM)
564 <  
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) {
564 > #define fortranGetEAMcut FC_FUNC(geteamcut, GETEAMCUT)
565  
566 <    return OpenMD::EAM::Instance()->calc_eam_prepair_rho(*atid1, *atid2, *d,
567 <                                                         *rij, *r2,
568 <                                                         *rho_i_at_j,  
569 <                                                         *rho_j_at_i);
566 >  
567 >  void fortranCalcDensity(int *atid1, int *atid2, RealType *rij,
568 >                          RealType *rho_i_at_j, RealType *rho_j_at_i) {
569 >    
570 >    return OpenMD::EAM::Instance()->calc_eam_prepair_rho(atid1, atid2, rij,
571 >                                                         rho_i_at_j,  
572 >                                                         rho_j_at_i);
573    }
574 <  RealType fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho,
575 <                                 RealType *dfrhodrho) {  
576 <
577 <    return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(*atid1,
578 <                                                           *rho,
579 <                                                           *frho,
564 <                                                           *dfrhodrho);
565 <
574 >  void fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho,
575 >                             RealType *dfrhodrho) {  
576 >    
577 >    return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(atid1, rho, frho,
578 >                                                           dfrhodrho);
579 >    
580    }
581 <  void fortranSetEAMCutoff(RealType *rcut) {
581 >  void fortranSetCutoffEAM(RealType *rcut) {
582      return OpenMD::EAM::Instance()->setCutoffEAM(rcut);
583    }
584 <  void fortranDoEAMPair(int *atid1, int *atid2, RealType *d, RealType *rij,
584 >  void fortranCalcForce(int *atid1, int *atid2, RealType *d, RealType *rij,
585                          RealType *r2, RealType *sw, RealType *vpair,
586                          RealType *pot, RealType *f1, RealType *rho1,
587                          RealType *rho2, RealType *dfrho1, RealType *dfrho2,
588                          RealType *fshift1, RealType *fshift2){
589      
590 <    return OpenMD::EAM::Instance()->do_eam_pair(*atid1, *atid2, *d, *rij,
591 <                                                *r2, *sw,  *vpair,
592 <                                                *pot, *f1,  *rho1,
593 <                                                *rho2,  *dfrho1,  *dfrho2,
594 <                                                *fshift1,  *fshift2);
590 >    return OpenMD::EAM::Instance()->do_eam_pair(atid1, atid2, d, rij,
591 >                                                r2, sw, vpair,
592 >                                                pot, f1, rho1,
593 >                                                rho2, dfrho1, dfrho2,
594 >                                                fshift1,  fshift2);
595    }
596 +  RealType fortranGetEAMcut(int* atid) {
597 +    return OpenMD::EAM::Instance()->getEAMcut(atid);
598 +  }
599 +
600   }

Diff Legend

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