diff --git a/doc/src/fix_move.rst b/doc/src/fix_move.rst index 5d3a9de47f..ff1f8403df 100644 --- a/doc/src/fix_move.rst +++ b/doc/src/fix_move.rst @@ -12,7 +12,7 @@ Syntax * ID, group-ID are documented in :doc:`fix ` command * move = style name of this fix command -* style = *linear* or *wiggle* or *rotate* or *variable* +* style = *linear* or *wiggle* or *rotate* or *transrot* or *variable* .. parsed-literal:: @@ -25,6 +25,11 @@ Syntax Px,Py,Pz = origin point of axis of rotation (distance units) Rx,Ry,Rz = axis of rotation vector period = period of rotation (time units) + *transrot* args = Vx Vy Vz Px Py Pz Rx Ry Rz period + Vx,Vy,Vz = components of velocity vector (velocity units) + Px,Py,Pz = origin point of axis of rotation (distance units) + Rx,Ry,Rz = axis of rotation vector + period = period of rotation (time units) *variable* args = v_dx v_dy v_dz v_vx v_vy v_vz v_dx,v_dy,v_dz = 3 variable names that calculate x,y,z displacement as function of time, any component can be specified as NULL 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 @@ -44,6 +49,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 Description """"""""""" @@ -55,15 +61,17 @@ whose movement can influence nearby atoms. .. note:: - The atoms affected by this fix should not normally be time - integrated by other fixes (e.g. :doc:`fix nve `, :doc:`fix nvt `), since that will change their positions and - velocities twice. + The atoms affected by this fix should not normally be time integrated + by other fixes (e.g. :doc:`fix nve `, :doc:`fix nvt + `), since that will change their positions and velocities + twice. .. note:: As atoms move due to this fix, they will pass through periodic boundaries and be remapped to the other side of the simulation box, - just as they would during normal time integration (e.g. via the :doc:`fix nve ` command). It is up to you to decide whether + just as they would during normal time integration (e.g. via the + :doc:`fix nve ` command). It is up to you to decide whether periodic boundaries are appropriate with the kind of atom motion you are prescribing with this fix. @@ -73,11 +81,11 @@ whose movement can influence nearby atoms. position at the time the fix is specified. These initial coordinates are stored by the fix in "unwrapped" form, by using the image flags associated with each atom. See the :doc:`dump custom ` command - for a discussion of "unwrapped" coordinates. See the Atoms section of - the :doc:`read_data ` command for a discussion of image flags - and how they are set for each atom. You can reset the image flags - (e.g. to 0) before invoking this fix by using the :doc:`set image ` - command. + for a discussion of "unwrapped" coordinates. See the Atoms section + of the :doc:`read_data ` command for a discussion of image + flags and how they are set for each atom. You can reset the image + flags (e.g. to 0) before invoking this fix by using the :doc:`set + image ` command. ---------- @@ -118,13 +126,13 @@ notation as where *X0* = (x0,y0,z0) is their position at the time the fix is specified, *A* is the specified amplitude vector with components -(Ax,Ay,Az), *omega* is 2 PI / *period*, and *delta* is the time -elapsed since the fix was specified. This style also sets the -velocity of each atom to the time derivative of this expression. If -any of the amplitude components is specified as NULL, then the -position and velocity of that component is time integrated the same as -the :doc:`fix nve ` command would perform, using the -corresponding force component on the atom. +(Ax,Ay,Az), *omega* is 2 PI / *period*, and *delta* is the time elapsed +since the fix was specified. This style also sets the velocity of each +atom to the time derivative of this expression. If any of the amplitude +components is specified as NULL, then the position and velocity of that +component is time integrated the same as the :doc:`fix nve ` +command would perform, using the corresponding force component on the +atom. Note that the *wiggle* style is identical to using the *variable* style with :doc:`equal-style variables ` that use the @@ -139,21 +147,29 @@ swiggle() and cwiggle() functions. E.g. variable v equal v_omega*($A-cwiggle(0.0,$A,$T)) fix 1 boundary move variable v_x NULL NULL v_v NULL NULL -The *rotate* style rotates atoms around a rotation axis *R* = -(Rx,Ry,Rz) that goes through a point *P* = (Px,Py,Pz). The *period* of -the rotation is also specified. The direction of rotation for the -atoms around the rotation axis is consistent with the right-hand rule: -if your right-hand thumb points along *R*, then your fingers wrap -around the axis in the direction of rotation. +The *rotate* style rotates atoms around a rotation axis *R* = (Rx,Ry,Rz) +that goes through a point *P* = (Px,Py,Pz). The *period* of the +rotation is also specified. The direction of rotation for the atoms +around the rotation axis is consistent with the right-hand rule: if your +right-hand thumb points along *R*, then your fingers wrap around the +axis in the direction of rotation. This style also sets the velocity of each atom to (omega cross Rperp) where omega is its angular velocity around the rotation axis and Rperp is a perpendicular vector from the rotation axis to the atom. If the defined :doc:`atom_style ` assigns an angular velocity or -angular momentum or orientation to each atom (:doc:`atom styles ` sphere, ellipsoid, line, tri, body), then +angular momentum or orientation to each atom (:doc:`atom styles +` sphere, ellipsoid, line, tri, body), then those properties are also updated appropriately to correspond to the atom's motion and rotation over time. +The *transrot* style combines the effects of *rotate* and *linear* so +that it is possible to prescribe a rotating group of atoms that also +moves at a constant velocity. The arguments are for the translation +first and then for the rotation. Since the rotation affects all +coordinate components, it is not possible to set any of the +translation vector components to NULL. + The *variable* style allows the position and velocity components of each atom to be set by formulas specified via the :doc:`variable ` command. Each of the 6 variables is @@ -165,10 +181,10 @@ Each variable must be of either the *equal* or *atom* style. a function of the timestep as well as of other simulation values. *Atom*\ -style variables compute a numeric quantity for each atom, that can be a function per-atom quantities, such as the atom's position, as -well as of the timestep and other simulation values. Note that this -fix stores the original coordinates of each atom (see note below) so -that per-atom quantity can be used in an atom-style variable formula. -See the :doc:`variable ` command for details. +well as of the timestep and other simulation values. Note that this fix +stores the original coordinates of each atom (see note below) so that +per-atom quantity can be used in an atom-style variable formula. See +the :doc:`variable ` command for details. The first 3 variables (v_dx,v_dy,v_dz) specified for the *variable* style are used to calculate a displacement from the atom's original @@ -206,8 +222,9 @@ spacings can be different in x,y,z. Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" -This fix writes the original coordinates of moving atoms to :doc:`binary restart files `, as well as the initial timestep, so that -the motion can be continuous in a restarted simulation. See the +This fix writes the original coordinates of moving atoms to :doc:`binary +restart files `, as well as the initial timestep, so that the +motion can be continuous in a restarted simulation. See the :doc:`read_restart ` command for info on how to re-specify a fix in an input script that reads a restart file, so that the operation of the fix continues in an uninterrupted fashion. @@ -224,11 +241,12 @@ fix. This fix produces a per-atom array which can be accessed by various :doc:`output commands `. The number of columns for each -atom is 3, and the columns store the original unwrapped x,y,z coords -of each atom. The per-atom values can be accessed on any timestep. +atom is 3, and the columns store the original unwrapped x,y,z coords of +each atom. The per-atom values can be accessed on any timestep. No parameter of this fix can be used with the *start/stop* keywords of -the :doc:`run ` command. This fix is not invoked during :doc:`energy minimization `. +the :doc:`run ` command. This fix is not invoked during +:doc:`energy minimization `. For :doc:`rRESPA time integration `, this fix adjusts the position and velocity of atoms on the outermost rRESPA level. diff --git a/src/fix_move.cpp b/src/fix_move.cpp index f7bc4d3640..756292d06b 100644 --- a/src/fix_move.cpp +++ b/src/fix_move.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -40,21 +39,19 @@ using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; -enum{LINEAR,WIGGLE,ROTATE,VARIABLE}; -enum{EQUAL,ATOM}; +enum { LINEAR, WIGGLE, ROTATE, VARIABLE, TRANSROT }; +enum { EQUAL, ATOM }; -#define INERTIA 0.2 // moment of inertia prefactor for ellipsoid +#define INERTIA 0.2 // moment of inertia prefactor for ellipsoid /* ---------------------------------------------------------------------- */ 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) + 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) { - if (narg < 4) error->all(FLERR,"Illegal fix move command"); + if (narg < 4) error->all(FLERR, "Illegal fix move command"); restart_global = 1; restart_peratom = 1; @@ -71,125 +68,162 @@ 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 (strcmp(arg[3], "linear") == 0) { + if (narg < 7) error->all(FLERR, "Illegal fix move command"); iarg = 7; mstyle = LINEAR; - if (strcmp(arg[4],"NULL") == 0) vxflag = 0; + if (strcmp(arg[4], "NULL") == 0) + vxflag = 0; else { vxflag = 1; - vx = utils::numeric(FLERR,arg[4],false,lmp); + vx = utils::numeric(FLERR, arg[4], false, lmp); } - if (strcmp(arg[5],"NULL") == 0) vyflag = 0; + if (strcmp(arg[5], "NULL") == 0) + vyflag = 0; else { vyflag = 1; - vy = utils::numeric(FLERR,arg[5],false,lmp); + vy = utils::numeric(FLERR, arg[5], false, lmp); } - if (strcmp(arg[6],"NULL") == 0) vzflag = 0; + if (strcmp(arg[6], "NULL") == 0) + vzflag = 0; else { vzflag = 1; - vz = utils::numeric(FLERR,arg[6],false,lmp); + vz = utils::numeric(FLERR, arg[6], false, lmp); } - } else if (strcmp(arg[3],"wiggle") == 0) { - if (narg < 8) error->all(FLERR,"Illegal fix move command"); + } else if (strcmp(arg[3], "wiggle") == 0) { + if (narg < 8) error->all(FLERR, "Illegal fix move command"); iarg = 8; mstyle = WIGGLE; - if (strcmp(arg[4],"NULL") == 0) axflag = 0; + if (strcmp(arg[4], "NULL") == 0) + axflag = 0; else { axflag = 1; - ax = utils::numeric(FLERR,arg[4],false,lmp); + ax = utils::numeric(FLERR, arg[4], false, lmp); } - if (strcmp(arg[5],"NULL") == 0) ayflag = 0; + if (strcmp(arg[5], "NULL") == 0) + ayflag = 0; else { ayflag = 1; - ay = utils::numeric(FLERR,arg[5],false,lmp); + ay = utils::numeric(FLERR, arg[5], false, lmp); } - if (strcmp(arg[6],"NULL") == 0) azflag = 0; + if (strcmp(arg[6], "NULL") == 0) + azflag = 0; else { azflag = 1; - az = utils::numeric(FLERR,arg[6],false,lmp); + 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"); + period = utils::numeric(FLERR, arg[7], false, lmp); + if (period <= 0.0) error->all(FLERR, "Illegal fix move command"); - } else if (strcmp(arg[3],"rotate") == 0) { - if (narg < 11) error->all(FLERR,"Illegal fix move command"); + } else if (strcmp(arg[3], "rotate") == 0) { + if (narg < 11) error->all(FLERR, "Illegal fix move command"); iarg = 11; mstyle = ROTATE; - point[0] = utils::numeric(FLERR,arg[4],false,lmp); - point[1] = utils::numeric(FLERR,arg[5],false,lmp); - point[2] = utils::numeric(FLERR,arg[6],false,lmp); - axis[0] = utils::numeric(FLERR,arg[7],false,lmp); - 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"); + point[0] = utils::numeric(FLERR, arg[4], false, lmp); + point[1] = utils::numeric(FLERR, arg[5], false, lmp); + point[2] = utils::numeric(FLERR, arg[6], false, lmp); + axis[0] = utils::numeric(FLERR, arg[7], false, lmp); + 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"); - } else if (strcmp(arg[3],"variable") == 0) { - if (narg < 10) error->all(FLERR,"Illegal fix move command"); + } else if (strcmp(arg[3], "transrot") == 0) { + if (narg < 11) error->all(FLERR, "Illegal fix move command"); + iarg = 14; + mstyle = TRANSROT; + vxflag = vyflag = vzflag = 1; + vx = utils::numeric(FLERR, arg[4], false, lmp); + vy = utils::numeric(FLERR, arg[5], false, lmp); + vz = utils::numeric(FLERR, arg[6], false, lmp); + point[0] = utils::numeric(FLERR, arg[7], false, lmp); + point[1] = utils::numeric(FLERR, arg[8], false, lmp); + point[2] = utils::numeric(FLERR, arg[9], false, lmp); + axis[0] = utils::numeric(FLERR, arg[10], false, lmp); + 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"); + + } else if (strcmp(arg[3], "variable") == 0) { + if (narg < 10) error->all(FLERR, "Illegal fix move command"); iarg = 10; mstyle = VARIABLE; - if (strcmp(arg[4],"NULL") == 0) xvarstr = nullptr; - else if (utils::strmatch(arg[4],"^v_")) { - xvarstr = utils::strdup(arg[4]+2); - } else error->all(FLERR,"Illegal fix move command"); - 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"); - 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"); - 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"); - 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"); - 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"); + if (strcmp(arg[4], "NULL") == 0) + xvarstr = nullptr; + else if (utils::strmatch(arg[4], "^v_")) { + xvarstr = utils::strdup(arg[4] + 2); + } else + error->all(FLERR, "Illegal fix move command"); + 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"); + 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"); + 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"); + 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"); + 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"); - } else error->all(FLERR,"Illegal fix move command"); + } else + error->all(FLERR, "Illegal fix move command"); // optional args int scaleflag = 1; while (iarg < narg) { - if (strcmp(arg[iarg],"units") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix move command"); - 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"); + if (strcmp(arg[iarg], "units") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix move command"); + 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"); iarg += 2; - } else error->all(FLERR,"Illegal fix move command"); + } else + error->all(FLERR, "Illegal fix move command"); } // error checks and warnings if (domain->dimension == 2) { - if (mstyle == LINEAR && vzflag && vz != 0.0) - error->all(FLERR,"Fix move cannot set linear z motion for 2d problem"); + if (((mstyle == LINEAR) || (mstyle == TRANSROT)) && vzflag && (vz != 0.0)) + error->all(FLERR, "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"); - if (mstyle == ROTATE && (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, "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"); if (mstyle == VARIABLE && (zvarstr || vzvarstr)) - error->all(FLERR, - "Fix move cannot define z or vz variable for 2d problem"); + error->all(FLERR, "Fix move cannot define z or vz variable for 2d problem"); } // setup scaling and apply scaling factors to velocity & amplitude - if ((mstyle == LINEAR || mstyle == WIGGLE || mstyle == ROTATE) && - scaleflag) { + if ((mstyle != VARIABLE) && scaleflag) { double xscale = domain->lattice->xlattice; double yscale = domain->lattice->ylattice; double zscale = domain->lattice->zlattice; @@ -206,22 +240,29 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) : point[0] *= xscale; point[1] *= yscale; point[2] *= zscale; + } else if (mstyle == TRANSROT) { + vx *= xscale; + vy *= yscale; + vz *= zscale; + point[0] *= xscale; + point[1] *= yscale; + point[2] *= zscale; } } // set omega_rotate from period - if (mstyle == WIGGLE || mstyle == ROTATE) omega_rotate = MY_2PI / period; + if ((mstyle == WIGGLE) || (mstyle == ROTATE) || (mstyle == TRANSROT)) + omega_rotate = MY_2PI / period; // runit = unit vector along rotation axis - if (mstyle == ROTATE) { - double len = sqrt(axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2]); - if (len == 0.0) - error->all(FLERR,"Zero length rotation vector with fix move"); - runit[0] = axis[0]/len; - runit[1] = axis[1]/len; - runit[2] = axis[2]/len; + if ((mstyle == ROTATE) || (mstyle == TRANSROT)) { + double len = sqrt(axis[0] * axis[0] + axis[1] * axis[1] + axis[2] * axis[2]); + if (len == 0.0) error->all(FLERR, "Zero length rotation vector with fix move"); + runit[0] = axis[0] / len; + runit[1] = axis[1] / len; + runit[2] = axis[2] / len; } // set flags for extra attributes particles may store @@ -273,15 +314,18 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) : int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]); - else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0; + if (mask[i] & groupbit) + domain->unmap(x[i], image[i], xoriginal[i]); + else + xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0; } if (theta_flag) { for (int i = 0; i < nlocal; i++) { if ((mask[i] & groupbit) && line[i] >= 0) toriginal[i] = avec_line->bonus[line[i]].theta; - else toriginal[i] = 0.0; + else + toriginal[i] = 0.0; } } @@ -302,8 +346,8 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) : qoriginal[i][1] = quat[1]; qoriginal[i][2] = quat[2]; qoriginal[i][3] = quat[3]; - } else qoriginal[i][0] = qoriginal[i][1] = - qoriginal[i][2] = qoriginal[i][3] = 0.0; + } else + qoriginal[i][0] = qoriginal[i][1] = qoriginal[i][2] = qoriginal[i][3] = 0.0; } } @@ -325,8 +369,8 @@ FixMove::~FixMove() { // unregister callbacks to this fix from Atom class - atom->delete_callback(id,Atom::GROW); - atom->delete_callback(id,Atom::RESTART); + atom->delete_callback(id, Atom::GROW); + atom->delete_callback(id, Atom::RESTART); // delete locally stored arrays @@ -336,12 +380,12 @@ FixMove::~FixMove() memory->destroy(displace); memory->destroy(velocity); - delete [] xvarstr; - delete [] yvarstr; - delete [] zvarstr; - delete [] vxvarstr; - delete [] vyvarstr; - delete [] vzvarstr; + delete[] xvarstr; + delete[] yvarstr; + delete[] zvarstr; + delete[] vxvarstr; + delete[] vyvarstr; + delete[] vzvarstr; } /* ---------------------------------------------------------------------- */ @@ -371,51 +415,63 @@ 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 (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"); + if (xvar < 0) error->all(FLERR, "Variable name for fix move does not exist"); + 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"); } if (yvarstr) { yvar = input->variable->find(yvarstr); - if (yvar < 0) error->all(FLERR, - "Variable name for fix move does not exist"); - 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"); + if (yvar < 0) error->all(FLERR, "Variable name for fix move does not exist"); + 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"); } if (zvarstr) { zvar = input->variable->find(zvarstr); - if (zvar < 0) error->all(FLERR, - "Variable name for fix move does not exist"); - 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"); + if (zvar < 0) error->all(FLERR, "Variable name for fix move does not exist"); + 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"); } if (vxvarstr) { vxvar = input->variable->find(vxvarstr); - if (vxvar < 0) error->all(FLERR, - "Variable name for fix move does not exist"); - 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"); + if (vxvar < 0) error->all(FLERR, "Variable name for fix move does not exist"); + 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"); } if (vyvarstr) { vyvar = input->variable->find(vyvarstr); - if (vyvar < 0) error->all(FLERR, - "Variable name for fix move does not exist"); - 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"); + if (vyvar < 0) error->all(FLERR, "Variable name for fix move does not exist"); + 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"); } if (vzvarstr) { vzvar = input->variable->find(vzvarstr); - if (vzvar < 0) error->all(FLERR, - "Variable name for fix move does not exist"); - 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"); + if (vzvar < 0) error->all(FLERR, "Variable name for fix move does not exist"); + 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"); } if (xvarstr && xvarstyle == ATOM) displaceflag = 1; @@ -429,12 +485,16 @@ void FixMove::init() maxatom = atom->nmax; memory->destroy(displace); memory->destroy(velocity); - if (displaceflag) memory->create(displace,maxatom,3,"move:displace"); - else displace = nullptr; - if (velocityflag) memory->create(velocity,maxatom,3,"move:velocity"); - else velocity = nullptr; + if (displaceflag) + memory->create(displace, maxatom, 3, "move:displace"); + else + displace = nullptr; + if (velocityflag) + memory->create(velocity, maxatom, 3, "move:velocity"); + else + velocity = nullptr; - if (utils::strmatch(update->integrate_style,"^respa")) + if (utils::strmatch(update->integrate_style, "^respa")) nlevels_respa = ((Respa *) update->integrate)->nlevels; } @@ -445,11 +505,11 @@ void FixMove::init() void FixMove::initial_integrate(int /*vflag*/) { int flag; - double ddotr,dx,dy,dz; - double dtfm,theta_new; - double xold[3],a[3],b[3],c[3],d[3],disp[3],w[3],ex[3],ey[3],ez[3]; - double inertia_ellipsoid[3],qrotate[4]; - double *quat,*inertia,*shape; + double ddotr, dx, dy, dz; + double dtfm, theta_new; + double xold[3], a[3], b[3], c[3], d[3], disp[3], w[3], ex[3], ey[3], ez[3]; + double inertia_ellipsoid[3], qrotate[4]; + double *quat, *inertia, *shape; double delta = (update->ntimestep - time_origin) * dt; @@ -481,7 +541,7 @@ void FixMove::initial_integrate(int /*vflag*/) if (vxflag) { v[i][0] = vx; - x[i][0] = xoriginal[i][0] + vx*delta; + x[i][0] = xoriginal[i][0] + vx * delta; } else if (rmass) { dtfm = dtf / rmass[i]; v[i][0] += dtfm * f[i][0]; @@ -494,7 +554,7 @@ void FixMove::initial_integrate(int /*vflag*/) if (vyflag) { v[i][1] = vy; - x[i][1] = xoriginal[i][1] + vy*delta; + x[i][1] = xoriginal[i][1] + vy * delta; } else if (rmass) { dtfm = dtf / rmass[i]; v[i][1] += dtfm * f[i][1]; @@ -507,7 +567,7 @@ void FixMove::initial_integrate(int /*vflag*/) if (vzflag) { v[i][2] = vz; - x[i][2] = xoriginal[i][2] + vz*delta; + x[i][2] = xoriginal[i][2] + vz * delta; } else if (rmass) { dtfm = dtf / rmass[i]; v[i][2] += dtfm * f[i][2]; @@ -518,7 +578,7 @@ void FixMove::initial_integrate(int /*vflag*/) x[i][2] += dtv * v[i][2]; } - domain->remap_near(x[i],xold); + domain->remap_near(x[i], xold); } } @@ -536,8 +596,8 @@ void FixMove::initial_integrate(int /*vflag*/) xold[2] = x[i][2]; if (axflag) { - v[i][0] = ax*omega_rotate*cosine; - x[i][0] = xoriginal[i][0] + ax*sine; + v[i][0] = ax * omega_rotate * cosine; + x[i][0] = xoriginal[i][0] + ax * sine; } else if (rmass) { dtfm = dtf / rmass[i]; v[i][0] += dtfm * f[i][0]; @@ -549,8 +609,8 @@ void FixMove::initial_integrate(int /*vflag*/) } if (ayflag) { - v[i][1] = ay*omega_rotate*cosine; - x[i][1] = xoriginal[i][1] + ay*sine; + v[i][1] = ay * omega_rotate * cosine; + x[i][1] = xoriginal[i][1] + ay * sine; } else if (rmass) { dtfm = dtf / rmass[i]; v[i][1] += dtfm * f[i][1]; @@ -562,8 +622,8 @@ void FixMove::initial_integrate(int /*vflag*/) } if (azflag) { - v[i][2] = az*omega_rotate*cosine; - x[i][2] = xoriginal[i][2] + az*sine; + v[i][2] = az * omega_rotate * cosine; + x[i][2] = xoriginal[i][2] + az * sine; } else if (rmass) { dtfm = dtf / rmass[i]; v[i][2] += dtfm * f[i][2]; @@ -574,7 +634,7 @@ void FixMove::initial_integrate(int /*vflag*/) x[i][2] += dtv * v[i][2]; } - domain->remap_near(x[i],xold); + domain->remap_near(x[i], xold); } } @@ -597,12 +657,12 @@ void FixMove::initial_integrate(int /*vflag*/) double cosine = cos(arg); double sine = sin(arg); - double qcosine = cos(0.5*arg); - double qsine = sin(0.5*arg); + double qcosine = cos(0.5 * arg); + double qsine = sin(0.5 * arg); qrotate[0] = qcosine; - qrotate[1] = runit[0]*qsine; - qrotate[2] = runit[1]*qsine; - qrotate[3] = runit[2]*qsine; + qrotate[1] = runit[0] * qsine; + qrotate[2] = runit[1] * qsine; + qrotate[3] = runit[2] * qsine; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -613,26 +673,26 @@ void FixMove::initial_integrate(int /*vflag*/) d[0] = xoriginal[i][0] - point[0]; d[1] = xoriginal[i][1] - point[1]; d[2] = xoriginal[i][2] - point[2]; - ddotr = d[0]*runit[0] + d[1]*runit[1] + d[2]*runit[2]; - c[0] = ddotr*runit[0]; - c[1] = ddotr*runit[1]; - c[2] = ddotr*runit[2]; + ddotr = d[0] * runit[0] + d[1] * runit[1] + d[2] * runit[2]; + c[0] = ddotr * runit[0]; + c[1] = ddotr * runit[1]; + c[2] = ddotr * runit[2]; a[0] = d[0] - c[0]; a[1] = d[1] - c[1]; a[2] = d[2] - c[2]; - b[0] = runit[1]*a[2] - runit[2]*a[1]; - b[1] = runit[2]*a[0] - runit[0]*a[2]; - b[2] = runit[0]*a[1] - runit[1]*a[0]; - disp[0] = a[0]*cosine + b[0]*sine; - disp[1] = a[1]*cosine + b[1]*sine; - disp[2] = a[2]*cosine + b[2]*sine; + b[0] = runit[1] * a[2] - runit[2] * a[1]; + b[1] = runit[2] * a[0] - runit[0] * a[2]; + b[2] = runit[0] * a[1] - runit[1] * a[0]; + disp[0] = a[0] * cosine + b[0] * sine; + disp[1] = a[1] * cosine + b[1] * sine; + disp[2] = a[2] * cosine + b[2] * sine; x[i][0] = point[0] + c[0] + disp[0]; x[i][1] = point[1] + c[1] + disp[1]; x[i][2] = point[2] + c[2] + disp[2]; - v[i][0] = omega_rotate * (runit[1]*disp[2] - runit[2]*disp[1]); - v[i][1] = omega_rotate * (runit[2]*disp[0] - runit[0]*disp[2]); - v[i][2] = omega_rotate * (runit[0]*disp[1] - runit[1]*disp[0]); + v[i][0] = omega_rotate * (runit[1] * disp[2] - runit[2] * disp[1]); + v[i][1] = omega_rotate * (runit[2] * disp[0] - runit[0] * disp[2]); + v[i][2] = omega_rotate * (runit[0] * disp[1] - runit[1] * disp[0]); // set any extra attributes affected by rotation @@ -646,9 +706,9 @@ void FixMove::initial_integrate(int /*vflag*/) if (line_flag && line[i] >= 0.0) flag = 1; if (tri_flag && tri[i] >= 0.0) flag = 1; if (flag) { - omega[i][0] = omega_rotate*runit[0]; - omega[i][1] = omega_rotate*runit[1]; - omega[i][2] = omega_rotate*runit[2]; + omega[i][0] = omega_rotate * runit[0]; + omega[i][1] = omega_rotate * runit[1]; + omega[i][2] = omega_rotate * runit[2]; } } @@ -660,11 +720,11 @@ void FixMove::initial_integrate(int /*vflag*/) quat = avec_ellipsoid->bonus[ellipsoid[i]].quat; shape = avec_ellipsoid->bonus[ellipsoid[i]].shape; inertia_ellipsoid[0] = - INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]); + INERTIA * rmass[i] * (shape[1] * shape[1] + shape[2] * shape[2]); inertia_ellipsoid[1] = - INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]); + INERTIA * rmass[i] * (shape[0] * shape[0] + shape[2] * shape[2]); inertia_ellipsoid[2] = - INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]); + INERTIA * rmass[i] * (shape[0] * shape[0] + shape[1] * shape[1]); inertia = inertia_ellipsoid; } else if (tri_flag && tri[i] >= 0) { quat = avec_tri->bonus[tri[i]].quat; @@ -674,18 +734,18 @@ void FixMove::initial_integrate(int /*vflag*/) inertia = avec_body->bonus[body[i]].inertia; } if (quat) { - w[0] = omega_rotate*runit[0]; - w[1] = omega_rotate*runit[1]; - w[2] = omega_rotate*runit[2]; - MathExtra::q_to_exyz(quat,ex,ey,ez); - MathExtra::omega_to_angmom(w,ex,ey,ez,inertia,angmom[i]); + w[0] = omega_rotate * runit[0]; + w[1] = omega_rotate * runit[1]; + w[2] = omega_rotate * runit[2]; + MathExtra::q_to_exyz(quat, ex, ey, ez); + MathExtra::omega_to_angmom(w, ex, ey, ez, inertia, angmom[i]); } } // theta for lines if (theta_flag && line[i] >= 0.0) { - theta_new = fmod(toriginal[i]+arg,MY_2PI); + theta_new = fmod(toriginal[i] + arg, MY_2PI); avec_line->bonus[atom->line[i]].theta = theta_new; } @@ -699,11 +759,149 @@ void FixMove::initial_integrate(int /*vflag*/) quat = avec_tri->bonus[tri[i]].quat; else if (body_flag && body[i] >= 0) quat = avec_body->bonus[body[i]].quat; - if (quat) MathExtra::quatquat(qrotate,qoriginal[i],quat); + if (quat) MathExtra::quatquat(qrotate, qoriginal[i], quat); } } - domain->remap_near(x[i],xold); + domain->remap_near(x[i], xold); + } + } + + // for rotate by right-hand rule around omega: + // P = point = vector = point of rotation + // R = vector = axis of rotation + // w = omega of rotation (from period) + // X0 = xoriginal = initial coord of atom + // R0 = runit = unit vector for R + // D = X0 - P = vector from P to X0 + // C = (D dot R0) R0 = projection of atom coord onto R line + // A = D - C = vector from R line to X0 + // B = R0 cross A = vector perp to A in plane of rotation + // A,B define plane of circular rotation around R line + // X = P + C + A cos(w*dt) + B sin(w*dt) + // V = w R0 cross (A cos(w*dt) + B sin(w*dt)) + // + // add translation after rotation + + } else if (mstyle == TRANSROT) { + double arg = omega_rotate * delta; + double cosine = cos(arg); + double sine = sin(arg); + + double qcosine = cos(0.5 * arg); + double qsine = sin(0.5 * arg); + qrotate[0] = qcosine; + qrotate[1] = runit[0] * qsine; + qrotate[2] = runit[1] * qsine; + qrotate[3] = runit[2] * qsine; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + xold[0] = x[i][0]; + xold[1] = x[i][1]; + xold[2] = x[i][2]; + + d[0] = xoriginal[i][0] - point[0]; + d[1] = xoriginal[i][1] - point[1]; + d[2] = xoriginal[i][2] - point[2]; + ddotr = d[0] * runit[0] + d[1] * runit[1] + d[2] * runit[2]; + c[0] = ddotr * runit[0]; + c[1] = ddotr * runit[1]; + c[2] = ddotr * runit[2]; + a[0] = d[0] - c[0]; + a[1] = d[1] - c[1]; + a[2] = d[2] - c[2]; + b[0] = runit[1] * a[2] - runit[2] * a[1]; + b[1] = runit[2] * a[0] - runit[0] * a[2]; + b[2] = runit[0] * a[1] - runit[1] * a[0]; + disp[0] = a[0] * cosine + b[0] * sine; + disp[1] = a[1] * cosine + b[1] * sine; + disp[2] = a[2] * cosine + b[2] * sine; + + x[i][0] = point[0] + c[0] + disp[0]; + x[i][1] = point[1] + c[1] + disp[1]; + x[i][2] = point[2] + c[2] + disp[2]; + v[i][0] = omega_rotate * (runit[1] * disp[2] - runit[2] * disp[1]); + v[i][1] = omega_rotate * (runit[2] * disp[0] - runit[0] * disp[2]); + v[i][2] = omega_rotate * (runit[0] * disp[1] - runit[1] * disp[0]); + + // set any extra attributes affected by rotation + + if (extra_flag) { + + // omega for spheres, lines, tris + + if (omega_flag) { + flag = 0; + if (radius_flag && radius[i] > 0.0) flag = 1; + if (line_flag && line[i] >= 0.0) flag = 1; + if (tri_flag && tri[i] >= 0.0) flag = 1; + if (flag) { + omega[i][0] = omega_rotate * runit[0]; + omega[i][1] = omega_rotate * runit[1]; + omega[i][2] = omega_rotate * runit[2]; + } + } + + // angmom for ellipsoids, tris, and bodies + + if (angmom_flag) { + quat = inertia = nullptr; + if (ellipsoid_flag && ellipsoid[i] >= 0) { + quat = avec_ellipsoid->bonus[ellipsoid[i]].quat; + shape = avec_ellipsoid->bonus[ellipsoid[i]].shape; + inertia_ellipsoid[0] = + INERTIA * rmass[i] * (shape[1] * shape[1] + shape[2] * shape[2]); + inertia_ellipsoid[1] = + INERTIA * rmass[i] * (shape[0] * shape[0] + shape[2] * shape[2]); + inertia_ellipsoid[2] = + INERTIA * rmass[i] * (shape[0] * shape[0] + shape[1] * shape[1]); + inertia = inertia_ellipsoid; + } else if (tri_flag && tri[i] >= 0) { + quat = avec_tri->bonus[tri[i]].quat; + inertia = avec_tri->bonus[tri[i]].inertia; + } else if (body_flag && body[i] >= 0) { + quat = avec_body->bonus[body[i]].quat; + inertia = avec_body->bonus[body[i]].inertia; + } + if (quat) { + w[0] = omega_rotate * runit[0]; + w[1] = omega_rotate * runit[1]; + w[2] = omega_rotate * runit[2]; + MathExtra::q_to_exyz(quat, ex, ey, ez); + MathExtra::omega_to_angmom(w, ex, ey, ez, inertia, angmom[i]); + } + } + + // theta for lines + + if (theta_flag && line[i] >= 0.0) { + theta_new = fmod(toriginal[i] + arg, MY_2PI); + avec_line->bonus[atom->line[i]].theta = theta_new; + } + + // quats for ellipsoids, tris, and bodies + + if (quat_flag) { + quat = nullptr; + if (ellipsoid_flag && ellipsoid[i] >= 0) + quat = avec_ellipsoid->bonus[ellipsoid[i]].quat; + else if (tri_flag && tri[i] >= 0) + quat = avec_tri->bonus[tri[i]].quat; + else if (body_flag && body[i] >= 0) + quat = avec_body->bonus[body[i]].quat; + if (quat) MathExtra::quatquat(qrotate, qoriginal[i], quat); + } + } + + v[i][0] += vx; + x[i][0] += vx * delta; + v[i][1] += vy; + x[i][1] += vy * delta; + v[i][2] += vz; + x[i][2] += vz * delta; + + domain->remap_near(x[i], xold); } } @@ -720,11 +918,11 @@ void FixMove::initial_integrate(int /*vflag*/) maxatom = atom->nmax; if (displaceflag) { memory->destroy(displace); - memory->create(displace,maxatom,3,"move:displace"); + memory->create(displace, maxatom, 3, "move:displace"); } if (velocityflag) { memory->destroy(velocity); - memory->create(velocity,maxatom,3,"move:velocity"); + memory->create(velocity, maxatom, 3, "move:velocity"); } } @@ -733,28 +931,40 @@ void FixMove::initial_integrate(int /*vflag*/) modify->clearstep_compute(); if (xvarstr) { - if (xvarstyle == EQUAL) dx = input->variable->compute_equal(xvar); - else input->variable->compute_atom(xvar,igroup,&displace[0][0],3,0); + if (xvarstyle == EQUAL) + dx = input->variable->compute_equal(xvar); + else + input->variable->compute_atom(xvar, igroup, &displace[0][0], 3, 0); } if (yvarstr) { - if (yvarstyle == EQUAL) dy = input->variable->compute_equal(yvar); - else input->variable->compute_atom(yvar,igroup,&displace[0][1],3,0); + if (yvarstyle == EQUAL) + dy = input->variable->compute_equal(yvar); + else + input->variable->compute_atom(yvar, igroup, &displace[0][1], 3, 0); } if (zvarstr) { - if (zvarstyle == EQUAL) dz = input->variable->compute_equal(zvar); - else input->variable->compute_atom(zvar,igroup,&displace[0][2],3,0); + if (zvarstyle == EQUAL) + dz = input->variable->compute_equal(zvar); + else + input->variable->compute_atom(zvar, igroup, &displace[0][2], 3, 0); } if (vxvarstr) { - if (vxvarstyle == EQUAL) vx = input->variable->compute_equal(vxvar); - else input->variable->compute_atom(vxvar,igroup,&velocity[0][0],3,0); + if (vxvarstyle == EQUAL) + vx = input->variable->compute_equal(vxvar); + else + input->variable->compute_atom(vxvar, igroup, &velocity[0][0], 3, 0); } if (vyvarstr) { - if (vyvarstyle == EQUAL) vy = input->variable->compute_equal(vyvar); - else input->variable->compute_atom(vyvar,igroup,&velocity[0][1],3,0); + if (vyvarstyle == EQUAL) + vy = input->variable->compute_equal(vyvar); + else + input->variable->compute_atom(vyvar, igroup, &velocity[0][1], 3, 0); } if (vzvarstr) { - if (vzvarstyle == EQUAL) vz = input->variable->compute_equal(vzvar); - else input->variable->compute_atom(vzvar,igroup,&velocity[0][2],3,0); + if (vzvarstyle == EQUAL) + vz = input->variable->compute_equal(vzvar); + else + input->variable->compute_atom(vzvar, igroup, &velocity[0][2], 3, 0); } modify->addstep_compute(update->ntimestep + 1); @@ -768,16 +978,24 @@ void FixMove::initial_integrate(int /*vflag*/) xold[2] = x[i][2]; if (xvarstr && vxvarstr) { - if (vxvarstyle == EQUAL) v[i][0] = vx; - else v[i][0] = velocity[i][0]; - if (xvarstyle == EQUAL) x[i][0] = xoriginal[i][0] + dx; - else x[i][0] = xoriginal[i][0] + displace[i][0]; + if (vxvarstyle == EQUAL) + v[i][0] = vx; + else + v[i][0] = velocity[i][0]; + if (xvarstyle == EQUAL) + x[i][0] = xoriginal[i][0] + dx; + else + x[i][0] = xoriginal[i][0] + displace[i][0]; } else if (xvarstr) { - if (xvarstyle == EQUAL) x[i][0] = xoriginal[i][0] + dx; - else x[i][0] = xoriginal[i][0] + displace[i][0]; + if (xvarstyle == EQUAL) + x[i][0] = xoriginal[i][0] + dx; + else + x[i][0] = xoriginal[i][0] + displace[i][0]; } else if (vxvarstr) { - if (vxvarstyle == EQUAL) v[i][0] = vx; - else v[i][0] = velocity[i][0]; + if (vxvarstyle == EQUAL) + v[i][0] = vx; + else + v[i][0] = velocity[i][0]; x[i][0] += dtv * v[i][0]; } else { if (rmass) { @@ -791,16 +1009,24 @@ void FixMove::initial_integrate(int /*vflag*/) } if (yvarstr && vyvarstr) { - if (vyvarstyle == EQUAL) v[i][1] = vy; - else v[i][1] = velocity[i][1]; - if (yvarstyle == EQUAL) x[i][1] = xoriginal[i][1] + dy; - else x[i][1] = xoriginal[i][1] + displace[i][1]; + if (vyvarstyle == EQUAL) + v[i][1] = vy; + else + v[i][1] = velocity[i][1]; + if (yvarstyle == EQUAL) + x[i][1] = xoriginal[i][1] + dy; + else + x[i][1] = xoriginal[i][1] + displace[i][1]; } else if (yvarstr) { - if (yvarstyle == EQUAL) x[i][1] = xoriginal[i][1] + dy; - else x[i][1] = xoriginal[i][1] + displace[i][1]; + if (yvarstyle == EQUAL) + x[i][1] = xoriginal[i][1] + dy; + else + x[i][1] = xoriginal[i][1] + displace[i][1]; } else if (vyvarstr) { - if (vyvarstyle == EQUAL) v[i][1] = vy; - else v[i][1] = velocity[i][1]; + if (vyvarstyle == EQUAL) + v[i][1] = vy; + else + v[i][1] = velocity[i][1]; x[i][1] += dtv * v[i][1]; } else { if (rmass) { @@ -814,16 +1040,24 @@ void FixMove::initial_integrate(int /*vflag*/) } if (zvarstr && vzvarstr) { - if (vzvarstyle == EQUAL) v[i][2] = vz; - else v[i][2] = velocity[i][2]; - if (zvarstyle == EQUAL) x[i][2] = xoriginal[i][2] + dz; - else x[i][2] = xoriginal[i][2] + displace[i][2]; + if (vzvarstyle == EQUAL) + v[i][2] = vz; + else + v[i][2] = velocity[i][2]; + if (zvarstyle == EQUAL) + x[i][2] = xoriginal[i][2] + dz; + else + x[i][2] = xoriginal[i][2] + displace[i][2]; } else if (zvarstr) { - if (zvarstyle == EQUAL) x[i][2] = xoriginal[i][2] + dz; - else x[i][2] = xoriginal[i][2] + displace[i][2]; + if (zvarstyle == EQUAL) + x[i][2] = xoriginal[i][2] + dz; + else + x[i][2] = xoriginal[i][2] + displace[i][2]; } else if (vzvarstr) { - if (vzvarstyle == EQUAL) v[i][2] = vz; - else v[i][2] = velocity[i][2]; + if (vzvarstyle == EQUAL) + v[i][2] = vz; + else + v[i][2] = velocity[i][2]; x[i][2] += dtv * v[i][2]; } else { if (rmass) { @@ -836,7 +1070,7 @@ void FixMove::initial_integrate(int /*vflag*/) x[i][2] += dtv * v[i][2]; } - domain->remap_near(x[i],xold); + domain->remap_near(x[i], xold); } } } @@ -851,22 +1085,40 @@ void FixMove::final_integrate() double dtfm; int xflag = 1; - if (mstyle == LINEAR && vxflag) xflag = 0; - else if (mstyle == WIGGLE && axflag) xflag = 0; - else if (mstyle == ROTATE) xflag = 0; - else if (mstyle == VARIABLE && (xvarstr || vxvarstr)) xflag = 0; + if (mstyle == LINEAR && vxflag) + xflag = 0; + else if (mstyle == WIGGLE && axflag) + xflag = 0; + else if (mstyle == ROTATE) + xflag = 0; + else if (mstyle == TRANSROT) + xflag = 0; + else if (mstyle == VARIABLE && (xvarstr || vxvarstr)) + xflag = 0; int yflag = 1; - if (mstyle == LINEAR && vyflag) yflag = 0; - else if (mstyle == WIGGLE && ayflag) yflag = 0; - else if (mstyle == ROTATE) yflag = 0; - else if (mstyle == VARIABLE && (yvarstr || vyvarstr)) yflag = 0; + if (mstyle == LINEAR && vyflag) + yflag = 0; + else if (mstyle == WIGGLE && ayflag) + yflag = 0; + else if (mstyle == ROTATE) + yflag = 0; + else if (mstyle == TRANSROT) + yflag = 0; + else if (mstyle == VARIABLE && (yvarstr || vyvarstr)) + yflag = 0; int zflag = 1; - if (mstyle == LINEAR && vzflag) zflag = 0; - else if (mstyle == WIGGLE && azflag) zflag = 0; - else if (mstyle == ROTATE) zflag = 0; - else if (mstyle == VARIABLE && (zvarstr || vzvarstr)) zflag = 0; + if (mstyle == LINEAR && vzflag) + zflag = 0; + else if (mstyle == WIGGLE && azflag) + zflag = 0; + else if (mstyle == ROTATE) + zflag = 0; + else if (mstyle == TRANSROT) + zflag = 0; + else if (mstyle == VARIABLE && (zvarstr || vzvarstr)) + zflag = 0; double **v = atom->v; double **f = atom->f; @@ -918,14 +1170,14 @@ void FixMove::initial_integrate_respa(int vflag, int ilevel, int /*iloop*/) // outermost level - update v and x // all other levels - nothing - if (ilevel == nlevels_respa-1) initial_integrate(vflag); + if (ilevel == nlevels_respa - 1) initial_integrate(vflag); } /* ---------------------------------------------------------------------- */ void FixMove::final_integrate_respa(int ilevel, int /*iloop*/) { - if (ilevel == nlevels_respa-1) final_integrate(); + if (ilevel == nlevels_respa - 1) final_integrate(); } /* ---------------------------------------------------------------------- @@ -934,11 +1186,11 @@ void FixMove::final_integrate_respa(int ilevel, int /*iloop*/) double FixMove::memory_usage() { - double bytes = (double)atom->nmax*3 * sizeof(double); - if (theta_flag) bytes += (double)atom->nmax * sizeof(double); - if (quat_flag) bytes += (double)atom->nmax*4 * sizeof(double); - if (displaceflag) bytes += (double)atom->nmax*3 * sizeof(double); - if (velocityflag) bytes += (double)atom->nmax*3 * sizeof(double); + double bytes = (double) atom->nmax * 3 * sizeof(double); + if (theta_flag) bytes += (double) atom->nmax * sizeof(double); + if (quat_flag) bytes += (double) atom->nmax * 4 * sizeof(double); + if (displaceflag) bytes += (double) atom->nmax * 3 * sizeof(double); + if (velocityflag) bytes += (double) atom->nmax * 3 * sizeof(double); return bytes; } @@ -954,8 +1206,8 @@ void FixMove::write_restart(FILE *fp) if (comm->me == 0) { int size = n * sizeof(double); - fwrite(&size,sizeof(int),1,fp); - fwrite(list,sizeof(double),n,fp); + fwrite(&size, sizeof(int), 1, fp); + fwrite(list, sizeof(double), n, fp); } } @@ -968,7 +1220,7 @@ void FixMove::restart(char *buf) int n = 0; double *list = (double *) buf; - time_origin = static_cast (list[n++]); + time_origin = static_cast(list[n++]); } /* ---------------------------------------------------------------------- @@ -977,9 +1229,9 @@ void FixMove::restart(char *buf) 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"); + 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"); array_atom = xoriginal; } @@ -1028,50 +1280,51 @@ void FixMove::set_arrays(int i) // current time still equal fix creation time if (update->ntimestep == time_origin) { - domain->unmap(x[i],image[i],xoriginal[i]); + domain->unmap(x[i], image[i], xoriginal[i]); return; } // backup particle to time_origin - if (mstyle == VARIABLE) - error->all(FLERR,"Cannot add atoms to fix move variable"); + if (mstyle == VARIABLE) error->all(FLERR, "Cannot add atoms to fix move variable"); - domain->unmap(x[i],image[i],xoriginal[i]); + domain->unmap(x[i], image[i], xoriginal[i]); double delta = (update->ntimestep - time_origin) * update->dt; if (mstyle == LINEAR) { if (vxflag) xoriginal[i][0] -= vx * delta; if (vyflag) xoriginal[i][1] -= vy * delta; if (vzflag) xoriginal[i][2] -= vz * delta; + } else if (mstyle == WIGGLE) { double arg = omega_rotate * delta; double sine = sin(arg); - if (axflag) xoriginal[i][0] -= ax*sine; - if (ayflag) xoriginal[i][1] -= ay*sine; - if (azflag) xoriginal[i][2] -= az*sine; + if (axflag) xoriginal[i][0] -= ax * sine; + if (ayflag) xoriginal[i][1] -= ay * sine; + if (azflag) xoriginal[i][2] -= az * sine; + } else if (mstyle == ROTATE) { - double a[3],b[3],c[3],d[3],disp[3],ddotr; - double arg = - omega_rotate * delta; + double a[3], b[3], c[3], d[3], disp[3], ddotr; + double arg = -omega_rotate * delta; double sine = sin(arg); double cosine = cos(arg); d[0] = x[i][0] - point[0]; d[1] = x[i][1] - point[1]; d[2] = x[i][2] - point[2]; - ddotr = d[0]*runit[0] + d[1]*runit[1] + d[2]*runit[2]; - c[0] = ddotr*runit[0]; - c[1] = ddotr*runit[1]; - c[2] = ddotr*runit[2]; + ddotr = d[0] * runit[0] + d[1] * runit[1] + d[2] * runit[2]; + c[0] = ddotr * runit[0]; + c[1] = ddotr * runit[1]; + c[2] = ddotr * runit[2]; a[0] = d[0] - c[0]; a[1] = d[1] - c[1]; a[2] = d[2] - c[2]; - b[0] = runit[1]*a[2] - runit[2]*a[1]; - b[1] = runit[2]*a[0] - runit[0]*a[2]; - b[2] = runit[0]*a[1] - runit[1]*a[0]; - disp[0] = a[0]*cosine + b[0]*sine; - disp[1] = a[1]*cosine + b[1]*sine; - disp[2] = a[2]*cosine + b[2]*sine; + b[0] = runit[1] * a[2] - runit[2] * a[1]; + b[1] = runit[2] * a[0] - runit[0] * a[2]; + b[2] = runit[0] * a[1] - runit[1] * a[0]; + disp[0] = a[0] * cosine + b[0] * sine; + disp[1] = a[1] * cosine + b[1] * sine; + disp[2] = a[2] * cosine + b[2] * sine; xoriginal[i][0] = point[0] + c[0] + disp[0]; xoriginal[i][1] = point[1] + c[1] + disp[1]; @@ -1085,7 +1338,64 @@ void FixMove::set_arrays(int i) if (theta_flag && line[i] >= 0.0) { theta = avec_line->bonus[atom->line[i]].theta; - toriginal[i] = theta - 0.0; // NOTE: edit this line + toriginal[i] = theta - 0.0; // NOTE: edit this line + } + + // quats for ellipsoids, tris, and bodies + + if (quat_flag) { + quat = nullptr; + if (ellipsoid_flag && ellipsoid[i] >= 0) + quat = avec_ellipsoid->bonus[ellipsoid[i]].quat; + else if (tri_flag && tri[i] >= 0) + quat = avec_tri->bonus[tri[i]].quat; + else if (body_flag && body[i] >= 0) + quat = avec_body->bonus[body[i]].quat; + if (quat) { + // qoriginal = f(quat,-delta); // NOTE: edit this line + } + } + } + xoriginal[i][0] -= vx * delta; + xoriginal[i][1] -= vy * delta; + xoriginal[i][2] -= vz * delta; + + } else if (mstyle == TRANSROT) { + double a[3], b[3], c[3], d[3], disp[3], ddotr; + double arg = -omega_rotate * delta; + double sine = sin(arg); + double cosine = cos(arg); + d[0] = x[i][0] - point[0]; + d[1] = x[i][1] - point[1]; + d[2] = x[i][2] - point[2]; + ddotr = d[0] * runit[0] + d[1] * runit[1] + d[2] * runit[2]; + c[0] = ddotr * runit[0]; + c[1] = ddotr * runit[1]; + c[2] = ddotr * runit[2]; + + a[0] = d[0] - c[0]; + a[1] = d[1] - c[1]; + a[2] = d[2] - c[2]; + b[0] = runit[1] * a[2] - runit[2] * a[1]; + b[1] = runit[2] * a[0] - runit[0] * a[2]; + b[2] = runit[0] * a[1] - runit[1] * a[0]; + disp[0] = a[0] * cosine + b[0] * sine; + disp[1] = a[1] * cosine + b[1] * sine; + disp[2] = a[2] * cosine + b[2] * sine; + + xoriginal[i][0] = point[0] + c[0] + disp[0]; + xoriginal[i][1] = point[1] + c[1] + disp[1]; + xoriginal[i][2] = point[2] + c[2] + disp[2]; + + // set theta and quat extra attributes affected by rotation + + if (extra_flag) { + + // theta for lines + + if (theta_flag && line[i] >= 0.0) { + theta = avec_line->bonus[atom->line[i]].theta; + toriginal[i] = theta - 0.0; // NOTE: edit this line } // quats for ellipsoids, tris, and bodies @@ -1180,7 +1490,7 @@ void FixMove::unpack_restart(int nlocal, int nth) // unpack the Nth first values this way because other fixes pack them int m = 0; - for (int i = 0; i < nth; i++) m += static_cast (extra[nlocal][m]); + for (int i = 0; i < nth; i++) m += static_cast(extra[nlocal][m]); m++; xoriginal[nlocal][0] = extra[nlocal][m++]; @@ -1217,5 +1527,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, "Resetting timestep size is not allowed with fix move"); }