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 1479 by gezelter, Mon Jul 26 19:00:48 2010 UTC vs.
Revision 1710 by gezelter, Fri May 18 21:44:02 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 50 | Line 51 | namespace OpenMD {
51  
52   namespace OpenMD {
53  
54 <  bool EAM::initialized_ = false;
55 <  RealType EAM::eamRcut_ = 0.0;
55 <  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() : name_("EAM"), initialized_(false), forceField_(NULL),
55 >               mixMeth_(eamJohnson), eamRcut_(0.0), haveCutoffRadius_(false) {}
56    
57 <  EAM* EAM::_instance = NULL;
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 <  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 <  }
63 >    // make the r grid:
64  
108 –  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 –  }
65  
66 <  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;
66 >    // we need phi out to the largest value we'll encounter in the radial space;
67      
68 <    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
69 <      
70 <    CubicSpline* cs = new CubicSpline();
135 <    cs->addPoints(rvals, eamParam.rho);
136 <    return cs;
137 <  }
68 >    RealType rmax = 0.0;
69 >    rmax = max(rmax, ea1.getRcut());
70 >    rmax = max(rmax, ea1.getNr() * ea1.getDr());
71  
72 <  CubicSpline* EAM::getF(AtomType* atomType) {    
73 <    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);
72 >    rmax = max(rmax, ea2.getRcut());
73 >    rmax = max(rmax, ea2.getNr() * ea2.getDr());
74  
75 <    // make the r grid:
75 >    // use the smallest dr (finest grid) to build our grid:
76  
77 <    // set rcut to be the smaller of the two atomic rcuts
77 >    RealType dr = min(ea1.getDr(), ea2.getDr());
78  
79 <    RealType rcut = eamParam1.rcut < eamParam2.rcut ?
167 <      eamParam1.rcut : eamParam2.rcut;
79 >    int nr = int(rmax/dr + 0.5);
80  
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);
81      vector<RealType> rvals;
82 <    for (int i = 0; i < nr; i++) rvals.push_back(i*dr);
82 >    for (int i = 0; i < nr; i++) rvals.push_back(RealType(i*dr));
83  
84      // construct the pair potential:
85  
# Line 184 | Line 92 | namespace OpenMD {
92  
93      for (int i = 1; i < rvals.size(); i++ ) {
94        r = rvals[i];
187 –      zi = z1->getValueAt(r);
188 –      zj = z2->getValueAt(r);
95  
96 +      // only use z(r) if we're inside this atom's cutoff radius,
97 +      // otherwise, we'll use zero for the charge.  This effectively
98 +      // means that our phi grid goes out beyond the cutoff of the
99 +      // pair potential
100 +
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 +
106        phivals.push_back(phi);
107      }
108        
109      CubicSpline* cs = new CubicSpline();
110      cs->addPoints(rvals, phivals);
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:
122      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
123 <    string EAMMixMeth = toUpperCopy(fopts.getEAMMixingMethod());
124 <
123 >    string EAMMixMeth = fopts.getEAMMixingMethod();
124 >    toUpper(EAMMixMeth);
125 >  
126      if (EAMMixMeth == "JOHNSON")
127        mixMeth_ = eamJohnson;    
128      else if (EAMMixMeth == "DAW")
# Line 231 | Line 152 | namespace OpenMD {
152        
153        if (nbt->isEAM()) {
154          
155 <        std::pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
155 >        pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
156          
157          GenericData* data = nbt->getPropertyByName("EAM");
158          if (data == NULL) {
# 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 <    std::pair<std::map<int,AtomType*>::iterator,bool> ret;    
207 <    ret = EAMlist.insert( std::pair<int, AtomType*>(atp.ident, atomType) );
206 >    pair<map<int,AtomType*>::iterator,bool> ret;    
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 297 | Line 218 | namespace OpenMD {
218      
219      // Now, iterate over all known types and add to the mixing map:
220      
221 <    std::map<AtomType*, EAMAtomData>::iterator it;
221 >    map<AtomType*, EAMAtomData>::iterator it;
222      for( it = EAMMap.begin(); it != EAMMap.end(); ++it) {
223        
224        AtomType* atype2 = (*it).first;
# Line 306 | Line 227 | namespace OpenMD {
227        mixer.phi = getPhi(atomType, atype2);
228        mixer.explicitlySet = false;
229  
230 <      std::pair<AtomType*, AtomType*> key1, key2;
231 <      key1 = std::make_pair(atomType, atype2);
232 <      key2 = std::make_pair(atype2, atomType);
230 >      pair<AtomType*, AtomType*> key1, key2;
231 >      key1 = make_pair(atomType, atype2);
232 >      key2 = make_pair(atype2, atomType);
233        
234        MixingMap[key1] = mixer;
235        if (key2 != key1) {
# Line 336 | Line 257 | namespace OpenMD {
257      mixer.phi = cs;
258      mixer.explicitlySet = true;
259  
260 <    std::pair<AtomType*, AtomType*> key1, key2;
261 <    key1 = std::make_pair(atype1, atype2);
262 <    key2 = std::make_pair(atype2, atype1);
260 >    pair<AtomType*, AtomType*> key1, key2;
261 >    key1 = make_pair(atype1, atype2);
262 >    key2 = make_pair(atype2, atype1);
263      
264      MixingMap[key1] = mixer;
265      if (key2 != key1) {
# Line 347 | Line 268 | namespace OpenMD {
268      return;
269    }
270  
271 <  void EAM::calcDensity(AtomType* at1, AtomType* at2, const RealType rij,
351 <                        RealType &rho_i_at_j, RealType &rho_j_at_i) {
271 >  void EAM::calcDensity(InteractionData &idat) {
272      
273      if (!initialized_) initialize();
274      
275 <    EAMAtomData data1 = EAMMap[at1];
276 <    EAMAtomData data2 = EAMMap[at2];
277 <
278 <    if (rij < data1.rcut) rho_i_at_j = data1.rho->getValueAt(rij);
279 <    if (rij < data2.rcut) rho_j_at_i = data2.rho->getValueAt(rij);
280 <    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(AtomType* at1, RealType rho, RealType &frho,
292 <                           RealType &dfrhodrho) {
365 <
290 >  
291 >  void EAM::calcFunctional(SelfData &sdat) {
292 >    
293      if (!initialized_) initialize();
294  
295 <    EAMAtomData data1 = EAMMap[at1];
295 >    EAMAtomData data1 = EAMMap[ sdat.atype ];
296          
297 <    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt(rho);
297 >    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt( *(sdat.rho) );
298  
299 <    frho = result.first;
300 <    dfrhodrho = result.second;
299 >    *(sdat.frho) = result.first;
300 >    *(sdat.dfrhodrho) = result.second;
301 >
302 >    (*(sdat.pot))[METALLIC_FAMILY] += result.first;
303 >    *(sdat.particlePot) += result.first;
304 >
305      return;
306    }
307  
308  
309 <  void EAM::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
379 <                      RealType rij, RealType r2, RealType sw,
380 <                      RealType &vpair, RealType &pot, Vector3d &f1,
381 <                      RealType rho_i, RealType rho_j,
382 <                      RealType dfrhodrho_i, RealType dfrhodrho_j,
383 <                      RealType &fshift_i, RealType &fshift_j) {
309 >  void EAM::calcForce(InteractionData &idat) {
310  
311      if (!initialized_) initialize();
312 <    
312 >
313 >    if (haveCutoffRadius_)
314 >      if ( *(idat.rij) > eamRcut_) return;
315 >  
316      pair<RealType, RealType> res;
317      
318 <    if (rij < eamRcut_) {
319 <
320 <      EAMAtomData data1 = EAMMap[at1];
321 <      EAMAtomData data2 = EAMMap[at2];
322 <
323 <      // get type-specific cutoff radii
324 <
325 <      RealType rci = data1.rcut;
326 <      RealType rcj = data2.rcut;
318 >    EAMAtomData data1 = EAMMap[idat.atypes.first];
319 >    EAMAtomData data2 = EAMMap[idat.atypes.second];
320 >    
321 >    // get type-specific cutoff radii
322 >    
323 >    RealType rci = data1.rcut;
324 >    RealType rcj = data2.rcut;
325 >    
326 >    RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0);
327 >    RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0);
328 >    RealType phab(0.0), dvpdr(0.0);
329 >    RealType drhoidr, drhojdr, dudr;
330 >    
331 >    if ( *(idat.rij) < rci) {
332 >      res = data1.rho->getValueAndDerivativeAt( *(idat.rij));
333 >      rha = res.first;
334 >      drha = res.second;
335        
336 <      RealType rha, drha, rhb, drhb;
337 <      RealType pha, dpha, phb, dphb;
338 <      RealType phab, dvpdr;
339 <      RealType drhoidr, drhojdr, dudr;
340 <      
341 <      if (rij < rci) {
342 <        res = data1.rho->getValueAndDerivativeAt(rij);
343 <        rha = res.first;
344 <        drha = res.second;
336 >      res = MixingMap[make_pair(idat.atypes.first, idat.atypes.first)].phi->getValueAndDerivativeAt( *(idat.rij) );
337 >      pha = res.first;
338 >      dpha = res.second;
339 >    }
340 >    
341 >    if ( *(idat.rij) < rcj) {
342 >      res = data2.rho->getValueAndDerivativeAt( *(idat.rij) );
343 >      rhb = res.first;
344 >      drhb = res.second;
345 >      
346 >      res = MixingMap[make_pair(idat.atypes.second, idat.atypes.second)].phi->getValueAndDerivativeAt( *(idat.rij) );
347 >      phb = res.first;
348 >      dphb = res.second;
349 >    }
350  
351 <        res = MixingMap[make_pair(at1, at1)].phi->getValueAndDerivativeAt(rij);
352 <        pha = res.first;
353 <        dpha = res.second;
351 >    switch(mixMeth_) {
352 >    case eamJohnson:
353 >      
354 >      if ( *(idat.rij) < rci) {
355 >        phab = phab + 0.5 * (rhb / rha) * pha;
356 >        dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
357 >                             pha*((drhb/rha) - (rhb*drha/rha/rha)));
358        }
359 <
360 <      if (rij < rcj) {
361 <        res = data2.rho->getValueAndDerivativeAt(rij);
362 <        rhb = res.first;
363 <        drhb = res.second;
364 <
365 <        res = MixingMap[make_pair(at2, at2)].phi->getValueAndDerivativeAt(rij);
420 <        phb = res.first;
421 <        dphb = res.second;
359 >      
360 >      
361 >      
362 >      if ( *(idat.rij) < rcj) {
363 >        phab = phab + 0.5 * (rha / rhb) * phb;
364 >        dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
365 >                               phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
366        }
423 –
424 –      phab = 0.0;
425 –      dvpdr = 0.0;
426 –
427 –      switch(mixMeth_) {
428 –      case eamJohnson:
429 –      
430 –        if (rij < rci) {
431 –          phab = phab + 0.5 * (rhb / rha) * pha;
432 –          dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
433 –                               pha*((drhb/rha) - (rhb*drha/rha/rha)));
434 –        }
435 –
436 –        if (rij < rcj) {
437 –          phab = phab + 0.5 * (rha / rhb) * phb;
438 –          dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
439 –                                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
440 –        }
441 –
442 –        break;
443 –
444 –      case eamDaw:
445 –                
446 –        res = MixingMap[make_pair(at1,at2)].phi->getValueAndDerivativeAt(rij);
447 –        phab = res.first;
448 –        dvpdr = res.second;
449 –
450 –        break;
451 –      case eamUnknown:
452 –      default:
453 –
454 –        sprintf(painCave.errMsg,
455 –                "EAM::calcForce hit a mixing method it doesn't know about!\n"
456 –                );
457 –        painCave.severity = OPENMD_ERROR;
458 –        painCave.isFatal = 1;
459 –        simError();        
460 –          
461 –      }
367        
368 <      drhoidr = drha;
369 <      drhojdr = drhb;
370 <
371 <      dudr = drhojdr*dfrhodrho_i + drhoidr*dfrhodrho_j + dvpdr;
372 <
373 <      f1 = d * dudr / rij;
374 <        
375 <      // particle_pot is the difference between the full potential
376 <      // and the full potential without the presence of a particular
377 <      // particle (atom1).
378 <      //
379 <      // This reduces the density at other particle locations, so
380 <      // we need to recompute the density at atom2 assuming atom1
381 <      // didn't contribute.  This then requires recomputing the
382 <      // density functional for atom2 as well.
383 <      //
384 <      // Most of the particle_pot heavy lifting comes from the
385 <      // pair interaction, and will be handled by vpair.
481 <    
482 <      fshift_i = data1.F->getValueAt( rho_i - rhb );
483 <      fshift_j = data1.F->getValueAt( rho_j - rha );
484 <
485 <      pot += phab;
486 <
487 <      vpair += phab;
368 >      break;
369 >      
370 >    case eamDaw:
371 >      res = MixingMap[idat.atypes].phi->getValueAndDerivativeAt( *(idat.rij));
372 >      phab = res.first;
373 >      dvpdr = res.second;
374 >      
375 >      break;
376 >    case eamUnknown:
377 >    default:
378 >      
379 >      sprintf(painCave.errMsg,
380 >              "EAM::calcForce hit a mixing method it doesn't know about!\n"
381 >              );
382 >      painCave.severity = OPENMD_ERROR;
383 >      painCave.isFatal = 1;
384 >      simError();        
385 >      
386      }
489 –
490 –    return;
387      
388 <  }
389 <
494 <
495 <  void EAM::calc_eam_prepair_rho(int *atid1, int *atid2, RealType *rij,
496 <                                 RealType* rho_i_at_j, RealType* rho_j_at_i){
497 <
498 <    if (!initialized_) initialize();
388 >    drhoidr = drha;
389 >    drhojdr = drhb;
390      
391 <    AtomType* atype1 = EAMlist[*atid1];
501 <    AtomType* atype2 = EAMlist[*atid2];
391 >    dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
392      
393 <    calcDensity(atype1, atype2, *rij, *rho_i_at_j, *rho_j_at_i);
394 <
395 <    return;    
396 <  }
397 <
398 <  void EAM::calc_eam_preforce_Frho(int *atid1, RealType *rho, RealType *frho,
399 <                                   RealType *dfrhodrho) {
400 <
401 <    if (!initialized_) initialize();
402 <
513 <    AtomType* atype1 = EAMlist[*atid1];  
514 <
515 <    calcFunctional(atype1, *rho, *frho, *dfrhodrho);
393 >    *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
394 >        
395 >    // particlePot is the difference between the full potential and
396 >    // the full potential without the presence of a particular
397 >    // particle (atom1).
398 >    //
399 >    // This reduces the density at other particle locations, so we
400 >    // need to recompute the density at atom2 assuming atom1 didn't
401 >    // contribute.  This then requires recomputing the density
402 >    // functional for atom2 as well.
403      
404 <    return;    
405 <  }
519 <  RealType EAM::getEAMcut(int *atid1) {
520 <
521 <    if (!initialized_) initialize();
404 >    *(idat.particlePot1) += data2.F->getValueAt( *(idat.rho2) - rha )
405 >      - *(idat.frho2);
406      
407 <    AtomType* atype1 = EAMlist[*atid1];  
408 <      
525 <    return getRcut(atype1);
526 <  }
527 <
528 <  void EAM::do_eam_pair(int *atid1, int *atid2, RealType *d, RealType *rij,
529 <                        RealType *r2, RealType *sw, RealType *vpair,
530 <                        RealType *pot, RealType *f1, RealType *rho1,
531 <                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
532 <                        RealType *fshift1, RealType *fshift2) {
533 <
534 <    if (!initialized_) initialize();
407 >    *(idat.particlePot2) += data1.F->getValueAt( *(idat.rho1) - rhb)
408 >      - *(idat.frho1);
409      
410 <    AtomType* atype1 = EAMlist[*atid1];
537 <    AtomType* atype2 = EAMlist[*atid2];
410 >    (*(idat.pot))[METALLIC_FAMILY] += phab;
411      
412 <    Vector3d disp(d[0], d[1], d[2]);
413 <    Vector3d frc(f1[0], f1[1], f1[2]);
412 >    *(idat.vpair) += phab;
413 >  
414 >    return;
415      
542 –    calcForce(atype1, atype2, disp, *rij, *r2, *sw, *vpair,  *pot, frc,
543 –              *rho1, *rho2, *dfrho1, *dfrho2, *fshift1, *fshift2);
544 –      
545 –    f1[0] = frc.x();
546 –    f1[1] = frc.y();
547 –    f1[2] = frc.z();
548 –
549 –    return;    
416    }
551 –  
552 –  void EAM::setCutoffEAM(RealType *thisRcut) {
553 –    eamRcut_ = *thisRcut;
554 –  }
555 – }
417  
418 < extern "C" {
419 <  
559 < #define fortranCalcDensity FC_FUNC(calc_eam_prepair_rho, CALC_EAM_PREPAIR_RHO)
560 < #define fortranCalcFunctional FC_FUNC(calc_eam_preforce_frho, CALC_EAM_PREFORCE_FRHO)
561 < #define fortranCalcForce FC_FUNC(do_eam_pair, DO_EAM_PAIR)
562 < #define fortranSetCutoffEAM FC_FUNC(setcutoffeam, SETCUTOFFEAM)
563 < #define fortranGetEAMcut FC_FUNC(geteamcut, GETEAMCUT)
418 >  RealType EAM::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
419 >    if (!initialized_) initialize();  
420  
421 <  
566 <  void fortranCalcDensity(int *atid1, int *atid2, RealType *rij,
567 <                          RealType *rho_i_at_j, RealType *rho_j_at_i) {
568 <    
569 <    return OpenMD::EAM::Instance()->calc_eam_prepair_rho(atid1, atid2, rij,
570 <                                                         rho_i_at_j,  
571 <                                                         rho_j_at_i);
572 <  }
573 <  void fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho,
574 <                             RealType *dfrhodrho) {  
575 <    
576 <    return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(atid1, rho, frho,
577 <                                                           dfrhodrho);
578 <    
579 <  }
580 <  void fortranSetCutoffEAM(RealType *rcut) {
581 <    return OpenMD::EAM::Instance()->setCutoffEAM(rcut);
582 <  }
583 <  void fortranCalcForce(int *atid1, int *atid2, RealType *d, RealType *rij,
584 <                        RealType *r2, RealType *sw, RealType *vpair,
585 <                        RealType *pot, RealType *f1, RealType *rho1,
586 <                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
587 <                        RealType *fshift1, RealType *fshift2){
588 <    
589 <    return OpenMD::EAM::Instance()->do_eam_pair(atid1, atid2, d, rij,
590 <                                                r2, sw, vpair,
591 <                                                pot, f1, rho1,
592 <                                                rho2, dfrho1, dfrho2,
593 <                                                fshift1,  fshift2);
594 <  }
595 <  RealType fortranGetEAMcut(int* atid) {
596 <    return OpenMD::EAM::Instance()->getEAMcut(atid);
597 <  }
421 >    RealType cut = 0.0;
422  
423 +    map<AtomType*, EAMAtomData>::iterator it;
424 +
425 +    it = EAMMap.find(atypes.first);
426 +    if (it != EAMMap.end()) {
427 +      EAMAtomData data1 = (*it).second;
428 +      cut = data1.rcut;
429 +    }
430 +
431 +    it = EAMMap.find(atypes.second);
432 +    if (it != EAMMap.end()) {
433 +      EAMAtomData data2 = (*it).second;
434 +      if (data2.rcut > cut)
435 +        cut = data2.rcut;
436 +    }
437 +
438 +    return cut;
439 +  }
440   }
441 +

Diff Legend

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