diff --git a/examples/cij_nvt/in.ljfcc_both b/examples/cij_nvt/in.ljfcc_both
new file mode 100644
index 0000000000..a332c4d31b
--- /dev/null
+++ b/examples/cij_nvt/in.ljfcc_both
@@ -0,0 +1,70 @@
+# Numerical difference calculation
+# of Born matrix
+
+# Note that because of cubic symmetry and central forces, we have:
+# C11, pure axial == positive mean value: 1,2,3
+# C44==C23, pure shear == positive mean value, (exactly) match in pairs: (4,12),(5,8),(6,7)
+# C14==C56, shear/axial(normal) == zero mean, (exactly) match in pairs: (9,21),(14,20),(18,19)
+# C15, shear/axial(in-plane) == zero mean: 10,11,13,15,16,17
+# the "(exactly)" above is because these symmetries seem to be sensitive to details
+# of how the finite shear deformations are applied
+
+# adjustable parameters
+
+variable    nsteps index 50     # length of run
+variable    nthermo index 10    # thermo output interval 
+variable    nlat index 3        # size of box
+variable    delta index 1.0e-6  # strain size 
+variable    temp index 0.3      # temperature
+variable    rho index 0.8442    # reduced density
+
+units		lj
+atom_style	atomic
+
+lattice		fcc ${rho}
+region		box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
+create_box	1 box
+create_atoms	1 box
+mass		1 1.0
+
+velocity	all create ${temp} 87287 loop geom
+
+pair_style	lj/cut 2.5
+pair_coeff	1 1 1.0 1.0
+
+neighbor     	0.0 bin
+neigh_modify 	every 1 delay 0 check no
+
+fix		1 all nve
+
+compute         myvirial all pressure NULL virial
+compute 	bornnum all born/matrix numdiff ${delta} myvirial
+compute 	born all born/matrix
+variable	bornerr01 equal c_born[1]-c_bornnum[1]
+variable	bornerr02 equal c_born[2]-c_bornnum[2]
+variable	bornerr03 equal c_born[3]-c_bornnum[3]
+variable	bornerr04 equal c_born[4]-c_bornnum[4]
+variable	bornerr05 equal c_born[5]-c_bornnum[5]
+variable	bornerr06 equal c_born[6]-c_bornnum[6]
+variable	bornerr07 equal c_born[7]-c_bornnum[7]
+variable	bornerr08 equal c_born[8]-c_bornnum[8]
+variable	bornerr09 equal c_born[9]-c_bornnum[9]
+variable	bornerr10 equal c_born[10]-c_bornnum[10]
+variable	bornerr11 equal c_born[11]-c_bornnum[11]
+variable	bornerr12 equal c_born[12]-c_bornnum[12]
+variable	bornerr13 equal c_born[13]-c_bornnum[13]
+variable	bornerr14 equal c_born[14]-c_bornnum[14]
+variable	bornerr15 equal c_born[15]-c_bornnum[15]
+variable	bornerr16 equal c_born[16]-c_bornnum[16]
+variable	bornerr17 equal c_born[17]-c_bornnum[17]
+variable	bornerr18 equal c_born[18]-c_bornnum[18]
+variable	bornerr19 equal c_born[19]-c_bornnum[19]
+variable	bornerr20 equal c_born[20]-c_bornnum[20]
+variable	bornerr21 equal c_born[21]-c_bornnum[21]
+thermo		${nthermo}
+thermo_style 	custom step temp pe press v_bornerr01 v_bornerr02 v_bornerr03 v_bornerr04 v_bornerr05 v_bornerr06 v_bornerr07 v_bornerr08 &
+v_bornerr09 v_bornerr10 v_bornerr11 v_bornerr12 v_bornerr13 v_bornerr14 v_bornerr15 v_bornerr16 v_bornerr17 v_bornerr18 v_bornerr19 v_bornerr20 v_bornerr21
+
+run		${nsteps}
+
+
diff --git a/src/EXTRA-COMPUTE/compute_born_matrix.cpp b/src/EXTRA-COMPUTE/compute_born_matrix.cpp
index e943b1c3e8..10dfaeae35 100644
--- a/src/EXTRA-COMPUTE/compute_born_matrix.cpp
+++ b/src/EXTRA-COMPUTE/compute_born_matrix.cpp
@@ -462,8 +462,7 @@ void ComputeBornMatrix::compute_pairs()
 	// Add contribution to Born tensor
 
 	pair->born_matrix(i, j, itype, jtype, rsq, factor_coul, factor_lj, dupair, du2pair);
-	//	pair_pref = du2pair - dupair * rinv;
-	pair_pref = du2pair;
+	pair_pref = du2pair - dupair * rinv;
 
 	// See albemunu in compute_born_matrix.h for indices order.
 	  
@@ -544,7 +543,7 @@ void ComputeBornMatrix::compute_numdiff()
   
   update_virial();
 
-  //  virial_addon();
+  virial_addon();
 
   // restore original forces for owned and ghost atoms
 
@@ -635,10 +634,10 @@ void ComputeBornMatrix::virial_addon()
     int ijvgt = idir; // this is it.
     double addon;
 
-    // extract the two indicies composing the voigt reprensentation
+    // extract the two indices composing the voigt reprensentation
 
-    id  = voigt3VtoM[idir][0];
-    jd  = voigt3VtoM[idir][1];
+    id  = voigt3VtoM[ijvgt][0];
+    jd  = voigt3VtoM[ijvgt][1];
 
     int SHEAR = 0;
     if( id != jd) SHEAR = 1;
@@ -646,12 +645,11 @@ void ComputeBornMatrix::virial_addon()
     for (int knvgt=ijvgt; knvgt < NDIR_VIRIAL; knvgt++) {
       kd = voigt3VtoM[knvgt][0];
       nd = voigt3VtoM[knvgt][1];
-      addon =  kronecker[id][nd]*sigv[virialMtoV[jd][kd]] + 
+      addon = kronecker[id][nd]*sigv[virialMtoV[jd][kd]] + 
 	kronecker[id][kd]*sigv[virialMtoV[jd][nd]]; 
       if(SHEAR)
 	addon += kronecker[jd][nd]*sigv[virialMtoV[id][kd]] + 
 	  kronecker[jd][kd]*sigv[virialMtoV[id][nd]];
-      //      values_global[voigt6MtoV[ijvgt][knvgt]] += addon;
       values_global[revalbe[ijvgt][knvgt]] += addon;
     }
   }