Merge pull request #489 from akohlmey/thread-safe-biasing

port thread-safe temperature biasing from LAMMPS-ICMS

src/ASPHERE/compute_temp_asphere.cpp

@@ -396,6 +396,15 @@ void ComputeTempAsphere::remove_bias(int i, double *v)
   if (tbias) tbias->remove_bias(i,v);
 }
 
+/* ----------------------------------------------------------------------
+   remove velocity bias from atom I to leave thermal velocity
+------------------------------------------------------------------------- */
+
+void ComputeTempAsphere::remove_bias_thr(int i, double *v, double *b)
+{
+  if (tbias) tbias->remove_bias_thr(i,v,b);
+}
+
 /* ----------------------------------------------------------------------
    add back in velocity bias to atom I removed by remove_bias()
    assume remove_bias() was previously called
@@ -405,3 +414,13 @@ void ComputeTempAsphere::restore_bias(int i, double *v)
 {
   if (tbias) tbias->restore_bias(i,v);
 }
+
+/* ----------------------------------------------------------------------
+   add back in velocity bias to atom I removed by remove_bias_thr()
+   assume remove_bias_thr() was previously called with the same buffer b
+------------------------------------------------------------------------- */
+
+void ComputeTempAsphere::restore_bias_thr(int i, double *v, double *b)
+{
+  if (tbias) tbias->restore_bias_thr(i,v,b);
+}

src/ASPHERE/compute_temp_asphere.h

@@ -35,6 +35,8 @@ class ComputeTempAsphere : public Compute {
 
   void remove_bias(int, double *);
   void restore_bias(int, double *);
+  void remove_bias_thr(int, double *, double *);
+  void restore_bias_thr(int, double *, double *);
 
  private:
   int mode;

src/USER-MISC/compute_temp_rotate.cpp

@@ -221,6 +221,17 @@ void ComputeTempRotate::remove_bias(int i, double *v)
   v[2] -= vbiasall[i][2];
 }
 
+/* ----------------------------------------------------------------------
+   remove velocity bias from atom I to leave thermal velocity
+------------------------------------------------------------------------- */
+
+void ComputeTempRotate::remove_bias_thr(int i, double *v, double *)
+{
+  v[0] -= vbiasall[i][0];
+  v[1] -= vbiasall[i][1];
+  v[2] -= vbiasall[i][2];
+}
+
 /* ----------------------------------------------------------------------
    remove velocity bias from all atoms to leave thermal velocity
 ------------------------------------------------------------------------- */
@@ -251,6 +262,18 @@ void ComputeTempRotate::restore_bias(int i, double *v)
   v[2] += vbiasall[i][2];
 }
 
