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

enable fix rigid to support 2d enforcement internally

Browse Source
This commit is contained in:
Steve Plimpton
2023-04-03 17:21:11 -06:00
parent c9b6fad6a0
commit 8d523c0604
8 changed files with 361 additions and 369 deletions
Download Patch File Download Diff File Expand all files Collapse all files

314
src/RIGID/fix_rigid.cpp
View File

@ -68,7 +68,6 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
thermo_virial = 1; thermo_virial = 1;
create_attribute = 1; create_attribute = 1;
dof_flag = 1; dof_flag = 1;
enforce2d_flag = 1;
centroidstressflag = CENTROID_NOTAVAIL; centroidstressflag = CENTROID_NOTAVAIL;
MPI_Comm_rank(world,&me); MPI_Comm_rank(world,&me);
@ -77,6 +76,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
// perform initial allocation of atom-based arrays // perform initial allocation of atom-based arrays
// register with Atom class // register with Atom class
dimension = domain->dimension;
extended = orientflag = dorientflag = 0; extended = orientflag = dorientflag = 0;
body = nullptr; body = nullptr;
xcmimage = nullptr; xcmimage = nullptr;
@ -301,7 +301,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
for (i = 0; i < nbody; i++) { for (i = 0; i < nbody; i++) {
fflag[i][0] = fflag[i][1] = fflag[i][2] = 1.0; fflag[i][0] = fflag[i][1] = fflag[i][2] = 1.0;
tflag[i][0] = tflag[i][1] = tflag[i][2] = 1.0; tflag[i][0] = tflag[i][1] = tflag[i][2] = 1.0;
if (domain->dimension == 2) fflag[i][2] = tflag[i][0] = tflag[i][1] = 0.0; if (dimension == 2) fflag[i][2] = tflag[i][0] = tflag[i][1] = 0.0;
} }
// number of linear rigid bodies is counted later // number of linear rigid bodies is counted later
@ -328,7 +328,6 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
pcouple = NONE; pcouple = NONE;
pstyle = ANISO; pstyle = ANISO;
dimension = domain->dimension;
for (i = 0; i < 3; i++) { for (i = 0; i < 3; i++) {
p_start[i] = p_stop[i] = p_period[i] = 0.0; p_start[i] = p_stop[i] = p_period[i] = 0.0;
@ -353,7 +352,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
else if (strcmp(arg[iarg+4],"on") == 0) zflag = 1.0; else if (strcmp(arg[iarg+4],"on") == 0) zflag = 1.0;
else error->all(FLERR,"Illegal fix rigid command"); else error->all(FLERR,"Illegal fix rigid command");
if (domain->dimension == 2 && zflag == 1.0) if (dimension == 2 && zflag == 1.0)
error->all(FLERR,"Fix rigid z force cannot be on for 2d simulation"); error->all(FLERR,"Fix rigid z force cannot be on for 2d simulation");
int count = 0; int count = 0;
@ -384,7 +383,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
else if (strcmp(arg[iarg+4],"on") == 0) zflag = 1.0; else if (strcmp(arg[iarg+4],"on") == 0) zflag = 1.0;
else error->all(FLERR,"Illegal fix rigid command"); else error->all(FLERR,"Illegal fix rigid command");
if (domain->dimension == 2 && (xflag == 1.0 || yflag == 1.0)) if (dimension == 2 && (xflag == 1.0 || yflag == 1.0))
error->all(FLERR,"Fix rigid xy torque cannot be on for 2d simulation"); error->all(FLERR,"Fix rigid xy torque cannot be on for 2d simulation");
int count = 0; int count = 0;
@ -434,7 +433,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
utils::numeric(FLERR,arg[iarg+3],false,lmp); utils::numeric(FLERR,arg[iarg+3],false,lmp);
p_flag[0] = p_flag[1] = p_flag[2] = 1; p_flag[0] = p_flag[1] = p_flag[2] = 1;
if (dimension == 2) { if (dimension == 2) {
p_start[2] = p_stop[2] = p_period[2] = 0.0; p_start[2] = p_stop[2] = p_period[2] = 0.0;
p_flag[2] = 0; p_flag[2] = 0;
} }
iarg += 4; iarg += 4;
@ -450,7 +449,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
p_flag[0] = p_flag[1] = p_flag[2] = 1; p_flag[0] = p_flag[1] = p_flag[2] = 1;
if (dimension == 2) { if (dimension == 2) {
p_start[2] = p_stop[2] = p_period[2] = 0.0; p_start[2] = p_stop[2] = p_period[2] = 0.0;
p_flag[2] = 0; p_flag[2] = 0;
} }
iarg += 4; iarg += 4;
@ -856,6 +855,10 @@ void FixRigid::setup(int vflag)
torque[ibody][2] = all[ibody][2]; torque[ibody][2] = all[ibody][2];
} }
// enforce 2d body forces and torques
if (dimension == 2) enforce2d();
// zero langextra in case Langevin thermostat not used // zero langextra in case Langevin thermostat not used
// no point to calling post_force() here since langextra // no point to calling post_force() here since langextra
// is only added to fcm/torque in final_integrate() // is only added to fcm/torque in final_integrate()
@ -936,6 +939,131 @@ void FixRigid::initial_integrate(int vflag)
set_xv(); set_xv();
} }
/* ----------------------------------------------------------------------
remap xcm of each rigid body back into periodic simulation box
done during pre_neighbor so will be after call to pbc()
and after fix_deform::pre_exchange() may have flipped box
use domain->remap() in case xcm is far away from box
due to first-time definition of rigid body in setup_bodies_static()
or due to box flip
also adjust imagebody = rigid body image flags, due to xcm remap
also reset body xcmimage flags of all atoms in bodies
xcmimage flags are relative to xcm so that body can be unwrapped
if don't do this, would need xcm to move with true image flags
then a body could end up very far away from box
set_xv() will then compute huge displacements every step to
reset coords of all body atoms to be back inside the box,
ditto for triclinic box flip, which causes numeric problems
------------------------------------------------------------------------- */
void FixRigid::pre_neighbor()
{
for (int ibody = 0; ibody < nbody; ibody++)
domain->remap(xcm[ibody],imagebody[ibody]);
image_shift();
}
/* ---------------------------------------------------------------------- */
void FixRigid::post_force(int /*vflag*/)
{
if (langflag) apply_langevin_thermostat();
if (earlyflag) compute_forces_and_torques();
}
/* ---------------------------------------------------------------------- */
void FixRigid::final_integrate()
{
int ibody;
double dtfm;
// compute forces and torques (after all post_force contributions)
// if 2d model, enforce2d() on body forces/torques
if (!earlyflag) compute_forces_and_torques();
if (dimension == 2) enforce2d();
// update vcm and angmom
// fflag,tflag = 0 for some dimensions in 2d
for (ibody = 0; ibody < nbody; ibody++) {
// update vcm by 1/2 step
dtfm = dtf / masstotal[ibody];
vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0];
vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1];
vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2];
// update angular momentum by 1/2 step
angmom[ibody][0] += dtf * torque[ibody][0] * tflag[ibody][0];
angmom[ibody][1] += dtf * torque[ibody][1] * tflag[ibody][1];
angmom[ibody][2] += dtf * torque[ibody][2] * tflag[ibody][2];
MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
ez_space[ibody],inertia[ibody],omega[ibody]);
}
// set velocity/rotation of atoms in rigid bodies
// virial is already setup from initial_integrate
set_v();
}
/* ---------------------------------------------------------------------- */
void FixRigid::initial_integrate_respa(int vflag, int ilevel, int /*iloop*/)
{
dtv = step_respa[ilevel];
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
dtq = 0.5 * step_respa[ilevel];
if (ilevel == 0) initial_integrate(vflag);
else final_integrate();
}
/* ---------------------------------------------------------------------- */
void FixRigid::final_integrate_respa(int ilevel, int /*iloop*/)
{
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
final_integrate();
}
/* ----------------------------------------------------------------------
reset body xcmimage flags of atoms in bodies
xcmimage flags are relative to xcm so that body can be unwrapped
xcmimage = true image flag - imagebody flag
------------------------------------------------------------------------- */
void FixRigid::image_shift()
{
int ibody;
imageint tdim,bdim,xdim[3];
imageint *image = atom->image;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
tdim = image[i] & IMGMASK;
bdim = imagebody[ibody] & IMGMASK;
xdim[0] = IMGMAX + tdim - bdim;
tdim = (image[i] >> IMGBITS) & IMGMASK;
bdim = (imagebody[ibody] >> IMGBITS) & IMGMASK;
xdim[1] = IMGMAX + tdim - bdim;
tdim = image[i] >> IMG2BITS;
bdim = imagebody[ibody] >> IMG2BITS;
xdim[2] = IMGMAX + tdim - bdim;
xcmimage[i] = (xdim[2] << IMG2BITS) | (xdim[1] << IMGBITS) | xdim[0];
}
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
apply Langevin thermostat to all 6 DOF of rigid bodies apply Langevin thermostat to all 6 DOF of rigid bodies
computed by proc 0, broadcast to other procs computed by proc 0, broadcast to other procs
@ -991,31 +1119,6 @@ void FixRigid::apply_langevin_thermostat()
MPI_Bcast(&langextra[0][0],6*nbody,MPI_DOUBLE,0,world); MPI_Bcast(&langextra[0][0],6*nbody,MPI_DOUBLE,0,world);
} }
/* ----------------------------------------------------------------------
called from FixEnforce2d post_force() for 2d problems
zero all body values that should be zero for 2d model
------------------------------------------------------------------------- */
void FixRigid::enforce2d()
{
for (int ibody = 0; ibody < nbody; ibody++) {
xcm[ibody][2] = 0.0;
vcm[ibody][2] = 0.0;
fcm[ibody][2] = 0.0;
torque[ibody][0] = 0.0;
torque[ibody][1] = 0.0;
angmom[ibody][0] = 0.0;
angmom[ibody][1] = 0.0;
omega[ibody][0] = 0.0;
omega[ibody][1] = 0.0;
if (langflag && langextra) {
langextra[ibody][2] = 0.0;
langextra[ibody][3] = 0.0;
langextra[ibody][4] = 0.0;
}
}
}
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixRigid::compute_forces_and_torques() void FixRigid::compute_forces_and_torques()
@ -1093,124 +1196,28 @@ void FixRigid::compute_forces_and_torques()
} }
} }
/* ---------------------------------------------------------------------- */
void FixRigid::post_force(int /*vflag*/)
{
if (langflag) apply_langevin_thermostat();
if (earlyflag) compute_forces_and_torques();
}
/* ---------------------------------------------------------------------- */
void FixRigid::final_integrate()
{
int ibody;
double dtfm;
if (!earlyflag) compute_forces_and_torques();
// update vcm and angmom
// fflag,tflag = 0 for some dimensions in 2d
for (ibody = 0; ibody < nbody; ibody++) {
// update vcm by 1/2 step
dtfm = dtf / masstotal[ibody];
vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0];
vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1];
vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2];
// update angular momentum by 1/2 step
angmom[ibody][0] += dtf * torque[ibody][0] * tflag[ibody][0];
angmom[ibody][1] += dtf * torque[ibody][1] * tflag[ibody][1];
angmom[ibody][2] += dtf * torque[ibody][2] * tflag[ibody][2];
MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
ez_space[ibody],inertia[ibody],omega[ibody]);
}
// set velocity/rotation of atoms in rigid bodies
// virial is already setup from initial_integrate
set_v();
}
/* ---------------------------------------------------------------------- */
void FixRigid::initial_integrate_respa(int vflag, int ilevel, int /*iloop*/)
{
dtv = step_respa[ilevel];
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
dtq = 0.5 * step_respa[ilevel];
if (ilevel == 0) initial_integrate(vflag);
else final_integrate();
}
/* ---------------------------------------------------------------------- */
void FixRigid::final_integrate_respa(int ilevel, int /*iloop*/)
{
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
final_integrate();
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
remap xcm of each rigid body back into periodic simulation box called from FixEnforce2d post_force() for 2d problems
done during pre_neighbor so will be after call to pbc() zero all body values that should be zero for 2d model
and after fix_deform::pre_exchange() may have flipped box
use domain->remap() in case xcm is far away from box
due to first-time definition of rigid body in setup_bodies_static()
or due to box flip
also adjust imagebody = rigid body image flags, due to xcm remap
also reset body xcmimage flags of all atoms in bodies
xcmimage flags are relative to xcm so that body can be unwrapped
if don't do this, would need xcm to move with true image flags
then a body could end up very far away from box
set_xv() will then compute huge displacements every step to
reset coords of all body atoms to be back inside the box,
ditto for triclinic box flip, which causes numeric problems
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void FixRigid::pre_neighbor() void FixRigid::enforce2d()
{ {
for (int ibody = 0; ibody < nbody; ibody++) for (int ibody = 0; ibody < nbody; ibody++) {
domain->remap(xcm[ibody],imagebody[ibody]); xcm[ibody][2] = 0.0;
image_shift(); vcm[ibody][2] = 0.0;
} fcm[ibody][2] = 0.0;
torque[ibody][0] = 0.0;
/* ---------------------------------------------------------------------- torque[ibody][1] = 0.0;
reset body xcmimage flags of atoms in bodies angmom[ibody][0] = 0.0;
xcmimage flags are relative to xcm so that body can be unwrapped angmom[ibody][1] = 0.0;
xcmimage = true image flag - imagebody flag omega[ibody][0] = 0.0;
------------------------------------------------------------------------- */ omega[ibody][1] = 0.0;
if (langflag && langextra) {
void FixRigid::image_shift() langextra[ibody][2] = 0.0;
{ langextra[ibody][3] = 0.0;
int ibody; langextra[ibody][4] = 0.0;
imageint tdim,bdim,xdim[3]; }
imageint *image = atom->image;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
tdim = image[i] & IMGMASK;
bdim = imagebody[ibody] & IMGMASK;
xdim[0] = IMGMAX + tdim - bdim;
tdim = (image[i] >> IMGBITS) & IMGMASK;
bdim = (imagebody[ibody] >> IMGBITS) & IMGMASK;
xdim[1] = IMGMAX + tdim - bdim;
tdim = image[i] >> IMG2BITS;
bdim = imagebody[ibody] >> IMG2BITS;
xdim[2] = IMGMAX + tdim - bdim;
xcmimage[i] = (xdim[2] << IMG2BITS) | (xdim[1] << IMGBITS) | xdim[0];
} }
} }
@ -1280,7 +1287,7 @@ int FixRigid::dof(int tgroup)
int n = 0; int n = 0;
nlinear = 0; nlinear = 0;
if (domain->dimension == 3) { if (dimension == 3) {
for (int ibody = 0; ibody < nbody; ibody++) for (int ibody = 0; ibody < nbody; ibody++)
if (nall[ibody]+mall[ibody] == nrigid[ibody]) { if (nall[ibody]+mall[ibody] == nrigid[ibody]) {
n += 3*nall[ibody] + 6*mall[ibody] - 6; n += 3*nall[ibody] + 6*mall[ibody] - 6;
@ -1290,7 +1297,7 @@ int FixRigid::dof(int tgroup)
nlinear++; nlinear++;
} }
} }
} else if (domain->dimension == 2) { } else if (dimension == 2) {
for (int ibody = 0; ibody < nbody; ibody++) for (int ibody = 0; ibody < nbody; ibody++)
if (nall[ibody]+mall[ibody] == nrigid[ibody]) if (nall[ibody]+mall[ibody] == nrigid[ibody])
n += 2*nall[ibody] + 3*mall[ibody] - 3; n += 2*nall[ibody] + 3*mall[ibody] - 3;
@ -1382,6 +1389,7 @@ void FixRigid::set_xv()
// x = displacement from center-of-mass, based on body orientation // x = displacement from center-of-mass, based on body orientation
// v = vcm + omega around center-of-mass // v = vcm + omega around center-of-mass
// enforce 2d x and v
MathExtra::matvec(ex_space[ibody],ey_space[ibody], MathExtra::matvec(ex_space[ibody],ey_space[ibody],
ez_space[ibody],displace[i],x[i]); ez_space[ibody],displace[i],x[i]);
@ -1390,6 +1398,11 @@ void FixRigid::set_xv()
v[i][1] = omega[ibody][2]*x[i][0] - omega[ibody][0]*x[i][2] + vcm[ibody][1]; v[i][1] = omega[ibody][2]*x[i][0] - omega[ibody][0]*x[i][2] + vcm[ibody][1];
v[i][2] = omega[ibody][0]*x[i][1] - omega[ibody][1]*x[i][0] + vcm[ibody][2]; v[i][2] = omega[ibody][0]*x[i][1] - omega[ibody][1]*x[i][0] + vcm[ibody][2];
if (dimension == 2) {
x[i][2] = 0.0;
v[i][2] = 0.0;
}
// add center of mass to displacement // add center of mass to displacement
// map back into periodic box via xbox,ybox,zbox // map back into periodic box via xbox,ybox,zbox
// for triclinic, add in box tilt factors as well // for triclinic, add in box tilt factors as well
@ -1541,10 +1554,15 @@ void FixRigid::set_v()
v2 = v[i][2]; v2 = v[i][2];
} }
// compute new v
// enforce 2d v
v[i][0] = omega[ibody][1]*delta[2] - omega[ibody][2]*delta[1] + vcm[ibody][0]; v[i][0] = omega[ibody][1]*delta[2] - omega[ibody][2]*delta[1] + vcm[ibody][0];
v[i][1] = omega[ibody][2]*delta[0] - omega[ibody][0]*delta[2] + vcm[ibody][1]; v[i][1] = omega[ibody][2]*delta[0] - omega[ibody][0]*delta[2] + vcm[ibody][1];
v[i][2] = omega[ibody][0]*delta[1] - omega[ibody][1]*delta[0] + vcm[ibody][2]; v[i][2] = omega[ibody][0]*delta[1] - omega[ibody][1]*delta[0] + vcm[ibody][2];
if (dimension == 2) v[i][2] = 0.0;
// virial = unwrapped coords dotted into body constraint force // virial = unwrapped coords dotted into body constraint force
// body constraint force = implied force due to v change minus f external // body constraint force = implied force due to v change minus f external
// assume f does not include forces internal to body // assume f does not include forces internal to body

