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 1586 by gezelter, Tue Jun 21 06:34:35 2011 UTC vs.
trunk/src/nonbonded/EAM.cpp (file contents), Revision 2033 by gezelter, Sat Nov 1 14:12:16 2014 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 53 | Line 54 | namespace OpenMD {
54    EAM::EAM() : name_("EAM"), initialized_(false), forceField_(NULL),
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  
63 <  CubicSpline* EAM::getZ(AtomType* atomType) {    
64 <    EAMParam eamParam = getEAMParam(atomType);
65 <    int nr = eamParam.nr;
66 <    RealType dr = eamParam.dr;
67 <    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);
63 >    // Thise prefactors convert the charge-charge interactions into
64 >    // kcal / mol all were computed assuming distances are measured in
65 >    // angstroms Charge-Charge, assuming charges are measured in
66 >    // electrons.  Matches value in Electrostatics.cpp
67 >    pre11_ = 332.0637778;
68  
69      // make the r grid:
70  
149 –
71      // we need phi out to the largest value we'll encounter in the radial space;
72      
73      RealType rmax = 0.0;
74 <    rmax = max(rmax, eamParam1.rcut);
75 <    rmax = max(rmax, eamParam1.nr * eamParam1.dr);
74 >    rmax = max(rmax, ea1.getRcut());
75 >    rmax = max(rmax, ea1.getNr() * ea1.getDr());
76  
77 <    rmax = max(rmax, eamParam2.rcut);
78 <    rmax = max(rmax, eamParam2.nr * eamParam2.dr);
77 >    rmax = max(rmax, ea2.getRcut());
78 >    rmax = max(rmax, ea2.getNr() * ea2.getDr());
79  
80      // use the smallest dr (finest grid) to build our grid:
81  
82 <    RealType dr = min(eamParam1.dr, eamParam2.dr);
82 >    RealType dr = min(ea1.getDr(), ea2.getDr());
83  
84      int nr = int(rmax/dr + 0.5);
85  
# Line 174 | Line 95 | namespace OpenMD {
95  
96      phivals.push_back(0.0);
97  
98 <    for (int i = 1; i < rvals.size(); i++ ) {
98 >    for (unsigned int i = 1; i < rvals.size(); i++ ) {
99        r = rvals[i];
100  
101        // only use z(r) if we're inside this atom's cutoff radius,
# Line 182 | Line 103 | namespace OpenMD {
103        // means that our phi grid goes out beyond the cutoff of the
104        // pair potential
105  
106 <      zi = r <= eamParam1.rcut ? z1->getValueAt(r) : 0.0;
107 <      zj = r <= eamParam2.rcut ? z2->getValueAt(r) : 0.0;
106 >      zi = r <= ea1.getRcut() ? z1->getValueAt(r) : 0.0;
107 >      zj = r <= ea2.getRcut() ? z2->getValueAt(r) : 0.0;
108  
109 <      phi = 331.999296 * (zi * zj) / r;
110 <
109 >      phi = pre11_ * (zi * zj) / r;
110 >      
111        phivals.push_back(phi);
112      }
192 –      
113      CubicSpline* cs = new CubicSpline();
114      cs->addPoints(rvals, phivals);
115      return cs;
# Line 201 | Line 121 | namespace OpenMD {
121    }
122  
123    void EAM::initialize() {
204 –
124      // set up the mixing method:
125      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
126      string EAMMixMeth = fopts.getEAMMixingMethod();
# Line 215 | Line 134 | namespace OpenMD {
134        mixMeth_ = eamUnknown;
135        
136      // find all of the EAM atom Types:
137 <    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
138 <    ForceField::AtomTypeContainer::MapTypeIterator i;
139 <    AtomType* at;
137 >    EAMtypes.clear();
138 >    EAMtids.clear();
139 >    EAMdata.clear();
140 >    MixingMap.clear();
141 >    nEAM_ = 0;
142 >    
143 >    EAMtids.resize( forceField_->getNAtomType(), -1);
144  
145 <    for (at = atomTypes->beginType(i); at != NULL;
146 <         at = atomTypes->nextType(i)) {
147 <      
225 <      if (at->isEAM())
226 <        addType(at);
145 >    set<AtomType*>::iterator at;
146 >    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
147 >      if ((*at)->isEAM()) nEAM_++;
148      }
149 +    EAMdata.resize(nEAM_);
150 +    MixingMap.resize(nEAM_);
151 +
152 +    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
153 +      if ((*at)->isEAM()) addType(*at);
154 +    }
155      
156      // find all of the explicit EAM interactions (setfl):
157      ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
# Line 277 | Line 204 | namespace OpenMD {
204  
205    void EAM::addType(AtomType* atomType){
206  
207 +    EAMAdapter ea = EAMAdapter(atomType);
208      EAMAtomData eamAtomData;
281 –    
282 –    eamAtomData.rho = getRho(atomType);
283 –    eamAtomData.F = getF(atomType);
284 –    eamAtomData.Z = getZ(atomType);
285 –    eamAtomData.rcut = getRcut(atomType);
209  
210 +    eamAtomData.rho = ea.getRho();
211 +    eamAtomData.F = ea.getF();
212 +    eamAtomData.Z = ea.getZ();
213 +    eamAtomData.rcut = ea.getRcut();
214 +      
215      // add it to the map:
216 <    AtomTypeProperties atp = atomType->getATP();    
216 >    int atid = atomType->getIdent();
217 >    int eamtid = EAMtypes.size();
218  
219 <    pair<map<int,AtomType*>::iterator,bool> ret;    
220 <    ret = EAMlist.insert( pair<int, AtomType*>(atp.ident, atomType) );
219 >    pair<set<int>::iterator,bool> ret;    
220 >    ret = EAMtypes.insert( atid );
221      if (ret.second == false) {
222        sprintf( painCave.errMsg,
223                 "EAM already had a previous entry with ident %d\n",
224 <               atp.ident);
224 >               atid);
225        painCave.severity = OPENMD_INFO;
226        painCave.isFatal = 0;
227        simError();        
228      }
229  
230 <    EAMMap[atomType] = eamAtomData;
230 >
231 >    EAMtids[atid] = eamtid;
232 >    EAMdata[eamtid] = eamAtomData;
233 >    MixingMap[eamtid].resize(nEAM_);
234      
235      // Now, iterate over all known types and add to the mixing map:
236      
237 <    map<AtomType*, EAMAtomData>::iterator it;
238 <    for( it = EAMMap.begin(); it != EAMMap.end(); ++it) {
237 >    std::set<int>::iterator it;
238 >    for( it = EAMtypes.begin(); it != EAMtypes.end(); ++it) {
239        
240 <      AtomType* atype2 = (*it).first;
240 >      int eamtid2 = EAMtids[ (*it) ];
241 >      AtomType* atype2 = forceField_->getAtomType( (*it) );
242  
243        EAMInteractionData mixer;
244        mixer.phi = getPhi(atomType, atype2);
245 +      mixer.rcut = mixer.phi->getLimits().second;
246        mixer.explicitlySet = false;
247  
248 <      pair<AtomType*, AtomType*> key1, key2;
315 <      key1 = make_pair(atomType, atype2);
316 <      key2 = make_pair(atype2, atomType);
248 >      MixingMap[eamtid2].resize( nEAM_ );
249        
250 <      MixingMap[key1] = mixer;
251 <      if (key2 != key1) {
252 <        MixingMap[key2] = mixer;
250 >      MixingMap[eamtid][eamtid2] = mixer;
251 >      if (eamtid2 != eamtid) {
252 >        MixingMap[eamtid2][eamtid] = mixer;
253        }
254      }      
255      return;
# Line 339 | Line 271 | namespace OpenMD {
271  
272      cs->addPoints(rVals, phiVals);
273      mixer.phi = cs;
274 +    mixer.rcut = mixer.phi->getLimits().second;
275      mixer.explicitlySet = true;
276  
277 <    pair<AtomType*, AtomType*> key1, key2;
278 <    key1 = make_pair(atype1, atype2);
346 <    key2 = make_pair(atype2, atype1);
277 >    int eamtid1 = EAMtids[ atype1->getIdent() ];
278 >    int eamtid2 = EAMtids[ atype2->getIdent() ];
279      
280 <    MixingMap[key1] = mixer;
281 <    if (key2 != key1) {
282 <      MixingMap[key2] = mixer;
280 >    MixingMap[eamtid1][eamtid2] = mixer;
281 >    if (eamtid2 != eamtid1) {
282 >      MixingMap[eamtid2][eamtid1] = mixer;
283      }    
284      return;
285    }
# Line 356 | Line 288 | namespace OpenMD {
288      
289      if (!initialized_) initialize();
290      
291 <    EAMAtomData data1 = EAMMap[idat.atypes.first];
292 <    EAMAtomData data2 = EAMMap[idat.atypes.second];
293 <    
291 >    EAMAtomData &data1 = EAMdata[EAMtids[idat.atid1]];
292 >    EAMAtomData &data2 = EAMdata[EAMtids[idat.atid2]];
293 >
294      if (haveCutoffRadius_)
295        if ( *(idat.rij) > eamRcut_) return;
296      
297      if ( *(idat.rij) < data1.rcut) {
298 <      *(idat.rho1) += data1.rho->getValueAt( *(idat.rij));
298 >      *(idat.rho2) += data1.rho->getValueAt( *(idat.rij));
299 >    }
300        
301 <      
302 <      if ( *(idat.rij) < data2.rcut)
370 <        *(idat.rho2) += data2.rho->getValueAt( *(idat.rij));
371 <      
372 <      return;
301 >    if ( *(idat.rij) < data2.rcut) {
302 >      *(idat.rho1) += data2.rho->getValueAt( *(idat.rij));
303      }
304 +    
305 +    return;  
306    }
307    
308    void EAM::calcFunctional(SelfData &sdat) {
309      
310      if (!initialized_) initialize();
311 +    EAMAtomData &data1 = EAMdata[ EAMtids[sdat.atid] ];
312 +            
313 +    data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho),
314 +                                      *(sdat.dfrhodrho) );
315  
316 <    EAMAtomData data1 = EAMMap[ sdat.atype ];
317 <        
318 <    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt( *(sdat.rho) );
316 >    (*(sdat.pot))[METALLIC_FAMILY] += *(sdat.frho);
317 >    if (sdat.doParticlePot) {
318 >      *(sdat.particlePot) += *(sdat.frho);
319 >    }
320  
384 –    *(sdat.frho) = result.first;
385 –    *(sdat.dfrhodrho) = result.second;
386 –
387 –    (*(sdat.pot))[METALLIC_FAMILY] += result.first;
388 –    *(sdat.particlePot) += result.first;
389 –
321      return;
322    }
323  
# Line 395 | Line 326 | namespace OpenMD {
326  
327      if (!initialized_) initialize();
328  
398 –
399 –
329      if (haveCutoffRadius_)
330        if ( *(idat.rij) > eamRcut_) return;
331    
332 <    pair<RealType, RealType> res;
332 >
333 >    int eamtid1 = EAMtids[idat.atid1];
334 >    int eamtid2 = EAMtids[idat.atid2];
335 >    EAMAtomData &data1 = EAMdata[eamtid1];
336 >    EAMAtomData &data2 = EAMdata[eamtid2];
337      
405 –    
406 –    EAMAtomData data1 = EAMMap[idat.atypes.first];
407 –    EAMAtomData data2 = EAMMap[idat.atypes.second];
408 –    
338      // get type-specific cutoff radii
339      
340      RealType rci = data1.rcut;
341      RealType rcj = data2.rcut;
342 +
343      
344      RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0);
345      RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0);
# Line 417 | Line 347 | namespace OpenMD {
347      RealType drhoidr, drhojdr, dudr;
348      
349      if ( *(idat.rij) < rci) {
350 <      res = data1.rho->getValueAndDerivativeAt( *(idat.rij));
351 <      rha = res.first;
352 <      drha = res.second;
423 <      
424 <      res = MixingMap[make_pair(idat.atypes.first, idat.atypes.first)].phi->getValueAndDerivativeAt( *(idat.rij) );
425 <      pha = res.first;
426 <      dpha = res.second;
350 >      data1.rho->getValueAndDerivativeAt( *(idat.rij), rha, drha);
351 >      CubicSpline* phi = MixingMap[eamtid1][eamtid1].phi;
352 >      phi->getValueAndDerivativeAt( *(idat.rij), pha, dpha);
353      }
354      
355      if ( *(idat.rij) < rcj) {
356 <      res = data2.rho->getValueAndDerivativeAt( *(idat.rij) );
357 <      rhb = res.first;
358 <      drhb = res.second;
433 <      
434 <      res = MixingMap[make_pair(idat.atypes.second, idat.atypes.second)].phi->getValueAndDerivativeAt( *(idat.rij) );
435 <      phb = res.first;
436 <      dphb = res.second;
356 >      data2.rho->getValueAndDerivativeAt( *(idat.rij), rhb, drhb );
357 >      CubicSpline* phi = MixingMap[eamtid2][eamtid2].phi;
358 >      phi->getValueAndDerivativeAt( *(idat.rij), phb, dphb);
359      }
438 –
360      switch(mixMeth_) {
361      case eamJohnson:
441 –      
362        if ( *(idat.rij) < rci) {
363          phab = phab + 0.5 * (rhb / rha) * pha;
364          dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
365                               pha*((drhb/rha) - (rhb*drha/rha/rha)));
366        }
367 <      
448 <      
449 <      
367 >                
368        if ( *(idat.rij) < rcj) {
369          phab = phab + 0.5 * (rha / rhb) * phb;
370          dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
371                                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
372        }
455 –      
373        break;
457 –      
374      case eamDaw:
375 <      res = MixingMap[idat.atypes].phi->getValueAndDerivativeAt( *(idat.rij));
376 <      phab = res.first;
377 <      dvpdr = res.second;
378 <      
375 >      if ( *(idat.rij) <  MixingMap[eamtid1][eamtid2].rcut) {
376 >        MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij),
377 >                                                                  phab, dvpdr);
378 >      }
379        break;
380      case eamUnknown:
381      default:
# Line 479 | Line 395 | namespace OpenMD {
395      dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
396      
397      *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
398 <        
399 <    // particlePot is the difference between the full potential and
400 <    // the full potential without the presence of a particular
401 <    // particle (atom1).
402 <    //
403 <    // This reduces the density at other particle locations, so we
404 <    // need to recompute the density at atom2 assuming atom1 didn't
405 <    // contribute.  This then requires recomputing the density
406 <    // functional for atom2 as well.
398 >
399 >    if (idat.doParticlePot) {
400 >      // particlePot is the difference between the full potential and
401 >      // the full potential without the presence of a particular
402 >      // particle (atom1).
403 >      //
404 >      // This reduces the density at other particle locations, so we
405 >      // need to recompute the density at atom2 assuming atom1 didn't
406 >      // contribute.  This then requires recomputing the density
407 >      // functional for atom2 as well.
408 >      
409 >      *(idat.particlePot1) += data2.F->getValueAt( *(idat.rho2) - rha )
410 >        - *(idat.frho2);
411 >      
412 >      *(idat.particlePot2) += data1.F->getValueAt( *(idat.rho1) - rhb)
413 >        - *(idat.frho1);
414 >    }
415      
416 <    *(idat.particlePot1) += data2.F->getValueAt( *(idat.rho2) - rha )
417 <      - *(idat.frho2);
418 <    
495 <    *(idat.particlePot2) += data1.F->getValueAt( *(idat.rho1) - rhb)
496 <      - *(idat.frho1);
497 <    
498 <    (*(idat.pot))[METALLIC_FAMILY] += phab;
499 <    
500 <    *(idat.vpair) += phab;
501 <  
502 <    return;
503 <    
416 >    (*(idat.pot))[METALLIC_FAMILY] += phab;    
417 >    *(idat.vpair) += phab;  
418 >    return;    
419    }
420  
421    RealType EAM::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
# Line 508 | Line 423 | namespace OpenMD {
423  
424      RealType cut = 0.0;
425  
426 <    map<AtomType*, EAMAtomData>::iterator it;
427 <
428 <    it = EAMMap.find(atypes.first);
429 <    if (it != EAMMap.end()) {
430 <      EAMAtomData data1 = (*it).second;
426 >    int atid1 = atypes.first->getIdent();
427 >    int atid2 = atypes.second->getIdent();
428 >    int eamtid1 = EAMtids[atid1];
429 >    int eamtid2 = EAMtids[atid2];
430 >    
431 >    if (eamtid1 != -1) {
432 >      EAMAtomData data1 = EAMdata[eamtid1];
433        cut = data1.rcut;
434      }
435  
436 <    it = EAMMap.find(atypes.second);
437 <    if (it != EAMMap.end()) {
521 <      EAMAtomData data2 = (*it).second;
436 >    if (eamtid2 != -1) {
437 >      EAMAtomData data2 = EAMdata[eamtid2];
438        if (data2.rcut > cut)
439          cut = data2.rcut;
440      }
441 <
441 >    
442      return cut;
443    }
444   }

Diff Legend

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