+/* ----------------------------------------------------------------------
+   add back in velocity bias to atom I removed by remove_bias_thr()
+   assume remove_bias_thr() was previously called
+------------------------------------------------------------------------- */
+
+void ComputeTempRotate::restore_bias_thr(int i, double *v, double *)
+{
+  v[0] += vbiasall[i][0];
+  v[1] += vbiasall[i][1];
+  v[2] += vbiasall[i][2];
+}
+
 /* ----------------------------------------------------------------------
    add back in velocity bias to all atoms removed by remove_bias_all()
    assume remove_bias_all() was previously called

src/USER-MISC/compute_temp_rotate.h

@@ -34,9 +34,12 @@ class ComputeTempRotate : public Compute {
   void compute_vector();
 
   void remove_bias(int, double *);
+  void remove_bias_thr(int, double *, double *);
   void remove_bias_all();
   void restore_bias(int, double *);
   void restore_bias_all();
+  void restore_bias_thr(int, double *, double *);
+
   double memory_usage();
 
  private:

src/USER-OMP/fix_nh_asphere_omp.cpp

@@ -183,12 +183,13 @@ void FixNHAsphereOMP::nh_v_temp()
 #pragma omp parallel for default(none) private(i) schedule(static)
 #endif
     for (i = 0; i < nlocal; i++) {
+      double buf[3];
       if (mask[i] & groupbit) {
-        temperature->remove_bias(i,&v[i].x);
+        temperature->remove_bias_thr(i,&v[i].x,buf);
         v[i].x *= factor_eta;
         v[i].y *= factor_eta;
         v[i].z *= factor_eta;
-        temperature->restore_bias(i,&v[i].x);
+        temperature->restore_bias_thr(i,&v[i].x,buf);
         angmom[i].x *= factor_eta;
         angmom[i].y *= factor_eta;
         angmom[i].z *= factor_eta;

src/USER-OMP/fix_nh_omp.cpp

@@ -261,8 +261,9 @@ void FixNHOMP::nh_v_press()
 #pragma omp parallel for default(none) private(i) schedule(static)
 #endif
     for (i = 0; i < nlocal; i++) {
+      double buf[3];
       if (mask[i] & groupbit) {
-        temperature->remove_bias(i,&v[i].x);
+        temperature->remove_bias_thr(i,&v[i].x,buf);
         v[i].x *= factor0;
         v[i].y *= factor1;
         v[i].z *= factor2;
@@ -273,7 +274,7 @@ void FixNHOMP::nh_v_press()
         v[i].x *= factor0;
         v[i].y *= factor1;
         v[i].z *= factor2;
-        temperature->restore_bias(i,&v[i].x);
+        temperature->restore_bias_thr(i,&v[i].x,buf);
       }
     }
   }
@@ -373,12 +374,13 @@ void FixNHOMP::nh_v_temp()
 #pragma omp parallel for default(none) private(i) schedule(static)
 #endif
     for (i = 0; i < nlocal; i++) {
+      double buf[3];
       if (mask[i] & groupbit) {
-        temperature->remove_bias(i,&v[i].x);
+        temperature->remove_bias_thr(i,&v[i].x,buf);
         v[i].x *= factor_eta;
         v[i].y *= factor_eta;
         v[i].z *= factor_eta;
-        temperature->restore_bias(i,&v[i].x);
+        temperature->restore_bias_thr(i,&v[i].x,buf);
       }
     }
   }

src/USER-OMP/fix_nh_sphere_omp.cpp

@@ -137,12 +137,13 @@ void FixNHSphereOMP::nh_v_temp()
 #pragma omp parallel for default(none) private(i) schedule(static)
 #endif
     for (i = 0; i < nlocal; i++) {
+      double buf[3];
       if (mask[i] & groupbit) {
-        temperature->remove_bias(i,&v[i].x);
+        temperature->remove_bias_thr(i,&v[i].x,buf);
         v[i].x *= factor_eta;
         v[i].y *= factor_eta;
         v[i].z *= factor_eta;
-        temperature->restore_bias(i,&v[i].x);
+        temperature->restore_bias_thr(i,&v[i].x,buf);
         omega[i].x *= factor_eta;
         omega[i].y *= factor_eta;
         omega[i].z *= factor_eta;

src/USER-OMP/fix_nvt_sllod_omp.cpp

@@ -121,16 +121,16 @@ void FixNVTSllodOMP::nh_v_temp()
 #pragma omp parallel for default(none) private(i) shared(h_two) schedule(static)
 #endif
   for (i = 0; i < nlocal; i++) {
-    double vdelu0,vdelu1,vdelu2;
+    double vdelu0,vdelu1,vdelu2,buf[3];
     if (mask[i] & groupbit) {
       vdelu0 = h_two[0]*v[i].x + h_two[5]*v[i].y + h_two[4]*v[i].z;
       vdelu1 = h_two[1]*v[i].y + h_two[3]*v[i].z;
       vdelu2 = h_two[2]*v[i].z;
-      temperature->remove_bias(i,&v[i].x);
+      temperature->remove_bias_thr(i,&v[i].x,buf);
       v[i].x = v[i].x*factor_eta - dthalf*vdelu0;
       v[i].y = v[i].y*factor_eta - dthalf*vdelu1;
       v[i].z = v[i].z*factor_eta - dthalf*vdelu2;
-      temperature->restore_bias(i,&v[i].x);
+      temperature->restore_bias_thr(i,&v[i].x,buf);
     }
   }
 }

src/compute.h

@@ -114,9 +114,11 @@ class Compute : protected Pointers {
   virtual void dof_remove_pre() {}
   virtual int dof_remove(int) {return 0;}
   virtual void remove_bias(int, double *) {}
+  virtual void remove_bias_thr(int, double *, double *) {}
   virtual void remove_bias_all() {}
   virtual void reapply_bias_all() {}
   virtual void restore_bias(int, double *) {}
+  virtual void restore_bias_thr(int, double *, double *) {}
   virtual void restore_bias_all() {}
 
   virtual void reset_extra_compute_fix(const char *);

src/compute_temp_com.cpp

@@ -170,6 +170,17 @@ void ComputeTempCOM::remove_bias(int i, double *v)
   v[2] -= vbias[2];
 }
 
+/* ----------------------------------------------------------------------
+   remove velocity bias from atom I to leave thermal velocity
+------------------------------------------------------------------------- */
+
+void ComputeTempCOM::remove_bias_thr(int, double *v, double *)
+{
+  v[0] -= vbias[0];
+  v[1] -= vbias[1];
+  v[2] -= vbias[2];
+}
+
 /* ----------------------------------------------------------------------
    remove velocity bias from all atoms to leave thermal velocity
 ------------------------------------------------------------------------- */