2
src/RIGID/fix_rigid.h
View File

@ -50,7 +50,6 @@ class FixRigid : public Fix {
void pre_neighbor() override; void pre_neighbor() override;
int dof(int) override; int dof(int) override;
void deform(int) override; void deform(int) override;
void enforce2d() override;
void reset_dt() override; void reset_dt() override;
void zero_momentum() override; void zero_momentum() override;
void zero_rotation() override; void zero_rotation() override;
@ -147,6 +146,7 @@ class FixRigid : public Fix {
void setup_bodies_dynamic(); void setup_bodies_dynamic();
void apply_langevin_thermostat(); void apply_langevin_thermostat();
void compute_forces_and_torques(); void compute_forces_and_torques();
void enforce2d();
void readfile(int, double *, double **, double **, double **, imageint *, int *); void readfile(int, double *, double **, double **, double **, imageint *, int *);
}; };

368
src/RIGID/fix_rigid_small.cpp
View File

@ -72,7 +72,6 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
thermo_virial = 1; thermo_virial = 1;
create_attribute = 1; create_attribute = 1;
dof_flag = 1; dof_flag = 1;
enforce2d_flag = 1;
stores_ids = 1; stores_ids = 1;
centroidstressflag = CENTROID_AVAIL; centroidstressflag = CENTROID_AVAIL;
@ -82,6 +81,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
// perform initial allocation of atom-based arrays // perform initial allocation of atom-based arrays
// register with Atom class // register with Atom class
dimension = domain->dimension;
extended = orientflag = dorientflag = customflag = 0; extended = orientflag = dorientflag = customflag = 0;
bodyown = nullptr; bodyown = nullptr;
bodytag = nullptr; bodytag = nullptr;
@ -258,8 +258,8 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
p_period[0] = p_period[1] = p_period[2] = p_period[0] = p_period[1] = p_period[2] =
utils::numeric(FLERR,arg[iarg+3],false,lmp); utils::numeric(FLERR,arg[iarg+3],false,lmp);
p_flag[0] = p_flag[1] = p_flag[2] = 1; p_flag[0] = p_flag[1] = p_flag[2] = 1;
if (domain->dimension == 2) { if (dimension == 2) {
p_start[2] = p_stop[2] = p_period[2] = 0.0; p_start[2] = p_stop[2] = p_period[2] = 0.0;
p_flag[2] = 0; p_flag[2] = 0;
} }
iarg += 4; iarg += 4;
@ -273,9 +273,9 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
p_period[0] = p_period[1] = p_period[2] = p_period[0] = p_period[1] = p_period[2] =
utils::numeric(FLERR,arg[iarg+3],false,lmp); utils::numeric(FLERR,arg[iarg+3],false,lmp);
p_flag[0] = p_flag[1] = p_flag[2] = 1; p_flag[0] = p_flag[1] = p_flag[2] = 1;
if (domain->dimension == 2) { if (dimension == 2) {
p_start[2] = p_stop[2] = p_period[2] = 0.0; p_start[2] = p_stop[2] = p_period[2] = 0.0;
p_flag[2] = 0; p_flag[2] = 0;
} }
iarg += 4; iarg += 4;
@ -383,7 +383,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
for (i = 0; i < 3; i++) for (i = 0; i < 3; i++)
if (p_flag[i]) pstat_flag = 1; if (p_flag[i]) pstat_flag = 1;
if (pcouple == XYZ || (domain->dimension == 2 && pcouple == XY)) pstyle = ISO; if (pcouple == XYZ || (dimension == 2 && pcouple == XY)) pstyle = ISO;
else pstyle = ANISO; else pstyle = ANISO;
// create rigid bodies based on molecule or custom ID // create rigid bodies based on molecule or custom ID
@ -690,6 +690,10 @@ void FixRigidSmall::setup(int vflag)
} }
} }
// enforce 2d body forces and torques
if (dimension == 2) enforce2d();
// reverse communicate fcm, torque of all bodies // reverse communicate fcm, torque of all bodies
commflag = FORCE_TORQUE; commflag = FORCE_TORQUE;
@ -780,6 +784,151 @@ void FixRigidSmall::initial_integrate(int vflag)
set_xv(); set_xv();
} }
/* ----------------------------------------------------------------------
remap xcm of each rigid body back into periodic simulation box
done during pre_neighbor so will be after call to pbc()
and after fix_deform::pre_exchange() may have flipped box
use domain->remap() in case xcm is far away from box
due to first-time definition of rigid body in setup_bodies_static()
or due to box flip
also adjust imagebody = rigid body image flags, due to xcm remap
then communicate bodies so other procs will know of changes to body xcm
then adjust xcmimage flags of all atoms in bodies via image_shift()
for two effects
(1) change in true image flags due to pbc() call during exchange
(2) change in imagebody due to xcm remap
xcmimage flags are always -1,0,-1 so that body can be unwrapped
around in-box xcm and stay close to simulation box
if just inferred unwrapped from atom image flags,
then a body could end up very far away
when unwrapped by true image flags
then set_xv() will compute huge displacements every step to reset coords of
all the body atoms to be back inside the box, ditto for triclinic box flip
note: so just want to avoid that numeric problem?
------------------------------------------------------------------------- */
void FixRigidSmall::pre_neighbor()
{
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
domain->remap(b->xcm,b->image);
}
nghost_body = 0;
commflag = FULL_BODY;
comm->forward_comm(this);
reset_atom2body();
//check(4);
image_shift();
}
/* ---------------------------------------------------------------------- */
void FixRigidSmall::post_force(int /*vflag*/)
{
if (langflag) apply_langevin_thermostat();
if (earlyflag) compute_forces_and_torques();
}
/* ---------------------------------------------------------------------- */
void FixRigidSmall::final_integrate()
{
double dtfm;
//check(3);
// compute forces and torques (after all post_force contributions)
// if 2d model, enforce2d() on body forces/torques
if (!earlyflag) compute_forces_and_torques();
if (dimension == 2) enforce2d();
// update vcm and angmom, recompute omega
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
// update vcm by 1/2 step
dtfm = dtf / b->mass;
b->vcm[0] += dtfm * b->fcm[0];
b->vcm[1] += dtfm * b->fcm[1];
b->vcm[2] += dtfm * b->fcm[2];
// update angular momentum by 1/2 step
b->angmom[0] += dtf * b->torque[0];
b->angmom[1] += dtf * b->torque[1];
b->angmom[2] += dtf * b->torque[2];
MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space,
b->ez_space,b->inertia,b->omega);
}
// forward communicate updated info of all bodies
commflag = FINAL;
comm->forward_comm(this,10);
// set velocity/rotation of atoms in rigid bodies
// virial is already setup from initial_integrate
set_v();
}
/* ---------------------------------------------------------------------- */
void FixRigidSmall::initial_integrate_respa(int vflag, int ilevel, int /*iloop*/)
{
dtv = step_respa[ilevel];
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
dtq = 0.5 * step_respa[ilevel];
if (ilevel == 0) initial_integrate(vflag);
else final_integrate();
}
/* ---------------------------------------------------------------------- */
void FixRigidSmall::final_integrate_respa(int ilevel, int /*iloop*/)
{
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
final_integrate();
}
/* ----------------------------------------------------------------------
reset body xcmimage flags of atoms in bodies
xcmimage flags are relative to xcm so that body can be unwrapped
xcmimage = true image flag - imagebody flag
------------------------------------------------------------------------- */
void FixRigidSmall::image_shift()
{
imageint tdim,bdim,xdim[3];
imageint *image = atom->image;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]];
tdim = image[i] & IMGMASK;
bdim = b->image & IMGMASK;
xdim[0] = IMGMAX + tdim - bdim;
tdim = (image[i] >> IMGBITS) & IMGMASK;
bdim = (b->image >> IMGBITS) & IMGMASK;
xdim[1] = IMGMAX + tdim - bdim;
tdim = image[i] >> IMG2BITS;
bdim = b->image >> IMG2BITS;
xdim[2] = IMGMAX + tdim - bdim;
xcmimage[i] = (xdim[2] << IMG2BITS) | (xdim[1] << IMGBITS) | xdim[0];
}
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
apply Langevin thermostat to all 6 DOF of rigid bodies I own apply Langevin thermostat to all 6 DOF of rigid bodies I own
unlike fix langevin, this stores extra force in extra arrays, unlike fix langevin, this stores extra force in extra arrays,
@ -845,52 +994,9 @@ void FixRigidSmall::apply_langevin_thermostat()
// convert langevin torques from body frame back to space frame // convert langevin torques from body frame back to space frame
MathExtra::matvec(ex_space,ey_space,ez_space,tbody,&langextra[ibody][3]); MathExtra::matvec(ex_space,ey_space,ez_space,tbody,&langextra[ibody][3]);
// enforce 2d motion
if (domain->dimension == 2)
langextra[ibody][2] = langextra[ibody][3] = langextra[ibody][4] = 0.0;
} }
} }
/* ----------------------------------------------------------------------
called from FixEnforce post_force() for 2d problems
zero all body values that should be zero for 2d model
------------------------------------------------------------------------- */
void FixRigidSmall::enforce2d()
{
Body *b;
for (int ibody = 0; ibody < nlocal_body; ibody++) {
b = &body[ibody];
b->xcm[2] = 0.0;
b->vcm[2] = 0.0;
b->fcm[2] = 0.0;
b->xgc[2] = 0.0;
b->torque[0] = 0.0;
b->torque[1] = 0.0;
b->angmom[0] = 0.0;
b->angmom[1] = 0.0;
b->omega[0] = 0.0;
b->omega[1] = 0.0;
if (langflag && langextra) {
langextra[ibody][2] = 0.0;
langextra[ibody][3] = 0.0;
langextra[ibody][4] = 0.0;
}
}
}
/* ---------------------------------------------------------------------- */
void FixRigidSmall::post_force(int /*vflag*/)
{
if (langflag) apply_langevin_thermostat();
if (earlyflag) compute_forces_and_torques();
}
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixRigidSmall::compute_forces_and_torques() void FixRigidSmall::compute_forces_and_torques()
@ -988,136 +1094,32 @@ void FixRigidSmall::compute_forces_and_torques()
} }
} }
/* ---------------------------------------------------------------------- */
void FixRigidSmall::final_integrate()
{
double dtfm;
//check(3);
if (!earlyflag) compute_forces_and_torques();
// update vcm and angmom, recompute omega
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
// update vcm by 1/2 step
dtfm = dtf / b->mass;
b->vcm[0] += dtfm * b->fcm[0];
b->vcm[1] += dtfm * b->fcm[1];
b->vcm[2] += dtfm * b->fcm[2];
// update angular momentum by 1/2 step
b->angmom[0] += dtf * b->torque[0];
b->angmom[1] += dtf * b->torque[1];
b->angmom[2] += dtf * b->torque[2];
MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space,
b->ez_space,b->inertia,b->omega);
}
// forward communicate updated info of all bodies
commflag = FINAL;
comm->forward_comm(this,10);
// set velocity/rotation of atoms in rigid bodies
// virial is already setup from initial_integrate
set_v();
}
/* ---------------------------------------------------------------------- */
void FixRigidSmall::initial_integrate_respa(int vflag, int ilevel, int /*iloop*/)
{
dtv = step_respa[ilevel];
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
dtq = 0.5 * step_respa[ilevel];
if (ilevel == 0) initial_integrate(vflag);
else final_integrate();
}
/* ---------------------------------------------------------------------- */
void FixRigidSmall::final_integrate_respa(int ilevel, int /*iloop*/)
{
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
final_integrate();
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
remap xcm of each rigid body back into periodic simulation box called from FixEnforce post_force() for 2d problems
done during pre_neighbor so will be after call to pbc() zero all body values that should be zero for 2d model
and after fix_deform::pre_exchange() may have flipped box
use domain->remap() in case xcm is far away from box
due to first-time definition of rigid body in setup_bodies_static()
or due to box flip
also adjust imagebody = rigid body image flags, due to xcm remap
then communicate bodies so other procs will know of changes to body xcm
then adjust xcmimage flags of all atoms in bodies via image_shift()
for two effects
(1) change in true image flags due to pbc() call during exchange
(2) change in imagebody due to xcm remap
xcmimage flags are always -1,0,-1 so that body can be unwrapped
around in-box xcm and stay close to simulation box
if just inferred unwrapped from atom image flags,
then a body could end up very far away
when unwrapped by true image flags
then set_xv() will compute huge displacements every step to reset coords of
all the body atoms to be back inside the box, ditto for triclinic box flip
note: so just want to avoid that numeric problem?
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void FixRigidSmall::pre_neighbor() void FixRigidSmall::enforce2d()
{ {
Body *b;
for (int ibody = 0; ibody < nlocal_body; ibody++) { for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody]; b = &body[ibody];
domain->remap(b->xcm,b->image); b->xcm[2] = 0.0;
} b->vcm[2] = 0.0;
b->fcm[2] = 0.0;
nghost_body = 0; b->xgc[2] = 0.0;
commflag = FULL_BODY; b->torque[0] = 0.0;
comm->forward_comm(this); b->torque[1] = 0.0;
reset_atom2body(); b->angmom[0] = 0.0;
//check(4); b->angmom[1] = 0.0;
b->omega[0] = 0.0;
image_shift(); b->omega[1] = 0.0;
} if (langflag && langextra) {
langextra[ibody][2] = 0.0;
/* ---------------------------------------------------------------------- langextra[ibody][3] = 0.0;
reset body xcmimage flags of atoms in bodies langextra[ibody][4] = 0.0;
xcmimage flags are relative to xcm so that body can be unwrapped }
xcmimage = true image flag - imagebody flag
------------------------------------------------------------------------- */
void FixRigidSmall::image_shift()
{
imageint tdim,bdim,xdim[3];
imageint *image = atom->image;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]];
tdim = image[i] & IMGMASK;
bdim = b->image & IMGMASK;
xdim[0] = IMGMAX + tdim - bdim;
tdim = (image[i] >> IMGBITS) & IMGMASK;
bdim = (b->image >> IMGBITS) & IMGMASK;
xdim[1] = IMGMAX + tdim - bdim;
tdim = image[i] >> IMG2BITS;
bdim = b->image >> IMG2BITS;
xdim[2] = IMGMAX + tdim - bdim;
xcmimage[i] = (xdim[2] << IMG2BITS) | (xdim[1] << IMGBITS) | xdim[0];
} }
} }
@ -1200,7 +1202,7 @@ int FixRigidSmall::dof(int tgroup)
int n = 0; int n = 0;
nlinear = 0; nlinear = 0;
if (domain->dimension == 3) { if (dimension == 3) {
for (int ibody = 0; ibody < nlocal_body; ibody++) { for (int ibody = 0; ibody < nlocal_body; ibody++) {
if (counts[ibody][0]+counts[ibody][1] == counts[ibody][2]) { if (counts[ibody][0]+counts[ibody][1] == counts[ibody][2]) {
n += 3*counts[ibody][0] + 6*counts[ibody][1] - 6; n += 3*counts[ibody][0] + 6*counts[ibody][1] - 6;
@ -1211,7 +1213,7 @@ int FixRigidSmall::dof(int tgroup)
} }
} }
} }
} else if (domain->dimension == 2) { } else if (dimension == 2) {
for (int ibody = 0; ibody < nlocal_body; ibody++) for (int ibody = 0; ibody < nlocal_body; ibody++)
if (counts[ibody][0]+counts[ibody][1] == counts[ibody][2]) if (counts[ibody][0]+counts[ibody][1] == counts[ibody][2])
n += 2*counts[ibody][0] + 3*counts[ibody][1] - 3; n += 2*counts[ibody][0] + 3*counts[ibody][1] - 3;
@ -1297,13 +1299,19 @@ void FixRigidSmall::set_xv()
// x = displacement from center-of-mass, based on body orientation // x = displacement from center-of-mass, based on body orientation
// v = vcm + omega around center-of-mass // v = vcm + omega around center-of-mass
// enforce 2d x and v
MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space,displace[i],x[i]); MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space,displace[i],x[i]);
v[i][0] = b->omega[1]*x[i][2] - b->omega[2]*x[i][1] + b->vcm[0]; v[i][0] = b->omega[1]*x[i][2] - b->omega[2]*x[i][1] + b->vcm[0];
v[i][1] = b->omega[2]*x[i][0] - b->omega[0]*x[i][2] + b->vcm[1]; v[i][1] = b->omega[2]*x[i][0] - b->omega[0]*x[i][2] + b->vcm[1];
v[i][2] = b->omega[0]*x[i][1] - b->omega[1]*x[i][0] + b->vcm[2]; v[i][2] = b->omega[0]*x[i][1] - b->omega[1]*x[i][0] + b->vcm[2];
if (dimension == 2) {
x[i][2] = 0.0;
v[i][2] = 0.0;
}
// add center of mass to displacement // add center of mass to displacement
// map back into periodic box via xbox,ybox,zbox // map back into periodic box via xbox,ybox,zbox
// for triclinic, add in box tilt factors as well // for triclinic, add in box tilt factors as well
@ -1345,6 +1353,7 @@ void FixRigidSmall::set_xv()
} }
// update the position of geometric center // update the position of geometric center
for (int ibody = 0; ibody < nlocal_body + nghost_body; ibody++) { for (int ibody = 0; ibody < nlocal_body + nghost_body; ibody++) {
Body *b = &body[ibody]; Body *b = &body[ibody];
MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space, MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space,
@ -1462,10 +1471,15 @@ void FixRigidSmall::set_v()
v2 = v[i][2]; v2 = v[i][2];
} }
// compute new v
// enforce 2d v
v[i][0] = b->omega[1]*delta[2] - b->omega[2]*delta[1] + b->vcm[0]; v[i][0] = b->omega[1]*delta[2] - b->omega[2]*delta[1] + b->vcm[0];
v[i][1] = b->omega[2]*delta[0] - b->omega[0]*delta[2] + b->vcm[1]; v[i][1] = b->omega[2]*delta[0] - b->omega[0]*delta[2] + b->vcm[1];
v[i][2] = b->omega[0]*delta[1] - b->omega[1]*delta[0] + b->vcm[2]; v[i][2] = b->omega[0]*delta[1] - b->omega[1]*delta[0] + b->vcm[2];
if (dimension == 2) v[i][2] = 0.0;
// virial = unwrapped coords dotted into body constraint force // virial = unwrapped coords dotted into body constraint force
// body constraint force = implied force due to v change minus f external // body constraint force = implied force due to v change minus f external
// assume f does not include forces internal to body // assume f does not include forces internal to body

