--- branches/development/src/nonbonded/GB.cpp	2010/10/02 19:53:32	1502
+++ branches/development/src/nonbonded/GB.cpp	2011/01/05 14:49:05	1536
@@ -263,7 +263,7 @@ namespace OpenMD {
     }      
   }
    
-  void GB::calcForce(InteractionData idat) {
+  void GB::calcForce(InteractionData &idat) {
 
     if (!initialized_) initialize();
     
@@ -357,14 +357,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(AtomType* at1, AtomType* at2) {
+    if (!initialized_) initialize();   
+
+    RealType cut = 0.0;
+
+    if (at1->isGayBerne()) {
+      GayBerneParam gb1 = getGayBerneParam(at1);
+      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 (at1->isLennardJones()) {
+      cut = max(cut, 2.5 * getLJSigma(at1));
+    }
+
+    if (at2->isGayBerne()) {
+      GayBerneParam gb2 = getGayBerneParam(at2);
+      RealType d2 = gb2.GB_d;
+      RealType l2 = gb2.GB_l;
+      cut = max(cut, 2.5 * sqrt(2.0) * max(d2, l2));
+    } else if (at1->isLennardJones()) {
+      cut = max(cut, 2.5 * getLJSigma(at2));
+    }
+   
+    return cut;
+  }
 }