Merge pull request #4490 from gsalkuin/fix-move-rotate

Add option to rotate dipoles in fix move rotate or transrot

doc/src/fix_move.rst

@@ -35,11 +35,12 @@ Syntax
          v_vx,v_vy,v_vz = 3 variable names that calculate x,y,z velocity as function of time, any component can be specified as NULL
 
 * zero or more keyword/value pairs may be appended
-* keyword = *units*
+* keyword = *units* or *update*
 
   .. parsed-literal::
 
        *units* value = *box* or *lattice*
+       *update* value = *dipole*
 
 Examples
 """"""""
@@ -49,7 +50,7 @@ Examples
    fix 1 boundary move wiggle 3.0 0.0 0.0 1.0 units box
    fix 2 boundary move rotate 0.0 0.0 0.0 0.0 0.0 1.0 5.0
    fix 2 boundary move variable v_myx v_myy NULL v_VX v_VY NULL
-   fix 3 boundary move transrot 0.1 0.1 0.0 0.0 0.0 0.0 0.0 0.0 1.0 5.0 units box
+   fix 3 boundary move transrot 0.1 0.1 0.0 0.0 0.0 0.0 0.0 0.0 1.0 5.0 units box update dipole
 
 Description
 """""""""""
@@ -217,6 +218,15 @@ been previously used to define the lattice spacing.  Each of these 3
 quantities may be dependent on the x,y,z dimension, since the lattice
 spacings can be different in x,y,z.
 
+.. versionadded:: TBD
+
+If the *update dipole* keyword/value pair is used together with the
+*rotate* or *transrot* style, then the orientation of the dipole moment
+of each particle is also updated appropriately to correspond with the rotation.
+This option should be used for models where a dipole moment is assigned to
+finite-size particles, e.g. spheroids via use of the :doc:`atom_style hybrid
+sphere dipole <atom_style>` command.
+
 ----------
 
 Restart, fix_modify, output, run start/stop, minimize info

src/fix_move.cpp

@@ -49,9 +49,9 @@ static constexpr double INERTIA = 0.2;    // moment of inertia prefactor for ell
 FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
     Fix(lmp, narg, arg), xvarstr(nullptr), yvarstr(nullptr), zvarstr(nullptr), vxvarstr(nullptr),
     vyvarstr(nullptr), vzvarstr(nullptr), xoriginal(nullptr), toriginal(nullptr),
-    qoriginal(nullptr), displace(nullptr), velocity(nullptr)
+    qoriginal(nullptr), muoriginal(nullptr), displace(nullptr), velocity(nullptr)
 {
-  if (narg < 4) error->all(FLERR, "Illegal fix move command");
+  if (narg < 4) utils::missing_cmd_args(FLERR, "fix move", error);
 
   restart_global = 1;
   restart_peratom = 1;
@@ -69,7 +69,7 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
   int iarg = 0;
 
   if (strcmp(arg[3], "linear") == 0) {
-    if (narg < 7) error->all(FLERR, "Illegal fix move command");
+    if (narg < 7) utils::missing_cmd_args(FLERR, "fix move linear", error);
     iarg = 7;
     mstyle = LINEAR;
     if (strcmp(arg[4], "NULL") == 0)
@@ -92,7 +92,7 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
     }
 
   } else if (strcmp(arg[3], "wiggle") == 0) {
-    if (narg < 8) error->all(FLERR, "Illegal fix move command");
+    if (narg < 8) utils::missing_cmd_args(FLERR, "fix move wiggle", error);
     iarg = 8;
     mstyle = WIGGLE;
     if (strcmp(arg[4], "NULL") == 0)
@@ -114,10 +114,10 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
       az = utils::numeric(FLERR, arg[6], false, lmp);
     }
     period = utils::numeric(FLERR, arg[7], false, lmp);
-    if (period <= 0.0) error->all(FLERR, "Illegal fix move command");
+    if (period <= 0.0) error->all(FLERR, 7, "Illegal fix move wiggle period {}", period);
 
   } else if (strcmp(arg[3], "rotate") == 0) {
-    if (narg < 11) error->all(FLERR, "Illegal fix move command");
+    if (narg < 11) utils::missing_cmd_args(FLERR, "fix move rotate", error);
     iarg = 11;
     mstyle = ROTATE;
     point[0] = utils::numeric(FLERR, arg[4], false, lmp);
@@ -127,10 +127,10 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
     axis[1] = utils::numeric(FLERR, arg[8], false, lmp);
     axis[2] = utils::numeric(FLERR, arg[9], false, lmp);
     period = utils::numeric(FLERR, arg[10], false, lmp);
-    if (period <= 0.0) error->all(FLERR, "Illegal fix move command");
+    if (period <= 0.0) error->all(FLERR, 10, "Illegal fix move rotate period {}", period);
 
   } else if (strcmp(arg[3], "transrot") == 0) {
-    if (narg < 11) error->all(FLERR, "Illegal fix move command");
+    if (narg < 11) utils::missing_cmd_args(FLERR, "fix move transrot", error);
     iarg = 14;
     mstyle = TRANSROT;
     vxflag = vyflag = vzflag = 1;
@@ -144,10 +144,10 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
     axis[1] = utils::numeric(FLERR, arg[11], false, lmp);
     axis[2] = utils::numeric(FLERR, arg[12], false, lmp);
     period = utils::numeric(FLERR, arg[13], false, lmp);
-    if (period <= 0.0) error->all(FLERR, "Illegal fix move command");
+    if (period <= 0.0) error->all(FLERR, 13, "Illegal fix move rotate period {}", period);
 
   } else if (strcmp(arg[3], "variable") == 0) {
-    if (narg < 10) error->all(FLERR, "Illegal fix move command");
+    if (narg < 10) utils::missing_cmd_args(FLERR, "fix move variable", error);
     iarg = 10;
     mstyle = VARIABLE;
     if (strcmp(arg[4], "NULL") == 0)
@@ -155,70 +155,83 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
     else if (utils::strmatch(arg[4], "^v_")) {
       xvarstr = utils::strdup(arg[4] + 2);
     } else
-      error->all(FLERR, "Illegal fix move command");
+      error->all(FLERR, 4, "Argument must be either NULL or variable reference");
     if (strcmp(arg[5], "NULL") == 0)
       yvarstr = nullptr;
     else if (utils::strmatch(arg[5], "^v_")) {
       yvarstr = utils::strdup(arg[5] + 2);
     } else
-      error->all(FLERR, "Illegal fix move command");
+      error->all(FLERR, 5, "Argument must be either NULL or variable reference");
     if (strcmp(arg[6], "NULL") == 0)
       zvarstr = nullptr;
     else if (utils::strmatch(arg[6], "^v_")) {
       zvarstr = utils::strdup(arg[6] + 2);
     } else
-      error->all(FLERR, "Illegal fix move command");
+      error->all(FLERR, 6, "Argument must be either NULL or variable reference");
     if (strcmp(arg[7], "NULL") == 0)
       vxvarstr = nullptr;
     else if (utils::strmatch(arg[7], "^v_")) {
       vxvarstr = utils::strdup(arg[7] + 2);
     } else
-      error->all(FLERR, "Illegal fix move command");
+      error->all(FLERR, 7, "Argument must be either NULL or variable reference");
     if (strcmp(arg[8], "NULL") == 0)
       vyvarstr = nullptr;
     else if (utils::strmatch(arg[8], "^v_")) {
       vyvarstr = utils::strdup(arg[8] + 2);
     } else
-      error->all(FLERR, "Illegal fix move command");
+      error->all(FLERR, 8, "Argument must be either NULL or variable reference");
     if (strcmp(arg[9], "NULL") == 0)
       vzvarstr = nullptr;
     else if (utils::strmatch(arg[9], "^v_")) {
       vzvarstr = utils::strdup(arg[9] + 2);
     } else
-      error->all(FLERR, "Illegal fix move command");
+      error->all(FLERR, 9, "Argument must be either NULL or variable reference");
 
-  } else
+  } else {
-    error->all(FLERR, "Illegal fix move command");
+    error->all(FLERR, 3, "Unknown fix move mode {}", arg[3]);
+  }
 
   // optional args
 
   int scaleflag = 1;
+  update_mu_flag = 0;
 
   while (iarg < narg) {
     if (strcmp(arg[iarg], "units") == 0) {
-      if (iarg + 2 > narg) error->all(FLERR, "Illegal fix move command");
+      if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix move units", error);
       if (strcmp(arg[iarg + 1], "box") == 0)
         scaleflag = 0;
       else if (strcmp(arg[iarg + 1], "lattice") == 0)
         scaleflag = 1;
       else
-        error->all(FLERR, "Illegal fix move command");
+        error->all(FLERR, iarg + 1, "Unknown fix move units setting {}", arg[iarg + 1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg], "update") == 0) {
+      if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix move update", error);
+      if (strcmp(arg[iarg + 1], "dipole") == 0) {
+        if ((mstyle != ROTATE) && (mstyle != TRANSROT))
+          error->all(FLERR, 3, "Fix move keyword update dipole requires style rotate or transrot");
+        if (!atom->mu_flag)
+          error->all(FLERR, iarg + 1, "Fix move keyword update dipole requires atom attribute mu");
+        update_mu_flag = 1;
+      } else
+        error->all(FLERR, iarg + 1, "Unkown fix move update argument {}", arg[iarg + 1]);
       iarg += 2;
     } else
-      error->all(FLERR, "Illegal fix move command");
+      error->all(FLERR, iarg, "Unknown fix move keyword {}", arg[iarg]);
   }
 
   // error checks and warnings
 
   if (domain->dimension == 2) {
     if (((mstyle == LINEAR) || (mstyle == TRANSROT)) && vzflag && (vz != 0.0))
-      error->all(FLERR, "Fix move cannot set linear z motion for 2d problem");
+      error->all(FLERR, 3, "Fix move cannot set linear z motion for 2d problem");
     if (mstyle == WIGGLE && azflag && az != 0.0)
-      error->all(FLERR, "Fix move cannot set wiggle z motion for 2d problem");
+      error->all(FLERR, 3, "Fix move cannot set wiggle z motion for 2d problem");
     if (((mstyle == ROTATE) || (mstyle == TRANSROT)) && (axis[0] != 0.0 || axis[1] != 0.0))
-      error->all(FLERR, "Fix move cannot rotate around non z-axis for 2d problem");
+      error->all(FLERR, 3, "Fix move cannot rotate around non z-axis for 2d problem");
     if (mstyle == VARIABLE && (zvarstr || vzvarstr))
-      error->all(FLERR, "Fix move cannot define z or vz variable for 2d problem");
+      error->all(FLERR, 3, "Fix move cannot define z or vz variable for 2d problem");
   }
 
   // setup scaling and apply scaling factors to velocity & amplitude
@@ -283,7 +296,7 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
   if (ellipsoid_flag || tri_flag || body_flag || quat_atom_flag) quat_flag = 1;
 
   extra_flag = 0;
-  if (omega_flag || angmom_flag || theta_flag || quat_flag) extra_flag = 1;
+  if (omega_flag || angmom_flag || theta_flag || quat_flag || update_mu_flag) extra_flag = 1;
 
   // perform initial allocation of atom-based array
   // register with Atom class
@@ -304,6 +317,7 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
   // xoriginal = initial unwrapped positions of atoms
   // toriginal = initial theta of lines
   // qoriginal = initial quat of extended particles
+  // muoriginal = initial dipole moment of extended particles
 
   double **x = atom->x;
   imageint *image = atom->image;
@@ -362,6 +376,18 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
     }
   }
 
+  if (update_mu_flag) {
+    double **mu = atom->mu;
+    for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+        muoriginal[i][0] = mu[i][0];
+        muoriginal[i][1] = mu[i][1];
+        muoriginal[i][2] = mu[i][2];
+        muoriginal[i][3] = mu[i][3];
+      }
+    }
+  }
+
   // nrestart = size of per-atom restart data
   // nrestart = 1 + xorig + torig + qorig
 
@@ -388,6 +414,7 @@ FixMove::~FixMove()
   memory->destroy(xoriginal);
   memory->destroy(toriginal);
   memory->destroy(qoriginal);
+  memory->destroy(muoriginal);
   memory->destroy(displace);
   memory->destroy(velocity);
 
@@ -426,63 +453,83 @@ void FixMove::init()
   if (mstyle == VARIABLE) {
     if (xvarstr) {
       xvar = input->variable->find(xvarstr);
-      if (xvar < 0) error->all(FLERR, "Variable name for fix move does not exist");
+      if (xvar < 0)
+        error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix move does not exist",
+                   xvarstr);
       if (input->variable->equalstyle(xvar))
         xvarstyle = EQUAL;
       else if (input->variable->atomstyle(xvar))
         xvarstyle = ATOM;
       else
-        error->all(FLERR, "Variable for fix move is invalid style");
+        error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix move is invalid style",
+                   xvarstr);
     }
     if (yvarstr) {
       yvar = input->variable->find(yvarstr);
-      if (yvar < 0) error->all(FLERR, "Variable name for fix move does not exist");
+      if (yvar < 0)
+        error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix move does not exist",
+                   yvarstr);
+
       if (input->variable->equalstyle(yvar))
         yvarstyle = EQUAL;
       else if (input->variable->atomstyle(yvar))
         yvarstyle = ATOM;
       else
-        error->all(FLERR, "Variable for fix move is invalid style");
+        error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix move is invalid style",
+                   yvarstr);
     }
     if (zvarstr) {
       zvar = input->variable->find(zvarstr);
-      if (zvar < 0) error->all(FLERR, "Variable name for fix move does not exist");
+      if (zvar < 0)
+        error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix move does not exist",
+                   zvarstr);
+
       if (input->variable->equalstyle(zvar))
         zvarstyle = EQUAL;
       else if (input->variable->atomstyle(zvar))
         zvarstyle = ATOM;
       else
-        error->all(FLERR, "Variable for fix move is invalid style");
+        error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix move is invalid style",
+                   zvarstr);
     }
     if (vxvarstr) {
       vxvar = input->variable->find(vxvarstr);
-      if (vxvar < 0) error->all(FLERR, "Variable name for fix move does not exist");
+      if (vxvar < 0)
+        error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix move does not exist",
+                   vxvarstr);
       if (input->variable->equalstyle(vxvar))
         vxvarstyle = EQUAL;
       else if (input->variable->atomstyle(vxvar))
         vxvarstyle = ATOM;
       else
-        error->all(FLERR, "Variable for fix move is invalid style");
+        error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix move is invalid style",
+                   vxvarstr);
     }
     if (vyvarstr) {
       vyvar = input->variable->find(vyvarstr);
-      if (vyvar < 0) error->all(FLERR, "Variable name for fix move does not exist");
+      if (vyvar < 0)
+        error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix move does not exist",
+                   vyvarstr);
       if (input->variable->equalstyle(vyvar))
         vyvarstyle = EQUAL;
       else if (input->variable->atomstyle(vyvar))
         vyvarstyle = ATOM;
       else
-        error->all(FLERR, "Variable for fix move is invalid style");
+        error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix move is invalid style",
+                   vyvarstr);
     }
     if (vzvarstr) {
       vzvar = input->variable->find(vzvarstr);
-      if (vzvar < 0) error->all(FLERR, "Variable name for fix move does not exist");
+      if (vzvar < 0)
+        error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix move does not exist",
+                   vzvarstr);
       if (input->variable->equalstyle(vzvar))
         vzvarstyle = EQUAL;
       else if (input->variable->atomstyle(vzvar))
         vzvarstyle = ATOM;
       else
-        error->all(FLERR, "Variable for fix move is invalid style");
+        error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix move is invalid style",
+                   vzvarstr);
     }
 
     if (xvarstr && xvarstyle == ATOM) displaceflag = 1;
@@ -533,6 +580,7 @@ void FixMove::initial_integrate(int /*vflag*/)
   double *rmass = atom->rmass;
   double *mass = atom->mass;
   double **quat_atom = atom->quat;
+  double **mu = atom->mu;
   int *type = atom->type;
   int *ellipsoid = atom->ellipsoid;
   int *line = atom->line;
@@ -676,6 +724,8 @@ void FixMove::initial_integrate(int /*vflag*/)
     qrotate[2] = runit[1] * qsine;
     qrotate[3] = runit[2] * qsine;
 
+    double rotmat[3][3], g[3];    // for dipole rotation
+
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
         xold[0] = x[i][0];
@@ -775,6 +825,15 @@ void FixMove::initial_integrate(int /*vflag*/)
           } else if (quat_atom_flag) {
             MathExtra::quatquat(qrotate, qoriginal[i], quat_atom[i]);
           }
+
+          if (update_mu_flag && mu[i][3] > 0.0) {
+            MathExtra::quat_to_mat(qrotate, rotmat);
+            MathExtra::matvec(rotmat, muoriginal[i], g);
+            mu[i][0] = g[0];
+            mu[i][1] = g[1];
+            mu[i][2] = g[2];
+            // mu[i][3] = sqrt(g[0] * g[0] + g[1] * g[1] + g[2] * g[2]);
+          }
         }
 
         domain->remap_near(x[i], xold);
@@ -809,6 +868,8 @@ void FixMove::initial_integrate(int /*vflag*/)
     qrotate[2] = runit[1] * qsine;
     qrotate[3] = runit[2] * qsine;
 
+    double rotmat[3][3], g[3];    // for dipole rotation
+
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
         xold[0] = x[i][0];
@@ -908,6 +969,15 @@ void FixMove::initial_integrate(int /*vflag*/)
           } else if (quat_atom_flag) {
             MathExtra::quatquat(qrotate, qoriginal[i], quat_atom[i]);
           }
+
+          if (update_mu_flag && mu[i][3] > 0.0) {
+            MathExtra::quat_to_mat(qrotate, rotmat);
+            MathExtra::matvec(rotmat, muoriginal[i], g);
+            mu[i][0] = g[0];
+            mu[i][1] = g[1];
+            mu[i][2] = g[2];
+            // mu[i][3] = sqrt(g[0] * g[0] + g[1] * g[1] + g[2] * g[2]);
+          }
         }
 
         v[i][0] += vx;
@@ -1248,6 +1318,7 @@ void FixMove::grow_arrays(int nmax)
   memory->grow(xoriginal, nmax, 3, "move:xoriginal");
   if (theta_flag) memory->grow(toriginal, nmax, "move:toriginal");
   if (quat_flag) memory->grow(qoriginal, nmax, 4, "move:qoriginal");
+  if (update_mu_flag) memory->grow(muoriginal, nmax, 4, "move:muoriginal");
   array_atom = xoriginal;
 }
 
@@ -1267,6 +1338,12 @@ void FixMove::copy_arrays(int i, int j, int /*delflag*/)
     qoriginal[j][2] = qoriginal[i][2];
     qoriginal[j][3] = qoriginal[i][3];
   }
+  if (update_mu_flag) {
+    muoriginal[j][0] = muoriginal[i][0];
+    muoriginal[j][1] = muoriginal[i][1];
+    muoriginal[j][2] = muoriginal[i][2];
+    muoriginal[j][3] = muoriginal[i][3];
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -1302,7 +1379,8 @@ void FixMove::set_arrays(int i)
 
   // backup particle to time_origin
 
-  if (mstyle == VARIABLE) error->all(FLERR, "Cannot add atoms to fix move variable");
+  if (mstyle == VARIABLE)
+    error->all(FLERR, Error::NOLASTLINE, "Cannot add atoms with fix move style variable");
 
   domain->unmap(x[i], image[i], xoriginal[i]);
   double delta = (update->ntimestep - time_origin) * update->dt;
@@ -1461,6 +1539,12 @@ int FixMove::pack_exchange(int i, double *buf)
     buf[n++] = qoriginal[i][2];
     buf[n++] = qoriginal[i][3];
   }
+  if (update_mu_flag) {
+    buf[n++] = muoriginal[i][0];
+    buf[n++] = muoriginal[i][1];
+    buf[n++] = muoriginal[i][2];
+    buf[n++] = muoriginal[i][3];
+  }
   return n;
 }
 
@@ -1481,6 +1565,12 @@ int FixMove::unpack_exchange(int nlocal, double *buf)
     qoriginal[nlocal][2] = buf[n++];
     qoriginal[nlocal][3] = buf[n++];
   }
+  if (update_mu_flag) {
+    muoriginal[nlocal][0] = buf[n++];
+    muoriginal[nlocal][1] = buf[n++];
+    muoriginal[nlocal][2] = buf[n++];
+    muoriginal[nlocal][3] = buf[n++];
+  }
   return n;
 }
 
@@ -1501,6 +1591,12 @@ int FixMove::pack_restart(int i, double *buf)
     buf[n++] = qoriginal[i][2];
     buf[n++] = qoriginal[i][3];
   }
+  if (update_mu_flag) {
+    buf[n++] = muoriginal[i][0];
+    buf[n++] = muoriginal[i][1];
+    buf[n++] = muoriginal[i][2];
+    buf[n++] = muoriginal[i][3];
+  }
   // pack buf[0] this way because other fixes unpack it
   buf[0] = n;
   return n;
@@ -1531,6 +1627,12 @@ void FixMove::unpack_restart(int nlocal, int nth)
     qoriginal[nlocal][2] = extra[nlocal][m++];
     qoriginal[nlocal][3] = extra[nlocal][m++];
   }
+  if (update_mu_flag) {
+    muoriginal[nlocal][0] = extra[nlocal][m++];
+    muoriginal[nlocal][1] = extra[nlocal][m++];
+    muoriginal[nlocal][2] = extra[nlocal][m++];
+    muoriginal[nlocal][3] = extra[nlocal][m++];
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -1555,5 +1657,5 @@ int FixMove::size_restart(int /*nlocal*/)
 
 void FixMove::reset_dt()
 {
-  error->all(FLERR, "Resetting timestep size is not allowed with fix move");
+  error->all(FLERR, Error::NOLASTLINE, "Resetting timestep size is not allowed with fix move");
 }

src/fix_move.h

@@ -62,13 +62,14 @@ class FixMove : public Fix {
   int xvarstyle, yvarstyle, zvarstyle, vxvarstyle, vyvarstyle, vzvarstyle;
   int extra_flag, omega_flag, angmom_flag;
   int radius_flag, ellipsoid_flag, line_flag, tri_flag, body_flag, quat_atom_flag;
-  int theta_flag, quat_flag;
+  int theta_flag, quat_flag, update_mu_flag;
   int nlevels_respa, nrestart;
   int time_origin;
 
   double **xoriginal;    // original coords of atoms
   double *toriginal;     // original theta of atoms
   double **qoriginal;    // original quat of atoms
+  double **muoriginal;    // original mu of atoms
   int displaceflag, velocityflag;
   int maxatom;
   double **displace, **velocity;