diff --git a/src/ASPHERE/compute_temp_asphere.cpp b/src/ASPHERE/compute_temp_asphere.cpp
index 14a737742d..c975e72a69 100755
--- a/src/ASPHERE/compute_temp_asphere.cpp
+++ b/src/ASPHERE/compute_temp_asphere.cpp
@@ -293,9 +293,9 @@ void ComputeTempAsphere::remove_bias_all()
    assume remove_bias() was previously called
 ------------------------------------------------------------------------- */
 
-void ComputeTempAsphere::restore_bias(double *v)
+void ComputeTempAsphere::restore_bias(int i, double *v)
 {
-  if (tbias) tbias->restore_bias(v);
+  if (tbias) tbias->restore_bias(i,v);
 }
 
 /* ----------------------------------------------------------------------
@@ -307,4 +307,3 @@ void ComputeTempAsphere::restore_bias_all()
 {
   if (tbias) tbias->restore_bias_all();
 }
-
diff --git a/src/ASPHERE/compute_temp_asphere.h b/src/ASPHERE/compute_temp_asphere.h
index 42441face1..95f0a4f71e 100755
--- a/src/ASPHERE/compute_temp_asphere.h
+++ b/src/ASPHERE/compute_temp_asphere.h
@@ -28,7 +28,7 @@ class ComputeTempAsphere : public Compute {
 
   void remove_bias(int, double *);
   void remove_bias_all();
-  void restore_bias(double *);
+  void restore_bias(int, double *);
   void restore_bias_all();
 
  private:
diff --git a/src/ASPHERE/fix_npt_asphere.cpp b/src/ASPHERE/fix_npt_asphere.cpp
index edb2897378..621c262911 100755
--- a/src/ASPHERE/fix_npt_asphere.cpp
+++ b/src/ASPHERE/fix_npt_asphere.cpp
@@ -122,7 +122,7 @@ void FixNPTAsphere::initial_integrate(int vflag)
 	v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
 	v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
 	v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
-	temperature->restore_bias(v[i]);
+	temperature->restore_bias(i,v[i]);
       }
     }
   }
@@ -205,7 +205,7 @@ void FixNPTAsphere::final_integrate()
 	v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
 	v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
 	v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
-	temperature->restore_bias(v[i]);
+	temperature->restore_bias(i,v[i]);
 	angmom[i][0] = (angmom[i][0] + dtf * torque[i][0]) * factor_rotate;
 	angmom[i][1] = (angmom[i][1] + dtf * torque[i][1]) * factor_rotate;
 	angmom[i][2] = (angmom[i][2] + dtf * torque[i][2]) * factor_rotate;
diff --git a/src/ASPHERE/fix_nvt_asphere.cpp b/src/ASPHERE/fix_nvt_asphere.cpp
index 1c820c766d..0bc410c5aa 100755
--- a/src/ASPHERE/fix_nvt_asphere.cpp
+++ b/src/ASPHERE/fix_nvt_asphere.cpp
@@ -118,7 +118,7 @@ void FixNVTAsphere::initial_integrate(int vflag)
 	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
 	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
 	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
-	temperature->restore_bias(v[i]);
+	temperature->restore_bias(i,v[i]);
 	x[i][0] += dtv * v[i][0];
 	x[i][1] += dtv * v[i][1];
 	x[i][2] += dtv * v[i][2];
@@ -178,7 +178,7 @@ void FixNVTAsphere::final_integrate()
 	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
 	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
 	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
-	temperature->restore_bias(v[i]);
+	temperature->restore_bias(i,v[i]);
 	angmom[i][0] = (angmom[i][0] + dtf*torque[i][0]) * factor;
 	angmom[i][1] = (angmom[i][1] + dtf*torque[i][1]) * factor;
 	angmom[i][2] = (angmom[i][2] + dtf*torque[i][2]) * factor;
diff --git a/src/compute.cpp b/src/compute.cpp
index 828e07e7f2..b97b914243 100644
--- a/src/compute.cpp
+++ b/src/compute.cpp
@@ -15,9 +15,10 @@
 #include "stdlib.h"
 #include "string.h"
 #include "ctype.h"
+#include "comm.h"
 #include "compute.h"
 #include "group.h"
-#include "domain.h"
+#include "modify.h"
 #include "lattice.h"
 #include "memory.h"
 #include "error.h"
@@ -124,6 +125,11 @@ void Compute::modify_params(int narg, char **arg)
 	int n = strlen(arg[iarg+1]) + 1;
 	id_bias = new char[n];
 	strcpy(id_bias,arg[iarg+1]);
+	int icompute = modify->find_compute(id_bias);
+	if (icompute < 0) error->all("Could not find compute_modify bias ID");
+	Compute *temperature = modify->compute[icompute];
+	if (temperature->igroup != igroup && comm->me == 0)
+	  error->warning("Group for compute_modify bias != compute group");
       }
       iarg += 2;
     } else error->all("Illegal compute_modify command");
diff --git a/src/compute.h b/src/compute.h
index 9f887bd100..adfbcd88c7 100644
--- a/src/compute.h
+++ b/src/compute.h
@@ -79,7 +79,7 @@ class Compute : protected Pointers {
 
   virtual void remove_bias(int, double *) {}
   virtual void remove_bias_all() {}
-  virtual void restore_bias(double *) {}
+  virtual void restore_bias(int, double *) {}
   virtual void restore_bias_all() {}
 
   void addstep(int);
diff --git a/src/compute_temp.cpp b/src/compute_temp.cpp
index 451247f910..0132b67969 100644
--- a/src/compute_temp.cpp
+++ b/src/compute_temp.cpp
@@ -183,9 +183,9 @@ void ComputeTemp::remove_bias_all()
    assume remove_bias() was previously called
 ------------------------------------------------------------------------- */
 
