diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp
index c457e87fda..a40c434570 100644
--- a/src/fix_langevin.cpp
+++ b/src/fix_langevin.cpp
@@ -12,7 +12,8 @@
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
-   Contributing author: Carolyn Phillips (U Mich), reservoir energy tally
+   Contributing authors: Carolyn Phillips (U Mich), reservoir energy tally
+                         Aidan Thompson (SNL) GJF formulation
 ------------------------------------------------------------------------- */
 
 #include "mpi.h"
@@ -90,9 +91,10 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
   // optional args
 
   for (int i = 1; i <= atom->ntypes; i++) ratio[i] = 1.0;
-  oflag = 0;
   ascale = 0.0;
-  tally = 0;
+  gjfflag = 0;
+  oflag = 0;
+  tallyflag = 0;
   zeroflag = 0;
 
   int iarg = 7;
@@ -102,6 +104,12 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
       if (strcmp(arg[iarg+1],"no") == 0) ascale = 0.0;
       else ascale = force->numeric(arg[iarg+1]);
       iarg += 2;
+    } else if (strcmp(arg[iarg],"gjf") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
+      if (strcmp(arg[iarg+1],"no") == 0) gjfflag = 0;
+      else if (strcmp(arg[iarg+1],"yes") == 0) gjfflag = 1;
+      else error->all(FLERR,"Illegal fix langevin command");
+      iarg += 2;
     } else if (strcmp(arg[iarg],"omega") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
       if (strcmp(arg[iarg+1],"no") == 0) oflag = 0;
@@ -118,8 +126,8 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
       iarg += 3;
     } else if (strcmp(arg[iarg],"tally") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
-      if (strcmp(arg[iarg+1],"no") == 0) tally = 0;
-      else if (strcmp(arg[iarg+1],"yes") == 0) tally = 1;
+      if (strcmp(arg[iarg+1],"no") == 0) tallyflag = 0;
+      else if (strcmp(arg[iarg+1],"yes") == 0) tallyflag = 1;
       else error->all(FLERR,"Illegal fix langevin command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"zero") == 0) {
@@ -141,8 +149,30 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
 
   energy = 0.0;
   flangevin = NULL;
+  franprev = NULL;
   tforce = NULL;
   maxatom1 = maxatom2 = 0;
+
+  // Setup atom-based array for franprev
+  // register with Atom class
+  // No need to set peratom_flag
+  // as this data is for internal use only
+
+  if (gjfflag) {
+    nvalues = 3;
+    grow_arrays(atom->nlocal);
+    atom->add_callback(0);
+
+  // initialize franprev to zero
+
+    int nlocal = atom->nlocal;
+    for (int i = 0; i < nlocal; i++) {
+      franprev[i][0] = 0.0;
+      franprev[i][1] = 0.0;
+      franprev[i][2] = 0.0;
+    }
+  }
+
 }
 
 /* ---------------------------------------------------------------------- */
@@ -157,6 +187,11 @@ FixLangevin::~FixLangevin()
   delete [] id_temp;
   memory->destroy(flangevin);
   memory->destroy(tforce);
+
+  if (gjfflag) {
+    memory->destroy(franprev);
+    atom->delete_callback(id,0);
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -166,6 +201,7 @@ int FixLangevin::setmask()
   int mask = 0;
   mask |= POST_FORCE;
   mask |= POST_FORCE_RESPA;
+  mask |= POST_INTEGRATE;
   mask |= END_OF_STEP;
   mask |= THERMO_ENERGY;
   return mask;
@@ -232,11 +268,14 @@ void FixLangevin::init()
     }
   }
 
-  if (temperature && temperature->tempbias) which = BIAS;
-  else which = NOBIAS;
+  if (temperature && temperature->tempbias) tbiasflag = BIAS;
+  else tbiasflag = NOBIAS;
 
   if (strstr(update->integrate_style,"respa"))
     nlevels_respa = ((Respa *) update->integrate)->nlevels;
