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 1554 by gezelter, Sat Apr 30 02:54:02 2011 UTC

# Line 347 | Line 347 | namespace OpenMD {
347      return;
348    }
349  
350 <  void EAM::calcDensity(DensityData ddat) {
350 >  void EAM::calcDensity(InteractionData &idat) {
351      
352      if (!initialized_) initialize();
353      
354 <    EAMAtomData data1 = EAMMap[ddat.atype1];
355 <    EAMAtomData data2 = EAMMap[ddat.atype2];
354 >    EAMAtomData data1 = EAMMap[idat.atypes->first];
355 >    EAMAtomData data2 = EAMMap[idat.atypes->second];
356  
357 <    if (ddat.rij < data1.rcut)
358 <      ddat.rho_i_at_j = data1.rho->getValueAt(ddat.rij);
357 >    if ( *(idat.rij) < data1.rcut)
358 >      *(idat.rho_i_at_j) = data1.rho->getValueAt( *(idat.rij));
359  
360 <    if (ddat.rij < data2.rcut)
361 <      ddat.rho_j_at_i = data2.rho->getValueAt(ddat.rij);
360 >    if ( *(idat.rij) < data2.rcut)
361 >      *(idat.rho_j_at_i) = data2.rho->getValueAt( *(idat.rij));
362  
363      return;
364    }
365  
366 <  void EAM::calcFunctional(FunctionalData fdat) {
366 >  void EAM::calcFunctional(SelfData &sdat) {
367  
368      if (!initialized_) initialize();
369  
370 <    EAMAtomData data1 = EAMMap[fdat.atype];
370 >    EAMAtomData data1 = EAMMap[ sdat.atype ];
371          
372 <    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt(fdat.rho);
372 >    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt( *(sdat.rho) );
373  
374 <    fdat.frho = result.first;
375 <    fdat.dfrhodrho = result.second;
374 >    *(sdat.frho) = result.first;
375 >    *(sdat.dfrhodrho) = result.second;
376      return;
377    }
378  
379  
380 <  void EAM::calcForce(InteractionData idat) {
380 >  void EAM::calcForce(InteractionData &idat) {
381  
382      if (!initialized_) initialize();
383  
384      pair<RealType, RealType> res;
385      
386 <    if (idat.rij < eamRcut_) {
386 >    if ( *(idat.rij) < eamRcut_) {
387  
388 <      EAMAtomData data1 = EAMMap[idat.atype1];
389 <      EAMAtomData data2 = EAMMap[idat.atype2];
388 >      EAMAtomData data1 = EAMMap[idat.atypes->first];
389 >      EAMAtomData data2 = EAMMap[idat.atypes->second];
390  
391        // get type-specific cutoff radii
392  
# Line 398 | Line 398 | namespace OpenMD {
398        RealType phab, dvpdr;
399        RealType drhoidr, drhojdr, dudr;
400        
401 <      if (idat.rij < rci) {
402 <        res = data1.rho->getValueAndDerivativeAt(idat.rij);
401 >      if ( *(idat.rij) < rci) {
402 >        res = data1.rho->getValueAndDerivativeAt( *(idat.rij));
403          rha = res.first;
404          drha = res.second;
405  
406 <        res = MixingMap[make_pair(idat.atype1, idat.atype1)].phi->getValueAndDerivativeAt(idat.rij);
406 >        res = MixingMap[make_pair(idat.atypes->first, idat.atypes->first)].phi->getValueAndDerivativeAt( *(idat.rij) );
407          pha = res.first;
408          dpha = res.second;
409        }
410  
411 <      if (idat.rij < rcj) {
412 <        res = data2.rho->getValueAndDerivativeAt(idat.rij);
411 >      if ( *(idat.rij) < rcj) {
412 >        res = data2.rho->getValueAndDerivativeAt( *(idat.rij) );
413          rhb = res.first;
414          drhb = res.second;
415  
416 <        res = MixingMap[make_pair(idat.atype2, idat.atype2)].phi->getValueAndDerivativeAt(idat.rij);
416 >        res = MixingMap[make_pair(idat.atypes->second, idat.atypes->second)].phi->getValueAndDerivativeAt( *(idat.rij) );
417          phb = res.first;
418          dphb = res.second;
419        }
# Line 424 | Line 424 | namespace OpenMD {
424        switch(mixMeth_) {
425        case eamJohnson:
426        
427 <        if (idat.rij < rci) {
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) {
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)));
# Line 439 | Line 439 | namespace OpenMD {
439          break;
440  
441        case eamDaw:
442 <        res = MixingMap[make_pair(idat.atype1,idat.atype2)].phi->getValueAndDerivativeAt(idat.rij);
442 >        res = MixingMap[*(idat.atypes)].phi->getValueAndDerivativeAt( *(idat.rij));
443          phab = res.first;
444          dvpdr = res.second;
445  
# Line 459 | Line 459 | namespace OpenMD {
459        drhoidr = drha;
460        drhojdr = drhb;
461  
462 <      dudr = drhojdr*idat.dfrho1 + drhoidr*idat.dfrho2 + dvpdr;
462 >      dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
463  
464 <      idat.f1 = idat.d * dudr / idat.rij;
464 >      *(idat.f1) = *(idat.d) * dudr / *(idat.rij);
465          
466        // particle_pot is the difference between the full potential
467        // and the full potential without the presence of a particular
# Line 475 | Line 475 | namespace OpenMD {
475        // Most of the particle_pot heavy lifting comes from the
476        // pair interaction, and will be handled by vpair.
477      
478 <      idat.fshift1 = data1.F->getValueAt( idat.rho1 - rhb );
479 <      idat.fshift2 = data1.F->getValueAt( idat.rho2 - rha );
478 >      *(idat.fshift1) = data1.F->getValueAt( *(idat.rho1) - rhb );
479 >      *(idat.fshift2) = data1.F->getValueAt( *(idat.rho2) - rha );
480  
481 <      idat.pot += phab;
481 >      idat.pot[METALLIC_FAMILY] += phab;
482  
483 <      idat.vpair += phab;
483 >      *(idat.vpair) += phab;
484      }
485  
486      return;
487      
488    }
489  
490 <  RealType EAM::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
490 >  RealType EAM::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
491      if (!initialized_) initialize();  
492  
493      RealType cut = 0.0;
494  
495      map<AtomType*, EAMAtomData>::iterator it;
496  
497 <    it = EAMMap.find(at1);
497 >    it = EAMMap.find(atypes.first);
498      if (it != EAMMap.end()) {
499        EAMAtomData data1 = (*it).second;
500        cut = data1.rcut;
501      }
502  
503 <    it = EAMMap.find(at2);
503 >    it = EAMMap.find(atypes.second);
504      if (it != EAMMap.end()) {
505        EAMAtomData data2 = (*it).second;
506        if (data2.rcut > cut)

Diff Legend

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