@@ -200,6 +211,18 @@ void ComputeTempCOM::restore_bias(int i, double *v)
   v[2] += vbias[2];
 }
 
+/* ----------------------------------------------------------------------
+   add back in velocity bias to atom I removed by remove_bias_thr()
+   assume remove_bias_thr() was previously called
+------------------------------------------------------------------------- */
+
+void ComputeTempCOM::restore_bias_thr(int, double *v, double *)
+{
+  v[0] += vbias[0];
+  v[1] += vbias[1];
+  v[2] += vbias[2];
+}
+
 /* ----------------------------------------------------------------------
    add back in velocity bias to all atoms removed by remove_bias_all()
    assume remove_bias_all() was previously called

src/compute_temp_com.h

@@ -34,9 +34,11 @@ class ComputeTempCOM : public Compute {
   void compute_vector();
 
   void remove_bias(int, double *);
+  void remove_bias_thr(int, double *, double *);
   void remove_bias_all();
   void restore_bias(int, double *);
   void restore_bias_all();
+  void restore_bias_thr(int, double *, double *);
 
  private:
   double tfactor,masstotal;

src/compute_temp_deform.cpp

@@ -221,6 +221,26 @@ void ComputeTempDeform::remove_bias(int i, double *v)
   v[2] -= vbias[2];
 }
 
+/* ----------------------------------------------------------------------
+   remove velocity bias from atom I to leave thermal velocity
+------------------------------------------------------------------------- */
+
+void ComputeTempDeform::remove_bias_thr(int i, double *v, double *b)
+{
+  double lamda[3];
+  double *h_rate = domain->h_rate;
+  double *h_ratelo = domain->h_ratelo;
+
+  domain->x2lamda(atom->x[i],lamda);
+  b[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] +
+    h_rate[4]*lamda[2] + h_ratelo[0];
+  b[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
+  b[2] = h_rate[2]*lamda[2] + h_ratelo[2];
+  v[0] -= b[0];
+  v[1] -= b[1];
+  v[2] -= b[2];
+}
+
 /* ----------------------------------------------------------------------
    remove velocity bias from all atoms to leave thermal velocity
 ------------------------------------------------------------------------- */
@@ -266,6 +286,18 @@ void ComputeTempDeform::restore_bias(int i, double *v)
   v[2] += vbias[2];
 }
 
+/* ----------------------------------------------------------------------
+   add back in velocity bias to atom I removed by remove_bias_thr()
+   assume remove_bias_thr() was previously called with the same buffer b
+------------------------------------------------------------------------- */
+
+void ComputeTempDeform::restore_bias_thr(int i, double *v, double *b)
+{
+  v[0] += b[0];
+  v[1] += b[1];
+  v[2] += b[2];
+}
+
 /* ----------------------------------------------------------------------
    add back in velocity bias to all atoms removed by remove_bias_all()
    assume remove_bias_all() was previously called

src/compute_temp_deform.h

@@ -34,8 +34,10 @@ class ComputeTempDeform : public Compute {
   virtual void compute_vector();
 
   void remove_bias(int, double *);
+  void remove_bias_thr(int, double *, double *);
   void remove_bias_all();
   void restore_bias(int, double *);
+  void restore_bias_thr(int, double *, double *);
   void restore_bias_all();
   double memory_usage();
 

src/compute_temp_partial.cpp

@@ -185,6 +185,26 @@ void ComputeTempPartial::remove_bias(int i, double *v)
   }
 }
 
+/* ----------------------------------------------------------------------
+   remove velocity bias from atom I to leave thermal velocity
+------------------------------------------------------------------------- */
+
+void ComputeTempPartial::remove_bias_thr(int i, double *v, double *b)
+{
+  if (!xflag) {
+    b[0] = v[0];
+    v[0] = 0.0;
+  }
+  if (!yflag) {
+    b[1] = v[1];
+    v[1] = 0.0;
+  }
+  if (!zflag) {
+    b[2] = v[2];
+    v[2] = 0.0;
+  }
+}
+
 /* ----------------------------------------------------------------------
    remove velocity bias from all atoms to leave thermal velocity
 ------------------------------------------------------------------------- */
