ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

add "transrot" style to fix move that allows to do translation and rotation at the same time

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-04-02 18:34:17 -04:00
parent 584b166823
commit c9d0889f25
2 changed files with 632 additions and 304 deletions
Download Patch File Download Diff File Expand all files Collapse all files

84
doc/src/fix_move.rst
View File

@ -12,7 +12,7 @@ Syntax
* ID, group-ID are documented in :doc:`fix <fix>` command * ID, group-ID are documented in :doc:`fix <fix>` command
* move = style name of this 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:: .. parsed-literal::
@ -25,6 +25,11 @@ Syntax
Px,Py,Pz = origin point of axis of rotation (distance units) Px,Py,Pz = origin point of axis of rotation (distance units)
Rx,Ry,Rz = axis of rotation vector Rx,Ry,Rz = axis of rotation vector
period = period of rotation (time units) 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 *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_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 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 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 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 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 Description
""""""""""" """""""""""
@ -55,15 +61,17 @@ whose movement can influence nearby atoms.
.. note:: .. note::
The atoms affected by this fix should not normally be time The atoms affected by this fix should not normally be time integrated
integrated by other fixes (e.g. :doc:`fix nve <fix_nve>`, :doc:`fix nvt <fix_nh>`), since that will change their positions and by other fixes (e.g. :doc:`fix nve <fix_nve>`, :doc:`fix nvt
velocities twice. <fix_nh>`), since that will change their positions and velocities
twice.
.. note:: .. note::
As atoms move due to this fix, they will pass through periodic As atoms move due to this fix, they will pass through periodic
boundaries and be remapped to the other side of the simulation box, 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 <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 <fix_nve>` command). It is up to you to decide whether
periodic boundaries are appropriate with the kind of atom motion you periodic boundaries are appropriate with the kind of atom motion you
are prescribing with this fix. 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 position at the time the fix is specified. These initial coordinates
are stored by the fix in "unwrapped" form, by using the image flags are stored by the fix in "unwrapped" form, by using the image flags
associated with each atom. See the :doc:`dump custom <dump>` command associated with each atom. See the :doc:`dump custom <dump>` command
for a discussion of "unwrapped" coordinates. See the Atoms section of for a discussion of "unwrapped" coordinates. See the Atoms section
the :doc:`read_data <read_data>` command for a discussion of image flags of the :doc:`read_data <read_data>` command for a discussion of image
and how they are set for each atom. You can reset the image flags flags and how they are set for each atom. You can reset the image
(e.g. to 0) before invoking this fix by using the :doc:`set image <set>` flags (e.g. to 0) before invoking this fix by using the :doc:`set
command. image <set>` command.
---------- ----------
@ -118,13 +126,13 @@ notation as
where *X0* = (x0,y0,z0) is their position at the time the fix is where *X0* = (x0,y0,z0) is their position at the time the fix is
specified, *A* is the specified amplitude vector with components specified, *A* is the specified amplitude vector with components
(Ax,Ay,Az), *omega* is 2 PI / *period*, and *delta* is the time (Ax,Ay,Az), *omega* is 2 PI / *period*, and *delta* is the time elapsed
elapsed since the fix was specified. This style also sets the since the fix was specified. This style also sets the velocity of each
velocity of each atom to the time derivative of this expression. If atom to the time derivative of this expression. If any of the amplitude
any of the amplitude components is specified as NULL, then the components is specified as NULL, then the position and velocity of that
position and velocity of that component is time integrated the same as component is time integrated the same as the :doc:`fix nve <fix_nve>`
the :doc:`fix nve <fix_nve>` command would perform, using the command would perform, using the corresponding force component on the
corresponding force component on the atom. atom.
Note that the *wiggle* style is identical to using the *variable* Note that the *wiggle* style is identical to using the *variable*
style with :doc:`equal-style variables <variable>` that use the style with :doc:`equal-style variables <variable>` that use the
@ -139,21 +147,29 @@ swiggle() and cwiggle() functions. E.g.
variable v equal v_omega*($A-cwiggle(0.0,$A,$T)) variable v equal v_omega*($A-cwiggle(0.0,$A,$T))
fix 1 boundary move variable v_x NULL NULL v_v NULL NULL fix 1 boundary move variable v_x NULL NULL v_v NULL NULL
The *rotate* style rotates atoms around a rotation axis *R* = The *rotate* style rotates atoms around a rotation axis *R* = (Rx,Ry,Rz)
(Rx,Ry,Rz) that goes through a point *P* = (Px,Py,Pz). The *period* of that goes through a point *P* = (Px,Py,Pz). The *period* of the
the rotation is also specified. The direction of rotation for the rotation is also specified. The direction of rotation for the atoms
atoms around the rotation axis is consistent with the right-hand rule: around the rotation axis is consistent with the right-hand rule: if your
if your right-hand thumb points along *R*, then your fingers wrap right-hand thumb points along *R*, then your fingers wrap around the
around the axis in the direction of rotation. axis in the direction of rotation.
This style also sets the velocity of each atom to (omega cross Rperp) 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 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 is a perpendicular vector from the rotation axis to the atom. If the
defined :doc:`atom_style <atom_style>` assigns an angular velocity or defined :doc:`atom_style <atom_style>` assigns an angular velocity or
angular momentum or orientation to each atom (:doc:`atom styles <atom_style>` sphere, ellipsoid, line, tri, body), then angular momentum or orientation to each atom (:doc:`atom styles
<atom_style>` sphere, ellipsoid, line, tri, body), then
those properties are also updated appropriately to correspond to the those properties are also updated appropriately to correspond to the
atom's motion and rotation over time. 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 The *variable* style allows the position and velocity components of
each atom to be set by formulas specified via the each atom to be set by formulas specified via the
:doc:`variable <variable>` command. Each of the 6 variables is :doc:`variable <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. a function of the timestep as well as of other simulation values.
*Atom*\ -style variables compute a numeric quantity for each atom, that *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 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 well as of the timestep and other simulation values. Note that this fix
fix stores the original coordinates of each atom (see note below) so stores the original coordinates of each atom (see note below) so that
that per-atom quantity can be used in an atom-style variable formula. per-atom quantity can be used in an atom-style variable formula. See
See the :doc:`variable <variable>` command for details. the :doc:`variable <variable>` command for details.
The first 3 variables (v_dx,v_dy,v_dz) specified for the *variable* 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 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 Restart, fix_modify, output, run start/stop, minimize info
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This fix writes the original coordinates of moving atoms to :doc:`binary restart files <restart>`, as well as the initial timestep, so that This fix writes the original coordinates of moving atoms to :doc:`binary
the motion can be continuous in a restarted simulation. See the restart files <restart>`, as well as the initial timestep, so that the
motion can be continuous in a restarted simulation. See the
:doc:`read_restart <read_restart>` command for info on how to re-specify :doc:`read_restart <read_restart>` command for info on how to re-specify
a fix in an input script that reads a restart file, so that the a fix in an input script that reads a restart file, so that the
operation of the fix continues in an uninterrupted fashion. 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 This fix produces a per-atom array which can be accessed by various
:doc:`output commands <Howto_output>`. The number of columns for each :doc:`output commands <Howto_output>`. The number of columns for each
atom is 3, and the columns store the original unwrapped x,y,z coords atom is 3, and the columns store the original unwrapped x,y,z coords of
of each atom. The per-atom values can be accessed on any timestep. 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 No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is not invoked during :doc:`energy minimization <minimize>`. the :doc:`run <run>` command. This fix is not invoked during
:doc:`energy minimization <minimize>`.
For :doc:`rRESPA time integration <run_style>`, this fix adjusts the For :doc:`rRESPA time integration <run_style>`, this fix adjusts the
position and velocity of atoms on the outermost rRESPA level. position and velocity of atoms on the outermost rRESPA level.

