--- trunk/src/nonbonded/EAM.cpp	2013/08/17 13:03:17	1928
+++ trunk/src/nonbonded/EAM.cpp	2015/03/07 21:41:51	2071
@@ -51,8 +51,9 @@ namespace OpenMD {
 
 namespace OpenMD {
 
-  EAM::EAM() : name_("EAM"), initialized_(false), forceField_(NULL), 
-               mixMeth_(eamJohnson), eamRcut_(0.0), haveCutoffRadius_(false) {}
+  EAM::EAM() : initialized_(false), haveCutoffRadius_(false), 
+	       forceField_(NULL), eamRcut_(0.0), mixMeth_(eamJohnson), 
+	       name_("EAM") {}
   
   CubicSpline* EAM::getPhi(AtomType* atomType1, AtomType* atomType2) {   
     EAMAdapter ea1 = EAMAdapter(atomType1);
@@ -107,10 +108,9 @@ namespace OpenMD {
       zj = r <= ea2.getRcut() ? z2->getValueAt(r) : 0.0;
 
       phi = pre11_ * (zi * zj) / r;
-
+      
       phivals.push_back(phi);
     }
-      
     CubicSpline* cs = new CubicSpline();
     cs->addPoints(rvals, phivals);
     return cs;
@@ -212,7 +212,6 @@ namespace OpenMD {
     eamAtomData.F = ea.getF();
     eamAtomData.Z = ea.getZ();
     eamAtomData.rcut = ea.getRcut();
-    eamAtomData.isFluctuating = atomType->isFluctuatingCharge();
       
     // add it to the map:
     int atid = atomType->getIdent();
@@ -227,23 +226,6 @@ namespace OpenMD {
       painCave.severity = OPENMD_INFO;
       painCave.isFatal = 0;
       simError();         
-    }
-
-    if (eamAtomData.isFluctuating) {
-      // compute charge to rho scaling:
-      RealType z0 = eamAtomData.Z->getValueAt(0.0);
-      RealType dr = ea.getDr();
-      RealType rmax = max(eamAtomData.rcut, ea.getNr() * dr);
-      int nr = int(rmax/dr + 0.5);
-      RealType r;
-      RealType sum(0.0);
-
-      for (int i = 0; i < nr; i++) {
-        r = RealType(i*dr);
-        sum += r * r * eamAtomData.rho->getValueAt(r) * dr;       
-      } 
-      sum *= 4.0 * M_PI;
-      eamAtomData.qToRhoScaling = sum / z0;
     }
 
 
@@ -314,21 +296,11 @@ namespace OpenMD {
       if ( *(idat.rij) > eamRcut_) return;
     
     if ( *(idat.rij) < data1.rcut) {
-      if (data1.isFluctuating) {
-        *(idat.rho2) += (1.0 -  *(idat.flucQ1) * data1.qToRhoScaling ) * 
-          data1.rho->getValueAt( *(idat.rij) );
-      } else {
-        *(idat.rho2) += data1.rho->getValueAt( *(idat.rij));
-      }
+      *(idat.rho2) += data1.rho->getValueAt( *(idat.rij));
     }
       
     if ( *(idat.rij) < data2.rcut) {
-      if (data2.isFluctuating) {
-        *(idat.rho1) += (1.0 -  *(idat.flucQ2) * data2.qToRhoScaling ) *
-          data2.rho->getValueAt( *(idat.rij) );
-      } else {
-        *(idat.rho1) += data2.rho->getValueAt( *(idat.rij));
-      }
+      *(idat.rho1) += data2.rho->getValueAt( *(idat.rij));
     }
     
     return;  
@@ -337,10 +309,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) {
@@ -361,7 +333,6 @@ namespace OpenMD {
 
     int eamtid1 = EAMtids[idat.atid1];
     int eamtid2 = EAMtids[idat.atid2];
-    
     EAMAtomData &data1 = EAMdata[eamtid1];
     EAMAtomData &data2 = EAMdata[eamtid2];
     
@@ -369,6 +340,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);
@@ -379,46 +351,32 @@ namespace OpenMD {
       data1.rho->getValueAndDerivativeAt( *(idat.rij), rha, drha);
       CubicSpline* phi = MixingMap[eamtid1][eamtid1].phi;
       phi->getValueAndDerivativeAt( *(idat.rij), pha, dpha);
-      if (data1.isFluctuating) {
-        *(idat.dVdFQ1) -= *(idat.dfrho2) * rha * data1.qToRhoScaling;
-      }
     }
     
     if ( *(idat.rij) < rcj) {
       data2.rho->getValueAndDerivativeAt( *(idat.rij), rhb, drhb );
       CubicSpline* phi = MixingMap[eamtid2][eamtid2].phi;
       phi->getValueAndDerivativeAt( *(idat.rij), phb, dphb);
-      if (data2.isFluctuating) {
-        *(idat.dVdFQ2) -= *(idat.dfrho1) * rhb * data2.qToRhoScaling;
-      }
     }
-
     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:
-      
       if ( *(idat.rij) <  MixingMap[eamtid1][eamtid2].rcut) {
         MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij), 
                                                                   phab, dvpdr);
       }
-      
       break;
     case eamUnknown:
     default:
@@ -439,7 +397,6 @@ namespace OpenMD {
     
     *(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
@@ -457,12 +414,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) {