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

src/ASPHERE/compute_temp_asphere.cpp

@@ -32,13 +32,16 @@
 
 using namespace LAMMPS_NS;
 
+enum{ROTATE,ALL};
+
+#define INERTIA 0.2          // moment of inertia for ellipsoid
+
 /* ---------------------------------------------------------------------- */
 
 ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg)
 {
-  if (narg != 3 && narg != 4)
+  if (narg < 3) error->all("Illegal compute temp/asphere command");
-    error->all("Illegal compute temp/asphere command");
 
   scalar_flag = vector_flag = 1;
   size_vector = 6;
@@ -48,11 +51,24 @@ ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :
 
   tempbias = 0;
   id_bias = NULL;
-  if (narg == 4) {
+  mode = ALL;
-    tempbias = 1;
+
-    int n = strlen(arg[3]) + 1;
+  int iarg = 3;
-    id_bias = new char[n];
+  while (iarg < narg) {
-    strcpy(id_bias,arg[3]);
+    if (strcmp(arg[iarg],"bias") == 0) {
+      if (iarg+2 > narg) error->all("Illegal compute temp/asphere command");
+      tempbias = 1;
+      int n = strlen(arg[iarg+1]) + 1;
+      id_bias = new char[n];
+      strcpy(id_bias,arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"dof") == 0) {
+      if (iarg+2 > narg) error->all("Illegal compute temp/asphere command");
+      if (strcmp(arg[iarg+1],"rotate") == 0) mode = ROTATE;
+      else if (strcmp(arg[iarg+1],"all") == 0) mode = ALL;
+      else error->all("Illegal compute temp/asphere command");
+      iarg += 2;
+    } else error->all("Illegal compute temp/asphere command");
   }
 
   vector = new double[6];
@@ -76,8 +92,7 @@ ComputeTempAsphere::~ComputeTempAsphere()
 
 void ComputeTempAsphere::init()
 {
-  // check that all particles are finite-size
+  // check that all particles are finite-size, no point particles allowed
-  // no point particles allowed, spherical is OK
 
   int *ellipsoid = atom->ellipsoid;
   int *mask = atom->mask;
@@ -114,18 +129,26 @@ void ComputeTempAsphere::init()
 void ComputeTempAsphere::dof_compute()
 {
   // 6 dof for 3d, 3 dof for 2d
+  // which dof are included also depends on mode
   // assume full rotation of extended particles
   // user should correct this via compute_modify if needed
 
   double natoms = group->count(igroup);
-  int nper = 6;
+  int nper;
-  if (domain->dimension == 2) nper = 3;
+  if (domain->dimension == 3) {
+    if (mode == ALL) nper = 6;
+    else nper = 3;
+  } else {
+    if (mode == ALL) nper = 3;
+    else nper = 1;
+  }
   dof = nper*natoms;
 
   // additional adjustments to dof
 
-  if (tempbias == 1) dof -= tbias->dof_remove(-1) * natoms;
+  if (tempbias == 1) {
-  else if (tempbias == 2) {
+    if (mode == ALL) dof -= tbias->dof_remove(-1) * natoms;
+  } else if (tempbias == 2) {
     int *mask = atom->mask;
     int nlocal = atom->nlocal;
     int count = 0;
@@ -154,46 +177,73 @@ double ComputeTempAsphere::compute_scalar()
   }
 
   AtomVecEllipsoid::Bonus *bonus = avec->bonus;
-  int *ellipsoid = atom->ellipsoid;
   double **v = atom->v;
   double **angmom = atom->angmom;
   double *rmass = atom->rmass;
+  int *ellipsoid = atom->ellipsoid;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   double *shape,*quat;
   double wbody[3],inertia[3];
   double rot[3][3];
-  double t = 0.0;
   
-  // sum translationals and rotational energy for each particle
+  // sum translational and rotational energy for each particle
   // no point particles since divide by inertia
 
-  for (int i = 0; i < nlocal; i++)
+  double t = 0.0;
-    if (mask[i] & groupbit) {
 
-      shape = bonus[ellipsoid[i]].shape;
+  if (mode == ALL) {
-      quat = bonus[ellipsoid[i]].quat;
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+	t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
 
-      t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
+	// principal moments of inertia
 
-      // principal moments of inertia
+	shape = bonus[ellipsoid[i]].shape;
+	quat = bonus[ellipsoid[i]].quat;
 
-      inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
+	inertia[0] = INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
-      inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
+	inertia[1] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
-      inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
+	inertia[2] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
 
-      // wbody = angular velocity in body frame
+	// wbody = angular velocity in body frame
       
-      MathExtra::quat_to_mat(quat,rot);
+	MathExtra::quat_to_mat(quat,rot);
-      MathExtra::transpose_matvec(rot,angmom[i],wbody);
+	MathExtra::transpose_matvec(rot,angmom[i],wbody);
-      wbody[0] /= inertia[0];
+	wbody[0] /= inertia[0];
-      wbody[1] /= inertia[1];
+	wbody[1] /= inertia[1];
-      wbody[2] /= inertia[2];
+	wbody[2] /= inertia[2];
+	
+	t += inertia[0]*wbody[0]*wbody[0] +
+	  inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2];
+      }
 
-      t += inertia[0]*wbody[0]*wbody[0] +
+  } else {
-	inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2];
+    for (int i = 0; i < nlocal; i++)
-    }
+      if (mask[i] & groupbit) {
+
+	// principal moments of inertia
+
+	shape = bonus[ellipsoid[i]].shape;
+	quat = bonus[ellipsoid[i]].quat;
+
+	inertia[0] = INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
+	inertia[1] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
+	inertia[2] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
+
+	// wbody = angular velocity in body frame
+      
+	MathExtra::quat_to_mat(quat,rot);
+	MathExtra::transpose_matvec(rot,angmom[i],wbody);
+	wbody[0] /= inertia[0];
+	wbody[1] /= inertia[1];
+	wbody[2] /= inertia[2];
+	
+	t += inertia[0]*wbody[0]*wbody[0] +
+	  inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2];
+      }
+  }
 
   if (tempbias) tbias->restore_bias_all();
 
@@ -217,58 +267,93 @@ void ComputeTempAsphere::compute_vector()
   }
 
   AtomVecEllipsoid::Bonus *bonus = avec->bonus;