+
+  if (gjfflag) gjffac = 1.0/(1.0+update->dt/2.0/t_period);
+
 }
 
 /* ---------------------------------------------------------------------- */
@@ -256,8 +295,143 @@ void FixLangevin::setup(int vflag)
 
 void FixLangevin::post_force(int vflag)
 {
-  if (tally) post_force_tally();
-  else post_force_no_tally();
+  double *rmass = atom->rmass;
+
+  // enumerate all 2^6 possibilities for template parameters
+  // this avoids testing them inside inner loop: 
+  // TSTYLEATOM, GJF, TALLY, BIAS, RMASS, ZERO
+
+#ifdef TEMPLATED_FIX_LANGEVIN
+  if (tstyle == ATOM)
+    if (gjfflag)
+      if (tallyflag)
+	if (tbiasflag == BIAS)
+	  if (rmass)
+	    if (zeroflag) post_force_templated<1,1,1,1,1,1>();
+	    else          post_force_templated<1,1,1,1,1,0>();
+	  else
+	    if (zeroflag) post_force_templated<1,1,1,1,0,1>();
+	    else          post_force_templated<1,1,1,1,0,0>();
+	else
+	  if (rmass)
+	    if (zeroflag) post_force_templated<1,1,1,0,1,1>();
+	    else          post_force_templated<1,1,1,0,1,0>();
+	  else
+	    if (zeroflag) post_force_templated<1,1,1,0,0,1>();
+	    else          post_force_templated<1,1,1,0,0,0>();
+      else
+	if (tbiasflag == BIAS)
+	  if (rmass)
+	    if (zeroflag) post_force_templated<1,1,0,1,1,1>();
+	    else          post_force_templated<1,1,0,1,1,0>();
+	  else
+	    if (zeroflag) post_force_templated<1,1,0,1,0,1>();
+	    else          post_force_templated<1,1,0,1,0,0>();
+	else
+	  if (rmass)
+	    if (zeroflag) post_force_templated<1,1,0,0,1,1>();
+	    else          post_force_templated<1,1,0,0,1,0>();
+	  else
+	    if (zeroflag) post_force_templated<1,1,0,0,0,1>();
+	    else          post_force_templated<1,1,0,0,0,0>();
+    else
+      if (tallyflag)
+	if (tbiasflag == BIAS)
+	  if (rmass)
+	    if (zeroflag) post_force_templated<1,0,1,1,1,1>();
+	    else          post_force_templated<1,0,1,1,1,0>();
+	  else
+	    if (zeroflag) post_force_templated<1,0,1,1,0,1>();
+	    else          post_force_templated<1,0,1,1,0,0>();
+	else
+	  if (rmass)
+	    if (zeroflag) post_force_templated<1,0,1,0,1,1>();
+	    else          post_force_templated<1,0,1,0,1,0>();
+	  else
+	    if (zeroflag) post_force_templated<1,0,1,0,0,1>();
+	    else          post_force_templated<1,0,1,0,0,0>();
+      else
+	if (tbiasflag == BIAS)
+	  if (rmass)
+	    if (zeroflag) post_force_templated<1,0,0,1,1,1>();
+	    else          post_force_templated<1,0,0,1,1,0>();
+	  else
+	    if (zeroflag) post_force_templated<1,0,0,1,0,1>();
+	    else          post_force_templated<1,0,0,1,0,0>();
+	else
+	  if (rmass)
+	    if (zeroflag) post_force_templated<1,0,0,0,1,1>();
+	    else          post_force_templated<1,0,0,0,1,0>();
+	  else
+	    if (zeroflag) post_force_templated<1,0,0,0,0,1>();
+	    else          post_force_templated<1,0,0,0,0,0>();
+  else
+    if (gjfflag)
+      if (tallyflag)
+	if (tbiasflag == BIAS)
+	  if (rmass)
+	    if (zeroflag) post_force_templated<0,1,1,1,1,1>();
+	    else          post_force_templated<0,1,1,1,1,0>();
+	  else
+	    if (zeroflag) post_force_templated<0,1,1,1,0,1>();
+	    else          post_force_templated<0,1,1,1,0,0>();
+	else
+	  if (rmass)
+	    if (zeroflag) post_force_templated<0,1,1,0,1,1>();
+	    else          post_force_templated<0,1,1,0,1,0>();
+	  else
+	    if (zeroflag) post_force_templated<0,1,1,0,0,1>();
+	    else          post_force_templated<0,1,1,0,0,0>();
+      else
+	if (tbiasflag == BIAS)
+	  if (rmass)
+	    if (zeroflag) post_force_templated<0,1,0,1,1,1>();
+	    else          post_force_templated<0,1,0,1,1,0>();
+	  else
+	    if (zeroflag) post_force_templated<0,1,0,1,0,1>();
+	    else          post_force_templated<0,1,0,1,0,0>();
+	else
+	  if (rmass)
+	    if (zeroflag) post_force_templated<0,1,0,0,1,1>();
+	    else          post_force_templated<0,1,0,0,1,0>();
+	  else
+	    if (zeroflag) post_force_templated<0,1,0,0,0,1>();
+	    else          post_force_templated<0,1,0,0,0,0>();
+    else
+      if (tallyflag)
+	if (tbiasflag == BIAS)
+	  if (rmass)
+	    if (zeroflag) post_force_templated<0,0,1,1,1,1>();
+	    else          post_force_templated<0,0,1,1,1,0>();
+	  else
+	    if (zeroflag) post_force_templated<0,0,1,1,0,1>();
+	    else          post_force_templated<0,0,1,1,0,0>();
+	else
+	  if (rmass)
+	    if (zeroflag) post_force_templated<0,0,1,0,1,1>();
+	    else          post_force_templated<0,0,1,0,1,0>();
+	  else
+	    if (zeroflag) post_force_templated<0,0,1,0,0,1>();
+	    else          post_force_templated<0,0,1,0,0,0>();
+      else
+	if (tbiasflag == BIAS)
+	  if (rmass)
+	    if (zeroflag) post_force_templated<0,0,0,1,1,1>();
+	    else          post_force_templated<0,0,0,1,1,0>();
+	  else
+	    if (zeroflag) post_force_templated<0,0,0,1,0,1>();
+	    else          post_force_templated<0,0,0,1,0,0>();
+	else
+	  if (rmass)
+	    if (zeroflag) post_force_templated<0,0,0,0,1,1>();
+	    else          post_force_templated<0,0,0,0,1,0>();
+	  else
+	    if (zeroflag) post_force_templated<0,0,0,0,0,1>();
+	    else          post_force_templated<0,0,0,0,0,0>();
+#else
+  post_force_untemplated(int(tstyle==ATOM), gjfflag, tallyflag, 
+			 int(tbiasflag==BIAS), int(rmass!=NULL), zeroflag);
+#endif
 }
 
 /* ---------------------------------------------------------------------- */
