git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4770 f3b2605a-c512-4ea7-a41b-209d697bcdaa

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;
+}

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;

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]);