--- branches/development/src/nonbonded/EAM.cpp	2010/10/02 19:53:32	1502
+++ branches/development/src/nonbonded/EAM.cpp	2011/06/14 20:41:44	1582
@@ -347,90 +347,94 @@ namespace OpenMD {
     return;
   }
 
-  void EAM::calcDensity(DensityData ddat) {
+  void EAM::calcDensity(InteractionData &idat) {
     
     if (!initialized_) initialize();
     
-    EAMAtomData data1 = EAMMap[ddat.atype1];
-    EAMAtomData data2 = EAMMap[ddat.atype2];
+    EAMAtomData data1 = EAMMap[idat.atypes.first];
+    EAMAtomData data2 = EAMMap[idat.atypes.second];
 
-    if (ddat.rij < data1.rcut) 
-      ddat.rho_i_at_j = data1.rho->getValueAt(ddat.rij);
+    if ( *(idat.rij) < data1.rcut) 
+      *(idat.rho1) += data1.rho->getValueAt( *(idat.rij));
 
-    if (ddat.rij < data2.rcut) 
-      ddat.rho_j_at_i = data2.rho->getValueAt(ddat.rij);
 
+    if ( *(idat.rij) < data2.rcut) 
+      *(idat.rho2) += data2.rho->getValueAt( *(idat.rij));
+
     return;
   }
 
-  void EAM::calcFunctional(FunctionalData fdat) {
+  void EAM::calcFunctional(SelfData &sdat) {
 
     if (!initialized_) initialize();
 
-    EAMAtomData data1 = EAMMap[fdat.atype];
+    EAMAtomData data1 = EAMMap[ sdat.atype ];
         
-    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt(fdat.rho);
+    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt( *(sdat.rho) );
 
-    fdat.frho = result.first;
-    fdat.dfrhodrho = result.second;
+    *(sdat.frho) = result.first;
+    *(sdat.dfrhodrho) = result.second;
+
+    sdat.pot[METALLIC_FAMILY] += result.first;
+    *(sdat.particlePot) += result.first;
+
     return;
   }
 
  
-  void EAM::calcForce(InteractionData idat) {
+  void EAM::calcForce(InteractionData &idat) {
 
     if (!initialized_) initialize();
 
     pair<RealType, RealType> res;
     
-    if (idat.rij < eamRcut_) {
+    if ( *(idat.rij) < eamRcut_) {
 
-      EAMAtomData data1 = EAMMap[idat.atype1];
-      EAMAtomData data2 = EAMMap[idat.atype2];
+      EAMAtomData data1 = EAMMap[idat.atypes.first];
+      EAMAtomData data2 = EAMMap[idat.atypes.second];
 
       // get type-specific cutoff radii
 
       RealType rci = data1.rcut;
       RealType rcj = data2.rcut;
       
-      RealType rha, drha, rhb, drhb;
-      RealType pha, dpha, phb, dphb;
-      RealType phab, dvpdr;
+      RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0);
+      RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0);
+      RealType phab(0.0), dvpdr(0.0);
       RealType drhoidr, drhojdr, dudr;
       
-      if (idat.rij < rci) {
-        res = data1.rho->getValueAndDerivativeAt(idat.rij);
+      if ( *(idat.rij) < rci) {
+        res = data1.rho->getValueAndDerivativeAt( *(idat.rij));
         rha = res.first;
         drha = res.second;
 
-        res = MixingMap[make_pair(idat.atype1, idat.atype1)].phi->getValueAndDerivativeAt(idat.rij);
+        res = MixingMap[make_pair(idat.atypes.first, idat.atypes.first)].phi->getValueAndDerivativeAt( *(idat.rij) );
         pha = res.first;
         dpha = res.second;
       }
 
-      if (idat.rij < rcj) {
-        res = data2.rho->getValueAndDerivativeAt(idat.rij);
+      if ( *(idat.rij) < rcj) {
+        res = data2.rho->getValueAndDerivativeAt( *(idat.rij) );
         rhb = res.first;
         drhb = res.second;
 
-        res = MixingMap[make_pair(idat.atype2, idat.atype2)].phi->getValueAndDerivativeAt(idat.rij);
+        res = MixingMap[make_pair(idat.atypes.second, idat.atypes.second)].phi->getValueAndDerivativeAt( *(idat.rij) );
         phb = res.first;
         dphb = res.second;
       }
 
-      phab = 0.0;
-      dvpdr = 0.0;
-
       switch(mixMeth_) {
       case eamJohnson:
        
-        if (idat.rij < rci) {
+        if ( *(idat.rij) < rci) {
           phab = phab + 0.5 * (rhb / rha) * pha;
           dvpdr = dvpdr + 0.5*((rhb/rha)*dpha + 
                                pha*((drhb/rha) - (rhb*drha/rha/rha)));
         }
+          
+          
 
-        if (idat.rij < rcj) {
+        if ( *(idat.rij) < rcj) {
           phab = phab + 0.5 * (rha / rhb) * phb;
           dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb + 
                                  phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
@@ -439,7 +443,7 @@ namespace OpenMD {
         break;
 
       case eamDaw:
-        res = MixingMap[make_pair(idat.atype1,idat.atype2)].phi->getValueAndDerivativeAt(idat.rij);
+        res = MixingMap[idat.atypes].phi->getValueAndDerivativeAt( *(idat.rij));
         phab = res.first;
         dvpdr = res.second;
 
@@ -459,32 +463,55 @@ namespace OpenMD {
       drhoidr = drha;
       drhojdr = drhb;
 
-      dudr = drhojdr*idat.dfrho1 + drhoidr*idat.dfrho2 + dvpdr; 
+      dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr; 
 
-      idat.f1 = idat.d * dudr / idat.rij;
+      *(idat.f1) = *(idat.d) * dudr / *(idat.rij);
         
-      // particle_pot is the difference between the full potential 
-      // and the full potential without the presence of a particular
+      // particlePot is the difference between the full potential and
+      // the full potential without the presence of a particular
       // particle (atom1).
       //
-      // This reduces the density at other particle locations, so
-      // we need to recompute the density at atom2 assuming atom1
-      // didn't contribute.  This then requires recomputing the
-      // density functional for atom2 as well.
-      //
-      // Most of the particle_pot heavy lifting comes from the
-      // pair interaction, and will be handled by vpair.
-     
-      idat.fshift1 = data1.F->getValueAt( idat.rho1 - rhb );
-      idat.fshift2 = data1.F->getValueAt( idat.rho2 - rha );
+      // This reduces the density at other particle locations, so we
+      // need to recompute the density at atom2 assuming atom1 didn't
+      // contribute.  This then requires recomputing the density
+      // functional for atom2 as well.
 
-      idat.pot += phab;
+      *(idat.particlePot1) += data2.F->getValueAt( *(idat.rho2) - rha ) 
+        - *(idat.frho2);
 
-      idat.vpair += phab;
+      *(idat.particlePot2) += data1.F->getValueAt( *(idat.rho1) - rhb) 
+        - *(idat.frho1);
+
+      (*(idat.pot))[METALLIC_FAMILY] += phab;
+
+      *(idat.vpair) += phab;
     }
 
     return;
     
   }
+
+  RealType EAM::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
+    if (!initialized_) initialize();   
+
+    RealType cut = 0.0;
+
+    map<AtomType*, EAMAtomData>::iterator it;
+
+    it = EAMMap.find(atypes.first);
+    if (it != EAMMap.end()) {
+      EAMAtomData data1 = (*it).second;
+      cut = data1.rcut;
+    }
+
+    it = EAMMap.find(atypes.second);
+    if (it != EAMMap.end()) {
+      EAMAtomData data2 = (*it).second;
+      if (data2.rcut > cut)
+        cut = data2.rcut;
+    }
+
+    return cut;
+  }
 }