--- branches/development/src/nonbonded/GB.cpp	2011/06/14 20:41:44	1582
+++ branches/development/src/nonbonded/GB.cpp	2012/02/16 15:59:20	1674
@@ -36,7 +36,8 @@
  * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
  * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
  * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
- * [4]  Vardeman & Gezelter, in progress (2009).                        
+ * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
+ * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
 
 #include <stdio.h>
@@ -222,31 +223,45 @@ namespace OpenMD {
         dw2 = 1.0;
       } 
                        
-      GBInteractionData mixer;        
+      GBInteractionData mixer1, mixer2;     
       
       //  Cleaver paper uses sqrt of squares to get sigma0 for
       //  mixed interactions.
             
-      mixer.sigma0 = sqrt(d1*d1 + d2*d2);
-      mixer.xa2 = (l1*l1 - d1*d1)/(l1*l1 + d2*d2);
-      mixer.xai2 = (l2*l2 - d2*d2)/(l2*l2 + d1*d1);
-      mixer.x2 = (l1*l1 - d1*d1) * (l2*l2 - d2*d2) /
+      mixer1.sigma0 = sqrt(d1*d1 + d2*d2);
+      mixer1.xa2 = (l1*l1 - d1*d1)/(l1*l1 + d2*d2);
+      mixer1.xai2 = (l2*l2 - d2*d2)/(l2*l2 + d1*d1);
+      mixer1.x2 = (l1*l1 - d1*d1) * (l2*l2 - d2*d2) /
         ((l2*l2 + d1*d1) * (l1*l1 + d2*d2));
+
+      mixer2.sigma0 = mixer1.sigma0;
+      // xa2 and xai2 for j-i pairs are reversed from the same i-j pairing.
+      // Swapping the particles reverses the anisotropy parameters:
+      mixer2.xa2 = mixer1.xai2;
+      mixer2.xai2 = mixer1.xa2;
+      mixer2.x2 = mixer1.x2;
  
       // assumed LB mixing rules for now:
 
-      mixer.dw = 0.5 * (dw1 + dw2);
-      mixer.eps0 = sqrt(e1 * e2);
+      mixer1.dw = 0.5 * (dw1 + dw2);
+      mixer1.eps0 = sqrt(e1 * e2);
+
+      mixer2.dw = mixer1.dw;
+      mixer2.eps0 = mixer1.eps0;
       
       RealType er = sqrt(er1 * er2);
-      RealType ermu = pow(er,(1.0 / mu_));
+      RealType ermu = pow(er, (RealType(1.0) / mu_));
       RealType xp = (1.0 - ermu) / (1.0 + ermu);
       RealType ap2 = 1.0 / (1.0 + ermu);
       
-      mixer.xp2 = xp * xp;
-      mixer.xpap2 = xp * ap2;
-      mixer.xpapi2 = xp / ap2;
+      mixer1.xp2 = xp * xp;
+      mixer1.xpap2 = xp * ap2;
+      mixer1.xpapi2 = xp / ap2;
 
+      mixer2.xp2 = mixer1.xp2;
+      mixer2.xpap2 = mixer1.xpap2;
+      mixer2.xpapi2 = mixer1.xpapi2;
+
       // only add this pairing if at least one of the atoms is a Gay-Berne atom
 
       if (atomType->isGayBerne() || atype2->isGayBerne()) {
@@ -255,9 +270,9 @@ namespace OpenMD {
         key1 = make_pair(atomType, atype2);
         key2 = make_pair(atype2, atomType);
         
-        MixingMap[key1] = mixer;
+        MixingMap[key1] = mixer1;
         if (key2 != key1) {
-          MixingMap[key2] = mixer;
+          MixingMap[key2] = mixer2;
         }
       }
     }      
@@ -352,17 +367,16 @@ namespace OpenMD {
       + 8.0 * eps * mu_ * (R12 - R6) * (xp2*au*bu - Hp*xp2*g) / 
       (1.0 - xp2 * g2) / e2 + 8.0 * eps * s3 * (3.0 * R7 - 6.0 * R13) * 
       (x2 * au * bu - H * x2 * g) / (1.0 - x2 * g2) / (dw * s03);
-    
 
     Vector3d rhat = *(idat.d) / *(idat.rij);   
     Vector3d rxu1 = cross(*(idat.d), ul1);
     Vector3d rxu2 = cross(*(idat.d), ul2);
     Vector3d uxu = cross(ul1, ul2);
-    
+
     (*(idat.pot))[VANDERWAALS_FAMILY] += U *  *(idat.sw);
     *(idat.f1) += dUdr * rhat + dUda * ul1 + dUdb * ul2;    
     *(idat.t1) += dUda * rxu1 - dUdg * uxu;
-    *(idat.t2) += dUdb * rxu2 - dUdg * uxu;
+    *(idat.t2) += dUdb * rxu2 + dUdg * uxu;
     *(idat.vpair) += U * *(idat.sw);
 
     return;
@@ -379,18 +393,18 @@ namespace OpenMD {
       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));
+      cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d1, l1));
     } else if (atypes.first->isLennardJones()) {
-      cut = max(cut, 2.5 * getLJSigma(atypes.first));
+      cut = max(cut, RealType(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));
+      cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d2, l2));
     } else if (atypes.second->isLennardJones()) {
-      cut = max(cut, 2.5 * getLJSigma(atypes.second));
+      cut = max(cut, RealType(2.5) * getLJSigma(atypes.second));
     }
    
     return cut;