--- trunk/src/nonbonded/EAM.cpp	2013/07/01 21:09:37	1895
+++ trunk/src/nonbonded/EAM.cpp	2014/11/01 14:12:16	2033
@@ -60,9 +60,14 @@ namespace OpenMD {
     CubicSpline* z1 = ea1.getZ();
     CubicSpline* z2 = ea2.getZ();
 
+    // Thise prefactors convert the charge-charge interactions into
+    // kcal / mol all were computed assuming distances are measured in
+    // angstroms Charge-Charge, assuming charges are measured in
+    // electrons.  Matches value in Electrostatics.cpp
+    pre11_ = 332.0637778;
+
     // make the r grid:
 
-
     // we need phi out to the largest value we'll encounter in the radial space;
     
     RealType rmax = 0.0;
@@ -101,11 +106,10 @@ namespace OpenMD {
       zi = r <= ea1.getRcut() ? z1->getValueAt(r) : 0.0;
       zj = r <= ea2.getRcut() ? z2->getValueAt(r) : 0.0;
 
-      phi = 331.999296 * (zi * zj) / r;
-
+      phi = pre11_ * (zi * zj) / r;
+      
       phivals.push_back(phi);
     }
-      
     CubicSpline* cs = new CubicSpline();
     cs->addPoints(rvals, phivals);
     return cs;
@@ -207,7 +211,7 @@ namespace OpenMD {
     eamAtomData.F = ea.getF();
     eamAtomData.Z = ea.getZ();
     eamAtomData.rcut = ea.getRcut();
-
+      
     // add it to the map:
     int atid = atomType->getIdent();
     int eamtid = EAMtypes.size();
@@ -223,6 +227,7 @@ namespace OpenMD {
       simError();         
     }
 
+
     EAMtids[atid] = eamtid;
     EAMdata[eamtid] = eamAtomData;
     MixingMap[eamtid].resize(nEAM_);
@@ -237,6 +242,7 @@ namespace OpenMD {
 
       EAMInteractionData mixer;
       mixer.phi = getPhi(atomType, atype2);
+      mixer.rcut = mixer.phi->getLimits().second;
       mixer.explicitlySet = false;
 
       MixingMap[eamtid2].resize( nEAM_ );
@@ -265,6 +271,7 @@ namespace OpenMD {
 
     cs->addPoints(rVals, phiVals);
     mixer.phi = cs;
+    mixer.rcut = mixer.phi->getLimits().second;
     mixer.explicitlySet = true;
 
     int eamtid1 = EAMtids[ atype1->getIdent() ];
@@ -287,12 +294,13 @@ namespace OpenMD {
     if (haveCutoffRadius_) 
       if ( *(idat.rij) > eamRcut_) return;
     
-    if ( *(idat.rij) < data1.rcut) 
-      *(idat.rho1) += data1.rho->getValueAt( *(idat.rij));
-    
+    if ( *(idat.rij) < data1.rcut) {
+      *(idat.rho2) += data1.rho->getValueAt( *(idat.rij));
+    }
       
-    if ( *(idat.rij) < data2.rcut) 
-      *(idat.rho2) += data2.rho->getValueAt( *(idat.rij));
+    if ( *(idat.rij) < data2.rcut) {
+      *(idat.rho1) += data2.rho->getValueAt( *(idat.rij));
+    }
     
     return;  
   }
@@ -300,10 +308,10 @@ namespace OpenMD {
   void EAM::calcFunctional(SelfData &sdat) {
     
     if (!initialized_) initialize();
-
     EAMAtomData &data1 = EAMdata[ EAMtids[sdat.atid] ];
-        
-    data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho), *(sdat.dfrhodrho) );
+            
+    data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho),
+                                      *(sdat.dfrhodrho) );
 
     (*(sdat.pot))[METALLIC_FAMILY] += *(sdat.frho);
     if (sdat.doParticlePot) {
@@ -324,7 +332,6 @@ namespace OpenMD {
 
     int eamtid1 = EAMtids[idat.atid1];
     int eamtid2 = EAMtids[idat.atid2];
-    
     EAMAtomData &data1 = EAMdata[eamtid1];
     EAMAtomData &data2 = EAMdata[eamtid2];
     
@@ -332,6 +339,7 @@ namespace OpenMD {
     
     RealType rci = data1.rcut;
     RealType rcj = data2.rcut;
+
     
     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);
@@ -349,29 +357,25 @@ namespace OpenMD {
       CubicSpline* phi = MixingMap[eamtid2][eamtid2].phi;
       phi->getValueAndDerivativeAt( *(idat.rij), phb, dphb);
     }
-
     switch(mixMeth_) {
     case eamJohnson:
-      
       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) {
         phab = phab + 0.5 * (rha / rhb) * phb;
         dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb + 
                                phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
       }
-      
       break;
-      
     case eamDaw:
-      MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij), phab, dvpdr);
-      
+      if ( *(idat.rij) <  MixingMap[eamtid1][eamtid2].rcut) {
+        MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij), 
+                                                                  phab, dvpdr);
+      }
       break;
     case eamUnknown:
     default:
@@ -391,7 +395,7 @@ namespace OpenMD {
     dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr; 
     
     *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
-        
+
     if (idat.doParticlePot) {
       // particlePot is the difference between the full potential and
       // the full potential without the presence of a particular
@@ -409,12 +413,9 @@ namespace OpenMD {
         - *(idat.frho1);
     }
     
-    (*(idat.pot))[METALLIC_FAMILY] += phab;
-    
-    *(idat.vpair) += phab;
-  
-    return;
-    
+    (*(idat.pot))[METALLIC_FAMILY] += phab;    
+    *(idat.vpair) += phab;  
+    return;    
   }
 
   RealType EAM::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {