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

src/ASPHERE/fix_npt_asphere.cpp

@@ -27,9 +27,12 @@
 #include "kspace.h"
 #include "update.h"
 #include "domain.h"
+#include "error.h"
 
 using namespace LAMMPS_NS;
 
+enum{NOBIAS,BIAS};
+
 /* ---------------------------------------------------------------------- */
 
 FixNPTASphere::FixNPTASphere(LAMMPS *lmp, int narg, char **arg) :
@@ -41,6 +44,14 @@ FixNPTASphere::FixNPTASphere(LAMMPS *lmp, int narg, char **arg) :
 	       "quat, angmom, torque, shape");
 }
 
+/* ---------------------------------------------------------------------- */
+
+void FixNPTASphere::init()
+{
+  FixNPT::init();
+  dtq = 0.5 * update->dt;
+}
+
 /* ----------------------------------------------------------------------
    1st half of Verlet update 
 ------------------------------------------------------------------------- */
@@ -48,6 +59,7 @@ FixNPTASphere::FixNPTASphere(LAMMPS *lmp, int narg, char **arg) :
 void FixNPTASphere::initial_integrate(int vflag)
 {
   int i;
+  double dtfm;
 
   double delta = update->ntimestep - update->beginstep;
   delta /= update->endstep - update->beginstep;
@@ -77,9 +89,9 @@ void FixNPTASphere::initial_integrate(int vflag)
     factor[i] = exp(-dthalf*(eta_dot+omega_dot[i]));
     dilation[i] = exp(dthalf*omega_dot[i]);
   }
-  ang_factor = exp(-dthalf*eta_dot);
+  factor_rotate = exp(-dthalf*eta_dot);
 
-  // v update only for atoms in group
+  // update v of only atoms in group
 
   double **x = atom->x;
   double **v = atom->v;
@@ -93,13 +105,25 @@ void FixNPTASphere::initial_integrate(int vflag)
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
-  double dtfm;
+  if (which == NOBIAS) {
-  for (i = 0; i < nlocal; i++) {
+    for (i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) {
+      if (mask[i] & groupbit) {
-      dtfm = dtf / mass[type[i]];
+	dtfm = dtf / mass[type[i]];
-      v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+	v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
-      v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
-      v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+      }
+    }
+  } else if (which == BIAS) {
+    for (i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+	temperature->remove_bias(i,v[i]);
+	dtfm = dtf / mass[type[i]];
+	v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+	v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+	v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+	temperature->restore_bias(v[i]);
+      }
     }
   }
 
@@ -120,11 +144,12 @@ void FixNPTASphere::initial_integrate(int vflag)
   // update angular momentum by 1/2 step
   // update quaternion a full step via Richardson iteration
   // returns new normalized quaternion
+
   for (i = 0; i < nlocal; i++) {    
     if (mask[i] & groupbit) {
-      angmom[i][0] = angmom[i][0] * ang_factor + dtf * torque[i][0];
+      angmom[i][0] = angmom[i][0]*factor_rotate + dtf*torque[i][0];
-      angmom[i][1] = angmom[i][1] * ang_factor + dtf * torque[i][1];
+      angmom[i][1] = angmom[i][1]*factor_rotate + dtf*torque[i][1];
-      angmom[i][2] = angmom[i][2] * ang_factor + dtf * torque[i][2];
+      angmom[i][2] = angmom[i][2]*factor_rotate + dtf*torque[i][2];
 		
       double inertia[3];
       calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia);
@@ -146,8 +171,9 @@ void FixNPTASphere::initial_integrate(int vflag)
 void FixNPTASphere::final_integrate()
 {
   int i;
+  double dtfm;
 
-  // v update only for atoms in group
+  // update v of only atoms in group
 
   double **v = atom->v;
   double **f = atom->f;
@@ -159,16 +185,31 @@ void FixNPTASphere::final_integrate()
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
-  double dtfm;
+  if (which == NOBIAS) { 
-  for (i = 0; i < nlocal; i++) {
+    for (i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) {
+      if (mask[i] & groupbit) {
-      dtfm = dtf / mass[type[i]];
+	dtfm = dtf / mass[type[i]];
-      v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+	v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
-      v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+	v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
-      v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+	v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
-      angmom[i][0] = (angmom[i][0] + dtf * torque[i][0]) * ang_factor;
+	angmom[i][0] = (angmom[i][0] + dtf * torque[i][0]) * factor_rotate;
-      angmom[i][1] = (angmom[i][1] + dtf * torque[i][1]) * ang_factor;
+	angmom[i][1] = (angmom[i][1] + dtf * torque[i][1]) * factor_rotate;
-      angmom[i][2] = (angmom[i][2] + dtf * torque[i][2]) * ang_factor;
+	angmom[i][2] = (angmom[i][2] + dtf * torque[i][2]) * factor_rotate;
+      }
+    }
+  } else if (which == BIAS) {
+    for (i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+	temperature->remove_bias(i,v[i]);
+	dtfm = dtf / mass[type[i]];
+	v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+	v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+	v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+	temperature->restore_bias(v[i]);
+	angmom[i][0] = (angmom[i][0] + dtf * torque[i][0]) * factor_rotate;
+	angmom[i][1] = (angmom[i][1] + dtf * torque[i][1]) * factor_rotate;
+	angmom[i][2] = (angmom[i][2] + dtf * torque[i][2]) * factor_rotate;
+      }
     }
   }
 
@@ -208,6 +249,14 @@ void FixNPTASphere::final_integrate()
   }
 }
 
+/* ---------------------------------------------------------------------- */
+
+void FixNPTASphere::reset_dt()
+{
+  FixNPT::reset_dt();
+  dtq = 0.5 * update->dt;
+}
+
 /* ----------------------------------------------------------------------
    Richardson iteration to update quaternion accurately
 ------------------------------------------------------------------------- */

src/ASPHERE/fix_npt_asphere.h

@@ -21,12 +21,14 @@ namespace LAMMPS_NS {
 class FixNPTASphere : public FixNPT {
  public:
   FixNPTASphere(class LAMMPS *, int, char **);
-  ~FixNPTASphere() {}
+  void init();
   void initial_integrate(int);
   void final_integrate();
+  void reset_dt();
 
  private:
-  double ang_factor;
+  double dtq;
+  double factor_rotate;
 
   void richardson(double *, double *, double *);
   void omega_from_mq(double *, double *, double *, double *);

src/ASPHERE/fix_nvt_asphere.cpp

@@ -32,6 +32,8 @@
 
 using namespace LAMMPS_NS;
 
+enum{NOBIAS,BIAS};
+
 /* ---------------------------------------------------------------------- */
 
 FixNVTASphere::FixNVTASphere(LAMMPS *lmp, int narg, char **arg) :
@@ -45,6 +47,14 @@ FixNVTASphere::FixNVTASphere(LAMMPS *lmp, int narg, char **arg) :
 
 /* ---------------------------------------------------------------------- */
 
+void FixNVTASphere::init()
+{
+  FixNVT::init();
+  dtq = 0.5 * update->dt;
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixNVTASphere::initial_integrate(int vflag)
 {
   double dtfm;
@@ -75,27 +85,56 @@ void FixNVTASphere::initial_integrate(int vflag)
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
-  for (int i = 0; i < nlocal; i++) {
+  if (which == NOBIAS) {
-    if (mask[i] & groupbit) {
+    for (int i = 0; i < nlocal; i++) {
-      dtfm = dtf / mass[type[i]];
+      if (mask[i] & groupbit) {
-      v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	dtfm = dtf / mass[type[i]];
-      v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
-      v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
-      x[i][0] += dtv * v[i][0];
+	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
-      x[i][1] += dtv * v[i][1];
+	x[i][0] += dtv * v[i][0];
-      x[i][2] += dtv * v[i][2];
+	x[i][1] += dtv * v[i][1];
+	x[i][2] += dtv * v[i][2];
+	
+	// update angular momentum by 1/2 step
+	// update quaternion a full step via Richardson iteration
+	// returns new normalized quaternion
+	
+	angmom[i][0] = angmom[i][0]*factor + dtf*torque[i][0];
+	angmom[i][1] = angmom[i][1]*factor + dtf*torque[i][1];
+	angmom[i][2] = angmom[i][2]*factor + dtf*torque[i][2];
+	
+	double inertia[3];
+	calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia);
+	richardson(quat[i],angmom[i],inertia);
+      }
+    }
 
-      // update angular momentum by 1/2 step
+  } else if (which == BIAS) {
-      // update quaternion a full step via Richardson iteration
+    for (int i = 0; i < nlocal; i++) {
-      // returns new normalized quaternion
+      if (mask[i] & groupbit) {
-      
+	temperature->remove_bias(i,v[i]);
-      angmom[i][0] = angmom[i][0] * factor + dtf * torque[i][0];
+	dtfm = dtf / mass[type[i]];
-      angmom[i][1] = angmom[i][1] * factor + dtf * torque[i][1];
+	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
-      angmom[i][2] = angmom[i][2] * factor + dtf * torque[i][2];
+	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
-
+	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
-      double inertia[3];
+	temperature->restore_bias(v[i]);
-      calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia);
+	x[i][0] += dtv * v[i][0];
-      richardson(quat[i],angmom[i],inertia);
+	x[i][1] += dtv * v[i][1];
+	x[i][2] += dtv * v[i][2];
+	
+	// update angular momentum by 1/2 step
+	// update quaternion a full step via Richardson iteration
+	// returns new normalized quaternion
+	
+	angmom[i][0] = angmom[i][0]*factor + dtf*torque[i][0];
+	angmom[i][1] = angmom[i][1]*factor + dtf*torque[i][1];
+	angmom[i][2] = angmom[i][2]*factor + dtf*torque[i][2];
+	
+	double inertia[3];
+	calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia);
+	richardson(quat[i],angmom[i],inertia);
+      }
     }
   }
 }
@@ -118,16 +157,32 @@ void FixNVTASphere::final_integrate()
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
-  for (int i = 0; i < nlocal; i++) {
+  if (which == NOBIAS) {
-    if (mask[i] & groupbit) {
+    for (int i = 0; i < nlocal; i++) {
-      dtfm = dtf / mass[type[i]] * factor;
+      if (mask[i] & groupbit) {
-      v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	dtfm = dtf / mass[type[i]] * factor;
-      v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
-      v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	angmom[i][0] = (angmom[i][0] + dtf*torque[i][0]) * factor;
+	angmom[i][1] = (angmom[i][1] + dtf*torque[i][1]) * factor;
+	angmom[i][2] = (angmom[i][2] + dtf*torque[i][2]) * factor;
+      }
+    }
 
-      angmom[i][0] = angmom[i][0] * factor + dtf * torque[i][0];
+  } else if (which == BIAS) {
-      angmom[i][1] = angmom[i][1] * factor + dtf * torque[i][1];
+    for (int i = 0; i < nlocal; i++) {
-      angmom[i][2] = angmom[i][2] * factor + dtf * torque[i][2];
+      if (mask[i] & groupbit) {
+	temperature->remove_bias(i,v[i]);
+	dtfm = dtf / mass[type[i]] * factor;
+	v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+	v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+	v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+	temperature->restore_bias(v[i]);
+	angmom[i][0] = (angmom[i][0] + dtf*torque[i][0]) * factor;
+	angmom[i][1] = (angmom[i][1] + dtf*torque[i][1]) * factor;
+	angmom[i][2] = (angmom[i][2] + dtf*torque[i][2]) * factor;
+      }
     }
   }
 
@@ -142,6 +197,14 @@ void FixNVTASphere::final_integrate()
   eta_dot *= drag_factor;
 }
 
+/* ---------------------------------------------------------------------- */
+
+void FixNVTASphere::reset_dt()
+{
+  FixNVT::reset_dt();
+  dtq = 0.5 * update->dt;
+}
+
 /* ----------------------------------------------------------------------
    Richardson iteration to update quaternion accurately
 ------------------------------------------------------------------------- */

src/ASPHERE/fix_nvt_asphere.h

@@ -21,11 +21,14 @@ namespace LAMMPS_NS {
 class FixNVTASphere : public FixNVT {
  public:
   FixNVTASphere(class LAMMPS *, int, char **);
-  ~FixNVTASphere() {}
+  void init();
   void initial_integrate(int);
   void final_integrate();
+  void reset_dt();
 
  private:
+  double dtq;
+
   void richardson(double *, double *, double *);
   void omega_from_mq(double *, double *, double *, double *);
   void calculate_inertia(double mass, double *shape, double *inertia);

src/compute.cpp

@@ -175,4 +175,3 @@ void Compute::restore_bias(double *v)
   v[1] += vbias[1];
   v[2] += vbias[2];
 }
-

src/fix_npt.cpp

@@ -178,6 +178,8 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
   if (strcmp(style,"npt") == 0) newarg[2] = (char *) "temp";
   else if (strcmp(style,"npt/asphere") == 0)
     newarg[2] = (char *) "temp/asphere";
+  else if (strcmp(style,"npt/sphere") == 0)
+    newarg[2] = (char *) "temp/sphere";
   modify->add_compute(3,newarg);
   delete [] newarg;
   tflag = 1;
@@ -270,7 +272,6 @@ void FixNPT::init()
 
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
-  dtq = 0.5 * update->dt;
   dthalf = 0.5 * update->dt;
 
   double freq = MAX(p_freq[0],p_freq[1]);
@@ -392,7 +393,6 @@ void FixNPT::initial_integrate(int vflag)
 	v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
       }
     }
-
   } else if (which == BIAS) {
     for (i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
@@ -456,7 +456,6 @@ void FixNPT::final_integrate()
 	v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
       }
     }
-
   } else if (which == BIAS) {
     for (i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
@@ -601,7 +600,6 @@ void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag)
 	  v[i][2] += dtfm*f[i][2];
 	}
       }
-
     } else if (which == BIAS) {
       for (int i = 0; i < nlocal; i++) {
 	if (mask[i] & groupbit) {
@@ -667,7 +665,6 @@ void FixNPT::final_integrate_respa(int ilevel)
 	  v[i][2] += dtfm*f[i][2];
 	}
       }
-
     } else if (which == BIAS) {
       for (int i = 0; i < nlocal; i++) {
 	if (mask[i] & groupbit) {
@@ -894,7 +891,6 @@ void FixNPT::reset_dt()
 {
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
-  dtq = 0.5 * update->dt;
   dthalf = 0.5 * update->dt;
 
   double freq = MAX(p_freq[0],p_freq[1]);

src/fix_npt.h

@@ -23,7 +23,7 @@ class FixNPT : public Fix {
   FixNPT(class LAMMPS *, int, char **);
   virtual ~FixNPT();
   int setmask();
-  void init();
+  virtual void init();
   void setup(int);
   void initial_integrate(int);
   virtual void final_integrate();
@@ -33,11 +33,11 @@ class FixNPT : public Fix {
   void write_restart(FILE *);
   void restart(char *);
   int modify_param(int, char **);
-  void reset_dt();
+  virtual void reset_dt();
 
  protected:
   int dimension,which;
-  double dtv,dtf,dtq,dthalf;
+  double dtv,dtf,dthalf;
   double boltz,nktv2p;
   double vol0;
 

src/fix_nvt.cpp

@@ -80,6 +80,8 @@ FixNVT::FixNVT(LAMMPS *lmp, int narg, char **arg) :
     newarg[2] = (char *) "temp/asphere";
   else if (strcmp(style,"nvt/sllod") == 0) 
     newarg[2] = (char *) "temp/deform";
+  else if (strcmp(style,"nvt/sphere") == 0) 
+    newarg[2] = (char *) "temp/sphere";
   modify->add_compute(3,newarg);
   delete [] newarg;
   tflag = 1;
@@ -126,7 +128,6 @@ void FixNVT::init()
 
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
-  dtq = 0.5 * update->dt;
   dthalf = 0.5 * update->dt;
 
   drag_factor = 1.0 - (update->dt * t_freq * drag);
@@ -435,7 +436,6 @@ void FixNVT::reset_dt()
 {
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
-  dtq = 0.5 * update->dt;
   dthalf = 0.5 * update->dt;
 
   drag_factor = 1.0 - (update->dt * t_freq * drag);

src/fix_nvt.h

@@ -23,7 +23,7 @@ class FixNVT : public Fix {
   FixNVT(class LAMMPS *, int, char **);
   virtual ~FixNVT();
   int setmask();
-  void init();
+  virtual void init();
   void setup(int);
   virtual void initial_integrate(int);
   virtual void final_integrate();
@@ -34,7 +34,7 @@ class FixNVT : public Fix {
   void restart(char *);
   int modify_param(int, char **);
   void reset_target(double);
-  void reset_dt();
+  virtual void reset_dt();
 
  protected:
   int which;
@@ -42,7 +42,7 @@ class FixNVT : public Fix {
   double t_current,t_target;
   double t_freq,drag,drag_factor;
   double f_eta,eta_dot,eta,factor;
-  double dtv,dtf,dtq,dthalf;
+  double dtv,dtf,dthalf;
 
   int nlevels_respa;
   double *step_respa;

src/fix_nvt_sllod.cpp

@@ -35,6 +35,7 @@
 using namespace LAMMPS_NS;
 
 enum{NO_REMAP,X_REMAP,V_REMAP};                   // same as fix_deform.cpp
+enum{NOBIAS,BIAS};
 
 /* ---------------------------------------------------------------------- */
 
@@ -47,6 +48,9 @@ void FixNVTSlodd::init()
 {
   FixNVT::init();
 
+  if (!temperature->tempbias)
+    error->all("Temperature for fix nvt/sllod does not have a bias");
+
   // check fix deform remap settings
 
   int i;
@@ -78,11 +82,10 @@ void FixNVTSlodd::initial_integrate(int vflag)
   eta += dtv*eta_dot;
   factor = exp(-dthalf*eta_dot);
 
-  // update vthermal and x of only atoms in group
+  // update v and x of only atoms in group
-  // lamda = 0-1 triclinic lamda coords
+  // remove and restore bias = streaming velocity = Hrate*lamda + Hratelo
-  // vstream = streaming velocity = Hrate*lamda + Hratelo
+  // thermostat thermal velocity only
-  // vthermal = thermal velocity = v - vstream
+  // vdelu = SLLOD correction = Hrate*Hinv*vthermal
-  // vdelu = Hrate*Hinv*vthermal
 
   double **x = atom->x;
   double **v = atom->v;
@@ -93,34 +96,21 @@ void FixNVTSlodd::initial_integrate(int vflag)
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
-  double *h_rate = domain->h_rate;
+  double h_two[6],vdelu[3];
-  double *h_ratelo = domain->h_ratelo;
+  MathExtra::multiply_shape_shape(domain->h_rate,domain->h_inv,h_two);
-  double h_two[6],lamda[3],vstream[3],vthermal[3],vdelu[3];
-  MathExtra::multiply_shape_shape(h_rate,domain->h_inv,h_two);
 
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
+      temperature->remove_bias(i,v[i]);
+      vdelu[0] = h_two[0]*v[i][0] + h_two[5]*v[i][1] + h_two[4]*v[i][2];
+      vdelu[1] = h_two[1]*v[i][1] + h_two[3]*v[i][2];
+      vdelu[2] = h_two[2]*v[i][2];
       dtfm = dtf / mass[type[i]];
+      v[i][0] = v[i][0]*factor + dtfm*f[i][0] - dthalf*vdelu[0];
+      v[i][1] = v[i][1]*factor + dtfm*f[i][1] - dthalf*vdelu[1];
+      v[i][2] = v[i][2]*factor + dtfm*f[i][2] - dthalf*vdelu[2];
+      temperature->restore_bias(v[i]);
 
-      domain->x2lamda(x[i],lamda);
-      vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + 
-	h_rate[4]*lamda[2] + h_ratelo[0];
-      vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
-      vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2];
-      vthermal[0] = v[i][0] - vstream[0];
-      vthermal[1] = v[i][1] - vstream[1];
-      vthermal[2] = v[i][2] - vstream[2];
-      vdelu[0] = h_two[0]*vthermal[0] + h_two[5]*vthermal[1] + 
-	h_two[4]*vthermal[2];
-      vdelu[1] = h_two[1]*vthermal[1] + h_two[3]*vthermal[2];
-      vdelu[2] = h_two[2]*vthermal[2];
-
-      v[i][0] = vstream[0] + 
-	vthermal[0]*factor + dtfm*f[i][0] - dthalf*vdelu[0];
-      v[i][1] = vstream[1] + 
-	vthermal[1]*factor + dtfm*f[i][1] - dthalf*vdelu[1];
-      v[i][2] = vstream[2] + 
-	vthermal[2]*factor + dtfm*f[i][2] - dthalf*vdelu[2];
       x[i][0] += dtv * v[i][0];
       x[i][1] += dtv * v[i][1];
       x[i][2] += dtv * v[i][2];
@@ -134,11 +124,10 @@ void FixNVTSlodd::final_integrate()
 {
   double dtfm;
 
-  // update vthermal of only atoms in group
+  // update v of only atoms in group
-  // lamda = 0-1 triclinic lamda coords
+  // remove and restore bias = streaming velocity = Hrate*lamda + Hratelo
-  // vstream = streaming velocity = Hrate*lamda + Hratelo
+  // thermostat thermal velocity only
-  // vthermal = thermal velocity = v - vstream
+  // vdelu = SLLOD correction = Hrate*Hinv*vthermal
-  // vdelu = Hrate*Hinv*vthermal
 
   double **x = atom->x;
   double **v = atom->v;
@@ -149,34 +138,21 @@ void FixNVTSlodd::final_integrate()
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
-  double *h_rate = domain->h_rate;
+  double h_two[6],vdelu[3];
-  double *h_ratelo = domain->h_ratelo;
+  MathExtra::multiply_shape_shape(domain->h_rate,domain->h_inv,h_two);
-  double h_two[6],lamda[3],vstream[3],vthermal[3],vdelu[3];
-  MathExtra::multiply_shape_shape(h_rate,domain->h_inv,h_two);
 
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
+      temperature->remove_bias(i,v[i]);
+      vdelu[0] = h_two[0]*v[i][0] + h_two[5]*v[i][1] + h_two[4]*v[i][2];
+      vdelu[1] = h_two[1]*v[i][1] + h_two[3]*v[i][2];
+      vdelu[2] = h_two[2]*v[i][2];
+      //dtfm = dtf / mass[type[i]] * factor;
       dtfm = dtf / mass[type[i]];
-
+      v[i][0] = v[i][0]*factor + dtfm*f[i][0] - dthalf*vdelu[0];
-      domain->x2lamda(x[i],lamda);
+      v[i][1] = v[i][1]*factor + dtfm*f[i][1] - dthalf*vdelu[1];
-      vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + 
+      v[i][2] = v[i][2]*factor + dtfm*f[i][2] - dthalf*vdelu[2];
-	h_rate[4]*lamda[2] + h_ratelo[0];
+      temperature->restore_bias(v[i]);
-      vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
-      vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2];
-      vthermal[0] = v[i][0] - vstream[0];
-      vthermal[1] = v[i][1] - vstream[1];
-      vthermal[2] = v[i][2] - vstream[2];
-      vdelu[0] = h_two[0]*vthermal[0] + h_two[5]*vthermal[1] + 
-	h_two[4]*vthermal[2];
-      vdelu[1] = h_two[1]*vthermal[1] + h_two[3]*vthermal[2];
-      vdelu[2] = h_two[2]*vthermal[2];
-
-      v[i][0] = vstream[0] + 
-	vthermal[0]*factor + dtfm*f[i][0] - dthalf*vdelu[0];
-      v[i][1] = vstream[1] + 
-	vthermal[1]*factor + dtfm*f[i][1] - dthalf*vdelu[1];
-      v[i][2] = vstream[2] + 
-	vthermal[2]*factor + dtfm*f[i][2] - dthalf*vdelu[2];
     }
   }
 
@@ -234,34 +210,20 @@ void FixNVTSlodd::initial_integrate_respa(int vflag, int ilevel, int flag)
 
   // update v of only atoms in group
 
-  double *h_rate = domain->h_rate;
+  double h_two[6],vdelu[3];
-  double *h_ratelo = domain->h_ratelo;
+  MathExtra::multiply_shape_shape(domain->h_rate,domain->h_inv,h_two);
-  double h_two[6],lamda[3],vstream[3],vthermal[3],vdelu[3];
-  MathExtra::multiply_shape_shape(h_rate,domain->h_inv,h_two);
 
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
+      temperature->remove_bias(i,v[i]);
+      vdelu[0] = h_two[0]*v[i][0] + h_two[5]*v[i][1] + h_two[4]*v[i][2];
+      vdelu[1] = h_two[1]*v[i][1] + h_two[3]*v[i][2];
+      vdelu[2] = h_two[2]*v[i][2];
       dtfm = dtf / mass[type[i]];
-
+      v[i][0] = v[i][0]*factor + dtfm*f[i][0] - dthalf*vdelu[0];
-      domain->x2lamda(x[i],lamda);
+      v[i][1] = v[i][1]*factor + dtfm*f[i][1] - dthalf*vdelu[1];
-      vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + 
+      v[i][2] = v[i][2]*factor + dtfm*f[i][2] - dthalf*vdelu[2];
-	h_rate[4]*lamda[2] + h_ratelo[0];
+      temperature->restore_bias(v[i]);
-      vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
-      vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2];
-      vthermal[0] = v[i][0] - vstream[0];
-      vthermal[1] = v[i][1] - vstream[1];
-      vthermal[2] = v[i][2] - vstream[2];
-      vdelu[0] = h_two[0]*vthermal[0] + h_two[5]*vthermal[1] + 
-	h_two[4]*vthermal[2];
-      vdelu[1] = h_two[1]*vthermal[1] + h_two[3]*vthermal[2];
-      vdelu[2] = h_two[2]*vthermal[2];
-
-      v[i][0] = vstream[0] + 
-	vthermal[0]*factor + dtfm*f[i][0] - dthalf*vdelu[0];
-      v[i][1] = vstream[1] + 
-	vthermal[1]*factor + dtfm*f[i][1] - dthalf*vdelu[1];
-      v[i][2] = vstream[2] + 
-	vthermal[2]*factor + dtfm*f[i][2] - dthalf*vdelu[2];
     }
   }
 

src/fix_nvt_sllod.h

@@ -21,7 +21,6 @@ namespace LAMMPS_NS {
 class FixNVTSlodd : public FixNVT {
  public:
   FixNVTSlodd(class LAMMPS *, int, char **);
-  ~FixNVTSlodd() {}
   void init();
   void initial_integrate(int);
   void final_integrate();

src/style.h

@@ -164,12 +164,14 @@ DumpStyle(xyz,DumpXYZ)
 #include "fix_momentum.h"
 #include "fix_nph.h"
 #include "fix_npt.h"
+#include "fix_npt_sphere.h"
 #include "fix_nve.h"
 #include "fix_nve_limit.h"
 #include "fix_nve_noforce.h"
 #include "fix_nve_sphere.h"
 #include "fix_nvt.h"
 #include "fix_nvt_sllod.h"
+#include "fix_nvt_sphere.h"
 #include "fix_plane_force.h"
 #include "fix_press_berendsen.h"
 #include "fix_print.h"
@@ -220,12 +222,14 @@ FixStyle(momentum,FixMomentum)
 FixStyle(msd,FixMSD)
 FixStyle(nph,FixNPH)
 FixStyle(npt,FixNPT)
+FixStyle(npt/sphere,FixNPTSphere)
 FixStyle(nve,FixNVE)
 FixStyle(nve/limit,FixNVELimit)
 FixStyle(nve/noforce,FixNVENoforce)
 FixStyle(nve/sphere,FixNVESphere)
 FixStyle(nvt,FixNVT)
 FixStyle(nvt/sllod,FixNVTSlodd)
+FixStyle(nvt/sphere,FixNVTSphere)
 FixStyle(orient/fcc,FixOrientFCC)
 FixStyle(press/berendsen,FixPressBerendsen)
 FixStyle(planeforce,FixPlaneForce)