572
src/fix_move.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -40,7 +39,7 @@ using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
using namespace MathConst; using namespace MathConst;
enum{LINEAR,WIGGLE,ROTATE,VARIABLE}; enum { LINEAR, WIGGLE, ROTATE, VARIABLE, TRANSROT };
enum { EQUAL, ATOM }; enum { EQUAL, ATOM };
#define INERTIA 0.2 // moment of inertia prefactor for ellipsoid #define INERTIA 0.2 // moment of inertia prefactor for ellipsoid
@ -48,11 +47,9 @@ enum{EQUAL,ATOM};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) : FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), Fix(lmp, narg, arg), xvarstr(nullptr), yvarstr(nullptr), zvarstr(nullptr), vxvarstr(nullptr),
xvarstr(nullptr), yvarstr(nullptr), zvarstr(nullptr), vxvarstr(nullptr), vyvarstr(nullptr), vzvarstr(nullptr), xoriginal(nullptr), toriginal(nullptr),
vyvarstr(nullptr), vzvarstr(nullptr), qoriginal(nullptr), displace(nullptr), velocity(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");
@ -75,17 +72,20 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
if (narg < 7) error->all(FLERR, "Illegal fix move command"); if (narg < 7) error->all(FLERR, "Illegal fix move command");
iarg = 7; iarg = 7;
mstyle = LINEAR; mstyle = LINEAR;
if (strcmp(arg[4],"NULL") == 0) vxflag = 0; if (strcmp(arg[4], "NULL") == 0)
vxflag = 0;
else { else {
vxflag = 1; 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 { else {
vyflag = 1; 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 { else {
vzflag = 1; vzflag = 1;
vz = utils::numeric(FLERR, arg[6], false, lmp); vz = utils::numeric(FLERR, arg[6], false, lmp);
@ -95,17 +95,20 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
if (narg < 8) error->all(FLERR, "Illegal fix move command"); if (narg < 8) error->all(FLERR, "Illegal fix move command");
iarg = 8; iarg = 8;
mstyle = WIGGLE; mstyle = WIGGLE;
if (strcmp(arg[4],"NULL") == 0) axflag = 0; if (strcmp(arg[4], "NULL") == 0)
axflag = 0;
else { else {
axflag = 1; 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 { else {
ayflag = 1; 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 { else {
azflag = 1; azflag = 1;
az = utils::numeric(FLERR, arg[6], false, lmp); az = utils::numeric(FLERR, arg[6], false, lmp);
@ -126,36 +129,66 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
period = utils::numeric(FLERR, arg[10], 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, "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) { } else if (strcmp(arg[3], "variable") == 0) {
if (narg < 10) error->all(FLERR, "Illegal fix move command"); if (narg < 10) error->all(FLERR, "Illegal fix move command");
iarg = 10; iarg = 10;
mstyle = VARIABLE; mstyle = VARIABLE;
if (strcmp(arg[4],"NULL") == 0) xvarstr = nullptr; if (strcmp(arg[4], "NULL") == 0)
xvarstr = nullptr;
else if (utils::strmatch(arg[4], "^v_")) { else if (utils::strmatch(arg[4], "^v_")) {
xvarstr = utils::strdup(arg[4] + 2); xvarstr = utils::strdup(arg[4] + 2);
} else error->all(FLERR,"Illegal fix move command"); } else
if (strcmp(arg[5],"NULL") == 0) yvarstr = nullptr; error->all(FLERR, "Illegal fix move command");
if (strcmp(arg[5], "NULL") == 0)
yvarstr = nullptr;
else if (utils::strmatch(arg[5], "^v_")) { else if (utils::strmatch(arg[5], "^v_")) {
yvarstr = utils::strdup(arg[5] + 2); yvarstr = utils::strdup(arg[5] + 2);
} else error->all(FLERR,"Illegal fix move command"); } else
if (strcmp(arg[6],"NULL") == 0) zvarstr = nullptr; error->all(FLERR, "Illegal fix move command");
if (strcmp(arg[6], "NULL") == 0)
zvarstr = nullptr;
else if (utils::strmatch(arg[6], "^v_")) { else if (utils::strmatch(arg[6], "^v_")) {
zvarstr = utils::strdup(arg[6] + 2); zvarstr = utils::strdup(arg[6] + 2);
} else error->all(FLERR,"Illegal fix move command"); } else
if (strcmp(arg[7],"NULL") == 0) vxvarstr = nullptr; error->all(FLERR, "Illegal fix move command");
if (strcmp(arg[7], "NULL") == 0)
vxvarstr = nullptr;
else if (utils::strmatch(arg[7], "^v_")) { else if (utils::strmatch(arg[7], "^v_")) {
vxvarstr = utils::strdup(arg[7] + 2); vxvarstr = utils::strdup(arg[7] + 2);
} else error->all(FLERR,"Illegal fix move command"); } else
if (strcmp(arg[8],"NULL") == 0) vyvarstr = nullptr; error->all(FLERR, "Illegal fix move command");
if (strcmp(arg[8], "NULL") == 0)
vyvarstr = nullptr;
else if (utils::strmatch(arg[8], "^v_")) { else if (utils::strmatch(arg[8], "^v_")) {
vyvarstr = utils::strdup(arg[8] + 2); vyvarstr = utils::strdup(arg[8] + 2);
} else error->all(FLERR,"Illegal fix move command"); } else
if (strcmp(arg[9],"NULL") == 0) vzvarstr = nullptr; error->all(FLERR, "Illegal fix move command");
if (strcmp(arg[9], "NULL") == 0)
vzvarstr = nullptr;
else if (utils::strmatch(arg[9], "^v_")) { else if (utils::strmatch(arg[9], "^v_")) {
vzvarstr = utils::strdup(arg[9] + 2); 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"); } else
error->all(FLERR, "Illegal fix move command");
// optional args // optional args
@ -164,32 +197,33 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg], "units") == 0) { if (strcmp(arg[iarg], "units") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix move command"); if (iarg + 2 > narg) error->all(FLERR, "Illegal fix move command");
if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; if (strcmp(arg[iarg + 1], "box") == 0)
else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; scaleflag = 0;
else error->all(FLERR,"Illegal fix move command"); else if (strcmp(arg[iarg + 1], "lattice") == 0)
scaleflag = 1;
else
error->all(FLERR, "Illegal fix move command");
iarg += 2; iarg += 2;
} else error->all(FLERR,"Illegal fix move command"); } else
error->all(FLERR, "Illegal fix move command");
} }
// error checks and warnings // error checks and warnings
if (domain->dimension == 2) { if (domain->dimension == 2) {
if (mstyle == LINEAR && vzflag && vz != 0.0) 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, "Fix move cannot set linear z motion for 2d problem");
if (mstyle == WIGGLE && azflag && az != 0.0) if (mstyle == WIGGLE && azflag && az != 0.0)
error->all(FLERR, "Fix move cannot set wiggle z motion for 2d problem"); error->all(FLERR, "Fix move cannot set wiggle z motion for 2d problem");
if (mstyle == ROTATE && (axis[0] != 0.0 || axis[1] != 0.0)) if (((mstyle == ROTATE) || (mstyle == TRANSROT)) && (axis[0] != 0.0 || axis[1] != 0.0))
error->all(FLERR, error->all(FLERR, "Fix move cannot rotate around non z-axis for 2d problem");
"Fix move cannot rotate around non z-axis for 2d problem");
if (mstyle == VARIABLE && (zvarstr || vzvarstr)) if (mstyle == VARIABLE && (zvarstr || vzvarstr))
error->all(FLERR, error->all(FLERR, "Fix move cannot define z or vz variable for 2d problem");
"Fix move cannot define z or vz variable for 2d problem");
} }
// setup scaling and apply scaling factors to velocity & amplitude // setup scaling and apply scaling factors to velocity & amplitude
if ((mstyle == LINEAR || mstyle == WIGGLE || mstyle == ROTATE) && if ((mstyle != VARIABLE) && scaleflag) {
scaleflag) {
double xscale = domain->lattice->xlattice; double xscale = domain->lattice->xlattice;
double yscale = domain->lattice->ylattice; double yscale = domain->lattice->ylattice;
double zscale = domain->lattice->zlattice; double zscale = domain->lattice->zlattice;
@ -206,19 +240,26 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
point[0] *= xscale; point[0] *= xscale;
point[1] *= yscale; point[1] *= yscale;
point[2] *= zscale; 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 // 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 // runit = unit vector along rotation axis
if (mstyle == ROTATE) { if ((mstyle == ROTATE) || (mstyle == TRANSROT)) {
double len = sqrt(axis[0] * axis[0] + axis[1] * axis[1] + axis[2] * axis[2]); double len = sqrt(axis[0] * axis[0] + axis[1] * axis[1] + axis[2] * axis[2]);
if (len == 0.0) if (len == 0.0) error->all(FLERR, "Zero length rotation vector with fix move");
error->all(FLERR,"Zero length rotation vector with fix move");
runit[0] = axis[0] / len; runit[0] = axis[0] / len;
runit[1] = axis[1] / len; runit[1] = axis[1] / len;
runit[2] = axis[2] / len; runit[2] = axis[2] / len;
@ -273,15 +314,18 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]); if (mask[i] & groupbit)
else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0; domain->unmap(x[i], image[i], xoriginal[i]);
else
xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
} }
if (theta_flag) { if (theta_flag) {
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && line[i] >= 0) if ((mask[i] & groupbit) && line[i] >= 0)
toriginal[i] = avec_line->bonus[line[i]].theta; 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][1] = quat[1];
qoriginal[i][2] = quat[2]; qoriginal[i][2] = quat[2];
qoriginal[i][3] = quat[3]; qoriginal[i][3] = quat[3];
} else qoriginal[i][0] = qoriginal[i][1] = } else
qoriginal[i][2] = qoriginal[i][3] = 0.0; qoriginal[i][0] = qoriginal[i][1] = qoriginal[i][2] = qoriginal[i][3] = 0.0;
} }
} }
@ -371,51 +415,63 @@ void FixMove::init()
if (mstyle == VARIABLE) { if (mstyle == VARIABLE) {
if (xvarstr) { if (xvarstr) {
xvar = input->variable->find(xvarstr); xvar = input->variable->find(xvarstr);
if (xvar < 0) error->all(FLERR, if (xvar < 0) error->all(FLERR, "Variable name for fix move does not exist");
"Variable name for fix move does not exist"); if (input->variable->equalstyle(xvar))
if (input->variable->equalstyle(xvar)) xvarstyle = EQUAL; xvarstyle = EQUAL;
else if (input->variable->atomstyle(xvar)) xvarstyle = ATOM; else if (input->variable->atomstyle(xvar))
else error->all(FLERR,"Variable for fix move is invalid style"); xvarstyle = ATOM;
else
error->all(FLERR, "Variable for fix move is invalid style");
} }
if (yvarstr) { if (yvarstr) {
yvar = input->variable->find(yvarstr); yvar = input->variable->find(yvarstr);
if (yvar < 0) error->all(FLERR, if (yvar < 0) error->all(FLERR, "Variable name for fix move does not exist");
"Variable name for fix move does not exist"); if (input->variable->equalstyle(yvar))
if (input->variable->equalstyle(yvar)) yvarstyle = EQUAL; yvarstyle = EQUAL;
else if (input->variable->atomstyle(yvar)) yvarstyle = ATOM; else if (input->variable->atomstyle(yvar))
else error->all(FLERR,"Variable for fix move is invalid style"); yvarstyle = ATOM;
else
error->all(FLERR, "Variable for fix move is invalid style");
} }
if (zvarstr) { if (zvarstr) {
zvar = input->variable->find(zvarstr); zvar = input->variable->find(zvarstr);
if (zvar < 0) error->all(FLERR, if (zvar < 0) error->all(FLERR, "Variable name for fix move does not exist");
"Variable name for fix move does not exist"); if (input->variable->equalstyle(zvar))
if (input->variable->equalstyle(zvar)) zvarstyle = EQUAL; zvarstyle = EQUAL;
else if (input->variable->atomstyle(zvar)) zvarstyle = ATOM; else if (input->variable->atomstyle(zvar))
else error->all(FLERR,"Variable for fix move is invalid style"); zvarstyle = ATOM;
else
error->all(FLERR, "Variable for fix move is invalid style");
} }
if (vxvarstr) { if (vxvarstr) {
vxvar = input->variable->find(vxvarstr); vxvar = input->variable->find(vxvarstr);
if (vxvar < 0) error->all(FLERR, if (vxvar < 0) error->all(FLERR, "Variable name for fix move does not exist");
"Variable name for fix move does not exist"); if (input->variable->equalstyle(vxvar))
if (input->variable->equalstyle(vxvar)) vxvarstyle = EQUAL; vxvarstyle = EQUAL;
else if (input->variable->atomstyle(vxvar)) vxvarstyle = ATOM; else if (input->variable->atomstyle(vxvar))
else error->all(FLERR,"Variable for fix move is invalid style"); vxvarstyle = ATOM;
else
error->all(FLERR, "Variable for fix move is invalid style");
} }
if (vyvarstr) { if (vyvarstr) {
vyvar = input->variable->find(vyvarstr); vyvar = input->variable->find(vyvarstr);
if (vyvar < 0) error->all(FLERR, if (vyvar < 0) error->all(FLERR, "Variable name for fix move does not exist");
"Variable name for fix move does not exist"); if (input->variable->equalstyle(vyvar))
if (input->variable->equalstyle(vyvar)) vyvarstyle = EQUAL; vyvarstyle = EQUAL;
else if (input->variable->atomstyle(vyvar)) vyvarstyle = ATOM; else if (input->variable->atomstyle(vyvar))
else error->all(FLERR,"Variable for fix move is invalid style"); vyvarstyle = ATOM;
else
error->all(FLERR, "Variable for fix move is invalid style");
} }
if (vzvarstr) { if (vzvarstr) {
vzvar = input->variable->find(vzvarstr); vzvar = input->variable->find(vzvarstr);
if (vzvar < 0) error->all(FLERR, if (vzvar < 0) error->all(FLERR, "Variable name for fix move does not exist");
"Variable name for fix move does not exist"); if (input->variable->equalstyle(vzvar))
if (input->variable->equalstyle(vzvar)) vzvarstyle = EQUAL; vzvarstyle = EQUAL;
else if (input->variable->atomstyle(vzvar)) vzvarstyle = ATOM; else if (input->variable->atomstyle(vzvar))
else error->all(FLERR,"Variable for fix move is invalid style"); vzvarstyle = ATOM;
else
error->all(FLERR, "Variable for fix move is invalid style");
} }
if (xvarstr && xvarstyle == ATOM) displaceflag = 1; if (xvarstr && xvarstyle == ATOM) displaceflag = 1;
@ -429,10 +485,14 @@ void FixMove::init()
maxatom = atom->nmax; maxatom = atom->nmax;
memory->destroy(displace); memory->destroy(displace);
memory->destroy(velocity); memory->destroy(velocity);
if (displaceflag) memory->create(displace,maxatom,3,"move:displace"); if (displaceflag)
else displace = nullptr; memory->create(displace, maxatom, 3, "move:displace");
if (velocityflag) memory->create(velocity,maxatom,3,"move:velocity"); else
else velocity = nullptr; 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; nlevels_respa = ((Respa *) update->integrate)->nlevels;
@ -707,6 +767,144 @@ void FixMove::initial_integrate(int /*vflag*/)
} }
} }
// 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);
}
}
// for variable: compute x,v from variables // for variable: compute x,v from variables
// NOTE: also allow for changes to extra attributes? // NOTE: also allow for changes to extra attributes?
// omega, angmom, theta, quat // omega, angmom, theta, quat
@ -733,28 +931,40 @@ void FixMove::initial_integrate(int /*vflag*/)
modify->clearstep_compute(); modify->clearstep_compute();
if (xvarstr) { if (xvarstr) {
if (xvarstyle == EQUAL) dx = input->variable->compute_equal(xvar); if (xvarstyle == EQUAL)
else input->variable->compute_atom(xvar,igroup,&displace[0][0],3,0); dx = input->variable->compute_equal(xvar);
else
input->variable->compute_atom(xvar, igroup, &displace[0][0], 3, 0);
} }
if (yvarstr) { if (yvarstr) {
if (yvarstyle == EQUAL) dy = input->variable->compute_equal(yvar); if (yvarstyle == EQUAL)
else input->variable->compute_atom(yvar,igroup,&displace[0][1],3,0); dy = input->variable->compute_equal(yvar);
else
input->variable->compute_atom(yvar, igroup, &displace[0][1], 3, 0);
} }
if (zvarstr) { if (zvarstr) {
if (zvarstyle == EQUAL) dz = input->variable->compute_equal(zvar); if (zvarstyle == EQUAL)
else input->variable->compute_atom(zvar,igroup,&displace[0][2],3,0); dz = input->variable->compute_equal(zvar);
else
input->variable->compute_atom(zvar, igroup, &displace[0][2], 3, 0);
} }
if (vxvarstr) { if (vxvarstr) {
if (vxvarstyle == EQUAL) vx = input->variable->compute_equal(vxvar); if (vxvarstyle == EQUAL)
else input->variable->compute_atom(vxvar,igroup,&velocity[0][0],3,0); vx = input->variable->compute_equal(vxvar);
else
input->variable->compute_atom(vxvar, igroup, &velocity[0][0], 3, 0);
} }
if (vyvarstr) { if (vyvarstr) {
if (vyvarstyle == EQUAL) vy = input->variable->compute_equal(vyvar); if (vyvarstyle == EQUAL)
else input->variable->compute_atom(vyvar,igroup,&velocity[0][1],3,0); vy = input->variable->compute_equal(vyvar);
else
input->variable->compute_atom(vyvar, igroup, &velocity[0][1], 3, 0);
} }
if (vzvarstr) { if (vzvarstr) {
if (vzvarstyle == EQUAL) vz = input->variable->compute_equal(vzvar); if (vzvarstyle == EQUAL)
else input->variable->compute_atom(vzvar,igroup,&velocity[0][2],3,0); vz = input->variable->compute_equal(vzvar);
else
input->variable->compute_atom(vzvar, igroup, &velocity[0][2], 3, 0);
} }
modify->addstep_compute(update->ntimestep + 1); modify->addstep_compute(update->ntimestep + 1);
@ -768,16 +978,24 @@ void FixMove::initial_integrate(int /*vflag*/)
xold[2] = x[i][2]; xold[2] = x[i][2];
if (xvarstr && vxvarstr) { if (xvarstr && vxvarstr) {
if (vxvarstyle == EQUAL) v[i][0] = vx; if (vxvarstyle == EQUAL)
else v[i][0] = velocity[i][0]; v[i][0] = vx;
if (xvarstyle == EQUAL) x[i][0] = xoriginal[i][0] + dx; else
else x[i][0] = xoriginal[i][0] + displace[i][0]; 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) { } else if (xvarstr) {
if (xvarstyle == EQUAL) x[i][0] = xoriginal[i][0] + dx; if (xvarstyle == EQUAL)
else x[i][0] = xoriginal[i][0] + displace[i][0]; x[i][0] = xoriginal[i][0] + dx;
else
x[i][0] = xoriginal[i][0] + displace[i][0];
} else if (vxvarstr) { } else if (vxvarstr) {
if (vxvarstyle == EQUAL) v[i][0] = vx; if (vxvarstyle == EQUAL)
else v[i][0] = velocity[i][0]; v[i][0] = vx;
else
v[i][0] = velocity[i][0];
x[i][0] += dtv * v[i][0]; x[i][0] += dtv * v[i][0];
} else { } else {
if (rmass) { if (rmass) {
@ -791,16 +1009,24 @@ void FixMove::initial_integrate(int /*vflag*/)
} }
if (yvarstr && vyvarstr) { if (yvarstr && vyvarstr) {
if (vyvarstyle == EQUAL) v[i][1] = vy; if (vyvarstyle == EQUAL)
else v[i][1] = velocity[i][1]; v[i][1] = vy;
if (yvarstyle == EQUAL) x[i][1] = xoriginal[i][1] + dy; else
else x[i][1] = xoriginal[i][1] + displace[i][1]; 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) { } else if (yvarstr) {
if (yvarstyle == EQUAL) x[i][1] = xoriginal[i][1] + dy; if (yvarstyle == EQUAL)
else x[i][1] = xoriginal[i][1] + displace[i][1]; x[i][1] = xoriginal[i][1] + dy;
else
x[i][1] = xoriginal[i][1] + displace[i][1];
} else if (vyvarstr) { } else if (vyvarstr) {
if (vyvarstyle == EQUAL) v[i][1] = vy; if (vyvarstyle == EQUAL)
else v[i][1] = velocity[i][1]; v[i][1] = vy;
else
v[i][1] = velocity[i][1];
x[i][1] += dtv * v[i][1]; x[i][1] += dtv * v[i][1];
} else { } else {
if (rmass) { if (rmass) {
@ -814,16 +1040,24 @@ void FixMove::initial_integrate(int /*vflag*/)
} }
if (zvarstr && vzvarstr) { if (zvarstr && vzvarstr) {
if (vzvarstyle == EQUAL) v[i][2] = vz; if (vzvarstyle == EQUAL)
else v[i][2] = velocity[i][2]; v[i][2] = vz;
if (zvarstyle == EQUAL) x[i][2] = xoriginal[i][2] + dz; else
else x[i][2] = xoriginal[i][2] + displace[i][2]; 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) { } else if (zvarstr) {
if (zvarstyle == EQUAL) x[i][2] = xoriginal[i][2] + dz; if (zvarstyle == EQUAL)
else x[i][2] = xoriginal[i][2] + displace[i][2]; x[i][2] = xoriginal[i][2] + dz;
else
x[i][2] = xoriginal[i][2] + displace[i][2];
} else if (vzvarstr) { } else if (vzvarstr) {
if (vzvarstyle == EQUAL) v[i][2] = vz; if (vzvarstyle == EQUAL)
else v[i][2] = velocity[i][2]; v[i][2] = vz;
else
v[i][2] = velocity[i][2];
x[i][2] += dtv * v[i][2]; x[i][2] += dtv * v[i][2];
} else { } else {
if (rmass) { if (rmass) {
@ -851,22 +1085,40 @@ void FixMove::final_integrate()
double dtfm; double dtfm;
int xflag = 1; int xflag = 1;
if (mstyle == LINEAR && vxflag) xflag = 0; if (mstyle == LINEAR && vxflag)
else if (mstyle == WIGGLE && axflag) xflag = 0; xflag = 0;
else if (mstyle == ROTATE) xflag = 0; else if (mstyle == WIGGLE && axflag)
else if (mstyle == VARIABLE && (xvarstr || vxvarstr)) xflag = 0; 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; int yflag = 1;
if (mstyle == LINEAR && vyflag) yflag = 0; if (mstyle == LINEAR && vyflag)
else if (mstyle == WIGGLE && ayflag) yflag = 0; yflag = 0;
else if (mstyle == ROTATE) yflag = 0; else if (mstyle == WIGGLE && ayflag)
else if (mstyle == VARIABLE && (yvarstr || vyvarstr)) yflag = 0; 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; int zflag = 1;
if (mstyle == LINEAR && vzflag) zflag = 0; if (mstyle == LINEAR && vzflag)
else if (mstyle == WIGGLE && azflag) zflag = 0; zflag = 0;
else if (mstyle == ROTATE) zflag = 0; else if (mstyle == WIGGLE && azflag)
else if (mstyle == VARIABLE && (zvarstr || vzvarstr)) zflag = 0; 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 **v = atom->v;
double **f = atom->f; double **f = atom->f;
@ -1034,8 +1286,7 @@ void FixMove::set_arrays(int i)
// backup particle to time_origin // backup particle to time_origin
if (mstyle == VARIABLE) if (mstyle == VARIABLE) error->all(FLERR, "Cannot add atoms to fix move 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; double delta = (update->ntimestep - time_origin) * update->dt;
@ -1044,12 +1295,14 @@ void FixMove::set_arrays(int i)
if (vxflag) xoriginal[i][0] -= vx * delta; if (vxflag) xoriginal[i][0] -= vx * delta;
if (vyflag) xoriginal[i][1] -= vy * delta; if (vyflag) xoriginal[i][1] -= vy * delta;
if (vzflag) xoriginal[i][2] -= vz * delta; if (vzflag) xoriginal[i][2] -= vz * delta;
} else if (mstyle == WIGGLE) { } else if (mstyle == WIGGLE) {
double arg = omega_rotate * delta; double arg = omega_rotate * delta;
double sine = sin(arg); double sine = sin(arg);
if (axflag) xoriginal[i][0] -= ax * sine; if (axflag) xoriginal[i][0] -= ax * sine;
if (ayflag) xoriginal[i][1] -= ay * sine; if (ayflag) xoriginal[i][1] -= ay * sine;
if (azflag) xoriginal[i][2] -= az * sine; if (azflag) xoriginal[i][2] -= az * sine;
} else if (mstyle == ROTATE) { } else if (mstyle == ROTATE) {
double a[3], b[3], c[3], d[3], disp[3], ddotr; double a[3], b[3], c[3], d[3], disp[3], ddotr;
double arg = -omega_rotate * delta; double arg = -omega_rotate * delta;
@ -1079,6 +1332,63 @@ void FixMove::set_arrays(int i)
// set theta and quat extra attributes affected by rotation // 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
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) { if (extra_flag) {
// theta for lines // theta for lines