-void ComputeTemp::restore_bias(double *v)
+void ComputeTemp::restore_bias(int i, double *v)
 {
-  if (tbias) tbias->restore_bias(v);
+  if (tbias) tbias->restore_bias(i,v);
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/compute_temp.h b/src/compute_temp.h
index 2405abbd66..68dd470d66 100644
--- a/src/compute_temp.h
+++ b/src/compute_temp.h
@@ -28,7 +28,7 @@ class ComputeTemp : public Compute {
 
   void remove_bias(int, double *);
   void remove_bias_all();
-  void restore_bias(double *);
+  void restore_bias(int, double *);
   void restore_bias_all();
 
  private:
diff --git a/src/compute_temp_com.cpp b/src/compute_temp_com.cpp
index 2cfbf72164..be00aa9419 100644
--- a/src/compute_temp_com.cpp
+++ b/src/compute_temp_com.cpp
@@ -188,9 +188,11 @@ void ComputeTempCOM::compute_vector()
 void ComputeTempCOM::remove_bias(int i, double *v)
 {
   if (tbias) tbias->remove_bias(i,v);
-  v[0] -= vbias[0];
-  v[1] -= vbias[1];
-  v[2] -= vbias[2];
+  if (atom->mask[i] & groupbit) {
+    v[0] -= vbias[0];
+    v[1] -= vbias[1];
+    v[2] -= vbias[2];
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -218,12 +220,14 @@ void ComputeTempCOM::remove_bias_all()
    assume remove_bias() was previously called
 ------------------------------------------------------------------------- */
 
-void ComputeTempCOM::restore_bias(double *v)
+void ComputeTempCOM::restore_bias(int i, double *v)
 {
-  v[0] += vbias[0];
-  v[1] += vbias[1];
-  v[2] += vbias[2];
-  if (tbias) tbias->restore_bias(v);
+  if (atom->mask[i] & groupbit) {
+    v[0] += vbias[0];
+    v[1] += vbias[1];
+    v[2] += vbias[2];
+  }
+  if (tbias) tbias->restore_bias(i,v);
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/compute_temp_com.h b/src/compute_temp_com.h
index 4944b2986c..c46b750023 100644
--- a/src/compute_temp_com.h
+++ b/src/compute_temp_com.h
@@ -28,7 +28,7 @@ class ComputeTempCOM : public Compute {
 
   void remove_bias(int, double *);
   void remove_bias_all();
-  void restore_bias(double *);
+  void restore_bias(int, double *);
   void restore_bias_all();
 
  private:
diff --git a/src/compute_temp_deform.cpp b/src/compute_temp_deform.cpp
index 40ba16136e..1ae5fb29c3 100644
--- a/src/compute_temp_deform.cpp
+++ b/src/compute_temp_deform.cpp
@@ -222,18 +222,20 @@ void ComputeTempDeform::remove_bias(int i, double *v)
 {
   if (tbias) tbias->remove_bias(i,v);
 
-  double lamda[3];
-  double *h_rate = domain->h_rate;
-  double *h_ratelo = domain->h_ratelo;
-
-  domain->x2lamda(atom->x[i],lamda);
-  vbias[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + 
-    h_rate[4]*lamda[2] + h_ratelo[0];
-  vbias[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
-  vbias[2] = h_rate[2]*lamda[2] + h_ratelo[2];
-  v[0] -= vbias[0];
-  v[1] -= vbias[1];
-  v[2] -= vbias[2];
+  if (atom->mask[i] & groupbit) {
+    double lamda[3];
+    double *h_rate = domain->h_rate;
+    double *h_ratelo = domain->h_ratelo;
+    
+    domain->x2lamda(atom->x[i],lamda);
+    vbias[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + 
+      h_rate[4]*lamda[2] + h_ratelo[0];
+    vbias[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
+    vbias[2] = h_rate[2]*lamda[2] + h_ratelo[2];
+    v[0] -= vbias[0];
+    v[1] -= vbias[1];
+    v[2] -= vbias[2];
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -277,12 +279,14 @@ void ComputeTempDeform::remove_bias_all()
    assume remove_bias() was previously called
 ------------------------------------------------------------------------- */
 
-void ComputeTempDeform::restore_bias(double *v)
+void ComputeTempDeform::restore_bias(int i, double *v)
 {
-  v[0] += vbias[0];
-  v[1] += vbias[1];
-  v[2] += vbias[2];
-  if (tbias) tbias->restore_bias(v);
+  if (atom->mask[i] & groupbit) {
+    v[0] += vbias[0];
+    v[1] += vbias[1];
+    v[2] += vbias[2];
+  }
+  if (tbias) tbias->restore_bias(i,v);
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/compute_temp_deform.h b/src/compute_temp_deform.h
index f6ab938322..f60343600b 100644
--- a/src/compute_temp_deform.h
+++ b/src/compute_temp_deform.h
@@ -28,7 +28,7 @@ class ComputeTempDeform : public Compute {
 
   void remove_bias(int, double *);
   void remove_bias_all();
-  void restore_bias(double *);
+  void restore_bias(int, double *);
   void restore_bias_all();
   double memory_usage();
 
diff --git a/src/compute_temp_partial.cpp b/src/compute_temp_partial.cpp
index e9b17ed930..21199a1958 100644
--- a/src/compute_temp_partial.cpp
+++ b/src/compute_temp_partial.cpp
@@ -176,17 +176,19 @@ void ComputeTempPartial::remove_bias(int i, double *v)
 {
   if (tbias) tbias->remove_bias(i,v);
 
-  if (!xflag) {
-    vbias[0] = v[0];
-    v[0] = 0.0;
-  }
-  if (!yflag) {
-    vbias[1] = v[1];
-    v[1] = 0.0;
-  }
-  if (!zflag) {
-    vbias[2] = v[2];
-    v[2] = 0.0;
+  if (atom->mask[i] & groupbit) {
+    if (!xflag) {
+      vbias[0] = v[0];
+      v[0] = 0.0;
+    }
+    if (!yflag) {
+      vbias[1] = v[1];
+      v[1] = 0.0;
+    }
+    if (!zflag) {
+      vbias[2] = v[2];
+      v[2] = 0.0;
+    }
   }
 }
 
@@ -231,12 +233,14 @@ void ComputeTempPartial::remove_bias_all()
    assume remove_bias() was previously called
 ------------------------------------------------------------------------- */
 
-void ComputeTempPartial::restore_bias(double *v)
+void ComputeTempPartial::restore_bias(int i, double *v)
 {
-  v[0] += vbias[0];
-  v[1] += vbias[1];
-  v[2] += vbias[2];
-  if (tbias) tbias->restore_bias(v);
+  if (atom->mask[i] & groupbit) {
+    v[0] += vbias[0];
+    v[1] += vbias[1];
+    v[2] += vbias[2];
+  }
+  if (tbias) tbias->restore_bias(i,v);
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/compute_temp_partial.h b/src/compute_temp_partial.h
index a327ce6391..48f3c90592 100644
--- a/src/compute_temp_partial.h
+++ b/src/compute_temp_partial.h
@@ -28,7 +28,7 @@ class ComputeTempPartial : public Compute {
 
   void remove_bias(int, double *);
   void remove_bias_all();
-  void restore_bias(double *);
+  void restore_bias(int, double *);
   void restore_bias_all();
   double memory_usage();
 
diff --git a/src/compute_temp_ramp.cpp b/src/compute_temp_ramp.cpp
index f068f233bb..ed9ed0baf1 100644
--- a/src/compute_temp_ramp.cpp
+++ b/src/compute_temp_ramp.cpp
@@ -260,11 +260,14 @@ void ComputeTempRamp::remove_bias(int i, double *v)
 {
   if (tbias) tbias->remove_bias(i,v);
 
-  double fraction = (atom->x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo);
-  fraction = MAX(fraction,0.0);
-  fraction = MIN(fraction,1.0);
-  vbias[v_dim] = v_lo + fraction*(v_hi - v_lo);
-  v[v_dim] -= vbias[v_dim];
+  if (atom->mask[i] & groupbit) {
+    double fraction = 
+      (atom->x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo);
+    fraction = MAX(fraction,0.0);
+    fraction = MIN(fraction,1.0);
+    vbias[v_dim] = v_lo + fraction*(v_hi - v_lo);
+    v[v_dim] -= vbias[v_dim];
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -302,10 +305,10 @@ void ComputeTempRamp::remove_bias_all()
    assume remove_bias() was previously called
 ------------------------------------------------------------------------- */
 
-void ComputeTempRamp::restore_bias(double *v)
+void ComputeTempRamp::restore_bias(int i, double *v)
 {
-  v[v_dim] += vbias[v_dim];
-  if (tbias) tbias->restore_bias(v);
+  if (atom->mask[i] & groupbit) v[v_dim] += vbias[v_dim];
+  if (tbias) tbias->restore_bias(i,v);
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/compute_temp_ramp.h b/src/compute_temp_ramp.h
index 44753eeba8..47df64e4eb 100644
--- a/src/compute_temp_ramp.h
+++ b/src/compute_temp_ramp.h
@@ -28,7 +28,7 @@ class ComputeTempRamp : public Compute {
 
   void remove_bias(int, double *);
   void remove_bias_all();
-  void restore_bias(double *);
+  void restore_bias(int, double *);
   void restore_bias_all();
   double memory_usage();
 
diff --git a/src/compute_temp_region.cpp b/src/compute_temp_region.cpp
index 873d1ae2eb..3a6f46b199 100644
--- a/src/compute_temp_region.cpp
+++ b/src/compute_temp_region.cpp
@@ -174,15 +174,16 @@ void ComputeTempRegion::remove_bias(int i, double *v)
 {
   if (tbias) tbias->remove_bias(i,v);
 
-  double *x = atom->x[i];
-  if (atom->mask[i] & groupbit && 
-      domain->regions[iregion]->match(x[0],x[1],x[2]))
-    vbias[0] = vbias[1] = vbias[2] = 0.0;
-  else {
-    vbias[0] = v[0];
-    vbias[1] = v[1];
-    vbias[2] = v[2];
-    v[0] = v[1] = v[2] = 0.0;
+  if (atom->mask[i] & groupbit) {
+    double *x = atom->x[i];
+    if (domain->regions[iregion]->match(x[0],x[1],x[2]))
+      vbias[0] = vbias[1] = vbias[2] = 0.0;
+    else {
+      vbias[0] = v[0];
+      vbias[1] = v[1];
+      vbias[2] = v[2];
+      v[0] = v[1] = v[2] = 0.0;
+    }
   }
 }
 
@@ -208,8 +209,7 @@ void ComputeTempRegion::remove_bias_all()
   
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
-      if (atom->mask[i] & groupbit && 
-	  domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]))
+      if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]))
 	vbiasall[i][0] = vbiasall[i][1] = vbiasall[i][2] = 0.0;
       else {
 	vbiasall[i][0] = v[i][0];
@@ -225,12 +225,14 @@ void ComputeTempRegion::remove_bias_all()
    assume remove_bias() was previously called
 ------------------------------------------------------------------------- */
 
-void ComputeTempRegion::restore_bias(double *v)
+void ComputeTempRegion::restore_bias(int i, double *v)
 {
-  v[0] += vbias[0];
-  v[1] += vbias[1];
-  v[2] += vbias[2];
-  if (tbias) tbias->restore_bias(v);
+  if (atom->mask[i] & groupbit) {
+    v[0] += vbias[0];
+    v[1] += vbias[1];
+    v[2] += vbias[2];
+  }
+  if (tbias) tbias->restore_bias(i,v);
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/compute_temp_region.h b/src/compute_temp_region.h
index 5d434c077b..2dd8273f65 100644
--- a/src/compute_temp_region.h
+++ b/src/compute_temp_region.h
@@ -28,7 +28,7 @@ class ComputeTempRegion : public Compute {
 
   void remove_bias(int, double *);
   void remove_bias_all();
-  void restore_bias(double *);
+  void restore_bias(int, double *);
   void restore_bias_all();
   double memory_usage();
 
diff --git a/src/compute_temp_sphere.cpp b/src/compute_temp_sphere.cpp
index 4e381b11e4..8e7b7baf3d 100644
--- a/src/compute_temp_sphere.cpp
+++ b/src/compute_temp_sphere.cpp
@@ -237,9 +237,9 @@ void ComputeTempSphere::remove_bias_all()
    assume remove_bias() was previously called
 ------------------------------------------------------------------------- */
 
-void ComputeTempSphere::restore_bias(double *v)
+void ComputeTempSphere::restore_bias(int i, double *v)
 {
-  if (tbias) tbias->restore_bias(v);
+  if (tbias) tbias->restore_bias(i,v);
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/compute_temp_sphere.h b/src/compute_temp_sphere.h
index d056ee7b5d..3e7d576997 100644
--- a/src/compute_temp_sphere.h
+++ b/src/compute_temp_sphere.h
@@ -28,7 +28,7 @@ class ComputeTempSphere : public Compute {
 
   void remove_bias(int, double *);
   void remove_bias_all();
-  void restore_bias(double *);
+  void restore_bias(int, double *);
   void restore_bias_all();
 
  private:
diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp
index 897e8ab911..f598ea3ca2 100644
--- a/src/fix_langevin.cpp
+++ b/src/fix_langevin.cpp
@@ -191,7 +191,7 @@ void FixLangevin::post_force(int vflag)
 	    f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
 	  if (v[i][2] != 0.0)
 	    f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
-	  temperature->restore_bias(v[i]);
+	  temperature->restore_bias(i,v[i]);
 	}
       }
     }
@@ -235,7 +235,7 @@ void FixLangevin::post_force(int vflag)
 	    f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
 	  if (v[i][2] != 0.0)
 	    f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
-	  temperature->restore_bias(v[i]);
+	  temperature->restore_bias(i,v[i]);
 	}
       }
     }
diff --git a/src/fix_npt.cpp b/src/fix_npt.cpp
index a377fb6dae..2d455d2d06 100644
--- a/src/fix_npt.cpp
+++ b/src/fix_npt.cpp
@@ -401,7 +401,7 @@ void FixNPT::initial_integrate(int vflag)
 	v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
 	v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
 	v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
-	temperature->restore_bias(v[i]);
+	temperature->restore_bias(i,v[i]);
       }
     }
   }
@@ -464,7 +464,7 @@ void FixNPT::final_integrate()
 	v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
 	v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
 	v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
-	temperature->restore_bias(v[i]);
+	temperature->restore_bias(i,v[i]);
       }
     }
   }
@@ -608,7 +608,7 @@ void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag)
 	  v[i][0] += dtfm*f[i][0];
 	  v[i][1] += dtfm*f[i][1];
 	  v[i][2] += dtfm*f[i][2];
-	  temperature->restore_bias(v[i]);
+	  temperature->restore_bias(i,v[i]);
 	}
       }
     }
@@ -673,7 +673,7 @@ void FixNPT::final_integrate_respa(int ilevel)
 	  v[i][0] += dtfm*f[i][0];
 	  v[i][1] += dtfm*f[i][1];
 	  v[i][2] += dtfm*f[i][2];
-	  temperature->restore_bias(v[i]);
+	  temperature->restore_bias(i,v[i]);
 	}
       }
     }
diff --git a/src/fix_nvt.cpp b/src/fix_nvt.cpp
index 1ea854955e..632b765bdc 100644
--- a/src/fix_nvt.cpp
+++ b/src/fix_nvt.cpp
@@ -196,7 +196,7 @@ void FixNVT::initial_integrate(int vflag)
 	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
 	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
 	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
-	temperature->restore_bias(v[i]);
+	temperature->restore_bias(i,v[i]);
 	x[i][0] += dtv * v[i][0];
 	x[i][1] += dtv * v[i][1];
 	x[i][2] += dtv * v[i][2];
@@ -239,7 +239,7 @@ void FixNVT::final_integrate()
 	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
 	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
 	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
-	temperature->restore_bias(v[i]);
+	temperature->restore_bias(i,v[i]);
       }
     }
   }
