From 9a0acc958596575d3b53d0ff7771753800071b11 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 10 Sep 2010 17:28:15 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4770
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/fix_temp_berendsen.cpp | 17 ++++++++++++++++-
 src/fix_temp_berendsen.h   |  4 +++-
 src/fix_temp_rescale.cpp   |  4 ++--
 3 files changed, 21 insertions(+), 4 deletions(-)

diff --git a/src/fix_temp_berendsen.cpp b/src/fix_temp_berendsen.cpp
index adf260d3be..a74de7d36d 100644
--- a/src/fix_temp_berendsen.cpp
+++ b/src/fix_temp_berendsen.cpp
@@ -38,6 +38,9 @@ FixTempBerendsen::FixTempBerendsen(LAMMPS *lmp, int narg, char **arg) :
   // Berendsen thermostat should be applied every step
 
   nevery = 1;
+  scalar_flag = 1;
+  global_freq = nevery;
+  extscalar = 1;
 
   t_start = atof(arg[3]);
   t_stop = atof(arg[4]);
@@ -62,6 +65,8 @@ FixTempBerendsen::FixTempBerendsen(LAMMPS *lmp, int narg, char **arg) :
   modify->add_compute(3,newarg);
   delete [] newarg;
   tflag = 1;
+
+  energy = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -80,6 +85,7 @@ int FixTempBerendsen::setmask()
 {
   int mask = 0;
   mask |= END_OF_STEP;
+  mask |= THERMO_ENERGY;
   return mask;
 }
 
@@ -106,7 +112,7 @@ void FixTempBerendsen::end_of_step()
 
   double delta = update->ntimestep - update->beginstep;
   delta /= update->endstep - update->beginstep;
-  t_target = t_start + delta * (t_stop-t_start);
+  double t_target = t_start + delta * (t_stop-t_start);
 
   // rescale velocities by lamda
   // for BIAS:
@@ -114,6 +120,8 @@ void FixTempBerendsen::end_of_step()
   //   OK to not test returned v = 0, since lamda is multiplied by v
 
   double lamda = sqrt(1.0 + update->dt/t_period*(t_target/t_current - 1.0));
+  double efactor = 0.5 * force->boltz * temperature->dof;
+  energy += t_current * (1.0-lamda*lamda) * efactor;
 
   double **v = atom->v;
   int *mask = atom->mask;
@@ -174,3 +182,10 @@ void FixTempBerendsen::reset_target(double t_new)
 {
   t_start = t_stop = t_new;
 }
+
+/* ---------------------------------------------------------------------- */
+
+double FixTempBerendsen::compute_scalar()
+{
+  return energy;
+}
diff --git a/src/fix_temp_berendsen.h b/src/fix_temp_berendsen.h
index 3e698cc7b8..7187383674 100644
--- a/src/fix_temp_berendsen.h
+++ b/src/fix_temp_berendsen.h
@@ -33,10 +33,12 @@ class FixTempBerendsen : public Fix {
   void end_of_step();
   int modify_param(int, char **);
   void reset_target(double);
+  double compute_scalar();
 
  private:
   int which;
-  double t_start,t_stop,t_target,t_period;
+  double t_start,t_stop,t_period;
+  double energy;
 
   char *id_temp;
   class Compute *temperature;
diff --git a/src/fix_temp_rescale.cpp b/src/fix_temp_rescale.cpp
index 701dd89143..1b4c8f3f04 100644
--- a/src/fix_temp_rescale.cpp
+++ b/src/fix_temp_rescale.cpp
@@ -127,8 +127,9 @@ void FixTempRescale::end_of_step()
     int *mask = atom->mask;
     int nlocal = atom->nlocal;
 
+    energy += (t_current-t_target) * efactor;
+
     if (which == NOBIAS) {
-      energy += (t_current-t_target) * efactor;
       for (int i = 0; i < nlocal; i++) {
 	if (mask[i] & groupbit) {
 	  v[i][0] *= factor;
@@ -137,7 +138,6 @@ void FixTempRescale::end_of_step()
 	}
       }
     } else {
-      energy += (t_current-t_target) * efactor;
       for (int i = 0; i < nlocal; i++) {
 	if (mask[i] & groupbit) {
 	  temperature->remove_bias(i,v[i]);