-  int *ellipsoid = atom->ellipsoid;
   double **v = atom->v;
   double **angmom = atom->angmom;
   double *rmass = atom->rmass;
+  int *ellipsoid = atom->ellipsoid;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   double *shape,*quat;
-  double wbody[3],inertia[3];
+  double wbody[3],inertia[3],t[6];
   double rot[3][3];
-  double massone,t[6];
+  double massone;
+
+  // sum translational and rotational energy for each particle
+  // no point particles since divide by inertia
+
   for (i = 0; i < 6; i++) t[i] = 0.0;
 
-  for (i = 0; i < nlocal; i++)
+  if (mode == ALL) {
-    if (mask[i] & groupbit) {
+    for (i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+	massone = rmass[i];
+	t[0] += massone * v[i][0]*v[i][0];
+	t[1] += massone * v[i][1]*v[i][1];
+	t[2] += massone * v[i][2]*v[i][2];
+	t[3] += massone * v[i][0]*v[i][1];
+	t[4] += massone * v[i][0]*v[i][2];
+	t[5] += massone * v[i][1]*v[i][2];
+	
+	// principal moments of inertia
 
-      shape = bonus[ellipsoid[i]].shape;
+	shape = bonus[ellipsoid[i]].shape;
-      quat = bonus[ellipsoid[i]].quat;
+	quat = bonus[ellipsoid[i]].quat;
 
-      // translational kinetic energy
+	inertia[0] = INERTIA*massone * (shape[1]*shape[1]+shape[2]*shape[2]);
+	inertia[1] = INERTIA*massone * (shape[0]*shape[0]+shape[2]*shape[2]);
+	inertia[2] = INERTIA*massone * (shape[0]*shape[0]+shape[1]*shape[1]);
+	
+	// wbody = angular velocity in body frame
+	
+	MathExtra::quat_to_mat(quat,rot);
+	MathExtra::transpose_matvec(rot,angmom[i],wbody);
+	wbody[0] /= inertia[0];
+	wbody[1] /= inertia[1];
+	wbody[2] /= inertia[2];
+	
+	// rotational kinetic energy
+	
+	t[0] += inertia[0]*wbody[0]*wbody[0];
+	t[1] += inertia[1]*wbody[1]*wbody[1];
+	t[2] += inertia[2]*wbody[2]*wbody[2];
+	t[3] += inertia[0]*wbody[0]*wbody[1];
+	t[4] += inertia[1]*wbody[0]*wbody[2];
+	t[5] += inertia[2]*wbody[1]*wbody[2];
+      }
 
-      massone = rmass[i];
+  } else {
-      t[0] += massone * v[i][0]*v[i][0];
+    for (i = 0; i < nlocal; i++)
-      t[1] += massone * v[i][1]*v[i][1];
+      if (mask[i] & groupbit) {
-      t[2] += massone * v[i][2]*v[i][2];
+	
-      t[3] += massone * v[i][0]*v[i][1];
+	// principal moments of inertia
-      t[4] += massone * v[i][0]*v[i][2];
-      t[5] += massone * v[i][1]*v[i][2];
 
-      // principal moments of inertia
+	shape = bonus[ellipsoid[i]].shape;
+	quat = bonus[ellipsoid[i]].quat;
+	massone = rmass[i];
 
-      inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
+	inertia[0] = INERTIA*massone * (shape[1]*shape[1]+shape[2]*shape[2]);
-      inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
+	inertia[1] = INERTIA*massone * (shape[0]*shape[0]+shape[2]*shape[2]);
-      inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
+	inertia[2] = INERTIA*massone * (shape[0]*shape[0]+shape[1]*shape[1]);
-
+	
-      // wbody = angular velocity in body frame
+	// wbody = angular velocity in body frame
-
+	
-      MathExtra::quat_to_mat(quat,rot);
+	MathExtra::quat_to_mat(quat,rot);
-      MathExtra::transpose_matvec(rot,angmom[i],wbody);
+	MathExtra::transpose_matvec(rot,angmom[i],wbody);
-      wbody[0] /= inertia[0];
+	wbody[0] /= inertia[0];
-      wbody[1] /= inertia[1];
+	wbody[1] /= inertia[1];
-      wbody[2] /= inertia[2];
+	wbody[2] /= inertia[2];
-
+	
-      // rotational kinetic energy
+	// rotational kinetic energy
-
+	
-      t[0] += inertia[0]*wbody[0]*wbody[0];
+	t[0] += inertia[0]*wbody[0]*wbody[0];
-      t[1] += inertia[1]*wbody[1]*wbody[1];
+	t[1] += inertia[1]*wbody[1]*wbody[1];
-      t[2] += inertia[2]*wbody[2]*wbody[2];
+	t[2] += inertia[2]*wbody[2]*wbody[2];
-      t[3] += inertia[0]*wbody[0]*wbody[1];
+	t[3] += inertia[0]*wbody[0]*wbody[1];
-      t[4] += inertia[1]*wbody[0]*wbody[2];
+	t[4] += inertia[1]*wbody[0]*wbody[2];
-      t[5] += inertia[2]*wbody[1]*wbody[2];
+	t[5] += inertia[2]*wbody[1]*wbody[2];
-    }
+      }
+  }
 
   if (tempbias) tbias->restore_bias_all();
 

src/ASPHERE/compute_temp_asphere.h

