diff --git a/src/fix_box_relax.cpp b/src/fix_box_relax.cpp
index 19f95e5831..a50f6a0bd8 100644
--- a/src/fix_box_relax.cpp
+++ b/src/fix_box_relax.cpp
@@ -44,7 +44,7 @@ enum{ISO,ANISO,TRICLINIC};
 FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg)
 {
-  if (narg < 4) error->all("Illegal fix box/relax command");
+  if (narg < 5) error->all("Illegal fix box/relax command");
 
   scalar_flag = 1;
   extscalar = 1;
@@ -62,6 +62,8 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) :
 
   p_target[0] = p_target[1] = p_target[2] = 
     p_target[3] = p_target[4] = p_target[5] = 0.0;
+  p_flag[0] = p_flag[1] = p_flag[2] = 
+    p_flag[3] = p_flag[4] = p_flag[5] = 0;
 
   // process keywords
 
@@ -71,7 +73,7 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) :
 
   while (iarg < narg) {
     if (strcmp(arg[iarg],"iso") == 0) {
-      if (iarg+4 > narg) error->all("Illegal fix box/relax command");
+      if (iarg+2 > narg) error->all("Illegal fix box/relax command");
       pcouple = XYZ;
       p_target[0] = p_target[1] = p_target[2] = atof(arg[iarg+1]);
       p_flag[0] = p_flag[1] = p_flag[2] = 1;
@@ -79,9 +81,9 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) :
 	p_target[2] = 0.0;
 	p_flag[2] = 0;
       }
-      iarg += 4; 
+      iarg += 2; 
     } else if (strcmp(arg[iarg],"aniso") == 0) {
-      if (iarg+4 > narg) error->all("Illegal fix box/relax command");
+      if (iarg+2 > narg) error->all("Illegal fix box/relax command");
       pcouple = NONE;
       p_target[0] = p_target[1] = p_target[2] = atof(arg[iarg+1]);
       p_flag[0] = p_flag[1] = p_flag[2] = 1;
@@ -89,9 +91,9 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) :
 	p_target[2] = 0.0;
 	p_flag[2] = 0;
       }
-      iarg += 4;
+      iarg += 2;
     } else if (strcmp(arg[iarg],"tri") == 0) {
-      if (iarg+4 > narg) error->all("Illegal fix box/relax command");
+      if (iarg+2 > narg) error->all("Illegal fix box/relax command");
       pcouple = NONE;
       p_target[0] = p_target[1] = p_target[2] = atof(arg[iarg+1]);
       p_flag[0] = p_flag[1] = p_flag[2] = 1;
@@ -101,51 +103,51 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) :
 	p_target[2] = p_target[3] = p_target[4] = 0.0;
 	p_flag[2] = p_flag[3] = p_flag[4] = 0;
       }
-      iarg += 4;
+      iarg += 2;
 
     } else if (strcmp(arg[iarg],"x") == 0) {
-      if (iarg+4 > narg) error->all("Illegal fix box/relax command");
+      if (iarg+2 > narg) error->all("Illegal fix box/relax command");
       p_target[0] = atof(arg[iarg+1]);
       p_flag[0] = 1;
       deviatoric_flag = 1;
-      iarg += 4; 
+      iarg += 2; 
     } else if (strcmp(arg[iarg],"y") == 0) {
-      if (iarg+4 > narg) error->all("Illegal fix box/relax command");
+      if (iarg+2 > narg) error->all("Illegal fix box/relax command");
       p_target[1] = atof(arg[iarg+1]);
       p_flag[1] = 1;
       deviatoric_flag = 1;
-      iarg += 4; 
+      iarg += 2; 
     } else if (strcmp(arg[iarg],"z") == 0) {
-      if (iarg+4 > narg) error->all("Illegal fix box/relax command");
+      if (iarg+2 > narg) error->all("Illegal fix box/relax command");
       p_target[2] = atof(arg[iarg+1]);
       p_flag[2] = 1;
       deviatoric_flag = 1;
-      iarg += 4; 
+      iarg += 2; 
       if (dimension == 2)
 	error->all("Invalid fix box/relax command for a 2d simulation");
 
     } else if (strcmp(arg[iarg],"yz") == 0) {
-      if (iarg+4 > narg) error->all("Illegal fix box/relax command");
+      if (iarg+2 > narg) error->all("Illegal fix box/relax command");
       p_target[3] = atof(arg[iarg+1]);
       p_flag[3] = 1;
       deviatoric_flag = 1;
-      iarg += 4; 
+      iarg += 2; 
       if (dimension == 2)
 	error->all("Invalid fix box/relax command for a 2d simulation");
     } else if (strcmp(arg[iarg],"xz") == 0) {
-      if (iarg+4 > narg) error->all("Illegal fix box/relax command");
+      if (iarg+2 > narg) error->all("Illegal fix box/relax command");
       p_target[4] = atof(arg[iarg+1]);
       p_flag[4] = 1;
       deviatoric_flag = 1;
-      iarg += 4; 
+      iarg += 2; 
       if (dimension == 2)
 	error->all("Invalid fix box/relax command for a 2d simulation");
     } else if (strcmp(arg[iarg],"xy") == 0) {
-      if (iarg+4 > narg) error->all("Illegal fix box/relax command");
+      if (iarg+2 > narg) error->all("Illegal fix box/relax command");
       p_target[5] = atof(arg[iarg+1]);
       p_flag[5] = 1;
       deviatoric_flag = 1;
-      iarg += 4; 
+      iarg += 2; 
 
     } else if (strcmp(arg[iarg],"couple") == 0) {
       if (iarg+2 > narg) error->all("Illegal fix box/relax command");
diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp
index f1bd93a3a8..d755c12046 100644
--- a/src/fix_nh.cpp
+++ b/src/fix_nh.cpp
@@ -493,6 +493,7 @@ void FixNH::init()
   dthalf = 0.5 * update->dt;
   dt4 = 0.25 * update->dt;
   dt8 = 0.125 * update->dt;
+  dto = dthalf;
 
   p_freq_max = 0.0;
   if (pstat_flag) {
@@ -503,11 +504,11 @@ void FixNH::init()
       p_freq_max = MAX(p_freq_max,p_freq[4]);
       p_freq_max = MAX(p_freq_max,p_freq[5]);
     }
-    drag_factor = 1.0 - (update->dt * p_freq_max * drag);
+    pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain);
   }
 
   if (tstat_flag)
