ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/nonbonded/EAM.cpp
(Generate patch)

Comparing:
branches/development/src/nonbonded/EAM.cpp (file contents), Revision 1480 by gezelter, Mon Jul 26 19:50:53 2010 UTC vs.
trunk/src/nonbonded/EAM.cpp (file contents), Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

# Line 35 | Line 35
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include <stdio.h>
# Line 50 | Line 51 | namespace OpenMD {
51  
52   namespace OpenMD {
53  
54 <  bool EAM::initialized_ = false;
55 <  RealType EAM::eamRcut_ = 0.0;
56 <  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 <
54 >  EAM::EAM() : initialized_(false), haveCutoffRadius_(false),
55 >               forceField_(NULL), eamRcut_(0.0), mixMeth_(eamJohnson),
56 >               name_("EAM") {}
57    
58 <  EAM* EAM::_instance = NULL;
58 >  CubicSpline* EAM::getPhi(AtomType* atomType1, AtomType* atomType2) {  
59 >    EAMAdapter ea1 = EAMAdapter(atomType1);
60 >    EAMAdapter ea2 = EAMAdapter(atomType2);
61 >    CubicSpline* z1 = ea1.getZ();
62 >    CubicSpline* z2 = ea2.getZ();
63  
64 <  EAM* EAM::Instance() {
65 <    if (!_instance) {
66 <      _instance = new EAM();
67 <    }
68 <    return _instance;
69 <  }
70 <  
71 <  EAMParam EAM::getEAMParam(AtomType* atomType) {
72 <    
73 <    // Do sanity checking on the AtomType we were passed before
74 <    // building any data structures:
75 <    if (!atomType->isEAM()) {
76 <      sprintf( painCave.errMsg,
77 <               "EAM::getEAMParam was passed an atomType (%s) that does not\n"
78 <               "\tappear to be an embedded atom method (EAM) atom.\n",
79 <               atomType->getName().c_str());
80 <      painCave.severity = OPENMD_ERROR;
81 <      painCave.isFatal = 1;
82 <      simError();
83 <    }
84 <    
85 <    GenericData* data = atomType->getPropertyByName("EAM");
86 <    if (data == NULL) {
87 <      sprintf( painCave.errMsg, "EAM::getEAMParam could not find EAM\n"
88 <               "\tparameters for atomType %s.\n",
89 <               atomType->getName().c_str());
90 <      painCave.severity = OPENMD_ERROR;
91 <      painCave.isFatal = 1;
92 <      simError();
93 <    }
94 <    
95 <    EAMParamGenericData* eamData = dynamic_cast<EAMParamGenericData*>(data);
96 <    if (eamData == NULL) {
97 <      sprintf( painCave.errMsg,
98 <               "EAM::getEAMParam could not convert GenericData to EAMParam for\n"
99 <               "\tatom type %s\n", atomType->getName().c_str());
100 <      painCave.severity = OPENMD_ERROR;
101 <      painCave.isFatal = 1;
102 <      simError();          
103 <    }
104 <    
105 <    return eamData->getData();
106 <  }
64 >    // Thise prefactors convert the charge-charge interactions into
65 >    // kcal / mol all were computed assuming distances are measured in
66 >    // angstroms Charge-Charge, assuming charges are measured in
67 >    // electrons.  Matches value in Electrostatics.cpp
68 >    pre11_ = 332.0637778;
69  
70 <  CubicSpline* EAM::getZ(AtomType* atomType) {    
109 <    EAMParam eamParam = getEAMParam(atomType);
110 <    int nr = eamParam.nr;
111 <    RealType dr = eamParam.dr;
112 <    vector<RealType> rvals;
113 <    
114 <    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
115 <      
116 <    CubicSpline* cs = new CubicSpline();
117 <    cs->addPoints(rvals, eamParam.Z);
118 <    return cs;
119 <  }
70 >    // make the r grid:
71  
72 <  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;
72 >    // we need phi out to the largest value we'll encounter in the radial space;
73      
74 <    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
75 <      
76 <    CubicSpline* cs = new CubicSpline();
135 <    cs->addPoints(rvals, eamParam.rho);
136 <    return cs;
137 <  }
74 >    RealType rmax = 0.0;
75 >    rmax = max(rmax, ea1.getRcut());
76 >    rmax = max(rmax, ea1.getNr() * ea1.getDr());
77  
78 <  CubicSpline* EAM::getF(AtomType* atomType) {    
79 <    EAMParam eamParam = getEAMParam(atomType);
141 <    int nrho = eamParam.nrho;
142 <    RealType drho = eamParam.drho;
143 <    vector<RealType> rhovals;
144 <    vector<RealType> scaledF;
145 <    
146 <    for (int i = 0; i < nrho; i++) {
147 <      rhovals.push_back(i * drho);
148 <      scaledF.push_back( eamParam.F[i] * 23.06054 );
149 <    }
150 <      
151 <    CubicSpline* cs = new CubicSpline();
152 <    cs->addPoints(rhovals, eamParam.F);
153 <    return cs;
154 <  }
155 <  
156 <  CubicSpline* EAM::getPhi(AtomType* atomType1, AtomType* atomType2) {    
157 <    EAMParam eamParam1 = getEAMParam(atomType1);
158 <    EAMParam eamParam2 = getEAMParam(atomType2);
159 <    CubicSpline* z1 = getZ(atomType1);
160 <    CubicSpline* z2 = getZ(atomType2);
78 >    rmax = max(rmax, ea2.getRcut());
79 >    rmax = max(rmax, ea2.getNr() * ea2.getDr());
80  
81 <    // make the r grid:
81 >    // use the smallest dr (finest grid) to build our grid:
82  
83 <    // set rcut to be the smaller of the two atomic rcuts
83 >    RealType dr = min(ea1.getDr(), ea2.getDr());
84  
85 <    RealType rcut = eamParam1.rcut < eamParam2.rcut ?
167 <      eamParam1.rcut : eamParam2.rcut;
85 >    int nr = int(rmax/dr + 0.5);
86  
169 –    // use the smallest dr (finest grid) to build our grid:
170 –
171 –    RealType dr = eamParam1.dr < eamParam2.dr ? eamParam1.dr : eamParam2.dr;
172 –    int nr = int(rcut/dr);
87      vector<RealType> rvals;
88 <    for (int i = 0; i < nr; i++) rvals.push_back(i*dr);
88 >    for (int i = 0; i < nr; i++) rvals.push_back(RealType(i*dr));
89  
90      // construct the pair potential:
91  
# Line 182 | Line 96 | namespace OpenMD {
96  
97      phivals.push_back(0.0);
98  
99 <    for (int i = 1; i < rvals.size(); i++ ) {
99 >    for (unsigned int i = 1; i < rvals.size(); i++ ) {
100        r = rvals[i];
187 –      zi = z1->getValueAt(r);
188 –      zj = z2->getValueAt(r);
101  
102 <      phi = 331.999296 * (zi * zj) / r;
102 >      // only use z(r) if we're inside this atom's cutoff radius,
103 >      // otherwise, we'll use zero for the charge.  This effectively
104 >      // means that our phi grid goes out beyond the cutoff of the
105 >      // pair potential
106 >
107 >      zi = r <= ea1.getRcut() ? z1->getValueAt(r) : 0.0;
108 >      zj = r <= ea2.getRcut() ? z2->getValueAt(r) : 0.0;
109 >
110 >      phi = pre11_ * (zi * zj) / r;
111 >      
112        phivals.push_back(phi);
113      }
193 –      
114      CubicSpline* cs = new CubicSpline();
115      cs->addPoints(rvals, phivals);
116      return cs;
117    }
118  
119 <  void EAM::initialize() {
119 >  void EAM::setCutoffRadius( RealType rCut ) {
120 >    eamRcut_ = rCut;
121 >    haveCutoffRadius_ = true;
122 >  }
123  
124 +  void EAM::initialize() {
125      // set up the mixing method:
126      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
127 <    std::string EAMMixMeth = fopts.getEAMMixingMethod();
127 >    string EAMMixMeth = fopts.getEAMMixingMethod();
128      toUpper(EAMMixMeth);
129    
130      if (EAMMixMeth == "JOHNSON")
# Line 211 | Line 135 | namespace OpenMD {
135        mixMeth_ = eamUnknown;
136        
137      // find all of the EAM atom Types:
138 <    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
139 <    ForceField::AtomTypeContainer::MapTypeIterator i;
140 <    AtomType* at;
138 >    EAMtypes.clear();
139 >    EAMtids.clear();
140 >    EAMdata.clear();
141 >    MixingMap.clear();
142 >    nEAM_ = 0;
143 >    
144 >    EAMtids.resize( forceField_->getNAtomType(), -1);
145  
146 <    for (at = atomTypes->beginType(i); at != NULL;
147 <         at = atomTypes->nextType(i)) {
148 <      
221 <      if (at->isEAM())
222 <        addType(at);
146 >    set<AtomType*>::iterator at;
147 >    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
148 >      if ((*at)->isEAM()) nEAM_++;
149      }
150 +    EAMdata.resize(nEAM_);
151 +    MixingMap.resize(nEAM_);
152 +
153 +    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
154 +      if ((*at)->isEAM()) addType(*at);
155 +    }
156      
157      // find all of the explicit EAM interactions (setfl):
158      ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
# Line 232 | Line 164 | namespace OpenMD {
164        
165        if (nbt->isEAM()) {
166          
167 <        std::pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
167 >        pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
168          
169          GenericData* data = nbt->getPropertyByName("EAM");
170          if (data == NULL) {
# Line 273 | Line 205 | namespace OpenMD {
205  
206    void EAM::addType(AtomType* atomType){
207  
208 +    EAMAdapter ea = EAMAdapter(atomType);
209      EAMAtomData eamAtomData;
277 –    
278 –    eamAtomData.rho = getRho(atomType);
279 –    eamAtomData.F = getF(atomType);
280 –    eamAtomData.Z = getZ(atomType);
281 –    eamAtomData.rcut = getRcut(atomType);
210  
211 +    eamAtomData.rho = ea.getRho();
212 +    eamAtomData.F = ea.getF();
213 +    eamAtomData.Z = ea.getZ();
214 +    eamAtomData.rcut = ea.getRcut();
215 +      
216      // add it to the map:
217 <    AtomTypeProperties atp = atomType->getATP();    
217 >    int atid = atomType->getIdent();
218 >    int eamtid = EAMtypes.size();
219  
220 <    std::pair<std::map<int,AtomType*>::iterator,bool> ret;    
221 <    ret = EAMlist.insert( std::pair<int, AtomType*>(atp.ident, atomType) );
220 >    pair<set<int>::iterator,bool> ret;    
221 >    ret = EAMtypes.insert( atid );
222      if (ret.second == false) {
223        sprintf( painCave.errMsg,
224                 "EAM already had a previous entry with ident %d\n",
225 <               atp.ident);
225 >               atid);
226        painCave.severity = OPENMD_INFO;
227        painCave.isFatal = 0;
228        simError();        
229      }
230  
231 <    EAMMap[atomType] = eamAtomData;
231 >
232 >    EAMtids[atid] = eamtid;
233 >    EAMdata[eamtid] = eamAtomData;
234 >    MixingMap[eamtid].resize(nEAM_);
235      
236      // Now, iterate over all known types and add to the mixing map:
237      
238 <    std::map<AtomType*, EAMAtomData>::iterator it;
239 <    for( it = EAMMap.begin(); it != EAMMap.end(); ++it) {
238 >    std::set<int>::iterator it;
239 >    for( it = EAMtypes.begin(); it != EAMtypes.end(); ++it) {
240        
241 <      AtomType* atype2 = (*it).first;
241 >      int eamtid2 = EAMtids[ (*it) ];
242 >      AtomType* atype2 = forceField_->getAtomType( (*it) );
243  
244        EAMInteractionData mixer;
245        mixer.phi = getPhi(atomType, atype2);
246 +      mixer.rcut = mixer.phi->getLimits().second;
247        mixer.explicitlySet = false;
248  
249 <      std::pair<AtomType*, AtomType*> key1, key2;
311 <      key1 = std::make_pair(atomType, atype2);
312 <      key2 = std::make_pair(atype2, atomType);
249 >      MixingMap[eamtid2].resize( nEAM_ );
250        
251 <      MixingMap[key1] = mixer;
252 <      if (key2 != key1) {
253 <        MixingMap[key2] = mixer;
251 >      MixingMap[eamtid][eamtid2] = mixer;
252 >      if (eamtid2 != eamtid) {
253 >        MixingMap[eamtid2][eamtid] = mixer;
254        }
255      }      
256      return;
# Line 335 | Line 272 | namespace OpenMD {
272  
273      cs->addPoints(rVals, phiVals);
274      mixer.phi = cs;
275 +    mixer.rcut = mixer.phi->getLimits().second;
276      mixer.explicitlySet = true;
277  
278 <    std::pair<AtomType*, AtomType*> key1, key2;
279 <    key1 = std::make_pair(atype1, atype2);
342 <    key2 = std::make_pair(atype2, atype1);
278 >    int eamtid1 = EAMtids[ atype1->getIdent() ];
279 >    int eamtid2 = EAMtids[ atype2->getIdent() ];
280      
281 <    MixingMap[key1] = mixer;
282 <    if (key2 != key1) {
283 <      MixingMap[key2] = mixer;
281 >    MixingMap[eamtid1][eamtid2] = mixer;
282 >    if (eamtid2 != eamtid1) {
283 >      MixingMap[eamtid2][eamtid1] = mixer;
284      }    
285      return;
286    }
287  
288 <  void EAM::calcDensity(AtomType* at1, AtomType* at2, const RealType rij,
352 <                        RealType &rho_i_at_j, RealType &rho_j_at_i) {
288 >  void EAM::calcDensity(InteractionData &idat) {
289      
290      if (!initialized_) initialize();
291      
292 <    EAMAtomData data1 = EAMMap[at1];
293 <    EAMAtomData data2 = EAMMap[at2];
292 >    EAMAtomData &data1 = EAMdata[EAMtids[idat.atid1]];
293 >    EAMAtomData &data2 = EAMdata[EAMtids[idat.atid2]];
294  
295 <    if (rij < data1.rcut) rho_i_at_j = data1.rho->getValueAt(rij);
296 <    if (rij < data2.rcut) rho_j_at_i = data2.rho->getValueAt(rij);
361 <    return;
362 <  }
363 <
364 <  void EAM::calcFunctional(AtomType* at1, RealType rho, RealType &frho,
365 <                           RealType &dfrhodrho) {
366 <
367 <    if (!initialized_) initialize();
368 <
369 <    EAMAtomData data1 = EAMMap[at1];
370 <        
371 <    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt(rho);
372 <
373 <    frho = result.first;
374 <    dfrhodrho = result.second;
375 <    return;
376 <  }
377 <
378 <
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 rho_i, RealType rho_j,
383 <                      RealType dfrhodrho_i, RealType dfrhodrho_j,
384 <                      RealType &fshift_i, RealType &fshift_j) {
385 <
386 <    if (!initialized_) initialize();
295 >    if (haveCutoffRadius_)
296 >      if ( *(idat.rij) > eamRcut_) return;
297      
298 <    pair<RealType, RealType> res;
299 <    
300 <    if (rij < eamRcut_) {
391 <
392 <      EAMAtomData data1 = EAMMap[at1];
393 <      EAMAtomData data2 = EAMMap[at2];
394 <
395 <      // get type-specific cutoff radii
396 <
397 <      RealType rci = data1.rcut;
398 <      RealType rcj = data2.rcut;
298 >    if ( *(idat.rij) < data1.rcut) {
299 >      *(idat.rho2) += data1.rho->getValueAt( *(idat.rij));
300 >    }
301        
302 <      RealType rha, drha, rhb, drhb;
303 <      RealType pha, dpha, phb, dphb;
402 <      RealType phab, dvpdr;
403 <      RealType drhoidr, drhojdr, dudr;
404 <      
405 <      if (rij < rci) {
406 <        res = data1.rho->getValueAndDerivativeAt(rij);
407 <        rha = res.first;
408 <        drha = res.second;
409 <
410 <        res = MixingMap[make_pair(at1, at1)].phi->getValueAndDerivativeAt(rij);
411 <        pha = res.first;
412 <        dpha = res.second;
413 <      }
414 <
415 <      if (rij < rcj) {
416 <        res = data2.rho->getValueAndDerivativeAt(rij);
417 <        rhb = res.first;
418 <        drhb = res.second;
419 <
420 <        res = MixingMap[make_pair(at2, at2)].phi->getValueAndDerivativeAt(rij);
421 <        phb = res.first;
422 <        dphb = res.second;
423 <      }
424 <
425 <      phab = 0.0;
426 <      dvpdr = 0.0;
427 <
428 <      switch(mixMeth_) {
429 <      case eamJohnson:
430 <      
431 <        if (rij < rci) {
432 <          phab = phab + 0.5 * (rhb / rha) * pha;
433 <          dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
434 <                               pha*((drhb/rha) - (rhb*drha/rha/rha)));
435 <        }
436 <
437 <        if (rij < rcj) {
438 <          phab = phab + 0.5 * (rha / rhb) * phb;
439 <          dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
440 <                                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
441 <        }
442 <
443 <        break;
444 <
445 <      case eamDaw:
446 <                
447 <        res = MixingMap[make_pair(at1,at2)].phi->getValueAndDerivativeAt(rij);
448 <        phab = res.first;
449 <        dvpdr = res.second;
450 <
451 <        break;
452 <      case eamUnknown:
453 <      default:
454 <
455 <        sprintf(painCave.errMsg,
456 <                "EAM::calcForce hit a mixing method it doesn't know about!\n"
457 <                );
458 <        painCave.severity = OPENMD_ERROR;
459 <        painCave.isFatal = 1;
460 <        simError();        
461 <          
462 <      }
463 <      
464 <      drhoidr = drha;
465 <      drhojdr = drhb;
466 <
467 <      dudr = drhojdr*dfrhodrho_i + drhoidr*dfrhodrho_j + dvpdr;
468 <
469 <      f1 = d * dudr / rij;
470 <        
471 <      // particle_pot is the difference between the full potential
472 <      // and the full potential without the presence of a particular
473 <      // particle (atom1).
474 <      //
475 <      // This reduces the density at other particle locations, so
476 <      // we need to recompute the density at atom2 assuming atom1
477 <      // didn't contribute.  This then requires recomputing the
478 <      // density functional for atom2 as well.
479 <      //
480 <      // Most of the particle_pot heavy lifting comes from the
481 <      // pair interaction, and will be handled by vpair.
482 <    
483 <      fshift_i = data1.F->getValueAt( rho_i - rhb );
484 <      fshift_j = data1.F->getValueAt( rho_j - rha );
485 <
486 <      pot += phab;
487 <
488 <      vpair += phab;
302 >    if ( *(idat.rij) < data2.rcut) {
303 >      *(idat.rho1) += data2.rho->getValueAt( *(idat.rij));
304      }
490 –
491 –    return;
305      
306 +    return;  
307    }
308 <
309 <
496 <  void EAM::calc_eam_prepair_rho(int *atid1, int *atid2, RealType *rij,
497 <                                 RealType* rho_i_at_j, RealType* rho_j_at_i){
498 <
499 <    if (!initialized_) initialize();
308 >  
309 >  void EAM::calcFunctional(SelfData &sdat) {
310      
311 <    AtomType* atype1 = EAMlist[*atid1];
312 <    AtomType* atype2 = EAMlist[*atid2];
313 <    
314 <    calcDensity(atype1, atype2, *rij, *rho_i_at_j, *rho_j_at_i);
311 >    if (!initialized_) initialize();
312 >    EAMAtomData &data1 = EAMdata[ EAMtids[sdat.atid] ];
313 >            
314 >    data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho),
315 >                                      *(sdat.dfrhodrho) );
316  
317 <    return;    
317 >    (*(sdat.pot))[METALLIC_FAMILY] += *(sdat.frho);
318 >    if (sdat.doParticlePot) {
319 >      *(sdat.particlePot) += *(sdat.frho);
320 >    }
321 >
322 >    return;
323    }
324  
325 <  void EAM::calc_eam_preforce_Frho(int *atid1, RealType *rho, RealType *frho,
326 <                                   RealType *dfrhodrho) {
325 >
326 >  void EAM::calcForce(InteractionData &idat) {
327  
328      if (!initialized_) initialize();
329  
330 <    AtomType* atype1 = EAMlist[*atid1];  
330 >    if (haveCutoffRadius_)
331 >      if ( *(idat.rij) > eamRcut_) return;
332 >  
333  
334 <    calcFunctional(atype1, *rho, *frho, *dfrhodrho);
334 >    int eamtid1 = EAMtids[idat.atid1];
335 >    int eamtid2 = EAMtids[idat.atid2];
336 >    EAMAtomData &data1 = EAMdata[eamtid1];
337 >    EAMAtomData &data2 = EAMdata[eamtid2];
338      
339 <    return;    
519 <  }
520 <  RealType EAM::getEAMcut(int *atid1) {
521 <
522 <    if (!initialized_) initialize();
339 >    // get type-specific cutoff radii
340      
341 <    AtomType* atype1 = EAMlist[*atid1];  
342 <      
526 <    return getRcut(atype1);
527 <  }
341 >    RealType rci = data1.rcut;
342 >    RealType rcj = data2.rcut;
343  
529 –  void EAM::do_eam_pair(int *atid1, int *atid2, RealType *d, RealType *rij,
530 –                        RealType *r2, RealType *sw, RealType *vpair,
531 –                        RealType *pot, RealType *f1, RealType *rho1,
532 –                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
533 –                        RealType *fshift1, RealType *fshift2) {
534 –
535 –    if (!initialized_) initialize();
344      
345 <    AtomType* atype1 = EAMlist[*atid1];
346 <    AtomType* atype2 = EAMlist[*atid2];
345 >    RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0);
346 >    RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0);
347 >    RealType phab(0.0), dvpdr(0.0);
348 >    RealType drhoidr, drhojdr, dudr;
349      
350 <    Vector3d disp(d[0], d[1], d[2]);
351 <    Vector3d frc(f1[0], f1[1], f1[2]);
350 >    if ( *(idat.rij) < rci) {
351 >      data1.rho->getValueAndDerivativeAt( *(idat.rij), rha, drha);
352 >      CubicSpline* phi = MixingMap[eamtid1][eamtid1].phi;
353 >      phi->getValueAndDerivativeAt( *(idat.rij), pha, dpha);
354 >    }
355      
356 <    calcForce(atype1, atype2, disp, *rij, *r2, *sw, *vpair,  *pot, frc,
357 <              *rho1, *rho2, *dfrho1, *dfrho2, *fshift1, *fshift2);
356 >    if ( *(idat.rij) < rcj) {
357 >      data2.rho->getValueAndDerivativeAt( *(idat.rij), rhb, drhb );
358 >      CubicSpline* phi = MixingMap[eamtid2][eamtid2].phi;
359 >      phi->getValueAndDerivativeAt( *(idat.rij), phb, dphb);
360 >    }
361 >    switch(mixMeth_) {
362 >    case eamJohnson:
363 >      if ( *(idat.rij) < rci) {
364 >        phab = phab + 0.5 * (rhb / rha) * pha;
365 >        dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
366 >                             pha*((drhb/rha) - (rhb*drha/rha/rha)));
367 >      }
368 >                
369 >      if ( *(idat.rij) < rcj) {
370 >        phab = phab + 0.5 * (rha / rhb) * phb;
371 >        dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
372 >                               phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
373 >      }
374 >      break;
375 >    case eamDaw:
376 >      if ( *(idat.rij) <  MixingMap[eamtid1][eamtid2].rcut) {
377 >        MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij),
378 >                                                                  phab, dvpdr);
379 >      }
380 >      break;
381 >    case eamUnknown:
382 >    default:
383        
384 <    f1[0] = frc.x();
385 <    f1[1] = frc.y();
386 <    f1[2] = frc.z();
384 >      sprintf(painCave.errMsg,
385 >              "EAM::calcForce hit a mixing method it doesn't know about!\n"
386 >              );
387 >      painCave.severity = OPENMD_ERROR;
388 >      painCave.isFatal = 1;
389 >      simError();        
390 >      
391 >    }
392 >    
393 >    drhoidr = drha;
394 >    drhojdr = drhb;
395 >    
396 >    dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
397 >    
398 >    *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
399  
400 +    if (idat.doParticlePot) {
401 +      // particlePot is the difference between the full potential and
402 +      // the full potential without the presence of a particular
403 +      // particle (atom1).
404 +      //
405 +      // This reduces the density at other particle locations, so we
406 +      // need to recompute the density at atom2 assuming atom1 didn't
407 +      // contribute.  This then requires recomputing the density
408 +      // functional for atom2 as well.
409 +      
410 +      *(idat.particlePot1) += data2.F->getValueAt( *(idat.rho2) - rha )
411 +        - *(idat.frho2);
412 +      
413 +      *(idat.particlePot2) += data1.F->getValueAt( *(idat.rho1) - rhb)
414 +        - *(idat.frho1);
415 +    }
416 +    
417 +    (*(idat.pot))[METALLIC_FAMILY] += phab;    
418 +    *(idat.vpair) += phab;  
419      return;    
420    }
552 –  
553 –  void EAM::setCutoffEAM(RealType *thisRcut) {
554 –    eamRcut_ = *thisRcut;
555 –  }
556 – }
421  
422 < extern "C" {
423 <  
560 < #define fortranCalcDensity FC_FUNC(calc_eam_prepair_rho, CALC_EAM_PREPAIR_RHO)
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 < #define fortranGetEAMcut FC_FUNC(geteamcut, GETEAMCUT)
422 >  RealType EAM::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
423 >    if (!initialized_) initialize();  
424  
425 <  
426 <  void fortranCalcDensity(int *atid1, int *atid2, RealType *rij,
427 <                          RealType *rho_i_at_j, RealType *rho_j_at_i) {
425 >    RealType cut = 0.0;
426 >
427 >    int atid1 = atypes.first->getIdent();
428 >    int atid2 = atypes.second->getIdent();
429 >    int eamtid1 = EAMtids[atid1];
430 >    int eamtid2 = EAMtids[atid2];
431      
432 <    return OpenMD::EAM::Instance()->calc_eam_prepair_rho(atid1, atid2, rij,
433 <                                                         rho_i_at_j,  
434 <                                                         rho_j_at_i);
435 <  }
436 <  void fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho,
437 <                             RealType *dfrhodrho) {  
432 >    if (eamtid1 != -1) {
433 >      EAMAtomData data1 = EAMdata[eamtid1];
434 >      cut = data1.rcut;
435 >    }
436 >
437 >    if (eamtid2 != -1) {
438 >      EAMAtomData data2 = EAMdata[eamtid2];
439 >      if (data2.rcut > cut)
440 >        cut = data2.rcut;
441 >    }
442      
443 <    return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(atid1, rho, frho,
578 <                                                           dfrhodrho);
579 <    
443 >    return cut;
444    }
581 –  void fortranSetCutoffEAM(RealType *rcut) {
582 –    return OpenMD::EAM::Instance()->setCutoffEAM(rcut);
583 –  }
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);
595 –  }
596 –  RealType fortranGetEAMcut(int* atid) {
597 –    return OpenMD::EAM::Instance()->getEAMcut(atid);
598 –  }
599 –
445   }
446 +

Diff Legend

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