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 1629 by gezelter, Wed Sep 14 21:15:17 2011 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 <  ForceField* EAM::forceField_ = NULL;
54 <  std::map<int, AtomType*> EAM::EAMlist;
55 <  std::map<AtomType*, EAMAtomData> EAM::EAMMap;
56 <  std::map<std::pair<AtomType*, AtomType*>, EAMInteractionData> EAM::MixingMap;
53 >  EAM::EAM() : name_("EAM"), initialized_(false), forceField_(NULL),
54 >               mixMeth_(eamJohnson), eamRcut_(0.0), haveCutoffRadius_(false) {}
55    
58 –  EAM* EAM::_instance = NULL;
59 –
60 –  EAM* EAM::Instance() {
61 –    if (!_instance) {
62 –      _instance = new EAM();
63 –    }
64 –    return _instance;
65 –  }
66 –  
56    EAMParam EAM::getEAMParam(AtomType* atomType) {
57      
58      // Do sanity checking on the AtomType we were passed before
# Line 107 | Line 96 | namespace OpenMD {
96      RealType dr = eamParam.dr;
97      vector<RealType> rvals;
98      
99 <    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
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(i * dr);
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);
# Line 135 | Line 129 | namespace OpenMD {
129      vector<RealType> scaledF;
130      
131      for (int i = 0; i < nrho; i++) {
132 <      rhovals.push_back(i * drho);
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, eamParam.F);
137 >    cs->addPoints(rhovals, scaledF);
138      return cs;
139    }
140    
# Line 152 | Line 146 | namespace OpenMD {
146  
147      // make the r grid:
148  
155 –    // set rcut to be the smaller of the two atomic rcuts
149  
150 <    RealType rcut = eamParam1.rcut < eamParam2.rcut ?
151 <      eamParam1.rcut : eamParam2.rcut;
150 >    // we need phi out to the largest value we'll encounter in the radial space;
151 >    
152 >    RealType rmax = 0.0;
153 >    rmax = max(rmax, eamParam1.rcut);
154 >    rmax = max(rmax, eamParam1.nr * eamParam1.dr);
155  
156 +    rmax = max(rmax, eamParam2.rcut);
157 +    rmax = max(rmax, eamParam2.nr * eamParam2.dr);
158 +
159      // use the smallest dr (finest grid) to build our grid:
160  
161 <    RealType dr = eamParam1.dr < eamParam2.dr ? eamParam1.dr : eamParam2.dr;
162 <    int nr = int(rcut/dr);
161 >    RealType dr = min(eamParam1.dr, eamParam2.dr);
162 >
163 >    int nr = int(rmax/dr + 0.5);
164 >
165      vector<RealType> rvals;
166 <    for (int i = 0; i < nr; i++) rvals.push_back(i*dr);
166 >    for (int i = 0; i < nr; i++) rvals.push_back(RealType(i*dr));
167  
168      // construct the pair potential:
169  
# Line 175 | Line 176 | namespace OpenMD {
176  
177      for (int i = 1; i < rvals.size(); i++ ) {
178        r = rvals[i];
178 –      zi = z1->getValueAt(r);
179 –      zj = z2->getValueAt(r);
179  
180 +      // only use z(r) if we're inside this atom's cutoff radius,
181 +      // otherwise, we'll use zero for the charge.  This effectively
182 +      // means that our phi grid goes out beyond the cutoff of the
183 +      // pair potential
184 +
185 +      zi = r <= eamParam1.rcut ? z1->getValueAt(r) : 0.0;
186 +      zj = r <= eamParam2.rcut ? z2->getValueAt(r) : 0.0;
187 +
188        phi = 331.999296 * (zi * zj) / r;
189 +
190        phivals.push_back(phi);
191      }
192        
193      CubicSpline* cs = new CubicSpline();
194      cs->addPoints(rvals, phivals);
195      return cs;
196 +  }
197 +
198 +  void EAM::setCutoffRadius( RealType rCut ) {
199 +    eamRcut_ = rCut;
200 +    haveCutoffRadius_ = true;
201    }
202  
203    void EAM::initialize() {
204  
205      // set up the mixing method:
206 <    ForceFieldOptions ffo = forceField_->getForceFieldOptions();
207 <    string EAMMixMeth = toUpperCopy(ffo.getEAMMixingMethod());
208 <
206 >    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
207 >    string EAMMixMeth = fopts.getEAMMixingMethod();
208 >    toUpper(EAMMixMeth);
209 >  
210      if (EAMMixMeth == "JOHNSON")
211        mixMeth_ = eamJohnson;    
212      else if (EAMMixMeth == "DAW")
# Line 222 | Line 236 | namespace OpenMD {
236        
237        if (nbt->isEAM()) {
238          
239 <        std::pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
239 >        pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
240          
241          GenericData* data = nbt->getPropertyByName("EAM");
242          if (data == NULL) {
# Line 247 | Line 261 | namespace OpenMD {
261            simError();          
262          }
263          
264 <        EAMMix eamParam = eamData->getData();
264 >        EAMMixingParam eamParam = eamData->getData();
265  
266 <        vector<RealType> phiAB = eamParam.phiAB;
266 >        vector<RealType> phiAB = eamParam.phi;
267          RealType dr = eamParam.dr;
268          int nr = eamParam.nr;
269  
# Line 264 | Line 278 | namespace OpenMD {
278    void EAM::addType(AtomType* atomType){
279  
280      EAMAtomData eamAtomData;
281 <
281 >    
282      eamAtomData.rho = getRho(atomType);
283      eamAtomData.F = getF(atomType);
284      eamAtomData.Z = getZ(atomType);
# Line 273 | Line 287 | namespace OpenMD {
287      // add it to the map:
288      AtomTypeProperties atp = atomType->getATP();    
289  
290 <    std::pair<std::map<int,AtomType*>::iterator,bool> ret;    
291 <    ret = EAMlist.insert( std::pair<int, AtomType*>(atp.ident, atomType) );
290 >    pair<map<int,AtomType*>::iterator,bool> ret;    
291 >    ret = EAMlist.insert( pair<int, AtomType*>(atp.ident, atomType) );
292      if (ret.second == false) {
293        sprintf( painCave.errMsg,
294                 "EAM already had a previous entry with ident %d\n",
# Line 288 | Line 302 | namespace OpenMD {
302      
303      // Now, iterate over all known types and add to the mixing map:
304      
305 <    std::map<int, AtomType*>::iterator it;
305 >    map<AtomType*, EAMAtomData>::iterator it;
306      for( it = EAMMap.begin(); it != EAMMap.end(); ++it) {
307        
308 <      AtomType* atype2 = (*it).second;
308 >      AtomType* atype2 = (*it).first;
309  
310        EAMInteractionData mixer;
311        mixer.phi = getPhi(atomType, atype2);
312        mixer.explicitlySet = false;
313  
314 <      std::pair<AtomType*, AtomType*> key1, key2;
315 <      key1 = std::make_pair(atomType, atype2);
316 <      key2 = std::make_pair(atype2, atomType);
314 >      pair<AtomType*, AtomType*> key1, key2;
315 >      key1 = make_pair(atomType, atype2);
316 >      key2 = make_pair(atype2, atomType);
317        
318        MixingMap[key1] = mixer;
319        if (key2 != key1) {
# Line 319 | Line 333 | namespace OpenMD {
333  
334      EAMInteractionData mixer;
335      CubicSpline* cs = new CubicSpline();
336 <    vector<RealType> rvals;
336 >    vector<RealType> rVals;
337  
338 <    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
338 >    for (int i = 0; i < nr; i++) rVals.push_back(i * dr);
339  
340      cs->addPoints(rVals, phiVals);
341      mixer.phi = cs;
342      mixer.explicitlySet = true;
343  
344 <    std::pair<AtomType*, AtomType*> key1, key2;
345 <    key1 = std::make_pair(atype1, atype2);
346 <    key2 = std::make_pair(atype2, atype1);
344 >    pair<AtomType*, AtomType*> key1, key2;
345 >    key1 = make_pair(atype1, atype2);
346 >    key2 = make_pair(atype2, atype1);
347      
348      MixingMap[key1] = mixer;
349      if (key2 != key1) {
# Line 338 | Line 352 | namespace OpenMD {
352      return;
353    }
354  
355 <  void EAM::calcDensity(AtomType* at1, AtomType* at2, Vector3d d,
356 <                        RealType rij, RealType r2, RealType rho_i_at_j,
343 <                        RealType rho_j_at_i) {
344 <
355 >  void EAM::calcDensity(InteractionData &idat) {
356 >    
357      if (!initialized_) initialize();
358 <
359 <    EAMAtomData data1 = EAMMap[at1];
360 <    EAMAtomData data2 = EAMMap[at2];
361 <
362 <    if (rij < data1.rcut) rho_i_at_j = data1.rho->getValueAt(rij);
363 <    if (rij < data2.rcut) rho_j_at_i = data2.rho->getValueAt(rij);
364 <    return;
358 >    
359 >    EAMAtomData data1 = EAMMap[idat.atypes.first];
360 >    EAMAtomData data2 = EAMMap[idat.atypes.second];
361 >    
362 >    if (haveCutoffRadius_)
363 >      if ( *(idat.rij) > eamRcut_) return;
364 >    
365 >    if ( *(idat.rij) < data1.rcut)
366 >      *(idat.rho1) += data1.rho->getValueAt( *(idat.rij));
367 >    
368 >      
369 >    if ( *(idat.rij) < data2.rcut)
370 >      *(idat.rho2) += data2.rho->getValueAt( *(idat.rij));
371 >    
372 >    return;  
373    }
374 <
375 <  void EAM::calcFunctional(AtomType* at1, RealType rho, RealType frho,
376 <                           RealType dfrhodrho) {
357 <
374 >  
375 >  void EAM::calcFunctional(SelfData &sdat) {
376 >    
377      if (!initialized_) initialize();
378  
379 <    EAMAtomData data1 = EAMMap[at1];
379 >    EAMAtomData data1 = EAMMap[ sdat.atype ];
380          
381 <    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt(rho);
381 >    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt( *(sdat.rho) );
382  
383 <    frho = result.first;
384 <    dfrhodrho = result.second;
383 >    *(sdat.frho) = result.first;
384 >    *(sdat.dfrhodrho) = result.second;
385 >
386 >    (*(sdat.pot))[METALLIC_FAMILY] += result.first;
387 >    *(sdat.particlePot) += result.first;
388 >
389      return;
390    }
391  
392  
393 <  void EAM::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
371 <                      RealType rij, RealType r2, RealType sw,
372 <                      RealType &vpair, RealType &pot, Vector3d &f1,
373 <                      RealType rho1, RealType rho2, RealType dfrho1,
374 <                      RealType dfrho2, RealType fshift1, RealType fshift2) {
393 >  void EAM::calcForce(InteractionData &idat) {
394  
395      if (!initialized_) initialize();
396 <    
396 >
397 >    if (haveCutoffRadius_)
398 >      if ( *(idat.rij) > eamRcut_) return;
399 >  
400      pair<RealType, RealType> res;
401      
402 <    if (rij < eamRcut_) {
403 <
404 <      EAMAtomData data1 = EAMMap[at1];
405 <      EAMAtomData data2 = EAMMap[at2];
406 <
407 <      // get type-specific cutoff radii
408 <
409 <      RealType rci = data1.rcut;
410 <      RealType rcj = data2.rcut;
402 >    EAMAtomData data1 = EAMMap[idat.atypes.first];
403 >    EAMAtomData data2 = EAMMap[idat.atypes.second];
404 >    
405 >    // get type-specific cutoff radii
406 >    
407 >    RealType rci = data1.rcut;
408 >    RealType rcj = data2.rcut;
409 >    
410 >    RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0);
411 >    RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0);
412 >    RealType phab(0.0), dvpdr(0.0);
413 >    RealType drhoidr, drhojdr, dudr;
414 >    
415 >    if ( *(idat.rij) < rci) {
416 >      res = data1.rho->getValueAndDerivativeAt( *(idat.rij));
417 >      rha = res.first;
418 >      drha = res.second;
419        
420 <      RealType rha, drha, rhb, drhb;
421 <      RealType pha, dpha, phb, dphb;
422 <      RealType phab, dvpdr;
423 <      RealType drhoidr, drhojdr, dudr;
420 >      res = MixingMap[make_pair(idat.atypes.first, idat.atypes.first)].phi->getValueAndDerivativeAt( *(idat.rij) );
421 >      pha = res.first;
422 >      dpha = res.second;
423 >    }
424 >    
425 >    if ( *(idat.rij) < rcj) {
426 >      res = data2.rho->getValueAndDerivativeAt( *(idat.rij) );
427 >      rhb = res.first;
428 >      drhb = res.second;
429        
430 <      if (rij < rci) {
431 <        res = data1.rho->getValueAndDerivativeAt(rij);
432 <        rha = res.first;
433 <        drha = res.second;
430 >      res = MixingMap[make_pair(idat.atypes.second, idat.atypes.second)].phi->getValueAndDerivativeAt( *(idat.rij) );
431 >      phb = res.first;
432 >      dphb = res.second;
433 >    }
434  
435 <        res = MixingMap[make_pair(at1, at1)].phi->getValueAndDerivativeAt(rij);
436 <        pha = res.first;
437 <        dpha = res.second;
435 >    switch(mixMeth_) {
436 >    case eamJohnson:
437 >      
438 >      if ( *(idat.rij) < rci) {
439 >        phab = phab + 0.5 * (rhb / rha) * pha;
440 >        dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
441 >                             pha*((drhb/rha) - (rhb*drha/rha/rha)));
442        }
443 <
444 <      if (rij < rcj) {
445 <        res = data2.rho->getValueAndDerivativeAt(rij);
446 <        rhb = res.first;
447 <        drhb = res.second;
448 <
449 <        res = MixingMap[make_pair(at2, at2)].phi->getValueAndDerivativeAt(rij);
411 <        phb = res.first;
412 <        dphb = res.second;
413 <      }
414 <
415 <      phab = 0.0;
416 <      dvpdr = 0.0;
417 <
418 <      switch(mixMeth_) {
419 <      case eamJohnson:
420 <      
421 <        if (rij < rci) {
422 <          phab = phab + 0.5 * (rhb / rha) * pha;
423 <          dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
424 <                               pha*((drhb/rha) - (rhb*drha/rha/rha)));
425 <        }
426 <
427 <        if (rij < rcj) {
428 <          phab = phab + 0.5 * (rha / rhb) * phb;
429 <          dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
430 <                                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
431 <        }
432 <
433 <        break;
434 <
435 <      case eamDaw:
436 <                
437 <        res = MixingMap[make_pair(at1,at2)].phi->getValueAndDerivativeAt(rij);
438 <        phab = res.first;
439 <        dvpdr = res.second;
440 <
441 <        break;
442 <      case eamUnknown:
443 <      default:
444 <
445 <        sprintf(painCave.errMsg,
446 <                "EAM::calcForce hit a mixing method it doesn't know about!\n"
447 <                );
448 <        painCave.severity = OPENMD_ERROR;
449 <        painCave.isFatal = 1;
450 <        simError();        
451 <          
443 >      
444 >      
445 >      
446 >      if ( *(idat.rij) < rcj) {
447 >        phab = phab + 0.5 * (rha / rhb) * phb;
448 >        dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
449 >                               phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
450        }
451        
452 <      drhoidr = drha;
453 <      drhojdr = drhb;
454 <
455 <      dudr = drhojdr*dfrhodrho_i + drhoidr*dfrhodrho_j + dvpdr;
456 <
457 <      f1 = d * dudr / rij;
458 <        
459 <      // particle_pot is the difference between the full potential
460 <      // and the full potential without the presence of a particular
461 <      // particle (atom1).
462 <      //
463 <      // This reduces the density at other particle locations, so
464 <      // we need to recompute the density at atom2 assuming atom1
465 <      // didn't contribute.  This then requires recomputing the
466 <      // density functional for atom2 as well.
467 <      //
468 <      // Most of the particle_pot heavy lifting comes from the
469 <      // pair interaction, and will be handled by vpair.
472 <    
473 <      fshift_i = data1.F->getValueAt( rho_i - rhb );
474 <      fshift_j = data1.F->getValueAt( rho_j - rha );
475 <
476 <      pot += phab;
477 <
478 <      vpair += phab;
452 >      break;
453 >      
454 >    case eamDaw:
455 >      res = MixingMap[idat.atypes].phi->getValueAndDerivativeAt( *(idat.rij));
456 >      phab = res.first;
457 >      dvpdr = res.second;
458 >      
459 >      break;
460 >    case eamUnknown:
461 >    default:
462 >      
463 >      sprintf(painCave.errMsg,
464 >              "EAM::calcForce hit a mixing method it doesn't know about!\n"
465 >              );
466 >      painCave.severity = OPENMD_ERROR;
467 >      painCave.isFatal = 1;
468 >      simError();        
469 >      
470      }
480 –
481 –    return;
471      
472 <  }
473 <
485 <
486 <  void EAM::calc_eam_prepair_rho(int *atid1, int *atid2, RealType *d,
487 <                                 RealType *rij, RealType *r2,
488 <                                 RealType* rho_i_at_j, RealType* rho_j_at_i){
489 <    if (!initialized_) initialize();
490 <
491 <    AtomType* atype1 = EAMlist[*atid1];
492 <    AtomType* atype2 = EAMlist[*atid2];
472 >    drhoidr = drha;
473 >    drhojdr = drhb;
474      
475 <    Vector3d disp(d[0], d[1], d[2]);
495 <
496 <    calcDensity(atype1, atype2, disp, *rij, *r2, *rho_i_at_j, *rho_j_at_i);
497 <
498 <    return;    
499 <  }
500 <
501 <  void EAM::calc_eam_preforce_Frho(int *atid1, RealType *rho, RealType *frho,
502 <                                   RealType *dfrhodrho) {
503 <
504 <    if (!initialized_) initialize();
505 <
506 <    AtomType* atype1 = EAMlist[*atid1];  
507 <
508 <    calcFunctional(atype1, *rho, *frho, *dfrhodrho);
475 >    dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
476      
477 <    return;    
478 <  }
479 <
480 <  void EAM::do_eam_pair(int *atid1, int *atid2, RealType *d, RealType *rij,
481 <                        RealType *r2, RealType *sw, RealType *vpair,
482 <                        RealType *pot, RealType *f1, RealType *rho1,
483 <                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
484 <                        RealType *fshift1, RealType *fshift2) {
485 <
486 <    if (!initialized_) initialize();
477 >    *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
478 >        
479 >    // particlePot is the difference between the full potential and
480 >    // the full potential without the presence of a particular
481 >    // particle (atom1).
482 >    //
483 >    // This reduces the density at other particle locations, so we
484 >    // need to recompute the density at atom2 assuming atom1 didn't
485 >    // contribute.  This then requires recomputing the density
486 >    // functional for atom2 as well.
487      
488 <    AtomType* atype1 = EAMMap[*atid1];
489 <    AtomType* atype2 = EAMMap[*atid2];
488 >    *(idat.particlePot1) += data2.F->getValueAt( *(idat.rho2) - rha )
489 >      - *(idat.frho2);
490      
491 <    Vector3d disp(d[0], d[1], d[2]);
492 <    Vector3d frc(f1[0], f1[1], f1[2]);
491 >    *(idat.particlePot2) += data1.F->getValueAt( *(idat.rho1) - rhb)
492 >      - *(idat.frho1);
493      
494 <    calcForce(atype1, atype2, disp, *rij, *r2, *sw, *vpair,  *pot, frc,
495 <              *rho1, *rho2, *dfrho1, *dfrho2, *fshift1, *fshift2);
496 <      
530 <    f1[0] = frc.x();
531 <    f1[1] = frc.y();
532 <    f1[2] = frc.z();
533 <
534 <    return;    
535 <  }
494 >    (*(idat.pot))[METALLIC_FAMILY] += phab;
495 >    
496 >    *(idat.vpair) += phab;
497    
498 <  void EAM::setCutoffEAM(RealType *thisRcut) {
499 <    eamRcut_ = thisRcut;
498 >    return;
499 >    
500    }
540 – }
501  
502 < extern "C" {
503 <  
544 < #define fortranCalcDensity FC_FUNC(calc_eam_prepair_rho, CALC_EAM_PREPAIR_RHO)
545 < #define fortranCalcFunctional FC_FUNC(calc_eam_preforce_frho, CALC_EAM_PREFORCE_FRHO)
546 < #define fortranCalcForce FC_FUNC(do_eam_pair, DO_EAM_PAIR)
547 < #define fortranSetCutoffEAM FC_FUNC(setcutoffeam, SETCUTOFFEAM)
548 <  
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) {
502 >  RealType EAM::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
503 >    if (!initialized_) initialize();  
504  
505 <    return OpenMD::EAM::Instance()->calc_eam_prepair_rho(*atid1, *atid2, *d,
554 <                                                         *rij, *r2,
555 <                                                         *rho_i_at_j,  
556 <                                                         *rho_j_at_i);
557 <  }
558 <  RealType fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho,
559 <                                 RealType *dfrhodrho) {  
505 >    RealType cut = 0.0;
506  
507 <    return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(*atid1,
562 <                                                           *rho,
563 <                                                           *frho,
564 <                                                           *dfrhodrho);
507 >    map<AtomType*, EAMAtomData>::iterator it;
508  
509 +    it = EAMMap.find(atypes.first);
510 +    if (it != EAMMap.end()) {
511 +      EAMAtomData data1 = (*it).second;
512 +      cut = data1.rcut;
513 +    }
514 +
515 +    it = EAMMap.find(atypes.second);
516 +    if (it != EAMMap.end()) {
517 +      EAMAtomData data2 = (*it).second;
518 +      if (data2.rcut > cut)
519 +        cut = data2.rcut;
520 +    }
521 +
522 +    return cut;
523    }
567 –  void fortranSetEAMCutoff(RealType *rcut) {
568 –    return OpenMD::EAM::Instance()->setCutoffEAM(rcut);
569 –  }
570 –  void fortranDoEAMPair(int *atid1, int *atid2, RealType *d, RealType *rij,
571 –                        RealType *r2, RealType *sw, RealType *vpair,
572 –                        RealType *pot, RealType *f1, RealType *rho1,
573 –                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
574 –                        RealType *fshift1, RealType *fshift2){
575 –    
576 –    return OpenMD::EAM::Instance()->do_eam_pair(*atid1, *atid2, *d, *rij,
577 –                                                *r2, *sw,  *vpair,
578 –                                                *pot, *f1,  *rho1,
579 –                                                *rho2,  *dfrho1,  *dfrho2,
580 –                                                *fshift1,  *fshift2);
581 –  }
524   }
525 +

Diff Legend

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