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 1505 by gezelter, Sun Oct 3 22:18:59 2010 UTC vs.
trunk/src/nonbonded/EAM.cpp (file contents), Revision 1929 by gezelter, Mon Aug 19 13:12:00 2013 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 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  
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 <  }
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  
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 –
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;
109 >      phi = pre11_ * (zi * zj) / r;
110  
111        phivals.push_back(phi);
112      }
# Line 195 | Line 116 | namespace OpenMD {
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      string EAMMixMeth = fopts.getEAMMixingMethod();
# Line 210 | 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 <      
220 <      if (at->isEAM())
221 <        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 272 | Line 205 | namespace OpenMD {
205  
206    void EAM::addType(AtomType* atomType){
207  
208 +    EAMAdapter ea = EAMAdapter(atomType);
209      EAMAtomData eamAtomData;
276 –    
277 –    eamAtomData.rho = getRho(atomType);
278 –    eamAtomData.F = getF(atomType);
279 –    eamAtomData.Z = getZ(atomType);
280 –    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 <    pair<map<int,AtomType*>::iterator,bool> ret;    
221 <    ret = EAMlist.insert( 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 <    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 <      pair<AtomType*, AtomType*> key1, key2;
310 <      key1 = make_pair(atomType, atype2);
311 <      key2 = 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 334 | 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 <    pair<AtomType*, AtomType*> key1, key2;
279 <    key1 = make_pair(atype1, atype2);
341 <    key2 = 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(DensityData ddat) {
288 >  void EAM::calcDensity(InteractionData &idat) {
289      
290      if (!initialized_) initialize();
291      
292 <    EAMAtomData data1 = EAMMap[ddat.atype1];
293 <    EAMAtomData data2 = EAMMap[ddat.atype2];
292 >    EAMAtomData &data1 = EAMdata[EAMtids[idat.atid1]];
293 >    EAMAtomData &data2 = EAMdata[EAMtids[idat.atid2]];
294  
295 <    if (ddat.rij < data1.rcut)
296 <      ddat.rho_i_at_j = data1.rho->getValueAt(ddat.rij);
297 <
298 <    if (ddat.rij < data2.rcut)
299 <      ddat.rho_j_at_i = data2.rho->getValueAt(ddat.rij);
300 <
301 <    return;
302 <  }
303 <
304 <  void EAM::calcFunctional(FunctionalData fdat) {
305 <
295 >    if (haveCutoffRadius_)
296 >      if ( *(idat.rij) > eamRcut_) return;
297 >    
298 >    if ( *(idat.rij) < data1.rcut) {
299 >      *(idat.rho2) += data1.rho->getValueAt( *(idat.rij));
300 >    }
301 >      
302 >    if ( *(idat.rij) < data2.rcut) {
303 >      *(idat.rho1) += data2.rho->getValueAt( *(idat.rij));
304 >    }
305 >    
306 >    return;  
307 >  }
308 >  
309 >  void EAM::calcFunctional(SelfData &sdat) {
310 >    
311      if (!initialized_) initialize();
312  
313 <    EAMAtomData data1 = EAMMap[fdat.atype];
314 <        
315 <    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt(fdat.rho);
313 >    EAMAtomData &data1 = EAMdata[ EAMtids[sdat.atid] ];
314 >            
315 >    data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho), *(sdat.dfrhodrho) );
316  
317 <    fdat.frho = result.first;
318 <    fdat.dfrhodrho = result.second;
317 >    (*(sdat.pot))[METALLIC_FAMILY] += *(sdat.frho);
318 >    if (sdat.doParticlePot) {
319 >      *(sdat.particlePot) += *(sdat.frho);
320 >    }
321 >
322      return;
323    }
324  
325  
326 <  void EAM::calcForce(InteractionData idat) {
326 >  void EAM::calcForce(InteractionData &idat) {
327  
328      if (!initialized_) initialize();
329  
330 <    pair<RealType, RealType> res;
330 >    if (haveCutoffRadius_)
331 >      if ( *(idat.rij) > eamRcut_) return;
332 >  
333 >
334 >    int eamtid1 = EAMtids[idat.atid1];
335 >    int eamtid2 = EAMtids[idat.atid2];
336      
337 <    if (idat.rij < eamRcut_) {
337 >    EAMAtomData &data1 = EAMdata[eamtid1];
338 >    EAMAtomData &data2 = EAMdata[eamtid2];
339 >    
340 >    // get type-specific cutoff radii
341 >    
342 >    RealType rci = data1.rcut;
343 >    RealType rcj = data2.rcut;
344 >    
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 >    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 >    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  
362 <      EAMAtomData data1 = EAMMap[idat.atype1];
363 <      EAMAtomData data2 = EAMMap[idat.atype2];
390 <
391 <      // get type-specific cutoff radii
392 <
393 <      RealType rci = data1.rcut;
394 <      RealType rcj = data2.rcut;
362 >    switch(mixMeth_) {
363 >    case eamJohnson:
364        
365 <      RealType rha, drha, rhb, drhb;
366 <      RealType pha, dpha, phb, dphb;
367 <      RealType phab, dvpdr;
368 <      RealType drhoidr, drhojdr, dudr;
365 >      if ( *(idat.rij) < rci) {
366 >        phab = phab + 0.5 * (rhb / rha) * pha;
367 >        dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
368 >                             pha*((drhb/rha) - (rhb*drha/rha/rha)));
369 >      }
370        
371 <      if (idat.rij < rci) {
372 <        res = data1.rho->getValueAndDerivativeAt(idat.rij);
373 <        rha = res.first;
374 <        drha = res.second;
375 <
376 <        res = MixingMap[make_pair(idat.atype1, idat.atype1)].phi->getValueAndDerivativeAt(idat.rij);
407 <        pha = res.first;
408 <        dpha = res.second;
371 >      
372 >      
373 >      if ( *(idat.rij) < rcj) {
374 >        phab = phab + 0.5 * (rha / rhb) * phb;
375 >        dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
376 >                               phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
377        }
378 <
379 <      if (idat.rij < rcj) {
380 <        res = data2.rho->getValueAndDerivativeAt(idat.rij);
381 <        rhb = res.first;
382 <        drhb = res.second;
383 <
384 <        res = MixingMap[make_pair(idat.atype2, idat.atype2)].phi->getValueAndDerivativeAt(idat.rij);
385 <        phb = res.first;
418 <        dphb = res.second;
378 >      
379 >      break;
380 >      
381 >    case eamDaw:
382 >      
383 >      if ( *(idat.rij) <  MixingMap[eamtid1][eamtid2].rcut) {
384 >        MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij),
385 >                                                                  phab, dvpdr);
386        }
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 –      }
387        
388 <      drhoidr = drha;
389 <      drhojdr = drhb;
388 >      break;
389 >    case eamUnknown:
390 >    default:
391 >      
392 >      sprintf(painCave.errMsg,
393 >              "EAM::calcForce hit a mixing method it doesn't know about!\n"
394 >              );
395 >      painCave.severity = OPENMD_ERROR;
396 >      painCave.isFatal = 1;
397 >      simError();        
398 >      
399 >    }
400 >    
401 >    drhoidr = drha;
402 >    drhojdr = drhb;
403 >    
404 >    dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
405 >    
406 >    *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
407  
462 –      dudr = drhojdr*idat.dfrho1 + drhoidr*idat.dfrho2 + dvpdr;
463 –
464 –      idat.f1 = idat.d * dudr / idat.rij;
408          
409 <      // particle_pot is the difference between the full potential
410 <      // and the full potential without the presence of a particular
409 >    if (idat.doParticlePot) {
410 >      // particlePot is the difference between the full potential and
411 >      // the full potential without the presence of a particular
412        // particle (atom1).
413        //
414 <      // This reduces the density at other particle locations, so
415 <      // we need to recompute the density at atom2 assuming atom1
416 <      // didn't contribute.  This then requires recomputing the
417 <      // density functional for atom2 as well.
418 <      //
419 <      // Most of the particle_pot heavy lifting comes from the
420 <      // pair interaction, and will be handled by vpair.
421 <    
422 <      idat.fshift1 = data1.F->getValueAt( idat.rho1 - rhb );
423 <      idat.fshift2 = data1.F->getValueAt( idat.rho2 - rha );
480 <
481 <      idat.pot += phab;
482 <
483 <      idat.vpair += phab;
414 >      // This reduces the density at other particle locations, so we
415 >      // need to recompute the density at atom2 assuming atom1 didn't
416 >      // contribute.  This then requires recomputing the density
417 >      // functional for atom2 as well.
418 >      
419 >      *(idat.particlePot1) += data2.F->getValueAt( *(idat.rho2) - rha )
420 >        - *(idat.frho2);
421 >      
422 >      *(idat.particlePot2) += data1.F->getValueAt( *(idat.rho1) - rhb)
423 >        - *(idat.frho1);
424      }
425 <
425 >    
426 >    (*(idat.pot))[METALLIC_FAMILY] += phab;
427 >    
428 >    *(idat.vpair) += phab;
429 >  
430      return;
431      
432    }
433  
434 <  RealType EAM::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
434 >  RealType EAM::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
435      if (!initialized_) initialize();  
436  
437      RealType cut = 0.0;
438  
439 <    map<AtomType*, EAMAtomData>::iterator it;
440 <
441 <    it = EAMMap.find(at1);
442 <    if (it != EAMMap.end()) {
443 <      EAMAtomData data1 = (*it).second;
439 >    int atid1 = atypes.first->getIdent();
440 >    int atid2 = atypes.second->getIdent();
441 >    int eamtid1 = EAMtids[atid1];
442 >    int eamtid2 = EAMtids[atid2];
443 >    
444 >    if (eamtid1 != -1) {
445 >      EAMAtomData data1 = EAMdata[eamtid1];
446        cut = data1.rcut;
447      }
448  
449 <    it = EAMMap.find(at2);
450 <    if (it != EAMMap.end()) {
505 <      EAMAtomData data2 = (*it).second;
449 >    if (eamtid2 != -1) {
450 >      EAMAtomData data2 = EAMdata[eamtid2];
451        if (data2.rcut > cut)
452          cut = data2.rcut;
453      }
454 <
454 >    
455      return cut;
456    }
457   }

Diff Legend

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