@@ -262,6 +282,18 @@ void ComputeTempPartial::restore_bias(int i, double *v)
   if (!zflag) v[2] += vbias[2];
 }
 
+/* ----------------------------------------------------------------------
+   add back in velocity bias to atom I removed by remove_bias_thr()
+   assume remove_bias_thr() was previously called with the same buffer b
+------------------------------------------------------------------------- */
+
+void ComputeTempPartial::restore_bias_thr(int i, double *v, double *b)
+{
+  if (!xflag) v[0] += b[0];
+  if (!yflag) v[1] += b[1];
+  if (!zflag) v[2] += b[2];
+}
+
 /* ----------------------------------------------------------------------
    add back in velocity bias to all atoms removed by remove_bias_all()
    assume remove_bias_all() was previously called

src/compute_temp_partial.h

@@ -35,9 +35,11 @@ class ComputeTempPartial : public Compute {
 
   int dof_remove(int);
   void remove_bias(int, double *);
+  void remove_bias_thr(int, double *, double *);
   void remove_bias_all();
   void reapply_bias_all();
   void restore_bias(int, double *);
+  void restore_bias_thr(int, double *, double *);
   void restore_bias_all();
   double memory_usage();
 

src/compute_temp_profile.cpp

@@ -359,6 +359,15 @@ void ComputeTempProfile::remove_bias(int i, double *v)
   if (zflag) v[2] -= binave[ibin][ivz];
 }
 
+/* ----------------------------------------------------------------------
+   remove velocity bias from atom I to leave thermal velocity
+------------------------------------------------------------------------- */
+
+void ComputeTempProfile::remove_bias_thr(int i, double *v, double *)
+{
+  remove_bias(i,v);
+}
+
 /* ----------------------------------------------------------------------
    remove velocity bias from all atoms to leave thermal velocity
 ------------------------------------------------------------------------- */
@@ -392,6 +401,16 @@ void ComputeTempProfile::restore_bias(int i, double *v)
   if (zflag) v[2] += binave[ibin][ivz];
 }
 
+/* ----------------------------------------------------------------------
+   add back in velocity bias to atom I removed by remove_bias_thr()
+   assume remove_bias_thr() was previously called
+------------------------------------------------------------------------- */
+
+void ComputeTempProfile::restore_bias_thr(int i, double *v, double *)
+{
+  restore_bias(i,v);
+}
+
 /* ----------------------------------------------------------------------
    add back in velocity bias to all atoms removed by remove_bias_all()
    assume remove_bias_all() was previously called

src/compute_temp_profile.h

@@ -35,8 +35,10 @@ class ComputeTempProfile : public Compute {
   void compute_array();
 
   void remove_bias(int, double *);
+  void remove_bias_thr(int, double *, double *);
   void remove_bias_all();
   void restore_bias(int, double *);
+  void restore_bias_thr(int, double *, double *);
   void restore_bias_all();
   double memory_usage();
 

src/compute_temp_ramp.cpp

@@ -234,6 +234,19 @@ void ComputeTempRamp::remove_bias(int i, double *v)
   v[v_dim] -= vbias[v_dim];
 }
 
+/* ----------------------------------------------------------------------
+   remove velocity bias from atom I to leave thermal velocity
+------------------------------------------------------------------------- */
+
+void ComputeTempRamp::remove_bias_thr(int i, double *v, double *b)
+{
+  double fraction = (atom->x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo);
+  fraction = MAX(fraction,0.0);
+  fraction = MIN(fraction,1.0);
+  b[v_dim] = v_lo + fraction*(v_hi - v_lo);
+  v[v_dim] -= b[v_dim];
+}
+
 /* ----------------------------------------------------------------------
    remove velocity bias from all atoms to leave thermal velocity
 ------------------------------------------------------------------------- */
@@ -271,6 +284,16 @@ void ComputeTempRamp::restore_bias(int i, double *v)
   v[v_dim] += vbias[v_dim];
 }
 