@@ -315,7 +315,7 @@ void FixNVT::initial_integrate_respa(int vflag, int ilevel, int flag)
 	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
 	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
 	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
-	temperature->restore_bias(v[i]);
+	temperature->restore_bias(i,v[i]);
       }
     }
   }
diff --git a/src/fix_nvt_sllod.cpp b/src/fix_nvt_sllod.cpp
index 423e6468e9..7c3f303a2f 100644
--- a/src/fix_nvt_sllod.cpp
+++ b/src/fix_nvt_sllod.cpp
@@ -109,7 +109,7 @@ void FixNVTSlodd::initial_integrate(int vflag)
       v[i][0] = v[i][0]*factor + dtfm*f[i][0] - dthalf*vdelu[0];
       v[i][1] = v[i][1]*factor + dtfm*f[i][1] - dthalf*vdelu[1];
       v[i][2] = v[i][2]*factor + dtfm*f[i][2] - dthalf*vdelu[2];
-      temperature->restore_bias(v[i]);
+      temperature->restore_bias(i,v[i]);
 
       x[i][0] += dtv * v[i][0];
       x[i][1] += dtv * v[i][1];
@@ -152,7 +152,7 @@ void FixNVTSlodd::final_integrate()
       v[i][0] = v[i][0]*factor + dtfm*f[i][0] - dthalf*vdelu[0];
       v[i][1] = v[i][1]*factor + dtfm*f[i][1] - dthalf*vdelu[1];
       v[i][2] = v[i][2]*factor + dtfm*f[i][2] - dthalf*vdelu[2];