3
src/RIGID/fix_rigid_small.h
View File

@ -56,7 +56,6 @@ class FixRigidSmall : public Fix {
void pre_neighbor() override; void pre_neighbor() override;
int dof(int) override; int dof(int) override;
void deform(int) override; void deform(int) override;
void enforce2d() override;
void reset_dt() override; void reset_dt() override;
void zero_momentum() override; void zero_momentum() override;
void zero_rotation() override; void zero_rotation() override;
@ -72,6 +71,7 @@ class FixRigidSmall : public Fix {
double dtv, dtf, dtq; double dtv, dtf, dtq;
double *step_respa; double *step_respa;
int triclinic; int triclinic;
int dimension;
char *inpfile; // file to read rigid body attributes from char *inpfile; // file to read rigid body attributes from
int setupflag; // 1 if body properties are setup, else 0 int setupflag; // 1 if body properties are setup, else 0
@ -203,6 +203,7 @@ class FixRigidSmall : public Fix {
void setup_bodies_dynamic(); void setup_bodies_dynamic();
void apply_langevin_thermostat(); void apply_langevin_thermostat();
void compute_forces_and_torques(); void compute_forces_and_torques();
void enforce2d();
void readfile(int, double **, int *); void readfile(int, double **, int *);
void grow_body(); void grow_body();
void reset_atom2body(); void reset_atom2body();

1
src/fix.cpp
View File

@ -72,7 +72,6 @@ Fix::Fix(LAMMPS *lmp, int /*narg*/, char **arg) :
dynamic = 0; dynamic = 0;
dof_flag = 0; dof_flag = 0;
special_alter_flag = 0; special_alter_flag = 0;
enforce2d_flag = 0;
respa_level_support = 0; respa_level_support = 0;
respa_level = -1; respa_level = -1;
maxexchange = 0; maxexchange = 0;

2
src/fix.h
View File

@ -72,7 +72,6 @@ class Fix : protected Pointers {
int dynamic_group_allow; // 1 if can be used with dynamic group, else 0 int dynamic_group_allow; // 1 if can be used with dynamic group, else 0
int dof_flag; // 1 if has dof() method (not min_dof()) int dof_flag; // 1 if has dof() method (not min_dof())
int special_alter_flag; // 1 if has special_alter() meth for spec lists int special_alter_flag; // 1 if has special_alter() meth for spec lists
int enforce2d_flag; // 1 if has enforce2d method
int respa_level_support; // 1 if fix supports fix_modify respa int respa_level_support; // 1 if fix supports fix_modify respa
int respa_level; // which respa level to apply fix (1-Nrespa) int respa_level; // which respa level to apply fix (1-Nrespa)
int maxexchange; // max # of per-atom values for Comm::exchange() int maxexchange; // max # of per-atom values for Comm::exchange()
@ -236,7 +235,6 @@ class Fix : protected Pointers {
virtual void deform(int) {} virtual void deform(int) {}
virtual void reset_target(double) {} virtual void reset_target(double) {}
virtual void reset_dt() {} virtual void reset_dt() {}
virtual void enforce2d() {}
virtual void read_data_header(char *) {} virtual void read_data_header(char *) {}
virtual void read_data_section(char *, int, char *, tagint) {} virtual void read_data_section(char *, int, char *, tagint) {}

36
src/fix_enforce2d.cpp
View File

@ -27,12 +27,9 @@ using namespace FixConst;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixEnforce2D::FixEnforce2D(LAMMPS *lmp, int narg, char **arg) : FixEnforce2D::FixEnforce2D(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), Fix(lmp, narg, arg)
flist(nullptr)
{ {
if (narg != 3) error->all(FLERR,"Illegal fix enforce2d command"); if (narg != 3) error->all(FLERR,"Illegal fix enforce2d command");
nfixlist = 0;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -40,8 +37,6 @@ FixEnforce2D::FixEnforce2D(LAMMPS *lmp, int narg, char **arg) :
FixEnforce2D::~FixEnforce2D() FixEnforce2D::~FixEnforce2D()
{ {
if (copymode) return; if (copymode) return;
delete [] flist;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -61,29 +56,6 @@ void FixEnforce2D::init()
{ {
if (domain->dimension == 3) if (domain->dimension == 3)
error->all(FLERR,"Cannot use fix enforce2d with 3d simulation"); error->all(FLERR,"Cannot use fix enforce2d with 3d simulation");
// list of fixes with enforce2d methods
nfixlist = 0;
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->enforce2d_flag) nfixlist++;
if (nfixlist) {
int myindex = -1;
delete [] flist;
flist = new Fix*[nfixlist];
nfixlist = 0;
for (int i = 0; i < modify->nfix; i++) {
if (modify->fix[i]->enforce2d_flag) {
if (myindex < 0)
flist[nfixlist++] = modify->fix[i];
else
error->all(FLERR,"Fix enforce2d must be defined after fix {}",
modify->fix[i]->style);
}
if (modify->fix[i] == this) myindex = i;
}
}
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -153,12 +125,6 @@ void FixEnforce2D::post_force(int /*vflag*/)
torque[i][1] = 0.0; torque[i][1] = 0.0;
} }
} }
// invoke other fixes that enforce 2d
// fix rigid variants
for (int m = 0; m < nfixlist; m++)
flist[m]->enforce2d();
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

4
src/fix_enforce2d.h
View File

@ -35,10 +35,6 @@ class FixEnforce2D : public Fix {
void post_force(int) override; void post_force(int) override;
void post_force_respa(int, int, int) override; void post_force_respa(int, int, int) override;
void min_post_force(int) override; void min_post_force(int) override;
protected:
int nfixlist;
class Fix **flist;
}; };
} // namespace LAMMPS_NS } // namespace LAMMPS_NS