@@ -36,7 +36,7 @@ class ComputeTempAsphere : public Compute {
   void restore_bias(int, double *);
 
  private:
-  int fix_dof;
+  int fix_dof,mode;
   double tfactor;
   char *id_bias;
   class Compute *tbias;              // ptr to additional bias compute

src/ASPHERE/fix_nve_asphere.cpp

@@ -29,6 +29,8 @@
 
 using namespace LAMMPS_NS;
 
+#define INERTIA 0.2          // moment of inertia for ellipsoid
+
 /* ---------------------------------------------------------------------- */
 
 FixNVEAsphere::FixNVEAsphere(LAMMPS *lmp, int narg, char **arg) : 
@@ -103,9 +105,9 @@ void FixNVEAsphere::initial_integrate(int vflag)
       shape = bonus[ellipsoid[i]].shape;
       quat = bonus[ellipsoid[i]].quat;
 
-      inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
+      inertia[0] = INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
-      inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
+      inertia[1] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
-      inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
+      inertia[2] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
 
       // compute omega at 1/2 step from angmom at 1/2 step and current q
       // update quaternion a full step via Richardson iteration

src/USER-EFF/fix_langevin_eff.cpp

@@ -88,7 +88,9 @@ void FixLangevinEff::post_force_no_tally()
 	f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
 	f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
 	f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
-        if (abs(spin[i])==1) erforce[i] += 0.75*gamma1*ervel[i] + 0.866025404*gamma2*(random->uniform()-0.5);
+        if (abs(spin[i])==1) 
+	  erforce[i] += 0.75*gamma1*ervel[i] + 
+	    0.866025404*gamma2*(random->uniform()-0.5);
       }
     }
   } else if (which == BIAS) {
@@ -105,7 +107,8 @@ void FixLangevinEff::post_force_no_tally()
 	if (v[i][2] != 0.0)
 	  f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
         if (abs(spin[i])==1 && ervel[i] != 0.0)
-          erforce[i] += 0.75*gamma1*ervel[i] + 0.866025404*gamma2*(random->uniform()-0.5);
+          erforce[i] += 0.75*gamma1*ervel[i] + 
+	    0.866025404*gamma2*(random->uniform()-0.5);
 	temperature->restore_bias(i,v[i]);
       }
     }
@@ -158,7 +161,8 @@ void FixLangevinEff::post_force_tally()
 	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);
-        erforcelangevin[i] = 0.75*gamma1*ervel[i]+0.866025404*gamma2*(random->uniform()-0.5);
+        erforcelangevin[i] = 0.75*gamma1*ervel[i] + 
+	  0.866025404*gamma2*(random->uniform()-0.5);
 	f[i][0] += flangevin[i][0];
 	f[i][1] += flangevin[i][1];
 	f[i][2] += flangevin[i][2];
@@ -175,14 +179,16 @@ void FixLangevinEff::post_force_tally()
 	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);
-        erforcelangevin[i] = 0.75*gamma1*ervel[i]+0.866025404*gamma2*(random->uniform()-0.5);
+        erforcelangevin[i] = 0.75*gamma1*ervel[i] + 
+	  0.866025404*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;
-        if (abs(spin[i])==1 && ervel[i] != 0.0) erforce[i] += erforcelangevin[i];
+        if (abs(spin[i])==1 && ervel[i] != 0.0)
+	  erforce[i] += erforcelangevin[i];
 	temperature->restore_bias(i,v[i]);
       }
     }

src/compute_temp_sphere.cpp

@@ -26,6 +26,8 @@
 
 using namespace LAMMPS_NS;
 
+enum{ROTATE,ALL};
+
 #define INERTIA 0.4          // moment of inertia for sphere
 
 /* ---------------------------------------------------------------------- */
@@ -33,8 +35,7 @@ using namespace LAMMPS_NS;
 ComputeTempSphere::ComputeTempSphere(LAMMPS *lmp, int narg, char **arg) : 
   Compute(lmp, narg, arg)
 {
-  if (narg != 3 && narg != 4)
+  if (narg < 3) error->all("Illegal compute temp/sphere command");
-    error->all("Illegal compute temp/sphere command");
 
   scalar_flag = vector_flag = 1;
   size_vector = 6;
@@ -44,11 +45,24 @@ ComputeTempSphere::ComputeTempSphere(LAMMPS *lmp, int narg, char **arg) :
 
   tempbias = 0;
   id_bias = NULL;
-  if (narg == 4) {
+  mode = ALL;
-    tempbias = 1;
+  
-    int n = strlen(arg[3]) + 1;
+  int iarg = 3;
-    id_bias = new char[n];
+  while (iarg < narg) {
-    strcpy(id_bias,arg[3]);
+    if (strcmp(arg[iarg],"bias") == 0) {
+      if (iarg+2 > narg) error->all("Illegal compute temp/sphere command");
+      tempbias = 1;
+      int n = strlen(arg[iarg+1]) + 1;
+      id_bias = new char[n];
+      strcpy(id_bias,arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"dof") == 0) {
+      if (iarg+2 > narg) error->all("Illegal compute temp/sphere command");
+      if (strcmp(arg[iarg+1],"rotate") == 0) mode = ROTATE;
+      else if (strcmp(arg[iarg+1],"all") == 0) mode = ALL;
+      else error->all("Illegal compute temp/sphere command");
+      iarg += 2;
+    } else error->all("Illegal compute temp/sphere command");
   }
 
   vector = new double[6];
@@ -100,27 +114,34 @@ void ComputeTempSphere::dof_compute()
 
   // 6 or 3 dof for extended/point particles for 3d
   // 3 or 2 dof for extended/point particles for 2d
+  // which dof are included also depends on mode
   // assume full rotation of extended particles
   // user should correct this via compute_modify if needed
 
-  int dimension = domain->dimension;
-
   double *radius = atom->radius;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   count = 0;
-  if (dimension == 3) {
+  if (domain->dimension == 3) {
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit) {
-	if (radius[i] == 0.0) count += 3;
+	if (radius[i] == 0.0) {
-	else count += 6;
+	  if (mode == ALL) count += 3;
+	} else {
+	  if (mode == ALL) count += 6;
+	  else count += 3;
+	}
       }
   } else {
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit) {
-	if (radius[i] == 0.0) count += 2;
+	if (radius[i] == 0.0) {
-	else count += 3;
+	  if (mode == ALL) count += 2;
+	} else {
+	  if (mode == ALL) count += 3;
+	  else count += 1;
+	}
       }
   }
 
@@ -130,28 +151,38 @@ void ComputeTempSphere::dof_compute()
   // additional adjustments to dof
 
   if (tempbias == 1) {
-    double natoms = group->count(igroup);
+    if (mode == ALL) {
-    dof -= tbias->dof_remove(-1) * natoms;
+      double natoms = group->count(igroup);
+      dof -= tbias->dof_remove(-1) * natoms;
+    }
 
   } else if (tempbias == 2) {
     int *mask = atom->mask;
     int nlocal = atom->nlocal;
 
     count = 0;
-    if (dimension == 3) {
+    if (domain->dimension == 3) {
       for (int i = 0; i < nlocal; i++)
 	if (mask[i] & groupbit) {
 	  if (tbias->dof_remove(i)) {
-	    if (radius[i] == 0.0) count += 3;
+	    if (radius[i] == 0.0) {
-	    else count += 6;
+	      if (mode == ALL) count += 3;
+	    } else {
+	      if (mode == ALL) count += 6;
+	      else count += 3;
+	    }
 	  }
 	}
     } else {
       for (int i = 0; i < nlocal; i++)
 	if (mask[i] & groupbit) {
 	  if (tbias->dof_remove(i)) {
-	    if (radius[i] == 0.0) count += 2;
+	    if (radius[i] == 0.0) {
-	    else count += 3;
+	      if (mode == ALL) count += 2;
+	    } else {
+	      if (mode == ALL) count += 3;
+	      else count += 1;
+	    }
 	  }
 	}
     }
@@ -187,12 +218,19 @@ double ComputeTempSphere::compute_scalar()
 
   double t = 0.0;
 
-  for (int i = 0; i < nlocal; i++)
+  if (mode == ALL) {
-    if (mask[i] & groupbit) {
+    for (int i = 0; i < nlocal; i++)
-      t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
+      if (mask[i] & groupbit) {
-      t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
+	t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
-	    omega[i][2]*omega[i][2]) * INERTIA*radius[i]*radius[i]*rmass[i];
+	t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
-    }
+	      omega[i][2]*omega[i][2]) * INERTIA*rmass[i]*radius[i]*radius[i];
+      }
+  } else {
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit)
+	t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
+	      omega[i][2]*omega[i][2]) * INERTIA*rmass[i]*radius[i]*radius[i];
+  }
 
   if (tempbias) tbias->restore_bias_all();
 