-    drag_factor = 1.0 - (update->dt * MAX(p_freq_max,t_freq) * drag);
+    tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain);
 
   // tally the number of dimensions that are barostatted
   // also compute the initial volume and reference cell  
@@ -536,6 +537,7 @@ void FixNH::init()
   if (strcmp(update->integrate_style,"respa") == 0) {
     nlevels_respa = ((Respa *) update->integrate)->nlevels;
     step_respa = ((Respa *) update->integrate)->step;
+    dto = 0.5*step_respa[0];
   }
 
   // detect if any rigid fixes exist so rigid bodies move when box is remapped
@@ -709,6 +711,7 @@ void FixNH::initial_integrate(int vflag)
 
 void FixNH::final_integrate()
 {
+
   nve_v();
 
   if (pstat_flag) nh_v_press();
@@ -779,14 +782,14 @@ void FixNH::initial_integrate_respa(int vflag, int ilevel, int iloop)
 	temperature->compute_scalar();
 	double tmp = pressure->compute_scalar();
       } else {
-	temperature->compute_vector();
+       	temperature->compute_vector();
 	pressure->compute_vector();
       }
       couple();
       pressure->addstep(update->ntimestep+1);
       if (mtk_flag) couple_ke();
     }
-
+    
     if (pstat_flag) {
       compute_press_target();
       nh_omega_dot();
@@ -795,24 +798,22 @@ void FixNH::initial_integrate_respa(int vflag, int ilevel, int iloop)
 
     nve_v();
 
-    // remap simulation box by 1/2 step
-
-    if (pstat_flag) remap();
-
   } else nve_v();
 
   // innermost level - also update x only for atoms in group
-  // if iloop = 0, then is 1st call at innermost level from rRESPA
-  // perform 2nd half of box remap
-  // redo KSpace coeffs since volume has changed
+  // if barostat, perform 1/2 step remap before and after
 
   if (ilevel == 0) {
+    if (pstat_flag) remap();
     nve_x();
-    if (iloop == 0) {
-      remap();
-      if (kspace_flag) force->kspace->setup();
-    }
+    if (pstat_flag) remap();
   }
+
+  // if barostat, redo KSpace coeffs at outermost level, 
+  // since volume has changed
+
+  if (ilevel == nlevels_respa-1 && kspace_flag && pstat_flag) force->kspace->setup();
+
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1339,7 +1340,7 @@ double FixNH::compute_vector(int n)
     if (pstyle == ISO) {
       ilen = 1;
       if (n < ilen) 
-	return 0.5*omega_dot[n]*omega_dot[n]*omega_mass[n];
+	return pdim*0.5*omega_dot[n]*omega_dot[n]*omega_mass[n];
       n -= ilen;
     } else if (pstyle == ANISO) {
       ilen = 3;
@@ -1396,6 +1397,12 @@ void FixNH::reset_dt()
   dthalf = 0.5 * update->dt;
   dt4 = 0.25 * update->dt;
   dt8 = 0.125 * update->dt;
+  dto = dthalf;
+
+  // If using respa, then remap is performed in innermost level
+  
+  if (strcmp(update->integrate_style,"respa") == 0)
+    dto = 0.5*step_respa[0];
   
   p_freq_max = 0.0;
   if (pstat_flag) {
@@ -1406,11 +1413,11 @@ void FixNH::reset_dt()
       p_freq_max = MAX(p_freq_max,p_freq[4]);
       p_freq_max = MAX(p_freq_max,p_freq[5]);
     }
-    drag_factor = 1.0 - (update->dt * p_freq_max * drag);
+    pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain);
   }
 
   if (tstat_flag)