+/* ----------------------------------------------------------------------
+   add back in velocity bias to atom I removed by remove_bias_thr()
+   assume remove_bias_thr() was previously called with the same buffer b
+------------------------------------------------------------------------- */
+
+void ComputeTempRamp::restore_bias_thr(int i, double *v, double *b)
+{
+  v[v_dim] += b[v_dim];
+}
+
 /* ----------------------------------------------------------------------
    add back in velocity bias to all atoms removed by remove_bias_all()
    assume remove_bias_all() was previously called

src/compute_temp_ramp.h

@@ -35,7 +35,9 @@ class ComputeTempRamp : public Compute {
 
   void remove_bias(int, double *);
   void remove_bias_all();
+  void remove_bias_thr(int, double *, double *);
   void restore_bias(int, double *);
+  void restore_bias_thr(int, double *, double *);
   void restore_bias_all();
   double memory_usage();
 

src/compute_temp_region.cpp

@@ -199,6 +199,23 @@ void ComputeTempRegion::remove_bias(int i, double *v)
   }
 }
 
+/* ----------------------------------------------------------------------
+   remove velocity bias from atom I to leave thermal velocity
+------------------------------------------------------------------------- */
+
+void ComputeTempRegion::remove_bias_thr(int i, double *v, double *b)
+{
+  double *x = atom->x[i];
+  if (domain->regions[iregion]->match(x[0],x[1],x[2]))
+    b[0] = b[1] = b[2] = 0.0;
+  else {
+    b[0] = v[0];
+    b[1] = v[1];
+    b[2] = v[2];
+    v[0] = v[1] = v[2] = 0.0;
+  }
+}
+
 /* ----------------------------------------------------------------------
    remove velocity bias from all atoms to leave thermal velocity
 ------------------------------------------------------------------------- */
@@ -243,6 +260,18 @@ void ComputeTempRegion::restore_bias(int i, double *v)
   v[2] += vbias[2];
 }
 
+/* ----------------------------------------------------------------------
+   add back in velocity bias to atom I removed by remove_bias_thr()
+   assume remove_bias_thr() was previously called with the same buffer b
+------------------------------------------------------------------------- */
+
+void ComputeTempRegion::restore_bias_thr(int i, double *v, double *b)
+{
+  v[0] += b[0];
+  v[1] += b[1];
+  v[2] += b[2];
+}
+
 /* ----------------------------------------------------------------------
    add back in velocity bias to all atoms removed by remove_bias_all()
    assume remove_bias_all() was previously called

src/compute_temp_region.h

@@ -35,10 +35,13 @@ class ComputeTempRegion : public Compute {
 
   void dof_remove_pre();
   int dof_remove(int);
+
   void remove_bias(int, double *);
+  void remove_bias_thr(int, double *, double *);
   void remove_bias_all();
   void restore_bias(int, double *);
   void restore_bias_all();
+  void restore_bias_thr(int, double *, double *);
   double memory_usage();
 
  protected:

src/compute_temp_sphere.cpp

@@ -331,6 +331,15 @@ void ComputeTempSphere::remove_bias(int i, double *v)
   tbias->remove_bias(i,v);
 }
 
+/* ----------------------------------------------------------------------
+   remove velocity bias from atom I to leave thermal velocity
+------------------------------------------------------------------------- */
+
+void ComputeTempSphere::remove_bias_thr(int i, double *v, double *b)
+{
+  tbias->remove_bias_thr(i,v,b);
+}
+
 /* ----------------------------------------------------------------------
    add back in velocity bias to atom I removed by remove_bias()
    assume remove_bias() was previously called
@@ -340,3 +349,13 @@ void ComputeTempSphere::restore_bias(int i, double *v)
 {
   tbias->restore_bias(i,v);
 }
+
+/* ----------------------------------------------------------------------
+   add back in velocity bias to atom I removed by remove_bias_thr()
+   assume remove_bias_thr() was previously called with the same buffer b
+------------------------------------------------------------------------- */
+
+void ComputeTempSphere::restore_bias_thr(int i, double *v, double *b)
+{
+  tbias->restore_bias_thr(i,v,b);
+}

src/compute_temp_sphere.h

@@ -34,12 +34,13 @@ class ComputeTempSphere : public Compute {
   void compute_vector();
 
   void remove_bias(int, double *);
+  void remove_bias_thr(int, double *, double *);
   void restore_bias(int, double *);
+  void restore_bias_thr(int, double *, double *);
 
  private:
   int mode;
   double tfactor;
-  double *inertia;
   char *id_bias;
   Compute *tbias;     // ptr to additional bias compute