add update dipole flag to nose-hoover integrator for extended particles

doc/fix_nh.txt

@@ -26,7 +26,7 @@ fix ID group-ID style_name keyword value ... :pre
 ID, group-ID are documented in "fix"_fix.html command :ulb,l
 style_name = {nvt} or {npt} or {nph} :l
 one or more keyword/value pairs may be appended :l
-keyword = {temp} or {iso} or {aniso} or {tri} or {x} or {y} or {z} or {xy} or {yz} or {xz} or {couple} or {tchain} or {pchain} or {mtk} or {tloop} or {ploop} or {nreset} or {drag} or {dilate} or {scalexy} or {scaleyz} or {scalexz} or {flip} or {fixedpoint}
+keyword = {temp} or {iso} or {aniso} or {tri} or {x} or {y} or {z} or {xy} or {yz} or {xz} or {couple} or {tchain} or {pchain} or {mtk} or {tloop} or {ploop} or {nreset} or {drag} or {dilate} or {scalexy} or {scaleyz} or {scalexz} or {flip} or {fixedpoint} or {update}
   {temp} values = Tstart Tstop Tdamp
     Tstart,Tstop = external temperature at start/end of run
     Tdamp = temperature damping parameter (time units)
@@ -56,8 +56,9 @@ keyword = {temp} or {iso} or {aniso} or {tri} or {x} or {y} or {z} or {xy} or {y
   {scalexz} value = {yes} or {no} = scale xz with lz
   {flip} value = {yes} or {no} = allow or disallow box flips when it becomes highly skewed
   {fixedpoint} values = x y z
-    x,y,z = perform barostat dilation/contraction around this point (distance units) :pre
-    
+    x,y,z = perform barostat dilation/contraction around this point (distance units)
+  {update} value = {dipole} update dipole orientation (only for sphere variants) :pre
+
 :ule
 
 [Examples:]
@@ -327,6 +328,12 @@ far. In all cases, the particle trajectories are unaffected by the
 chosen value, except for a time-dependent constant translation of
 positions.
 
+If the {update} keyword is used with the {dipole} value, then the
+orientation of the dipole moment of each particle is also updated
+during the time integration.  This option should be used for models
+where a dipole moment is assigned to extended particles via use of
+the "atom_style hybrid sphere dipole"_atom_style.html command.
+
 :line
 
 NOTE: Using a barostat coupled to tilt dimensions {xy}, {xz}, {yz} can

src/fix_nh.cpp

@@ -78,6 +78,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
   omega_mass_flag = 0;
   etap_mass_flag = 0;
   flipflag = 1;
+  dipole_flag = 0;
 
   // turn on tilt factor scaling, whenever applicable
 
@@ -321,6 +322,11 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
       else if (strcmp(arg[iarg+1],"no") == 0) flipflag = 0;
       else error->all(FLERR,"Illegal fix nvt/npt/nph command");
       iarg += 2;
+    } else if (strcmp(arg[iarg],"update") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
+      if (strcmp(arg[iarg+1],"dipole") == 0) dipole_flag = 1;
+      else error->all(FLERR,"Illegal fix nvt/npt/nph command");
+      iarg += 2;
     } else if (strcmp(arg[iarg],"fixedpoint") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       fixedpoint[0] = force->numeric(FLERR,arg[iarg+1]);
@@ -417,6 +423,13 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
        p_period[0] != p_period[2]))
     error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings");
 
+  if (dipole_flag) {
+    if (!atom->sphere_flag)
+      error->all(FLERR,"Using update dipole flag requires atom style sphere");
+    if (!atom->mu_flag)
+      error->all(FLERR,"Using update dipole flag requires atom attribute mu");
+  }
+
   if ((tstat_flag && t_period <= 0.0) ||
       (p_flag[0] && p_period[0] <= 0.0) ||
       (p_flag[1] && p_period[1] <= 0.0) ||

src/fix_nh.h

@@ -109,6 +109,7 @@ class FixNH : public Fix {
   int eta_mass_flag;               // 1 if eta_mass updated, 0 if not.
   int omega_mass_flag;             // 1 if omega_mass updated, 0 if not.
   int etap_mass_flag;              // 1 if etap_mass updated, 0 if not.
+  int dipole_flag;                 // 1 if dipole is updated, 0 if not.
 
   int scaleyz;                     // 1 if yz scaled with lz
   int scalexz;                     // 1 if xz scaled with lz

src/fix_nh_sphere.cpp

@@ -91,6 +91,42 @@ void FixNHSphere::nve_v()
     }
 }
 
+/* ----------------------------------------------------------------------
+   perform full-step update of position with dipole orientation, if requested
+-----------------------------------------------------------------------*/
+
+void FixNHSphere::nve_x()
+{
+  // standard nve_x position update
+
+  FixNH::nve_x();
+
+  // update mu for dipoles
+  // d_mu/dt = omega cross mu
+  // renormalize mu to dipole length
+
+  if (dipole_flag) {
+    double msq,scale,g[3];
+    double **mu = atom->mu;
+    double **omega = atom->omega;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit)
+        if (mu[i][3] > 0.0) {
+          g[0] = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]);
+          g[1] = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]);
+          g[2] = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]);
+          msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
+          scale = mu[i][3]/sqrt(msq);
+          mu[i][0] = g[0]*scale;
+          mu[i][1] = g[1]*scale;
+          mu[i][2] = g[2]*scale;
+        }
+  }
+}
+
 /* ----------------------------------------------------------------------
    perform half-step scaling of rotatonal velocities
 -----------------------------------------------------------------------*/

src/fix_nh_sphere.h

@@ -26,6 +26,7 @@ class FixNHSphere : public FixNH {
 
  protected:
   void nve_v();
+  void nve_x();
   void nh_v_temp();
 };