--- branches/development/src/nonbonded/GB.cpp	2010/10/02 19:53:32	1502
+++ branches/development/src/nonbonded/GB.cpp	2011/04/08 21:25:19	1545
@@ -263,12 +263,11 @@ namespace OpenMD {
     }      
   }
    
-  void GB::calcForce(InteractionData idat) {
+  void GB::calcForce(InteractionData &idat) {
 
     if (!initialized_) initialize();
     
-    pair<AtomType*, AtomType*> key = make_pair(idat.atype1, idat.atype2);
-    GBInteractionData mixer = MixingMap[key];
+    GBInteractionData mixer = MixingMap[idat.atypes];
 
     RealType sigma0 = mixer.sigma0;
     RealType dw     = mixer.dw;
@@ -285,8 +284,8 @@ namespace OpenMD {
 
     RealType a, b, g;
 
-    bool i_is_LJ = idat.atype1->isLennardJones();
-    bool j_is_LJ = idat.atype2->isLennardJones();
+    bool i_is_LJ = idat.atypes.first->isLennardJones();
+    bool j_is_LJ = idat.atypes.second->isLennardJones();
 
     if (i_is_LJ) {
       a = 0.0;
@@ -357,14 +356,41 @@ namespace OpenMD {
     Vector3d rxu2 = cross(idat.d, ul2);
     Vector3d uxu = cross(ul1, ul2);
     
-    idat.pot += U*idat.sw;
+    idat.pot[0] += U*idat.sw;
     idat.f1 += dUdr * rhat + dUda * ul1 + dUdb * ul2;    
     idat.t1 += dUda * rxu1 - dUdg * uxu;
     idat.t2 += dUdb * rxu2 - dUdg * uxu;
-    idat.vpair += U*idat.sw;
+    idat.vpair[0] += U*idat.sw;
 
     return;
+
+  }
+
+  RealType GB::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
+    if (!initialized_) initialize();   
+
+    RealType cut = 0.0;
+
+    if (atypes.first->isGayBerne()) {
+      GayBerneParam gb1 = getGayBerneParam(atypes.first);
+      RealType d1 = gb1.GB_d;
+      RealType l1 = gb1.GB_l;
+      // sigma is actually sqrt(2)*l  for prolate ellipsoids 
+      cut = max(cut, 2.5 * sqrt(2.0) * max(d1, l1));
+    } else if (atypes.first->isLennardJones()) {
+      cut = max(cut, 2.5 * getLJSigma(atypes.first));
+    }
 
+    if (atypes.second->isGayBerne()) {
+      GayBerneParam gb2 = getGayBerneParam(atypes.second);
+      RealType d2 = gb2.GB_d;
+      RealType l2 = gb2.GB_l;
+      cut = max(cut, 2.5 * sqrt(2.0) * max(d2, l2));
+    } else if (atypes.second->isLennardJones()) {
+      cut = max(cut, 2.5 * getLJSigma(atypes.second));
+    }
+   
+    return cut;
   }
 }