-    drag_factor = 1.0 - (update->dt * MAX(p_freq_max,t_freq) * drag);
+    tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain);
 }
 
 /* ----------------------------------------------------------------------
@@ -1433,14 +1440,14 @@ void FixNH::nhc_temp_integrate()
       expfac = exp(-ncfac*dt8*eta_dot[ich+1]);
       eta_dot[ich] *= expfac;
       eta_dot[ich] += eta_dotdot[ich] * ncfac*dt4;
-      eta_dot[ich] *= drag_factor;
+      eta_dot[ich] *= tdrag_factor;
       eta_dot[ich] *= expfac;
     }
 
     expfac = exp(-ncfac*dt8*eta_dot[1]);
     eta_dot[0] *= expfac;
     eta_dot[0] += eta_dotdot[0] * ncfac*dt4;
-    eta_dot[0] *= drag_factor;
+    eta_dot[0] *= tdrag_factor;
     eta_dot[0] *= expfac;
     
     factor_eta = exp(-ncfac*dthalf*eta_dot[0]);
@@ -1501,14 +1508,14 @@ void FixNH::nhc_press_integrate()
       expfac = exp(-ncfac*dt8*etap_dot[ich+1]);
       etap_dot[ich] *= expfac;
       etap_dot[ich] += etap_dotdot[ich] * ncfac*dt4;
-      etap_dot[ich] *= drag_factor;
+      etap_dot[ich] *= pdrag_factor;
       etap_dot[ich] *= expfac;
     }
     
     expfac = exp(-ncfac*dt8*etap_dot[1]);
     etap_dot[0] *= expfac;
     etap_dot[0] += etap_dotdot[0] * ncfac*dt4;
-    etap_dot[0] *= drag_factor;
+    etap_dot[0] *= pdrag_factor;
     etap_dot[0] *= expfac;
     
     for (ich = 0; ich < mpchain; ich++)
@@ -1854,11 +1861,11 @@ void FixNH::nh_omega_dot()
 	(omega_mass[i] * nktv2p);
       if (deviatoric_flag) f_omega -= fdev[i]/(omega_mass[i] * nktv2p);
       omega_dot[i] += f_omega*dthalf;
-      omega_dot[i] *= drag_factor;
+      omega_dot[i] *= pdrag_factor;
     }
     omega[i] += dthalf*omega_dot[i];
     factor[i] = exp(-dthalf*omega_dot[i]*mtk_factor);
-    dilation[i] = exp(dthalf*omega_dot[i]);
+    dilation[i] = exp(dto*omega_dot[i]);
   }
   
   if (pstyle == TRICLINIC) {
@@ -1868,11 +1875,11 @@ void FixNH::nh_omega_dot()
 	if (deviatoric_flag) 
 	  f_omega -= fdev[i]/(omega_mass[i] * nktv2p);
 	omega_dot[i] += f_omega*dthalf;
-	omega_dot[i] *= drag_factor;
+	omega_dot[i] *= pdrag_factor;
       }
       omega[i] += dthalf*omega_dot[i];
       factor[i] = -dthalf*omega_dot[i];
-      dilation[i] = dthalf*omega_dot[i];
+      dilation[i] = dto*omega_dot[i];
     } 
   }
 }
diff --git a/src/fix_nh.h b/src/fix_nh.h
index 7b4966f905..aa6a6edd49 100644
--- a/src/fix_nh.h
+++ b/src/fix_nh.h
@@ -38,7 +38,7 @@ class FixNH : public Fix {
 
  protected:
   int dimension,which;
-  double dtv,dtf,dthalf,dt4,dt8;
+  double dtv,dtf,dthalf,dt4,dt8,dto;
   double boltz,nktv2p,tdof;
   double vol0,t0;
 
@@ -56,7 +56,8 @@ class FixNH : public Fix {
   double omega[6],omega_dot[6];
   double omega_mass[6];
   double p_current[6],dilation[6];
-  double drag,drag_factor;         // drag factor on particle thermostat
+  double drag,tdrag_factor;        // drag factor on particle thermostat
+  double pdrag_factor;             // drag factor on barostat
   double factor[6];                // velocity scaling due to barostat
   int kspace_flag;                 // 1 if KSpace invoked, 0 if not
   int nrigid;                      // number of rigid fixes