@@ -225,25 +263,38 @@ void ComputeTempSphere::compute_vector()
   double massone,inertiaone,t[6];
   for (int i = 0; i < 6; i++) t[i] = 0.0;
 
-  for (int i = 0; i < nlocal; i++)
+  if (mode == ALL) {
-    if (mask[i] & groupbit) {
+    for (int i = 0; i < nlocal; i++)
-      massone = rmass[i];
+      if (mask[i] & groupbit) {
-      t[0] += massone * v[i][0]*v[i][0];
+	massone = rmass[i];
-      t[1] += massone * v[i][1]*v[i][1];
+	t[0] += massone * v[i][0]*v[i][0];
-      t[2] += massone * v[i][2]*v[i][2];
+	t[1] += massone * v[i][1]*v[i][1];
-      t[3] += massone * v[i][0]*v[i][1];
+	t[2] += massone * v[i][2]*v[i][2];
-      t[4] += massone * v[i][0]*v[i][2];
+	t[3] += massone * v[i][0]*v[i][1];
-      t[5] += massone * v[i][1]*v[i][2];
+	t[4] += massone * v[i][0]*v[i][2];
-      
+	t[5] += massone * v[i][1]*v[i][2];
-      inertiaone = INERTIA*radius[i]*radius[i]*rmass[i];
+	
-      t[0] += inertiaone * omega[i][0]*omega[i][0];
+	inertiaone = INERTIA*rmass[i]*radius[i]*radius[i];
-      t[1] += inertiaone * omega[i][1]*omega[i][1];
+	t[0] += inertiaone * omega[i][0]*omega[i][0];
-      t[2] += inertiaone * omega[i][2]*omega[i][2];
+	t[1] += inertiaone * omega[i][1]*omega[i][1];
-      t[3] += inertiaone * omega[i][0]*omega[i][1];
+	t[2] += inertiaone * omega[i][2]*omega[i][2];
-      t[4] += inertiaone * omega[i][0]*omega[i][2];
+	t[3] += inertiaone * omega[i][0]*omega[i][1];
-      t[5] += inertiaone * omega[i][1]*omega[i][2];
+	t[4] += inertiaone * omega[i][0]*omega[i][2];
-    }
+	t[5] += inertiaone * omega[i][1]*omega[i][2];
-
+      }
+  } else {
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+	inertiaone = INERTIA*rmass[i]*radius[i]*radius[i];
+	t[0] += inertiaone * omega[i][0]*omega[i][0];
+	t[1] += inertiaone * omega[i][1]*omega[i][1];
+	t[2] += inertiaone * omega[i][2]*omega[i][2];
+	t[3] += inertiaone * omega[i][0]*omega[i][1];
+	t[4] += inertiaone * omega[i][0]*omega[i][2];
+	t[5] += inertiaone * omega[i][1]*omega[i][2];
+      }
+  }
+    
   if (tempbias) tbias->restore_bias_all();
 
   MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);

src/compute_temp_sphere.h

@@ -36,7 +36,7 @@ class ComputeTempSphere : public Compute {
   void restore_bias(int, double *);
 
  private:
-  int fix_dof;
+  int fix_dof,mode;
   double tfactor;
   double *inertia;
   char *id_bias;

src/fix_langevin.cpp

@@ -20,7 +20,9 @@
 #include "string.h"
 #include "stdlib.h"
 #include "fix_langevin.h"
+#include "math_extra.h"
 #include "atom.h"
+#include "atom_vec_ellipsoid.h"
 #include "force.h"
 #include "update.h"
 #include "modify.h"
@@ -38,6 +40,9 @@ using namespace LAMMPS_NS;
 
 enum{NOBIAS,BIAS};
 
+#define SINERTIA 0.4          // moment of inertia for sphere
+#define EINERTIA 0.2          // moment of inertia for ellipsoid
+
 /* ---------------------------------------------------------------------- */
 
 FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
@@ -71,6 +76,7 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
   // optional args
 
   for (int i = 1; i <= atom->ntypes; i++) ratio[i] = 1.0;
+  oflag = aflag = 0;
   tally = 0;
   zeroflag = 0;
 
@@ -96,9 +102,29 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
       else if (strcmp(arg[iarg+1],"yes") == 0) zeroflag = 1;
       else error->all("Illegal fix langevin command");
       iarg += 2;
+    } else if (strcmp(arg[iarg],"omega") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix langevin command");
+      if (strcmp(arg[iarg+1],"no") == 0) oflag = 0;
+      else if (strcmp(arg[iarg+1],"yes") == 0) oflag = 1;
+      else error->all("Illegal fix langevin command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"angmom") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix langevin command");
+      if (strcmp(arg[iarg+1],"no") == 0) aflag = 0;
+      else if (strcmp(arg[iarg+1],"yes") == 0) aflag = 1;
+      else error->all("Illegal fix langevin command");
+      iarg += 2;
     } else error->all("Illegal fix langevin command");
   }
 