@@ -267,9 +441,19 @@ void FixLangevin::post_force_respa(int vflag, int ilevel, int iloop)
   if (ilevel == nlevels_respa-1) post_force(vflag);
 }
 
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   modify forces using one of the many Langevin styles
+------------------------------------------------------------------------- */
 
-void FixLangevin::post_force_no_tally()
+#ifdef TEMPLATED_FIX_LANGEVIN
+template < int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY, 
+	   int Tp_BIAS, int Tp_RMASS, int Tp_ZERO >
+void FixLangevin::post_force_templated()
+#else
+void FixLangevin::post_force_untemplated(int Tp_TSTYLEATOM, 
+					 int Tp_GJF, int Tp_TALLY, 
+					 int Tp_BIAS, int Tp_RMASS, int Tp_ZERO)
+#endif
 {
   double gamma1,gamma2;
 
@@ -280,161 +464,133 @@ void FixLangevin::post_force_no_tally()
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  double delta = update->ntimestep - update->beginstep;
-  if (delta != 0.0) delta /= update->endstep - update->beginstep;
-
-  // set current t_target and t_sqrt
-  // if variable temp, evaluate variable, wrap with clear/add
-  // reallocate tforce array if necessary
-
-  if (tstyle == CONSTANT) {
-    t_target = t_start + delta * (t_stop-t_start);
-    tsqrt = sqrt(t_target);
-  } else {
-    modify->clearstep_compute();
-    if (tstyle == EQUAL) {
-      t_target = input->variable->compute_equal(tvar);
-      if (t_target < 0.0)
-        error->one(FLERR,"Fix langevin variable returned negative temperature");
-      tsqrt = sqrt(t_target);
-    } else {
-      if (nlocal > maxatom2) {
-        maxatom2 = atom->nmax;
-        memory->destroy(tforce);
-        memory->create(tforce,maxatom2,"langevin:tforce");
-      }
-      input->variable->compute_atom(tvar,igroup,tforce,1,0);
-      for (int i = 0; i < nlocal; i++)
-        if (mask[i] & groupbit)
-            if (tforce[i] < 0.0)
-              error->one(FLERR,
-                         "Fix langevin variable returned negative temperature");
-    }
-    modify->addstep_compute(update->ntimestep + 1);
-  }
-
   // apply damping and thermostat to atoms in group
-  // for BIAS:
+
+  // for Tp_TSTYLEATOM:
+  //   use per-atom per-coord target temperature
+  // for Tp_GJF:
+  //   use Gronbech-Jensen/Farago algorithm
+  //   else use regular algorithm
+  // for Tp_TALLY:
+  //   store drag plus random forces in flangevin[nlocal][3]
+  // for Tp_BIAS:
   //   calculate temperature since some computes require temp
   //   computed on current nlocal atoms to remove bias
   //   test v = 0 since some computes mask non-participating atoms via v = 0
   //   and added force has extra term not multiplied by v = 0
-  // for ZEROFLAG:
+  // for Tp_RMASS:
+  //   use per-atom masses
+  //   else use per-type masses
+  // for Tp_ZERO:
   //   sum random force over all atoms in group
   //   subtract sum/count from each atom in group
 
-  double fran[3],fsum[3],fsumall[3];
-  fsum[0] = fsum[1] = fsum[2] = 0.0;
+  double fdrag[3],fran[3],fsum[3],fsumall[3];
   bigint count;
+  double fswap;
 
   double boltz = force->boltz;
   double dt = update->dt;
   double mvv2e = force->mvv2e;
   double ftm2v = force->ftm2v;
 
-  if (zeroflag) {
+  compute_target();
+
+  if (Tp_ZERO) {
+    fsum[0] = fsum[1] = fsum[2] = 0.0;
     count = group->count(igroup);
     if (count == 0)
       error->all(FLERR,"Cannot zero Langevin force of 0 atoms");
   }
 
-  if (rmass) {
-    if (which == NOBIAS) {
-      for (int i = 0; i < nlocal; i++) {
-        if (mask[i] & groupbit) {
-          if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
-          gamma1 = -rmass[i] / t_period / ftm2v;
-          gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
-          gamma1 *= 1.0/ratio[type[i]];
-          gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
-          fran[0] = gamma2*(random->uniform()-0.5);
-          fran[1] = gamma2*(random->uniform()-0.5);
-          fran[2] = gamma2*(random->uniform()-0.5);
-          f[i][0] += gamma1*v[i][0] + fran[0];
-          f[i][1] += gamma1*v[i][1] + fran[1];
-          f[i][2] += gamma1*v[i][2] + fran[2];
-          fsum[0] += fran[0];
-          fsum[1] += fran[1];
-          fsum[2] += fran[2];
-        }
-      }
+  // reallocate flangevin if necessary
 
-    } else if (which == BIAS) {
-      temperature->compute_scalar();
-      for (int i = 0; i < nlocal; i++) {
-        if (mask[i] & groupbit) {
-          if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
-          gamma1 = -rmass[i] / t_period / ftm2v;
-          gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
-          gamma1 *= 1.0/ratio[type[i]];
-          gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
-          temperature->remove_bias(i,v[i]);
-          fran[0] = gamma2*(random->uniform()-0.5);
-          fran[1] = gamma2*(random->uniform()-0.5);
-          fran[2] = gamma2*(random->uniform()-0.5);
-          if (v[i][0] != 0.0)
-            f[i][0] += gamma1*v[i][0] + fran[0];
-          if (v[i][1] != 0.0)
-            f[i][1] += gamma1*v[i][1] + fran[1];
-          if (v[i][2] != 0.0)
-            f[i][2] += gamma1*v[i][2] + fran[2];
-          fsum[0] += fran[0];
-          fsum[1] += fran[1];
-          fsum[2] += fran[2];
-          temperature->restore_bias(i,v[i]);
-        }
-      }
+  if (Tp_TALLY) {
+    if (atom->nlocal > maxatom1) {
+      memory->destroy(flangevin);
+      maxatom1 = atom->nmax;
+      memory->create(flangevin,maxatom1,3,"langevin:flangevin");
     }
+  }
 
-  } else {
-
-    if (which == NOBIAS) {
-      for (int i = 0; i < nlocal; i++) {
-        if (mask[i] & groupbit) {
-          if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
-          gamma1 = gfactor1[type[i]];
-          gamma2 = gfactor2[type[i]] * tsqrt;
-          fran[0] = gamma2*(random->uniform()-0.5);
-          fran[1] = gamma2*(random->uniform()-0.5);
-          fran[2] = gamma2*(random->uniform()-0.5);
-          f[i][0] += gamma1*v[i][0] + fran[0];
-          f[i][1] += gamma1*v[i][1] + fran[1];
-          f[i][2] += gamma1*v[i][2] + fran[2];
-          fsum[0] += fran[0];
-          fsum[1] += fran[1];
-          fsum[2] += fran[2];
-        }
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      if (Tp_TSTYLEATOM) tsqrt = sqrt(tforce[i]);
+      if (Tp_RMASS) {
+	gamma1 = -rmass[i] / t_period / ftm2v;
+	gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
+	gamma1 *= 1.0/ratio[type[i]];
+	gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
+      } else {
+	gamma1 = gfactor1[type[i]];
+	gamma2 = gfactor2[type[i]] * tsqrt;
       }
 
-    } else if (which == BIAS) {
-      temperature->compute_scalar();
-      for (int i = 0; i < nlocal; i++) {
-        if (mask[i] & groupbit) {
-          if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
-          gamma1 = gfactor1[type[i]];
-          gamma2 = gfactor2[type[i]] * tsqrt;
-          temperature->remove_bias(i,v[i]);
-          fran[0] = gamma2*(random->uniform()-0.5);
-          fran[1] = gamma2*(random->uniform()-0.5);
-          fran[2] = gamma2*(random->uniform()-0.5);
-          if (v[i][0] != 0.0)
-            f[i][0] += gamma1*v[i][0] + fran[0];
-          if (v[i][1] != 0.0)
-            f[i][1] += gamma1*v[i][1] + fran[1];
-          if (v[i][2] != 0.0)
-            f[i][2] += gamma1*v[i][2] + fran[2];
-          fsum[0] += fran[0];
-          fsum[1] += fran[1];
-          fsum[2] += fran[2];
-          temperature->restore_bias(i,v[i]);
-        }
+      fran[0] = gamma2*(random->uniform()-0.5);
+      fran[1] = gamma2*(random->uniform()-0.5);
+      fran[2] = gamma2*(random->uniform()-0.5);
+
+      if (Tp_BIAS) {
+	temperature->remove_bias(i,v[i]);
+	fdrag[1] = gamma1*v[i][1];
+	fdrag[2] = gamma1*v[i][2];
+	fdrag[0] = gamma1*v[i][0];
+	if (v[i][0] != 0.0) f[i][0] += fdrag[0] + fran[0];
+	else fran[0] = 0.0;
+	if (v[i][1] != 0.0) f[i][1] += fdrag[1] + fran[1];
+	else fran[1] = 0.0;
+	if (v[i][2] != 0.0) f[i][2] += fdrag[2] + fran[2];
+	else fran[2] = 0.0;
+	temperature->restore_bias(i,v[i]);
+      } else {
+	fdrag[0] = gamma1*v[i][0];
+	fdrag[1] = gamma1*v[i][1];
+	fdrag[2] = gamma1*v[i][2];
+      }
+
+      if (Tp_GJF) {
+	fswap = 0.5*(fran[0]+franprev[i][0]);
+	franprev[i][0] = fran[0];
+	fran[0] = fswap;
+	fswap = 0.5*(fran[1]+franprev[i][1]);
+	franprev[i][1] = fran[1];
+	fran[1] = fswap;
+	fswap = 0.5*(fran[2]+franprev[i][2]);
+	franprev[i][2] = fran[2];
+	fran[2] = fswap;
+
+	fdrag[0] *= gjffac;
+	fdrag[1] *= gjffac;
+	fdrag[2] *= gjffac;
+	fran[0] *= gjffac;
+	fran[1] *= gjffac;
+	fran[2] *= gjffac;
+	f[i][0] *= gjffac;
+	f[i][1] *= gjffac;
+	f[i][2] *= gjffac;
+      }
+
+      f[i][0] += fdrag[0] + fran[0];
+      f[i][1] += fdrag[1] + fran[1];
+      f[i][2] += fdrag[2] + fran[2];
+
+      if (Tp_TALLY) {
+	flangevin[i][0] = fdrag[0] + fran[0];
+	flangevin[i][1] = fdrag[1] + fran[1];
+	flangevin[i][2] = fdrag[2] + fran[2];
+      }
+
+      if (Tp_ZERO) {
+	fsum[0] += fran[0];
+	fsum[1] += fran[1];
+	fsum[2] += fran[2];
       }
     }
   }
 
   // set total force to zero
 
-  if (zeroflag) {
+  if (Tp_ZERO) {
     MPI_Allreduce(fsum,fsumall,3,MPI_DOUBLE,MPI_SUM,world);
     fsumall[0] /= count;
     fsumall[1] /= count;
@@ -454,31 +610,18 @@ void FixLangevin::post_force_no_tally()
   if (ascale) angmom_thermostat();
 }
 
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   set current t_target and t_sqrt
+------------------------------------------------------------------------- */
 
-void FixLangevin::post_force_tally()
+void FixLangevin::compute_target()
 {
-  double gamma1,gamma2;
-
-  // reallocate flangevin if necessary
-
-  if (atom->nlocal > maxatom1) {
-    memory->destroy(flangevin);
-    maxatom1 = atom->nmax;
-    memory->create(flangevin,maxatom1,3,"langevin:flangevin");
-  }
-
-  double **v = atom->v;
-  double **f = atom->f;
-  double *rmass = atom->rmass;
-  int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   double delta = update->ntimestep - update->beginstep;
   if (delta != 0.0) delta /= update->endstep - update->beginstep;
 
-  // set current t_target and t_sqrt
   // if variable temp, evaluate variable, wrap with clear/add
   // reallocate tforce array if necessary
 
@@ -507,104 +650,6 @@ void FixLangevin::post_force_tally()
     }
     modify->addstep_compute(update->ntimestep + 1);
   }
-
-  // apply damping and thermostat to appropriate atoms
-  // for BIAS:
-  //   calculate temperature since some computes require temp
-  //   computed on current nlocal atoms to remove bias
-  //   test v = 0 since some computes mask non-participating atoms via v = 0
-  //   and added force has extra term not multiplied by v = 0
-
-  double boltz = force->boltz;
-  double dt = update->dt;
-  double mvv2e = force->mvv2e;
-  double ftm2v = force->ftm2v;
-
-  if (rmass) {
-    if (which == NOBIAS) {
-      for (int i = 0; i < nlocal; i++) {
-        if (mask[i] & groupbit) {
-          if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
-          gamma1 = -rmass[i] / t_period / ftm2v;
-          gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
-          gamma1 *= 1.0/ratio[type[i]];
-          gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
-          flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
-          flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
-          flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
-          f[i][0] += flangevin[i][0];
-          f[i][1] += flangevin[i][1];
-          f[i][2] += flangevin[i][2];
-        }
-      }
-
-    } else if (which == BIAS) {
-      temperature->compute_scalar();
-      for (int i = 0; i < nlocal; i++) {
-        if (mask[i] & groupbit) {
-          if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
-          gamma1 = -rmass[i] / t_period / ftm2v;
-          gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
-          gamma1 *= 1.0/ratio[type[i]];
-          gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
-          temperature->remove_bias(i,v[i]);
-          flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
-          flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
-          flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
-          if (v[i][0] != 0.0) f[i][0] += flangevin[i][0];
-          else flangevin[i][0] = 0;
-          if (v[i][1] != 0.0) f[i][1] += flangevin[i][1];
-          else flangevin[i][1] = 0;
-          if (v[i][2] != 0.0) f[i][2] += flangevin[i][2];
-          else flangevin[i][2] = 0;
-          temperature->restore_bias(i,v[i]);
-        }
-      }
-    }
-
-  } else {
-    if (which == NOBIAS) {
-      for (int i = 0; i < nlocal; i++) {
-        if (mask[i] & groupbit) {
-          if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
-          gamma1 = gfactor1[type[i]];
-          gamma2 = gfactor2[type[i]] * tsqrt;
-          flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
-          flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
-          flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
-          f[i][0] += flangevin[i][0];
-          f[i][1] += flangevin[i][1];
-          f[i][2] += flangevin[i][2];
-        }
-      }
-
-    } else if (which == BIAS) {
-      temperature->compute_scalar();
-      for (int i = 0; i < nlocal; i++) {
-        if (mask[i] & groupbit) {
-          if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
-          gamma1 = gfactor1[type[i]];
-          gamma2 = gfactor2[type[i]] * tsqrt;
-          temperature->remove_bias(i,v[i]);
-          flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
-          flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
-          flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
-          if (v[i][0] != 0.0) f[i][0] += flangevin[i][0];
-          else flangevin[i][0] = 0.0;
-          if (v[i][1] != 0.0) f[i][1] += flangevin[i][1];
-          else flangevin[i][1] = 0.0;
-          if (v[i][2] != 0.0) f[i][2] += flangevin[i][2];
-          else flangevin[i][2] = 0.0;
-          temperature->restore_bias(i,v[i]);
-        }
-      }
-    }
-  }
-
-  // thermostat omega and angmom
-
-  if (oflag) omega_thermostat();
-  if (ascale) angmom_thermostat();
 }
 
 /* ----------------------------------------------------------------------
@@ -714,7 +759,7 @@ void FixLangevin::angmom_thermostat()
 
 void FixLangevin::end_of_step()
 {
-  if (!tally) return;
+  if (!tallyflag) return;
 
   double **v = atom->v;
   int *mask = atom->mask;
@@ -781,7 +826,7 @@ int FixLangevin::modify_param(int narg, char **arg)
 
 double FixLangevin::compute_scalar()
 {
-  if (!tally || flangevin == NULL) return 0.0;
+  if (!tallyflag || flangevin == NULL) return 0.0;
 
   // capture the very first energy transfer to thermal reservoir
 
@@ -798,6 +843,8 @@ double FixLangevin::compute_scalar()
     energy = 0.5*energy_onestep*update->dt;
   }
 
+  // convert midstep energy back to previous fullstep energy
+
   double energy_me = energy - 0.5*energy_onestep*update->dt;
 
   double energy_all;
@@ -825,7 +872,50 @@ void *FixLangevin::extract(const char *str, int &dim)
 double FixLangevin::memory_usage()
 {
   double bytes = 0.0;
-  if (tally) double bytes = atom->nmax*3 * sizeof(double);
-  if (tforce) bytes = atom->nmax * sizeof(double);
+  if (gjfflag) bytes += atom->nmax*3 * sizeof(double);
+  if (tallyflag) bytes += atom->nmax*3 * sizeof(double);
+  if (tforce) bytes += atom->nmax * sizeof(double);
   return bytes;
 }
+
+/* ----------------------------------------------------------------------
+   allocate atom-based array for franprev 
+------------------------------------------------------------------------- */
+
+void FixLangevin::grow_arrays(int nmax)
+{
+  if (!gjfflag) return;
+  memory->grow(franprev,nmax,3,"fix_langevin:franprev");
+  array = franprev;
+}
+
+/* ----------------------------------------------------------------------
+   copy values within local atom-based array
+------------------------------------------------------------------------- */
+
+void FixLangevin::copy_arrays(int i, int j, int delflag)
+{
+  for (int m = 0; m < nvalues; m++)
+    array[j][m] = array[i][m];
+}
+
+/* ----------------------------------------------------------------------
+   pack values in local atom-based array for exchange with another proc
+------------------------------------------------------------------------- */
+
+int FixLangevin::pack_exchange(int i, double *buf)
+{
+  for (int m = 0; m < nvalues; m++) buf[m] = array[i][m];
+  return nvalues;
+}
+
+/* ----------------------------------------------------------------------
+   unpack values in local atom-based array from exchange with another proc
+------------------------------------------------------------------------- */
+
+int FixLangevin::unpack_exchange(int nlocal, double *buf)
+{
+  for (int m = 0; m < nvalues; m++) array[nlocal][m] = buf[m];
+  return nvalues;
+}
+
diff --git a/src/fix_langevin.h b/src/fix_langevin.h
index 313670d44d..79fe982d70 100644
--- a/src/fix_langevin.h
+++ b/src/fix_langevin.h
@@ -40,15 +40,20 @@ class FixLangevin : public Fix {
   virtual double compute_scalar();
   double memory_usage();
   virtual void *extract(const char *, int &);
+  void grow_arrays(int);
+  void copy_arrays(int, int, int);
+  int pack_exchange(int, double *);
+  int unpack_exchange(int, double *);
 
  protected:
-  int which,tally,zeroflag,oflag;
+  int gjfflag,oflag,tallyflag,zeroflag,tbiasflag;
   double ascale;
   double t_start,t_stop,t_period,t_target;
   double *gfactor1,*gfactor2,*ratio;
   double energy,energy_onestep;
   double tsqrt;
   int tstyle,tvar;
+  double gjffac;
   char *tstr;
 
   class AtomVecEllipsoid *avec;
@@ -56,6 +61,8 @@ class FixLangevin : public Fix {
   int maxatom1,maxatom2;
   double **flangevin;
   double *tforce;
+  double **franprev, **array;
+  int nvalues;
 
   char *id_temp;
   class Compute *temperature;
@@ -63,10 +70,19 @@ class FixLangevin : public Fix {
   int nlevels_respa;
   class RanMars *random;
 
-  virtual void post_force_no_tally();
-  virtual void post_force_tally();
+  // comment next line to turn off templating
+#define TEMPLATED_FIX_LANGEVIN
+#ifdef TEMPLATED_FIX_LANGEVIN
+  template < int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY, 
+	     int Tp_BIAS, int Tp_RMASS, int Tp_ZERO > 
+  void post_force_templated();
+#else
+  void post_force_untemplated(int, int, int, 
+			      int, int, int);
+#endif
   void omega_thermostat();
   void angmom_thermostat();
+  void compute_target();
 };
 
 }