From 39b2a7831c02572fa7a13d663a244a8801b574aa Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Tue, 23 Jun 2009 20:11:54 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@2896
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/temper.cpp | 31 ++++++++++++++++++++++++++-----
 src/temper.h   |  1 +
 2 files changed, 27 insertions(+), 5 deletions(-)

diff --git a/src/temper.cpp b/src/temper.cpp
index c1a9a4cb43..9a5034f92a 100644
--- a/src/temper.cpp
+++ b/src/temper.cpp
@@ -176,7 +176,7 @@ void Temper::command(int narg, char **arg)
   // setup tempering runs
 
   int i,which,partner,swap,partner_set_temp,partner_world;
-  double pe,pe_partner,boltz_factor,new_temp;
+  double pe,pe_partner,ke,boltz_factor,new_temp;
   MPI_Status status;
 
   if (me_universe == 0 && universe->uscreen) 
@@ -193,16 +193,15 @@ void Temper::command(int narg, char **arg)
   timer->barrier_start(TIME_LOOP);
 
   for (int iswap = 0; iswap < nswaps; iswap++) {
-
+	
     // run for nevery timesteps
 
     update->integrate->iterate(nevery);
-
-    // compute PE, normalize if thermo PE does
+	
+    // compute PE
     // notify compute it will be called at next swap
 
     pe = pe_compute->compute_scalar();
-    if (output->thermo->normflag) pe /= atom->natoms;
     pe_compute->addstep(update->ntimestep + nevery);
 
     // which = which of 2 kinds of swaps to do (0,1)
@@ -268,6 +267,10 @@ void Temper::command(int narg, char **arg)
 
     MPI_Bcast(&swap,1,MPI_INT,0,world);
 
+    // rescale kinetic energy via velocities if move is accepted
+
+    if (swap) scale_velocities(partner_set_temp,my_set_temp);
+	
     // if my world swapped, all procs in world reset temp target of Fix
 
     if (swap) {
@@ -304,6 +307,24 @@ void Temper::command(int narg, char **arg)
   update->beginstep = update->endstep = 0;
 }
 
+/* ----------------------------------------------------------------------
+   scale kinetic energy via velocities a la Sugita
+------------------------------------------------------------------------- */
+
+void Temper::scale_velocities(int t_partner, int t_me)
+{
+  double sfactor = sqrt(set_temp[t_partner]/set_temp[t_me]);
+
+  double **v = atom->v;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    v[i][0] = v[i][0]*sfactor;
+    v[i][1] = v[i][1]*sfactor;
+    v[i][2] = v[i][2]*sfactor;
+  }
+}
+
 /* ----------------------------------------------------------------------
    proc 0 prints current tempering status
 ------------------------------------------------------------------------- */
diff --git a/src/temper.h b/src/temper.h
index bf717e13d4..7b6d92c7ff 100644
--- a/src/temper.h
+++ b/src/temper.h
@@ -43,6 +43,7 @@ class Temper : protected Pointers {
   int *world2temp;             // world2temp[i] = temp simulated by world i
   int *world2root;             // world2root[i] = root proc of world i
 
+  void scale_velocities(int, int);
   void print_status();
 };