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 1505 by gezelter, Sun Oct 3 22:18:59 2010 UTC vs.
Revision 1711 by gezelter, Sat May 19 02:58:35 2012 UTC

# Line 36 | Line 36
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).                        
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 51 | Line 52 | namespace OpenMD {
52   namespace OpenMD {
53  
54    EAM::EAM() : name_("EAM"), initialized_(false), forceField_(NULL),
55 <               mixMeth_(eamJohnson), eamRcut_(0.0) {}
55 >               mixMeth_(eamJohnson), eamRcut_(0.0), haveCutoffRadius_(false) {}
56    
57 <  EAMParam EAM::getEAMParam(AtomType* atomType) {
58 <    
59 <    // Do sanity checking on the AtomType we were passed before
60 <    // building any data structures:
61 <    if (!atomType->isEAM()) {
61 <      sprintf( painCave.errMsg,
62 <               "EAM::getEAMParam was passed an atomType (%s) that does not\n"
63 <               "\tappear to be an embedded atom method (EAM) atom.\n",
64 <               atomType->getName().c_str());
65 <      painCave.severity = OPENMD_ERROR;
66 <      painCave.isFatal = 1;
67 <      simError();
68 <    }
69 <    
70 <    GenericData* data = atomType->getPropertyByName("EAM");
71 <    if (data == NULL) {
72 <      sprintf( painCave.errMsg, "EAM::getEAMParam could not find EAM\n"
73 <               "\tparameters for atomType %s.\n",
74 <               atomType->getName().c_str());
75 <      painCave.severity = OPENMD_ERROR;
76 <      painCave.isFatal = 1;
77 <      simError();
78 <    }
79 <    
80 <    EAMParamGenericData* eamData = dynamic_cast<EAMParamGenericData*>(data);
81 <    if (eamData == NULL) {
82 <      sprintf( painCave.errMsg,
83 <               "EAM::getEAMParam could not convert GenericData to EAMParam for\n"
84 <               "\tatom type %s\n", atomType->getName().c_str());
85 <      painCave.severity = OPENMD_ERROR;
86 <      painCave.isFatal = 1;
87 <      simError();          
88 <    }
89 <    
90 <    return eamData->getData();
91 <  }
57 >  CubicSpline* EAM::getPhi(AtomType* atomType1, AtomType* atomType2) {  
58 >    EAMAdapter ea1 = EAMAdapter(atomType1);
59 >    EAMAdapter ea2 = EAMAdapter(atomType2);
60 >    CubicSpline* z1 = ea1.getZ();
61 >    CubicSpline* z2 = ea2.getZ();
62  
93 –  CubicSpline* EAM::getZ(AtomType* atomType) {    
94 –    EAMParam eamParam = getEAMParam(atomType);
95 –    int nr = eamParam.nr;
96 –    RealType dr = eamParam.dr;
97 –    vector<RealType> rvals;
98 –    
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(RealType(i) * dr);
118 –      
119 –    CubicSpline* cs = new CubicSpline();
120 –    cs->addPoints(rvals, eamParam.rho);
121 –    return cs;
122 –  }
123 –
124 –  CubicSpline* EAM::getF(AtomType* atomType) {    
125 –    EAMParam eamParam = getEAMParam(atomType);
126 –    int nrho = eamParam.nrho;
127 –    RealType drho = eamParam.drho;
128 –    vector<RealType> rhovals;
129 –    vector<RealType> scaledF;
130 –    
131 –    for (int i = 0; i < nrho; i++) {
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, scaledF);
138 –    return cs;
139 –  }
140 –  
141 –  CubicSpline* EAM::getPhi(AtomType* atomType1, AtomType* atomType2) {    
142 –    EAMParam eamParam1 = getEAMParam(atomType1);
143 –    EAMParam eamParam2 = getEAMParam(atomType2);
144 –    CubicSpline* z1 = getZ(atomType1);
145 –    CubicSpline* z2 = getZ(atomType2);
146 –
63      // make the r grid:
64  
65  
66      // we need phi out to the largest value we'll encounter in the radial space;
67      
68      RealType rmax = 0.0;
69 <    rmax = max(rmax, eamParam1.rcut);
70 <    rmax = max(rmax, eamParam1.nr * eamParam1.dr);
69 >    rmax = max(rmax, ea1.getRcut());
70 >    rmax = max(rmax, ea1.getNr() * ea1.getDr());
71  
72 <    rmax = max(rmax, eamParam2.rcut);
73 <    rmax = max(rmax, eamParam2.nr * eamParam2.dr);
72 >    rmax = max(rmax, ea2.getRcut());
73 >    rmax = max(rmax, ea2.getNr() * ea2.getDr());
74  
75      // use the smallest dr (finest grid) to build our grid:
76  
77 <    RealType dr = min(eamParam1.dr, eamParam2.dr);
77 >    RealType dr = min(ea1.getDr(), ea2.getDr());
78  
79      int nr = int(rmax/dr + 0.5);
80  
# Line 182 | Line 98 | namespace OpenMD {
98        // means that our phi grid goes out beyond the cutoff of the
99        // pair potential
100  
101 <      zi = r <= eamParam1.rcut ? z1->getValueAt(r) : 0.0;
102 <      zj = r <= eamParam2.rcut ? z2->getValueAt(r) : 0.0;
101 >      zi = r <= ea1.getRcut() ? z1->getValueAt(r) : 0.0;
102 >      zj = r <= ea2.getRcut() ? z2->getValueAt(r) : 0.0;
103  
104        phi = 331.999296 * (zi * zj) / r;
105  
# Line 195 | Line 111 | namespace OpenMD {
111      return cs;
112    }
113  
114 +  void EAM::setCutoffRadius( RealType rCut ) {
115 +    eamRcut_ = rCut;
116 +    haveCutoffRadius_ = true;
117 +  }
118 +
119    void EAM::initialize() {
120  
121      // set up the mixing method:
# Line 272 | Line 193 | namespace OpenMD {
193  
194    void EAM::addType(AtomType* atomType){
195  
196 +    EAMAdapter ea = EAMAdapter(atomType);
197      EAMAtomData eamAtomData;
276 –    
277 –    eamAtomData.rho = getRho(atomType);
278 –    eamAtomData.F = getF(atomType);
279 –    eamAtomData.Z = getZ(atomType);
280 –    eamAtomData.rcut = getRcut(atomType);
198  
199 +    eamAtomData.rho = ea.getRho();
200 +    eamAtomData.F = ea.getF();
201 +    eamAtomData.Z = ea.getZ();
202 +    eamAtomData.rcut = ea.getRcut();
203 +
204      // add it to the map:
283 –    AtomTypeProperties atp = atomType->getATP();    
205  
206      pair<map<int,AtomType*>::iterator,bool> ret;    
207 <    ret = EAMlist.insert( pair<int, AtomType*>(atp.ident, atomType) );
207 >    ret = EAMlist.insert( pair<int, AtomType*>(atomType->getIdent(), atomType) );
208      if (ret.second == false) {
209        sprintf( painCave.errMsg,
210                 "EAM already had a previous entry with ident %d\n",
211 <               atp.ident);
211 >               atomType->getIdent());
212        painCave.severity = OPENMD_INFO;
213        painCave.isFatal = 0;
214        simError();        
# Line 347 | Line 268 | namespace OpenMD {
268      return;
269    }
270  
271 <  void EAM::calcDensity(DensityData ddat) {
271 >  void EAM::calcDensity(InteractionData &idat) {
272      
273      if (!initialized_) initialize();
274      
275 <    EAMAtomData data1 = EAMMap[ddat.atype1];
276 <    EAMAtomData data2 = EAMMap[ddat.atype2];
277 <
278 <    if (ddat.rij < data1.rcut)
279 <      ddat.rho_i_at_j = data1.rho->getValueAt(ddat.rij);
280 <
281 <    if (ddat.rij < data2.rcut)
282 <      ddat.rho_j_at_i = data2.rho->getValueAt(ddat.rij);
283 <
284 <    return;
275 >    EAMAtomData data1 = EAMMap[idat.atypes.first];
276 >    EAMAtomData data2 = EAMMap[idat.atypes.second];
277 >    
278 >    if (haveCutoffRadius_)
279 >      if ( *(idat.rij) > eamRcut_) return;
280 >    
281 >    if ( *(idat.rij) < data1.rcut)
282 >      *(idat.rho1) += data1.rho->getValueAt( *(idat.rij));
283 >    
284 >      
285 >    if ( *(idat.rij) < data2.rcut)
286 >      *(idat.rho2) += data2.rho->getValueAt( *(idat.rij));
287 >    
288 >    return;  
289    }
290 <
291 <  void EAM::calcFunctional(FunctionalData fdat) {
292 <
290 >  
291 >  void EAM::calcFunctional(SelfData &sdat) {
292 >    
293      if (!initialized_) initialize();
294  
295 <    EAMAtomData data1 = EAMMap[fdat.atype];
295 >    EAMAtomData data1 = EAMMap[ sdat.atype ];
296          
297 <    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt(fdat.rho);
297 >    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt( *(sdat.rho) );
298  
299 <    fdat.frho = result.first;
300 <    fdat.dfrhodrho = result.second;
299 >    *(sdat.frho) = result.first;
300 >    *(sdat.dfrhodrho) = result.second;
301 >
302 >    (*(sdat.pot))[METALLIC_FAMILY] += result.first;
303 >    if (sdat.doParticlePot) {
304 >      *(sdat.particlePot) += result.first;
305 >    }
306 >
307      return;
308    }
309  
310  
311 <  void EAM::calcForce(InteractionData idat) {
311 >  void EAM::calcForce(InteractionData &idat) {
312  
313      if (!initialized_) initialize();
314  
315 +    if (haveCutoffRadius_)
316 +      if ( *(idat.rij) > eamRcut_) return;
317 +  
318      pair<RealType, RealType> res;
319      
320 <    if (idat.rij < eamRcut_) {
321 <
322 <      EAMAtomData data1 = EAMMap[idat.atype1];
323 <      EAMAtomData data2 = EAMMap[idat.atype2];
324 <
325 <      // get type-specific cutoff radii
326 <
327 <      RealType rci = data1.rcut;
328 <      RealType rcj = data2.rcut;
320 >    EAMAtomData data1 = EAMMap[idat.atypes.first];
321 >    EAMAtomData data2 = EAMMap[idat.atypes.second];
322 >    
323 >    // get type-specific cutoff radii
324 >    
325 >    RealType rci = data1.rcut;
326 >    RealType rcj = data2.rcut;
327 >    
328 >    RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0);
329 >    RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0);
330 >    RealType phab(0.0), dvpdr(0.0);
331 >    RealType drhoidr, drhojdr, dudr;
332 >    
333 >    if ( *(idat.rij) < rci) {
334 >      res = data1.rho->getValueAndDerivativeAt( *(idat.rij));
335 >      rha = res.first;
336 >      drha = res.second;
337        
338 <      RealType rha, drha, rhb, drhb;
339 <      RealType pha, dpha, phb, dphb;
340 <      RealType phab, dvpdr;
341 <      RealType drhoidr, drhojdr, dudr;
338 >      res = MixingMap[make_pair(idat.atypes.first, idat.atypes.first)].phi->getValueAndDerivativeAt( *(idat.rij) );
339 >      pha = res.first;
340 >      dpha = res.second;
341 >    }
342 >    
343 >    if ( *(idat.rij) < rcj) {
344 >      res = data2.rho->getValueAndDerivativeAt( *(idat.rij) );
345 >      rhb = res.first;
346 >      drhb = res.second;
347        
348 <      if (idat.rij < rci) {
349 <        res = data1.rho->getValueAndDerivativeAt(idat.rij);
350 <        rha = res.first;
351 <        drha = res.second;
348 >      res = MixingMap[make_pair(idat.atypes.second, idat.atypes.second)].phi->getValueAndDerivativeAt( *(idat.rij) );
349 >      phb = res.first;
350 >      dphb = res.second;
351 >    }
352  
353 <        res = MixingMap[make_pair(idat.atype1, idat.atype1)].phi->getValueAndDerivativeAt(idat.rij);
354 <        pha = res.first;
355 <        dpha = res.second;
353 >    switch(mixMeth_) {
354 >    case eamJohnson:
355 >      
356 >      if ( *(idat.rij) < rci) {
357 >        phab = phab + 0.5 * (rhb / rha) * pha;
358 >        dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
359 >                             pha*((drhb/rha) - (rhb*drha/rha/rha)));
360        }
361 <
362 <      if (idat.rij < rcj) {
363 <        res = data2.rho->getValueAndDerivativeAt(idat.rij);
364 <        rhb = res.first;
365 <        drhb = res.second;
366 <
367 <        res = MixingMap[make_pair(idat.atype2, idat.atype2)].phi->getValueAndDerivativeAt(idat.rij);
417 <        phb = res.first;
418 <        dphb = res.second;
361 >      
362 >      
363 >      
364 >      if ( *(idat.rij) < rcj) {
365 >        phab = phab + 0.5 * (rha / rhb) * phb;
366 >        dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
367 >                               phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
368        }
420 –
421 –      phab = 0.0;
422 –      dvpdr = 0.0;
423 –
424 –      switch(mixMeth_) {
425 –      case eamJohnson:
426 –      
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 (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)));
437 –        }
438 –
439 –        break;
440 –
441 –      case eamDaw:
442 –        res = MixingMap[make_pair(idat.atype1,idat.atype2)].phi->getValueAndDerivativeAt(idat.rij);
443 –        phab = res.first;
444 –        dvpdr = res.second;
445 –
446 –        break;
447 –      case eamUnknown:
448 –      default:
449 –
450 –        sprintf(painCave.errMsg,
451 –                "EAM::calcForce hit a mixing method it doesn't know about!\n"
452 –                );
453 –        painCave.severity = OPENMD_ERROR;
454 –        painCave.isFatal = 1;
455 –        simError();        
456 –          
457 –      }
369        
370 <      drhoidr = drha;
371 <      drhojdr = drhb;
372 <
373 <      dudr = drhojdr*idat.dfrho1 + drhoidr*idat.dfrho2 + dvpdr;
374 <
375 <      idat.f1 = idat.d * dudr / idat.rij;
370 >      break;
371 >      
372 >    case eamDaw:
373 >      res = MixingMap[idat.atypes].phi->getValueAndDerivativeAt( *(idat.rij));
374 >      phab = res.first;
375 >      dvpdr = res.second;
376 >      
377 >      break;
378 >    case eamUnknown:
379 >    default:
380 >      
381 >      sprintf(painCave.errMsg,
382 >              "EAM::calcForce hit a mixing method it doesn't know about!\n"
383 >              );
384 >      painCave.severity = OPENMD_ERROR;
385 >      painCave.isFatal = 1;
386 >      simError();        
387 >      
388 >    }
389 >    
390 >    drhoidr = drha;
391 >    drhojdr = drhb;
392 >    
393 >    dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
394 >    
395 >    *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
396          
397 <      // particle_pot is the difference between the full potential
398 <      // and the full potential without the presence of a particular
397 >    if (idat.doParticlePot) {
398 >      // particlePot is the difference between the full potential and
399 >      // the full potential without the presence of a particular
400        // particle (atom1).
401        //
402 <      // This reduces the density at other particle locations, so
403 <      // we need to recompute the density at atom2 assuming atom1
404 <      // didn't contribute.  This then requires recomputing the
405 <      // density functional for atom2 as well.
406 <      //
407 <      // Most of the particle_pot heavy lifting comes from the
408 <      // pair interaction, and will be handled by vpair.
409 <    
410 <      idat.fshift1 = data1.F->getValueAt( idat.rho1 - rhb );
411 <      idat.fshift2 = data1.F->getValueAt( idat.rho2 - rha );
480 <
481 <      idat.pot += phab;
482 <
483 <      idat.vpair += phab;
402 >      // This reduces the density at other particle locations, so we
403 >      // need to recompute the density at atom2 assuming atom1 didn't
404 >      // contribute.  This then requires recomputing the density
405 >      // functional for atom2 as well.
406 >      
407 >      *(idat.particlePot1) += data2.F->getValueAt( *(idat.rho2) - rha )
408 >        - *(idat.frho2);
409 >      
410 >      *(idat.particlePot2) += data1.F->getValueAt( *(idat.rho1) - rhb)
411 >        - *(idat.frho1);
412      }
413 <
413 >    
414 >    (*(idat.pot))[METALLIC_FAMILY] += phab;
415 >    
416 >    *(idat.vpair) += phab;
417 >  
418      return;
419      
420    }
421  
422 <  RealType EAM::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
422 >  RealType EAM::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
423      if (!initialized_) initialize();  
424  
425      RealType cut = 0.0;
426  
427      map<AtomType*, EAMAtomData>::iterator it;
428  
429 <    it = EAMMap.find(at1);
429 >    it = EAMMap.find(atypes.first);
430      if (it != EAMMap.end()) {
431        EAMAtomData data1 = (*it).second;
432        cut = data1.rcut;
433      }
434  
435 <    it = EAMMap.find(at2);
435 >    it = EAMMap.find(atypes.second);
436      if (it != EAMMap.end()) {
437        EAMAtomData data2 = (*it).second;
438        if (data2.rcut > cut)

Diff Legend

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