-      temperature->restore_bias(v[i]);
+      temperature->restore_bias(i,v[i]);
     }
   }
 
@@ -223,7 +223,7 @@ void FixNVTSlodd::initial_integrate_respa(int vflag, int ilevel, int flag)
       v[i][0] = v[i][0]*factor + dtfm*f[i][0] - dthalf*vdelu[0];
       v[i][1] = v[i][1]*factor + dtfm*f[i][1] - dthalf*vdelu[1];
       v[i][2] = v[i][2]*factor + dtfm*f[i][2] - dthalf*vdelu[2];
-      temperature->restore_bias(v[i]);
+      temperature->restore_bias(i,v[i]);
     }
   }
 
diff --git a/src/fix_temp_berendsen.cpp b/src/fix_temp_berendsen.cpp
index 82d0a61a11..3f58d8a325 100644
--- a/src/fix_temp_berendsen.cpp
+++ b/src/fix_temp_berendsen.cpp
@@ -135,7 +135,7 @@ void FixTempBerendsen::end_of_step()
 	v[i][0] *= lamda;
 	v[i][1] *= lamda;
 	v[i][2] *= lamda;
-	temperature->restore_bias(v[i]);
+	temperature->restore_bias(i,v[i]);
       }
     }
   }
diff --git a/src/fix_temp_rescale.cpp b/src/fix_temp_rescale.cpp
index 9c2343244b..8e9e39aed6 100644
--- a/src/fix_temp_rescale.cpp
+++ b/src/fix_temp_rescale.cpp
@@ -144,7 +144,7 @@ void FixTempRescale::end_of_step()
 	  v[i][0] *= factor;
 	  v[i][1] *= factor;
 	  v[i][2] *= factor;
-	  temperature->restore_bias(v[i]);
+	  temperature->restore_bias(i,v[i]);
 	}
       }
     }