+  // error check
+
+  if (aflag) {
+    avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+    if (!avec) 
+      error->all("Fix langevin angmom requires atom style ellipsoid");
+  }
+
   // set temperature = NULL, user can override via fix_modify if wants bias
 
   id_temp = NULL;
@@ -140,6 +166,35 @@ int FixLangevin::setmask()
 
 void FixLangevin::init()
 {
+  if (oflag && !atom->sphere_flag)
+    error->all("Fix langevin omega require atom style sphere");
+  if (aflag && !atom->ellipsoid_flag)
+    error->all("Fix langevin angmom require atom style ellipsoid");
+
+  // if oflag or aflag set, check that all group particles are finite-size
+
+  if (oflag) {
+    double *radius = atom->radius;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit)
+	if (radius[i] == 0.0)
+	  error->one("Fix langevin omega requires extended particles");
+  }
+
+  if (aflag) {
+    int *ellipsoid = atom->ellipsoid;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit)
+	if (ellipsoid[i] < 0)
+	  error->one("Fix langevin angmom requires extended particles");
+  }
+
   // set force prefactors
 
   if (!atom->rmass) {
@@ -219,6 +274,11 @@ void FixLangevin::post_force_no_tally()
   double fran[3],fsum[3],fsumall[3];
   fsum[0] = fsum[1] = fsum[2] = 0.0;
   bigint count;
+
+  double boltz = force->boltz;
+  double dt = update->dt;
+  double mvv2e = force->mvv2e;
+  double ftm2v = force->ftm2v;
   
   if (zeroflag) {
     count = group->count(igroup);
@@ -227,11 +287,6 @@ void FixLangevin::post_force_no_tally()
   }
   
   if (rmass) {
-    double boltz = force->boltz;
-    double dt = update->dt;
-    double mvv2e = force->mvv2e;
-    double ftm2v = force->ftm2v;
-
     if (which == NOBIAS) {
       for (int i = 0; i < nlocal; i++) {
 	if (mask[i] & groupbit) {
@@ -280,7 +335,6 @@ void FixLangevin::post_force_no_tally()
   } else {
     
     if (which == NOBIAS) {
-
       for (int i = 0; i < nlocal; i++) {
 	if (mask[i] & groupbit) {
 	  gamma1 = gfactor1[type[i]];
@@ -295,7 +349,6 @@ void FixLangevin::post_force_no_tally()
 	  fsum[1] += fran[1];
 	  fsum[2] += fran[2];
 	}
-
       }
 
     } else if (which == BIAS) {
@@ -338,6 +391,11 @@ void FixLangevin::post_force_no_tally()
       }
     }
   }
+
+  // thermostat omega and angmom
+
+  if (oflag) omega_thermostat(tsqrt);
+  if (aflag) angmom_thermostat(tsqrt);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -373,12 +431,12 @@ void FixLangevin::post_force_tally()
   //   test v = 0 since some computes mask non-participating atoms via v = 0
   //   and added force has extra term not multiplied by v = 0
 
-  if (rmass) {
+  double boltz = force->boltz;
-    double boltz = force->boltz;
+  double dt = update->dt;
-    double dt = update->dt;
+  double mvv2e = force->mvv2e;
-    double mvv2e = force->mvv2e;
+  double ftm2v = force->ftm2v;
-    double ftm2v = force->ftm2v;
 
+  if (rmass) {
     if (which == NOBIAS) {
       for (int i = 0; i < nlocal; i++) {
 	if (mask[i] & groupbit) {
@@ -454,6 +512,100 @@ void FixLangevin::post_force_tally()
       }
     }
   }
+
+  // thermostat omega and angmom
+
+  if (oflag) omega_thermostat(tsqrt);
+  if (aflag) angmom_thermostat(tsqrt);
+}
+
+/* ----------------------------------------------------------------------
+   thermostat rotational dof via omega
+------------------------------------------------------------------------- */
+
+void FixLangevin::omega_thermostat(double tsqrt)
+{
+  double gamma1,gamma2;
+
+  double boltz = force->boltz;
+  double dt = update->dt;
+  double mvv2e = force->mvv2e;
+  double ftm2v = force->ftm2v;
+
+  double **torque = atom->torque;
+  double **omega = atom->omega;
+  double *radius = atom->radius;
+  double *rmass = atom->rmass;
+  int *mask = atom->mask;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+
+  double tran[3];
+  double inertiaone;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      inertiaone = SINERTIA*radius[i]*radius[i]*rmass[i];
+      gamma1 = -inertiaone / t_period / ftm2v;
+      gamma2 = sqrt(inertiaone) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
+      gamma1 *= 1.0/ratio[type[i]];
+      gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
+      tran[0] = gamma2*(random->uniform()-0.5);
+      tran[1] = gamma2*(random->uniform()-0.5);
+      tran[2] = gamma2*(random->uniform()-0.5);
+      torque[i][0] += gamma1*omega[i][0] + tran[0];
+      torque[i][1] += gamma1*omega[i][1] + tran[1];
+      torque[i][2] += gamma1*omega[i][2] + tran[2];
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   thermostat rotational dof via angmom
+------------------------------------------------------------------------- */
+
+void FixLangevin::angmom_thermostat(double tsqrt)
+{
+  double gamma1,gamma2;
+
+  double boltz = force->boltz;
+  double dt = update->dt;
+  double mvv2e = force->mvv2e;
+  double ftm2v = force->ftm2v;
+
+  AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+  double **torque = atom->torque;
+  double **angmom = atom->angmom;
+  double *rmass = atom->rmass;
+  int *ellipsoid = atom->ellipsoid;
+  int *mask = atom->mask;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+
+  double inertia[3],wbody[3],omega[3],tran[3],rot[3][3];
+  double *shape,*quat;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      shape = bonus[ellipsoid[i]].shape;
+      inertia[0] = EINERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
+      inertia[1] = EINERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
+      inertia[2] = EINERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
+      quat = bonus[ellipsoid[i]].quat;
+      MathExtra::mq_to_omega(angmom[i],quat,inertia,omega);
+      
+      gamma1 = -1.0 / t_period / ftm2v;
+      gamma2 = sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
+      gamma1 *= 1.0/ratio[type[i]];
+      gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
+      tran[0] = sqrt(inertia[0])*gamma2*(random->uniform()-0.5);
+      tran[1] = sqrt(inertia[1])*gamma2*(random->uniform()-0.5);
+      tran[2] = sqrt(inertia[2])*gamma2*(random->uniform()-0.5);
+      torque[i][0] += inertia[0]*gamma1*omega[0] + tran[0];
+      torque[i][1] += inertia[1]*gamma1*omega[1] + tran[1];
+      torque[i][2] += inertia[2]*gamma1*omega[2] + tran[2];
+    }
+  }
 }
 
 /* ----------------------------------------------------------------------

src/fix_langevin.h

@@ -41,11 +41,13 @@ class FixLangevin : public Fix {
   double memory_usage();
 
  protected:
-  int which,tally,zeroflag;
+  int which,tally,zeroflag,oflag,aflag;
   double t_start,t_stop,t_period;
   double *gfactor1,*gfactor2,*ratio;
   double energy,energy_onestep;
 
+  class AtomVecEllipsoid *avec;
+
   int nmax;
   double **flangevin;
 
@@ -57,6 +59,8 @@ class FixLangevin : public Fix {
 
   virtual void post_force_no_tally();
   virtual void post_force_tally();
+  void omega_thermostat(double);
+  void angmom_thermostat(double);
 };
 
 }

src/fix_rigid.cpp

@@ -25,6 +25,7 @@
 #include "modify.h"
 #include "group.h"
 #include "comm.h"
+#include "random_mars.h"
 #include "force.h"
 #include "output.h"
 #include "memory.h"
@@ -45,11 +46,16 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
 {
   int i,ibody;
 
+  scalar_flag = 1;
+  extscalar = 0;
   time_integrate = 1;
   rigid_flag = 1;
   virial_flag = 1;
   create_attribute = 1;
 
+  MPI_Comm_rank(world,&me);
+  MPI_Comm_size(world,&nprocs);
+
   // perform initial allocation of atom-based arrays
   // register with Atom class
 
@@ -193,12 +199,14 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
   memory->create(imagebody,nbody,"rigid:imagebody");
   memory->create(fflag,nbody,3,"rigid:fflag");
   memory->create(tflag,nbody,3,"rigid:tflag");
+  memory->create(langextra,nbody,6,"rigid:langextra");
 
   memory->create(sum,nbody,6,"rigid:sum");
   memory->create(all,nbody,6,"rigid:all");
   memory->create(remapflag,nbody,4,"rigid:remapflag");
 
   // initialize force/torque flags to default = 1.0
+  // for 2d: fz, tx, ty = 0.0
 
   array_flag = 1;
   size_array_rows = nbody;
@@ -209,10 +217,13 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
   for (i = 0; i < nbody; i++) {
     fflag[i][0] = fflag[i][1] = fflag[i][2] = 1.0;
     tflag[i][0] = tflag[i][1] = tflag[i][2] = 1.0;
+    if (domain->dimension == 2) fflag[i][2] = tflag[i][0] = tflag[i][1] = 0.0;
   }
 
   // parse optional args
 
+  int seed;
+  langflag = 0;
   tempflag = 0;
   pressflag = 0;
   t_chain = 10;
@@ -238,6 +249,9 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
       else if (strcmp(arg[iarg+4],"on") == 0) zflag = 1.0;
       else error->all("Illegal fix rigid command");
 
+      if (domain->dimension == 2 && zflag == 1.0)
+	error->all("Fix rigid z force cannot be on for 2d simulation");
+
       int count = 0;
       for (int m = mlo; m <= mhi; m++) {
 	fflag[m-1][0] = xflag;
@@ -266,6 +280,9 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
       else if (strcmp(arg[iarg+4],"on") == 0) zflag = 1.0;
       else error->all("Illegal fix rigid command");
 
+      if (domain->dimension == 2 && (xflag == 1.0 || yflag == 1.0))
+	  error->all("Fix rigid xy torque cannot be on for 2d simulation");
+
       int count = 0;
       for (int m = mlo; m <= mhi; m++) {
 	tflag[m-1][0] = xflag;
@@ -277,10 +294,24 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
 
       iarg += 5;
 
+    } else if (strcmp(arg[iarg],"langevin") == 0) {
+      if (iarg+5 > narg) error->all("Illegal fix rigid command");
+      if (strcmp(style,"rigid") != 0 && strcmp(style,"rigid/nve") != 0)
+	error->all("Illegal fix rigid command");
+      langflag = 1;
+      t_start = atof(arg[iarg+1]);
+      t_stop = atof(arg[iarg+2]);
+      t_period = atof(arg[iarg+3]);
+      seed = atoi(arg[iarg+4]);
+      if (t_period <= 0.0) 
+	error->all("Fix rigid langevin period must be > 0.0");
+      if (seed <= 0) error->all("Illegal fix rigid command");
+      iarg += 5;
+
     } else if (strcmp(arg[iarg],"temp") == 0) {
       if (iarg+4 > narg) error->all("Illegal fix rigid command");
       if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0)
-	error->all("Illegal fix/rigid command");
+	error->all("Illegal fix rigid command");
       tempflag = 1;
       t_start = atof(arg[iarg+1]);
       t_stop = atof(arg[iarg+2]);
@@ -290,7 +321,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
     } else if (strcmp(arg[iarg],"press") == 0) {
       if (iarg+4 > narg) error->all("Illegal fix rigid command");
       if (strcmp(style,"rigid/npt") != 0)
-	error->all("Illegal fix/rigid command");
+	error->all("Illegal fix rigid command");
       pressflag = 1;
       p_start = atof(arg[iarg+1]);
       p_stop = atof(arg[iarg+2]);
@@ -300,7 +331,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
     } else if (strcmp(arg[iarg],"tparam") == 0) {
       if (iarg+4 > narg) error->all("Illegal fix rigid command");
       if (strcmp(style,"rigid/nvt") != 0)
-	error->all("Illegal fix/rigid command");
+	error->all("Illegal fix rigid command");
       t_chain = atoi(arg[iarg+1]);
       t_iter = atoi(arg[iarg+2]);
       t_order = atoi(arg[iarg+3]);
@@ -309,13 +340,18 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
     } else if (strcmp(arg[iarg],"pparam") == 0) {
       if (iarg+2 > narg) error->all("Illegal fix rigid command");
       if (strcmp(style,"rigid/npt") != 0)
-	error->all("Illegal fix/rigid command");
+	error->all("Illegal fix rigid command");
       p_chain = atoi(arg[iarg+1]);
       iarg += 2;
 
     } else error->all("Illegal fix rigid command");
   }
 
+  // initialize Marsaglia RNG with processor-unique seed
+
+  if (langflag) random = new RanMars(lmp,seed + me);
+  else random = NULL;
+
   // initialize vector output quantities in case accessed before run
 
   for (i = 0; i < nbody; i++) {
@@ -369,7 +405,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
   int nsum = 0;
   for (ibody = 0; ibody < nbody; ibody++) nsum += nrigid[ibody];
   
-  if (comm->me == 0) {
+  if (me == 0) {
     if (screen) fprintf(screen,"%d rigid bodies with %d atoms\n",nbody,nsum);
     if (logfile) fprintf(logfile,"%d rigid bodies with %d atoms\n",nbody,nsum);
   }
@@ -383,6 +419,8 @@ FixRigid::~FixRigid()
 
   atom->delete_callback(id,0);
   
+  delete random;
+
   // delete locally stored arrays
   
   memory->destroy(body);
@@ -409,6 +447,7 @@ FixRigid::~FixRigid()
   memory->destroy(imagebody);
   memory->destroy(fflag);
   memory->destroy(tflag);
+  memory->destroy(langextra);
 
   memory->destroy(sum);
   memory->destroy(all);
@@ -422,6 +461,7 @@ int FixRigid::setmask()
   int mask = 0;
   mask |= INITIAL_INTEGRATE;
   mask |= FINAL_INTEGRATE;
+  if (langflag) mask |= POST_FORCE;
   mask |= PRE_NEIGHBOR;
   mask |= INITIAL_INTEGRATE_RESPA;
   mask |= FINAL_INTEGRATE_RESPA;
@@ -441,7 +481,7 @@ void FixRigid::init()
   int count = 0;
   for (int i = 0; i < modify->nfix; i++)
     if (strcmp(modify->fix[i]->style,"rigid") == 0) count++;
-  if (count > 1 && comm->me == 0) error->warning("More than one fix rigid");
+  if (count > 1 && me == 0) error->warning("More than one fix rigid");
 
   // error if npt,nph fix comes before rigid fix
 
@@ -855,6 +895,15 @@ void FixRigid::init()
 	fabs(all[ibody][5]/norm) > TOLERANCE)
       error->all("Fix rigid: Bad principal moments");
   }
+
+  // temperature scale factor
+
+  double ndof = 0.0;
+  for (ibody = 0; ibody < nbody; ibody++) {
+    ndof += fflag[ibody][0] + fflag[ibody][1] + fflag[ibody][2];
+    ndof += tflag[ibody][0] + tflag[ibody][1] + tflag[ibody][2];
+  }
+  tfactor = force->mvv2e / (ndof * force->boltz);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -998,6 +1047,13 @@ void FixRigid::setup(int vflag)
     torque[ibody][2] = all[ibody][5];
   }
 
+  // zero langextra in case Langevin thermostat not used
+  // no point to calling post_force() here since langextra
+  //   is only added to fcm/torque in final_integrate()
+
+  for (ibody = 0; ibody < nbody; ibody++)
+    for (i = 0; i < 6; i++) langextra[ibody][i] = 0.0;
+
   // virial setup before call to set_v
 
   if (vflag) v_setup(vflag);
@@ -1072,6 +1128,50 @@ void FixRigid::initial_integrate(int vflag)
   set_xv();
 }
 
+/* ----------------------------------------------------------------------
+   apply Langevin thermostat to all 6 DOF of rigid bodies
+   computed by proc 0, broadcast to other procs
+   unlike fix langevin, this stores extra force in extra arrays,
+     which are added in when final_integrate() calculates a new fcm/torque
+------------------------------------------------------------------------- */
+
+void FixRigid::post_force(int vflag)
+{
+  if (me == 0) {
+    double gamma1,gamma2;
+
+    double delta = update->ntimestep - update->beginstep;
+    delta /= update->endstep - update->beginstep;
+    double t_target = t_start + delta * (t_stop-t_start);
+    double tsqrt = sqrt(t_target);
+
+    double boltz = force->boltz;
+    double dt = update->dt;
+    double mvv2e = force->mvv2e;
+    double ftm2v = force->ftm2v;
+    
+    for (int i = 0; i < nbody; i++) {
+      gamma1 = -masstotal[i] / t_period / ftm2v;
+      gamma2 = sqrt(masstotal[i]) * tsqrt * 
+	sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
+      langextra[i][0] = gamma1*vcm[i][0] + gamma2*(random->uniform()-0.5);
+      langextra[i][1] = gamma1*vcm[i][1] + gamma2*(random->uniform()-0.5);
+      langextra[i][2] = gamma1*vcm[i][2] + gamma2*(random->uniform()-0.5);
+
+      gamma1 = -1.0 / t_period / ftm2v;
+      gamma2 = tsqrt * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
+      langextra[i][3] = inertia[i][0]*gamma1*omega[i][0] + 
+	sqrt(inertia[i][0])*gamma2*(random->uniform()-0.5);
+      langextra[i][4] = inertia[i][1]*gamma1*omega[i][1] + 
+	sqrt(inertia[i][1])*gamma2*(random->uniform()-0.5);
+      langextra[i][5] = inertia[i][2]*gamma1*omega[i][2] + 
+	sqrt(inertia[i][2])*gamma2*(random->uniform()-0.5);
+    }
+  }
+
+  MPI_Bcast(&langextra[0][0],6*nbody,MPI_DOUBLE,0,world);
+}
+
 /* ---------------------------------------------------------------------- */
 
 void FixRigid::final_integrate()
@@ -1150,13 +1250,17 @@ void FixRigid::final_integrate()
 
   MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
   
+  // update vcm and angmom
+  // include Langevin thermostat forces
+  // fflag,tflag = 0 for some dimensions in 2d
+
   for (ibody = 0; ibody < nbody; ibody++) {
-    fcm[ibody][0] = all[ibody][0];
+    fcm[ibody][0] = all[ibody][0] + langextra[ibody][0];
-    fcm[ibody][1] = all[ibody][1];
+    fcm[ibody][1] = all[ibody][1] + langextra[ibody][1];
-    fcm[ibody][2] = all[ibody][2];
+    fcm[ibody][2] = all[ibody][2] + langextra[ibody][2];
-    torque[ibody][0] = all[ibody][3];
+    torque[ibody][0] = all[ibody][3] + langextra[ibody][3];
-    torque[ibody][1] = all[ibody][4];
+    torque[ibody][1] = all[ibody][4] + langextra[ibody][4];
-    torque[ibody][2] = all[ibody][5];
+    torque[ibody][2] = all[ibody][5] + langextra[ibody][5];
 
     // update vcm by 1/2 step
   
@@ -1360,7 +1464,7 @@ int FixRigid::dof(int igroup)
     if (nall[ibody]+mall[ibody] > 0 && 
 	nall[ibody]+mall[ibody] != nrigid[ibody]) flag = 1;
   }
-  if (flag && comm->me == 0)
+  if (flag && me == 0)
     error->warning("Computing temperature of portions of rigid bodies");
 
   // remove appropriate DOFs for each rigid body wholly in temperature group
@@ -1834,6 +1938,42 @@ void FixRigid::reset_dt()
   dtq = 0.5 * update->dt;
 }
 
+/* ----------------------------------------------------------------------
+   return temperature of collection of rigid bodies
+   non-active DOF are removed by fflag/tflag and in tfactor
+------------------------------------------------------------------------- */
+
+double FixRigid::compute_scalar()
+{
+  double wbody[3],rot[3][3];
+
+  double t = 0.0;
+
+  for (int i = 0; i < nbody; i++) {
+    t += masstotal[i] * (fflag[i][0]*vcm[i][0]*vcm[i][0] + 
+    			 fflag[i][1]*vcm[i][1]*vcm[i][1] +	\
+    			 fflag[i][2]*vcm[i][2]*vcm[i][2]);
+    
+    // wbody = angular velocity in body frame
+      
+    MathExtra::quat_to_mat(quat[i],rot);
+    MathExtra::transpose_matvec(rot,angmom[i],wbody);
+    if (inertia[i][0] == 0.0) wbody[0] = 0.0;
+    else wbody[0] /= inertia[i][0];
+    if (inertia[i][1] == 0.0) wbody[1] = 0.0;
+    else wbody[1] /= inertia[i][1];
+    if (inertia[i][2] == 0.0) wbody[2] = 0.0;
+    else wbody[2] /= inertia[i][2];
+    
+    t += tflag[i][0]*inertia[i][0]*wbody[0]*wbody[0] +
+      tflag[i][1]*inertia[i][1]*wbody[1]*wbody[1] + 
+      tflag[i][2]*inertia[i][2]*wbody[2]*wbody[2];
+  }
+
+  t *= tfactor;
+  return t;
+}
+
 /* ----------------------------------------------------------------------
    return attributes of a rigid body
    15 values per body

src/fix_rigid.h

@@ -32,9 +32,11 @@ class FixRigid : public Fix {
   virtual void init();
   virtual void setup(int);
   virtual void initial_integrate(int);
+  void post_force(int);
   virtual void final_integrate();
   void initial_integrate_respa(int, int, int);
   void final_integrate_respa(int, int);
+  virtual double compute_scalar();
 
   double memory_usage();
   void grow_arrays(int);
@@ -50,6 +52,7 @@ class FixRigid : public Fix {
   double compute_array(int, int);
 
  protected:
+  int me,nprocs;
   double dtv,dtf,dtq;
   double *step_respa;
   int triclinic;
@@ -70,6 +73,7 @@ class FixRigid : public Fix {
   int *imagebody;           // image flags of xcm of each rigid body
   double **fflag;           // flag for on/off of center-of-mass force
   double **tflag;           // flag for on/off of center-of-mass torque
+  double **langextra;       // Langevin thermostat forces and torques
 
   int *body;                // which body each atom is part of (-1 if none)
   double **displace;        // displacement of each atom in body coords
@@ -85,6 +89,9 @@ class FixRigid : public Fix {
   double **qorient;         // rotation state of ext particle wrt rigid body
   double **dorient;         // orientation of dipole mu wrt rigid body
 
+  double tfactor;           // scale factor on temperature of rigid bodies
+  int langflag;             // 0/1 = no/yes Langevin thermostat
+
   int tempflag;             // NVT settings
   double t_start,t_stop;
   double t_period,t_freq;
@@ -95,6 +102,7 @@ class FixRigid : public Fix {
   double p_period,p_freq;
   int p_chain;
 
+  class RanMars *random;
   class AtomVecEllipsoid *avec_ellipsoid;
 
                             // bitmasks for eflags