--- branches/development/src/nonbonded/LJ.cpp	2010/10/02 19:53:32	1502
+++ branches/development/src/nonbonded/LJ.cpp	2011/06/13 22:13:12	1581
@@ -240,53 +240,62 @@ namespace OpenMD {
     }    
   }
  
-  void LJ::calcForce(InteractionData idat) {
+  void LJ::calcForce(InteractionData &idat) {
     
     if (!initialized_) initialize();
     
-    RealType ros;
-    RealType rcos;
-    RealType myPot = 0.0;
-    RealType myPotC = 0.0;
-    RealType myDeriv = 0.0;
-    RealType myDerivC = 0.0;
-
-    std::pair<AtomType*, AtomType*> key = std::make_pair(idat.atype1, 
-                                                         idat.atype2);
-    LJInteractionData mixer = MixingMap[key];
-
-    RealType sigmai = mixer.sigmai;
-    RealType epsilon = mixer.epsilon;
+    map<pair<AtomType*, AtomType*>, LJInteractionData>::iterator it;
+    it = MixingMap.find( idat.atypes );
     
-
-    ros = idat.rij * sigmai;
-
-    getLJfunc(ros, myPot, myDeriv);
+    if (it != MixingMap.end())  {
+      
+      LJInteractionData mixer = (*it).second;
+      
+      RealType sigmai = mixer.sigmai;
+      RealType epsilon = mixer.epsilon;
+      
+      RealType ros;
+      RealType rcos;
+      RealType myPot = 0.0;
+      RealType myPotC = 0.0;
+      RealType myDeriv = 0.0;
+      RealType myDerivC = 0.0;
+     
+      ros = *(idat.rij) * sigmai;
+      
+      cerr << "ros = " << ros << "\n";
+      
+      getLJfunc(ros, myPot, myDeriv);
+      
+      if (shiftedPot_) {
+        rcos = *(idat.rcut) * sigmai;
+        getLJfunc(rcos, myPotC, myDerivC);
+        myDerivC = 0.0;
+      } else if (LJ::shiftedFrc_) {
+        rcos = *(idat.rcut) * sigmai;
+        getLJfunc(rcos, myPotC, myDerivC);
+        myPotC = myPotC + myDerivC * (*(idat.rij) - *(idat.rcut)) * sigmai;
+      } else {
+        myPotC = 0.0;
+        myDerivC = 0.0;        
+      }
 
-    if (shiftedPot_) {
-      rcos = idat.rcut * sigmai;
-      getLJfunc(rcos, myPotC, myDerivC);
-      myDerivC = 0.0;
-    } else if (LJ::shiftedFrc_) {
-      rcos = idat.rcut * sigmai;
-      getLJfunc(rcos, myPotC, myDerivC);
-      myPotC = myPotC + myDerivC * (idat.rij - idat.rcut) * sigmai;
-    } else {
-      myPotC = 0.0;
-      myDerivC = 0.0;
+      cerr << "myPot =  " << myPot << "\n";
+      cerr << "myPotC =  " << myPotC << "\n";
+      cerr << "myDerivC = " << myDerivC << "\n";
+      cerr << "epsilon =  " << epsilon << "\n";
+      cerr << "vdwm =  " << *(idat.vdwMult) << "\n";
+      cerr << "sw = " << *(idat.sw) << "\n";
+      RealType pot_temp = *(idat.vdwMult) * epsilon * (myPot - myPotC);
+      *(idat.vpair) += pot_temp;
+      
+      RealType dudr = *(idat.sw) * *(idat.vdwMult) * epsilon * (myDeriv - 
+                                                                myDerivC)*sigmai;
+      
+      (idat.pot)[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
+      *(idat.f1) = *(idat.d) * dudr / *(idat.rij);
     }
-
-    RealType pot_temp = idat.vdwMult * epsilon * (myPot - myPotC);
-    idat.vpair += pot_temp;
-
-    RealType dudr = idat.sw * idat.vdwMult * epsilon * (myDeriv - 
-                                                        myDerivC)*sigmai;
-    
-    idat.pot += idat.sw * pot_temp;
-    idat.f1 = idat.d * dudr / idat.rij;
-
     return;
-
   }
   
   void LJ::getLJfunc(RealType r, RealType &pot, RealType &deriv) {
@@ -304,5 +313,16 @@ namespace OpenMD {
     return;
   }
   
+  RealType LJ::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
+    if (!initialized_) initialize();   
+    map<pair<AtomType*, AtomType*>, LJInteractionData>::iterator it;
+    it = MixingMap.find(atypes);
+    if (it == MixingMap.end()) 
+      return 0.0;
+    else  {
+      LJInteractionData mixer = (*it).second;
+      return 2.5 * mixer.sigma;
+    }
+  }
 
 }