jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13010 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2015-01-30 22:22:27 +00:00
parent 8d657c5391
commit cef3a9da47
6 changed files with 252 additions and 210 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

223
src/RIGID/fix_rigid.cpp
View File

@ -75,6 +75,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
extended = orientflag = dorientflag = 0; extended = orientflag = dorientflag = 0;
body = NULL; body = NULL;
xcmimage = NULL;
displace = NULL; displace = NULL;
eflags = NULL; eflags = NULL;
orient = NULL; orient = NULL;
@ -352,7 +353,8 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"temp") == 0) { } else if (strcmp(arg[iarg],"temp") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0 && if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0 &&
strcmp(style,"rigid/nvt/omp") != 0 && strcmp(style,"rigid/npt/omp") != 0) strcmp(style,"rigid/nvt/omp") != 0 &&
strcmp(style,"rigid/npt/omp") != 0)
error->all(FLERR,"Illegal fix rigid command"); error->all(FLERR,"Illegal fix rigid command");
tstat_flag = 1; tstat_flag = 1;
t_start = force->numeric(FLERR,arg[iarg+1]); t_start = force->numeric(FLERR,arg[iarg+1]);
@ -363,8 +365,9 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"iso") == 0) { } else if (strcmp(arg[iarg],"iso") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 && if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
strcmp(style,"rigid/npt/omp") != 0 && strcmp(style,"rigid/nph/omp") != 0) strcmp(style,"rigid/npt/omp") != 0 &&
error->all(FLERR,"Illegal fix rigid command"); strcmp(style,"rigid/nph/omp") != 0)
error->all(FLERR,"Illegal fix rigid command");
pcouple = XYZ; pcouple = XYZ;
p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]); p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
@ -380,8 +383,9 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"aniso") == 0) { } else if (strcmp(arg[iarg],"aniso") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 && if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
strcmp(style,"rigid/npt/omp") != 0 && strcmp(style,"rigid/nph/omp") != 0) strcmp(style,"rigid/npt/omp") != 0 &&
error->all(FLERR,"Illegal fix rigid command"); strcmp(style,"rigid/nph/omp") != 0)
error->all(FLERR,"Illegal fix rigid command");
p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]); p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
p_period[0] = p_period[1] = p_period[2] = p_period[0] = p_period[1] = p_period[2] =
@ -396,8 +400,9 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"x") == 0) { } else if (strcmp(arg[iarg],"x") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 && if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
strcmp(style,"rigid/npt/omp") != 0 && strcmp(style,"rigid/nph/omp") != 0) strcmp(style,"rigid/npt/omp") != 0 &&
error->all(FLERR,"Illegal fix rigid command"); strcmp(style,"rigid/nph/omp") != 0)
error->all(FLERR,"Illegal fix rigid command");
p_start[0] = force->numeric(FLERR,arg[iarg+1]); p_start[0] = force->numeric(FLERR,arg[iarg+1]);
p_stop[0] = force->numeric(FLERR,arg[iarg+2]); p_stop[0] = force->numeric(FLERR,arg[iarg+2]);
p_period[0] = force->numeric(FLERR,arg[iarg+3]); p_period[0] = force->numeric(FLERR,arg[iarg+3]);
@ -407,8 +412,9 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"y") == 0) { } else if (strcmp(arg[iarg],"y") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 && if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
strcmp(style,"rigid/npt/omp") != 0 && strcmp(style,"rigid/nph/omp") != 0) strcmp(style,"rigid/npt/omp") != 0 &&
error->all(FLERR,"Illegal fix rigid command"); strcmp(style,"rigid/nph/omp") != 0)
error->all(FLERR,"Illegal fix rigid command");
p_start[1] = force->numeric(FLERR,arg[iarg+1]); p_start[1] = force->numeric(FLERR,arg[iarg+1]);
p_stop[1] = force->numeric(FLERR,arg[iarg+2]); p_stop[1] = force->numeric(FLERR,arg[iarg+2]);
p_period[1] = force->numeric(FLERR,arg[iarg+3]); p_period[1] = force->numeric(FLERR,arg[iarg+3]);
@ -418,8 +424,9 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"z") == 0) { } else if (strcmp(arg[iarg],"z") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 && if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
strcmp(style,"rigid/npt/omp") != 0 && strcmp(style,"rigid/nph/omp") != 0) strcmp(style,"rigid/npt/omp") != 0 &&
error->all(FLERR,"Illegal fix rigid command"); strcmp(style,"rigid/nph/omp") != 0)
error->all(FLERR,"Illegal fix rigid command");
p_start[2] = force->numeric(FLERR,arg[iarg+1]); p_start[2] = force->numeric(FLERR,arg[iarg+1]);
p_stop[2] = force->numeric(FLERR,arg[iarg+2]); p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
p_period[2] = force->numeric(FLERR,arg[iarg+3]); p_period[2] = force->numeric(FLERR,arg[iarg+3]);
@ -456,7 +463,8 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"tparam") == 0) { } else if (strcmp(arg[iarg],"tparam") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0 && if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0 &&
strcmp(style,"rigid/nvt/omp") != 0 && strcmp(style,"rigid/npt/omp") != 0) strcmp(style,"rigid/nvt/omp") != 0 &&
strcmp(style,"rigid/npt/omp") != 0)
error->all(FLERR,"Illegal fix rigid command"); error->all(FLERR,"Illegal fix rigid command");
t_chain = force->inumeric(FLERR,arg[iarg+1]); t_chain = force->inumeric(FLERR,arg[iarg+1]);
t_iter = force->inumeric(FLERR,arg[iarg+2]); t_iter = force->inumeric(FLERR,arg[iarg+2]);
@ -466,7 +474,8 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"pchain") == 0) { } else if (strcmp(arg[iarg],"pchain") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command"); if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command");
if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 && if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
strcmp(style,"rigid/npt/omp") != 0 && strcmp(style,"rigid/nph/omp") != 0) strcmp(style,"rigid/npt/omp") != 0 &&
strcmp(style,"rigid/nph/omp") != 0)
error->all(FLERR,"Illegal fix rigid command"); error->all(FLERR,"Illegal fix rigid command");
p_chain = force->inumeric(FLERR,arg[iarg+1]); p_chain = force->inumeric(FLERR,arg[iarg+1]);
iarg += 2; iarg += 2;
@ -569,9 +578,10 @@ FixRigid::~FixRigid()
memory->destroy(mol2body); memory->destroy(mol2body);
memory->destroy(body2mol); memory->destroy(body2mol);
// delete locally stored arrays // delete locally stored per-atom arrays
memory->destroy(body); memory->destroy(body);
memory->destroy(xcmimage);
memory->destroy(displace); memory->destroy(displace);
memory->destroy(eflags); memory->destroy(eflags);
memory->destroy(orient); memory->destroy(orient);
@ -675,6 +685,16 @@ void FixRigid::init()
else tfactor = 0.0; else tfactor = 0.0;
} }
/* ----------------------------------------------------------------------
invoke pre_neighbor() to insure body xcmimage flags are reset
needed if Verlet::setup::pbc() has remapped/migrated atoms for 2nd run
------------------------------------------------------------------------- */
void FixRigid::setup_pre_neighbor()
{
pre_neighbor();
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
compute initial fcm and torque on bodies, also initial virial compute initial fcm and torque on bodies, also initial virial
reset all particle velocities to be consistent with vcm and omega reset all particle velocities to be consistent with vcm and omega
@ -716,7 +736,6 @@ void FixRigid::setup(int vflag)
// torque = torque on each rigid body // torque = torque on each rigid body
imageint *image = atom->image;
double **x = atom->x; double **x = atom->x;
double dx,dy,dz; double dx,dy,dz;
@ -729,7 +748,7 @@ void FixRigid::setup(int vflag)
if (body[i] < 0) continue; if (body[i] < 0) continue;
ibody = body[i]; ibody = body[i];
domain->unmap(x[i],image[i],unwrap); domain->unmap(x[i],xcmimage[i],unwrap);
dx = unwrap[0] - xcm[ibody][0]; dx = unwrap[0] - xcm[ibody][0];
dy = unwrap[1] - xcm[ibody][1]; dy = unwrap[1] - xcm[ibody][1];
dz = unwrap[2] - xcm[ibody][2]; dz = unwrap[2] - xcm[ibody][2];
@ -780,6 +799,7 @@ void FixRigid::setup(int vflag)
for (ibody = 0; ibody < nbody; ibody++) for (ibody = 0; ibody < nbody; ibody++)
MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody], MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
ez_space[ibody],inertia[ibody],omega[ibody]); ez_space[ibody],inertia[ibody],omega[ibody]);
set_v(); set_v();
// guesstimate virial as 2x the set_v contribution // guesstimate virial as 2x the set_v contribution
@ -897,7 +917,6 @@ void FixRigid::final_integrate()
// sum over atoms to get force and torque on rigid body // sum over atoms to get force and torque on rigid body
imageint *image = atom->image;
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
@ -916,7 +935,7 @@ void FixRigid::final_integrate()
sum[ibody][1] += f[i][1]; sum[ibody][1] += f[i][1];
sum[ibody][2] += f[i][2]; sum[ibody][2] += f[i][2];
domain->unmap(x[i],image[i],unwrap); domain->unmap(x[i],xcmimage[i],unwrap);
dx = unwrap[0] - xcm[ibody][0]; dx = unwrap[0] - xcm[ibody][0];
dy = unwrap[1] - xcm[ibody][1]; dy = unwrap[1] - xcm[ibody][1];
dz = unwrap[2] - xcm[ibody][2]; dz = unwrap[2] - xcm[ibody][2];
@ -1005,72 +1024,54 @@ void FixRigid::final_integrate_respa(int ilevel, int iloop)
done during pre_neighbor so will be after call to pbc() done during pre_neighbor so will be after call to pbc()
and after fix_deform::pre_exchange() may have flipped box and after fix_deform::pre_exchange() may have flipped box
use domain->remap() in case xcm is far away from box use domain->remap() in case xcm is far away from box
due to 1st definition of rigid body or due to box flip due to first-time definition of rigid body in setup_bodies_static()
if don't do this, then atoms of a body which drifts far away or due to box flip
from a triclinic box will be remapped back into box also adjust imagebody = rigid body image flags, due to xcm remap
with huge displacements when the box tilt changes via set_x() also reset body xcmimage flags of all atoms in bodies
adjust image flag of body and image flags of all atoms in body 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::pre_neighbor()
{ {
imageint original,oldimage,newimage; for (int ibody = 0; ibody < nbody; ibody++)
for (int ibody = 0; ibody < nbody; ibody++) {
original = imagebody[ibody];
domain->remap(xcm[ibody],imagebody[ibody]); domain->remap(xcm[ibody],imagebody[ibody]);
image_shift();
}
if (original == imagebody[ibody]) remapflag[ibody][3] = 0; /* ----------------------------------------------------------------------
else { reset body xcmimage flags of atoms in bodies
oldimage = original & IMGMASK; xcmimage flags are relative to xcm so that body can be unwrapped
newimage = imagebody[ibody] & IMGMASK; xcmimage = true image flag - imagebody flag
remapflag[ibody][0] = newimage - oldimage; ------------------------------------------------------------------------- */
oldimage = (original >> IMGBITS) & IMGMASK;
newimage = (imagebody[ibody] >> IMGBITS) & IMGMASK;
remapflag[ibody][1] = newimage - oldimage;
oldimage = original >> IMG2BITS;
newimage = imagebody[ibody] >> IMG2BITS;
remapflag[ibody][2] = newimage - oldimage;
remapflag[ibody][3] = 1;
}
}
// adjust image flags of any atom in a rigid body whose xcm was remapped void FixRigid::image_shift()
// subtracting remapflag = new-old keeps ix,iy,iz near 0 {
// so body is always in central simulation box int ibody;
imageint tdim,bdim,xdim[3];
imageint *image = atom->image; imageint *image = atom->image;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int ibody;
imageint idim,otherdims;
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (body[i] == -1) continue; if (body[i] < 0) continue;
if (remapflag[body[i]][3] == 0) continue;
ibody = body[i]; ibody = body[i];
if (remapflag[ibody][0]) { tdim = image[i] & IMGMASK;
idim = image[i] & IMGMASK; bdim = imagebody[ibody] & IMGMASK;
otherdims = image[i] ^ idim; xdim[0] = IMGMAX + tdim - bdim;
idim -= remapflag[ibody][0]; tdim = (image[i] >> IMGBITS) & IMGMASK;
idim &= IMGMASK; bdim = (imagebody[ibody] >> IMGBITS) & IMGMASK;
image[i] = otherdims | idim; xdim[1] = IMGMAX + tdim - bdim;
} tdim = image[i] >> IMG2BITS;
if (remapflag[ibody][1]) { bdim = imagebody[ibody] >> IMG2BITS;
idim = (image[i] >> IMGBITS) & IMGMASK; xdim[2] = IMGMAX + tdim - bdim;
otherdims = image[i] ^ (idim << IMGBITS);
idim -= remapflag[ibody][1]; xcmimage[i] = (xdim[2] << IMG2BITS) | (xdim[1] << IMGBITS) | xdim[0];
idim &= IMGMASK;
image[i] = otherdims | (idim << IMGBITS);
}
if (remapflag[ibody][2]) {
idim = image[i] >> IMG2BITS;
otherdims = image[i] ^ (idim << IMG2BITS);
idim -= remapflag[ibody][2];
idim &= IMGMASK;
image[i] = otherdims | (idim << IMG2BITS);
}
} }
} }
@ -1191,7 +1192,6 @@ void FixRigid::set_xv()
double xy,xz,yz; double xy,xz,yz;
double ione[3],exone[3],eyone[3],ezone[3],vr[6],p[3][3]; double ione[3],exone[3],eyone[3],ezone[3],vr[6],p[3][3];
imageint *image = atom->image;
double **x = atom->x; double **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **f = atom->f; double **f = atom->f;
@ -1216,9 +1216,9 @@ void FixRigid::set_xv()
if (body[i] < 0) continue; if (body[i] < 0) continue;
ibody = body[i]; ibody = body[i];
xbox = (image[i] & IMGMASK) - IMGMAX; xbox = (xcmimage[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX; zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX;
// save old positions and velocities for virial // save old positions and velocities for virial
@ -1373,7 +1373,6 @@ void FixRigid::set_v()
double *rmass = atom->rmass; double *rmass = atom->rmass;
double *mass = atom->mass; double *mass = atom->mass;
int *type = atom->type; int *type = atom->type;
imageint *image = atom->image;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double xprd = domain->xprd; double xprd = domain->xprd;
@ -1422,9 +1421,9 @@ void FixRigid::set_v()
fc1 = massone*(v[i][1] - v1)/dtf - f[i][1]; fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; fc2 = massone*(v[i][2] - v2)/dtf - f[i][2];
xbox = (image[i] & IMGMASK) - IMGMAX; xbox = (xcmimage[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX; zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX;
if (triclinic == 0) { if (triclinic == 0) {
x0 = x[i][0] + xbox*xprd; x0 = x[i][0] + xbox*xprd;
@ -1585,11 +1584,19 @@ void FixRigid::setup_bodies_static()
} }
} }
// first-time setting of body xcmimage flags = true image flags
if (!staticflag) {
imageint *image = atom->image;
for (i = 0; i < nlocal; i++)
if (body[i] >= 0) xcmimage[i] = image[i];
else xcmimage[i] = 0;
}
// compute masstotal & center-of-mass of each rigid body // compute masstotal & center-of-mass of each rigid body
// error if image flag is not 0 in a non-periodic dim // error if image flag is not 0 in a non-periodic dim
double **x = atom->x; double **x = atom->x;
imageint *image = atom->image;
int *periodicity = domain->periodicity; int *periodicity = domain->periodicity;
double xprd = domain->xprd; double xprd = domain->xprd;
@ -1608,9 +1615,9 @@ void FixRigid::setup_bodies_static()
if (body[i] < 0) continue; if (body[i] < 0) continue;
ibody = body[i]; ibody = body[i];
xbox = (image[i] & IMGMASK) - IMGMAX; xbox = (xcmimage[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX; zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX;
if (rmass) massone = rmass[i]; if (rmass) massone = rmass[i];
else massone = mass[type[i]]; else massone = mass[type[i]];
@ -1654,9 +1661,10 @@ void FixRigid::setup_bodies_static()
readfile(0,masstotal,xcm,inbody); readfile(0,masstotal,xcm,inbody);
} }
// set image flags for each rigid body to default values // one-time set of rigid body image flags to default values
// then remap the xcm of each body back into simulation box if needed // staticflag insures this is only done once, not on successive runs
// staticflag check insures this in only done once, not on successive runs // then remap the xcm of each body back into simulation box
// and reset body xcmimage flags via pre_neighbor()
if (!staticflag) { if (!staticflag) {
for (ibody = 0; ibody < nbody; ibody++) for (ibody = 0; ibody < nbody; ibody++)
@ -1679,9 +1687,9 @@ void FixRigid::setup_bodies_static()
if (body[i] < 0) continue; if (body[i] < 0) continue;
ibody = body[i]; ibody = body[i];
xbox = (image[i] & IMGMASK) - IMGMAX; xbox = (xcmimage[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX; zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX;
if (triclinic == 0) { if (triclinic == 0) {
xunwrap = x[i][0] + xbox*xprd; xunwrap = x[i][0] + xbox*xprd;
@ -1834,9 +1842,9 @@ void FixRigid::setup_bodies_static()
ibody = body[i]; ibody = body[i];
xbox = (image[i] & IMGMASK) - IMGMAX; xbox = (xcmimage[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX; zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX;
if (triclinic == 0) { if (triclinic == 0) {
xunwrap = x[i][0] + xbox*xprd; xunwrap = x[i][0] + xbox*xprd;
@ -2022,7 +2030,6 @@ void FixRigid::setup_bodies_dynamic()
double **v = atom->v; double **v = atom->v;
double *rmass = atom->rmass; double *rmass = atom->rmass;
double *mass = atom->mass; double *mass = atom->mass;
imageint *image = atom->image;
int *type = atom->type; int *type = atom->type;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
@ -2043,7 +2050,7 @@ void FixRigid::setup_bodies_dynamic()
sum[ibody][1] += v[i][1] * massone; sum[ibody][1] += v[i][1] * massone;
sum[ibody][2] += v[i][2] * massone; sum[ibody][2] += v[i][2] * massone;
domain->unmap(x[i],image[i],unwrap); domain->unmap(x[i],xcmimage[i],unwrap);
dx = unwrap[0] - xcm[ibody][0]; dx = unwrap[0] - xcm[ibody][0];
dy = unwrap[1] - xcm[ibody][1]; dy = unwrap[1] - xcm[ibody][1];
dz = unwrap[2] - xcm[ibody][2]; dz = unwrap[2] - xcm[ibody][2];
@ -2264,6 +2271,7 @@ double FixRigid::memory_usage()
{ {
int nmax = atom->nmax; int nmax = atom->nmax;
double bytes = nmax * sizeof(int); double bytes = nmax * sizeof(int);
bytes += nmax * sizeof(imageint);
bytes += nmax*3 * sizeof(double); bytes += nmax*3 * sizeof(double);
bytes += maxvatom*6 * sizeof(double); // vatom bytes += maxvatom*6 * sizeof(double); // vatom
if (extended) { if (extended) {
@ -2281,6 +2289,7 @@ double FixRigid::memory_usage()
void FixRigid::grow_arrays(int nmax) void FixRigid::grow_arrays(int nmax)
{ {
memory->grow(body,nmax,"rigid:body"); memory->grow(body,nmax,"rigid:body");
memory->grow(xcmimage,nmax,"rigid:xcmimage");
memory->grow(displace,nmax,3,"rigid:displace"); memory->grow(displace,nmax,3,"rigid:displace");
if (extended) { if (extended) {
memory->grow(eflags,nmax,"rigid:eflags"); memory->grow(eflags,nmax,"rigid:eflags");
@ -2296,6 +2305,7 @@ void FixRigid::grow_arrays(int nmax)
void FixRigid::copy_arrays(int i, int j, int delflag) void FixRigid::copy_arrays(int i, int j, int delflag)
{ {
body[j] = body[i]; body[j] = body[i];
xcmimage[j] = xcmimage[i];
displace[j][0] = displace[i][0]; displace[j][0] = displace[i][0];
displace[j][1] = displace[i][1]; displace[j][1] = displace[i][1];
displace[j][2] = displace[i][2]; displace[j][2] = displace[i][2];
@ -2318,6 +2328,7 @@ void FixRigid::copy_arrays(int i, int j, int delflag)
void FixRigid::set_arrays(int i) void FixRigid::set_arrays(int i)
{ {
body[i] = -1; body[i] = -1;
xcmimage[i] = 0;
displace[i][0] = 0.0; displace[i][0] = 0.0;
displace[i][1] = 0.0; displace[i][1] = 0.0;
displace[i][2] = 0.0; displace[i][2] = 0.0;
@ -2329,13 +2340,14 @@ void FixRigid::set_arrays(int i)
int FixRigid::pack_exchange(int i, double *buf) int FixRigid::pack_exchange(int i, double *buf)
{ {
buf[0] = body[i]; buf[0] = ubuf(body[i]).d;
buf[1] = displace[i][0]; buf[1] = ubuf(xcmimage[i]).d;
buf[2] = displace[i][1]; buf[2] = displace[i][0];
buf[3] = displace[i][2]; buf[3] = displace[i][1];
if (!extended) return 4; buf[4] = displace[i][2];
if (!extended) return 5;
int m = 4; int m = 5;
buf[m++] = eflags[i]; buf[m++] = eflags[i];
for (int j = 0; j < orientflag; j++) for (int j = 0; j < orientflag; j++)
buf[m++] = orient[i][j]; buf[m++] = orient[i][j];
@ -2353,13 +2365,14 @@ int FixRigid::pack_exchange(int i, double *buf)
int FixRigid::unpack_exchange(int nlocal, double *buf) int FixRigid::unpack_exchange(int nlocal, double *buf)
{ {
body[nlocal] = static_cast<int> (buf[0]); body[nlocal] = (int) ubuf(buf[0]).i;
displace[nlocal][0] = buf[1]; xcmimage[nlocal] = (imageint) ubuf(buf[1]).i;
displace[nlocal][1] = buf[2]; displace[nlocal][0] = buf[2];
displace[nlocal][2] = buf[3]; displace[nlocal][1] = buf[3];
if (!extended) return 4; displace[nlocal][2] = buf[4];
if (!extended) return 5;
int m = 4; int m = 5;
eflags[nlocal] = static_cast<int> (buf[m++]); eflags[nlocal] = static_cast<int> (buf[m++]);
for (int j = 0; j < orientflag; j++) for (int j = 0; j < orientflag; j++)
orient[nlocal][j] = buf[m++]; orient[nlocal][j] = buf[m++];

4
src/RIGID/fix_rigid.h
View File

@ -47,6 +47,7 @@ class FixRigid : public Fix {
int pack_exchange(int, double *); int pack_exchange(int, double *);
int unpack_exchange(int, double *); int unpack_exchange(int, double *);
void setup_pre_neighbor();
void pre_neighbor(); void pre_neighbor();
int dof(int); int dof(int);
void deform(int); void deform(int);
@ -102,6 +103,8 @@ class FixRigid : public Fix {
int orientflag; // 1 if particles store spatial orientation int orientflag; // 1 if particles store spatial orientation
int dorientflag; // 1 if particles store dipole orientation int dorientflag; // 1 if particles store dipole orientation
imageint *xcmimage; // internal image flags for atoms in rigid bodies
// set relative to in-box xcm of each body
int *eflags; // flags for extended particles int *eflags; // flags for extended particles
double **orient; // orientation vector of particle wrt rigid body double **orient; // orientation vector of particle wrt rigid body
double **dorient; // orientation of dipole mu wrt rigid body double **dorient; // orientation of dipole mu wrt rigid body
@ -133,6 +136,7 @@ class FixRigid : public Fix {
int POINT,SPHERE,ELLIPSOID,LINE,TRIANGLE,DIPOLE; // bitmasks for eflags int POINT,SPHERE,ELLIPSOID,LINE,TRIANGLE,DIPOLE; // bitmasks for eflags
int OMEGA,ANGMOM,TORQUE; int OMEGA,ANGMOM,TORQUE;
void image_shift();
void set_xv(); void set_xv();
void set_v(); void set_v();
void setup_bodies_static(); void setup_bodies_static();

7
src/RIGID/fix_rigid_nh.cpp
View File

@ -615,7 +615,6 @@ void FixRigidNH::final_integrate()
// sum over atoms to get force and torque on rigid body // sum over atoms to get force and torque on rigid body
imageint *image = atom->image;
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
@ -642,9 +641,9 @@ void FixRigidNH::final_integrate()
sum[ibody][1] += f[i][1]; sum[ibody][1] += f[i][1];
sum[ibody][2] += f[i][2]; sum[ibody][2] += f[i][2];
xbox = (image[i] & IMGMASK) - IMGMAX; xbox = (xcmimage[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX; zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX;
if (triclinic == 0) { if (triclinic == 0) {
xunwrap = x[i][0] + xbox*xprd; xunwrap = x[i][0] + xbox*xprd;

36
src/RIGID/fix_rigid_nh_small.cpp
View File

@ -60,32 +60,42 @@ FixRigidNHSmall::FixRigidNHSmall(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Fix rigid/small npt/nph period must be > 0.0"); error->all(FLERR,"Fix rigid/small npt/nph period must be > 0.0");
if (domain->dimension == 2 && p_flag[2]) if (domain->dimension == 2 && p_flag[2])
error->all(FLERR,"Invalid fix rigid/small npt/nph command for a 2d simulation"); error->all(FLERR,"Invalid fix rigid/small npt/nph command "
"for a 2d simulation");
if (domain->dimension == 2 && (pcouple == YZ || pcouple == XZ)) if (domain->dimension == 2 && (pcouple == YZ || pcouple == XZ))
error->all(FLERR,"Invalid fix rigid/small npt/nph command for a 2d simulation"); error->all(FLERR,"Invalid fix rigid/small npt/nph command "
"for a 2d simulation");
if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0)) if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0))
error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings"); error->all(FLERR,"Invalid fix rigid/small npt/nph command "
"pressure settings");
if (pcouple == XYZ && domain->dimension == 3 && p_flag[2] == 0) if (pcouple == XYZ && domain->dimension == 3 && p_flag[2] == 0)
error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings"); error->all(FLERR,"Invalid fix rigid/small npt/nph command "
"pressure settings");
if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0)) if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0))
error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings"); error->all(FLERR,"Invalid fix rigid/small npt/nph command "
"pressure settings");
if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0)) if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0))
error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings"); error->all(FLERR,"Invalid fix rigid/small npt/nph command "
"pressure settings");
if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0)) if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0))
error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings"); error->all(FLERR,"Invalid fix rigid/small npt/nph command "
"pressure settings");
// require periodicity in tensile dimension // require periodicity in tensile dimension
if (p_flag[0] && domain->xperiodic == 0) if (p_flag[0] && domain->xperiodic == 0)
error->all(FLERR, error->all(FLERR,
"Cannot use fix rigid/small npt/nph on a non-periodic dimension"); "Cannot use fix rigid/small npt/nph on a "
"non-periodic dimension");
if (p_flag[1] && domain->yperiodic == 0) if (p_flag[1] && domain->yperiodic == 0)
error->all(FLERR, error->all(FLERR,
"Cannot use fix rigid/small npt/nph on a non-periodic dimension"); "Cannot use fix rigid/small npt/nph on a "
"non-periodic dimension");
if (p_flag[2] && domain->zperiodic == 0) if (p_flag[2] && domain->zperiodic == 0)
error->all(FLERR, error->all(FLERR,
"Cannot use fix rigid/small npt/nph on a non-periodic dimension"); "Cannot use fix rigid/small npt/nph on a "
"non-periodic dimension");
if (pcouple == XYZ && domain->dimension == 3 && if (pcouple == XYZ && domain->dimension == 3 &&
(p_start[0] != p_start[1] || p_start[0] != p_start[2] || (p_start[0] != p_start[1] || p_start[0] != p_start[2] ||
@ -113,7 +123,8 @@ FixRigidNHSmall::FixRigidNHSmall(LAMMPS *lmp, int narg, char **arg) :
(p_flag[0] && p_period[0] <= 0.0) || (p_flag[0] && p_period[0] <= 0.0) ||
(p_flag[1] && p_period[1] <= 0.0) || (p_flag[1] && p_period[1] <= 0.0) ||
(p_flag[2] && p_period[2] <= 0.0)) (p_flag[2] && p_period[2] <= 0.0))
error->all(FLERR,"Fix rigid/small nvt/npt/nph damping parameters must be > 0.0"); error->all(FLERR,"Fix rigid/small nvt/npt/nph damping parameters "
"must be > 0.0");
// memory allocation and initialization // memory allocation and initialization
@ -653,7 +664,6 @@ void FixRigidNHSmall::final_integrate()
// sum over atoms to get force and torque on rigid body // sum over atoms to get force and torque on rigid body
imageint *image = atom->image;
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
@ -678,7 +688,7 @@ void FixRigidNHSmall::final_integrate()
fcm[1] += f[i][1]; fcm[1] += f[i][1];
fcm[2] += f[i][2]; fcm[2] += f[i][2];
domain->unmap(x[i],image[i],unwrap); domain->unmap(x[i],xcmimage[i],unwrap);
xcm = b->xcm; xcm = b->xcm;
dx = unwrap[0] - xcm[0]; dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1]; dy = unwrap[1] - xcm[1];

176
src/RIGID/fix_rigid_small.cpp
View File

@ -89,6 +89,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
bodyown = NULL; bodyown = NULL;
bodytag = NULL; bodytag = NULL;
atom2body = NULL; atom2body = NULL;
xcmimage = NULL;
displace = NULL; displace = NULL;
eflags = NULL; eflags = NULL;
orient = NULL; orient = NULL;
@ -439,6 +440,7 @@ FixRigidSmall::~FixRigidSmall()
memory->destroy(bodyown); memory->destroy(bodyown);
memory->destroy(bodytag); memory->destroy(bodytag);
memory->destroy(atom2body); memory->destroy(atom2body);
memory->destroy(xcmimage);
memory->destroy(displace); memory->destroy(displace);
memory->destroy(eflags); memory->destroy(eflags);
memory->destroy(orient); memory->destroy(orient);
@ -507,11 +509,15 @@ void FixRigidSmall::init()
allows resetting of atom properties like mass between runs allows resetting of atom properties like mass between runs
only do static initialization once if using an infile only do static initialization once if using an infile
cannot do this until now, b/c requires comm->setup() to have setup stencil cannot do this until now, b/c requires comm->setup() to have setup stencil
invoke pre_neighbor() to insure body xcmimage flags are reset
needed if Verlet::setup::pbc() has remapped/migrated atoms for 2nd run
setup_bodies_static() invokes pre_neighbor itself
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void FixRigidSmall::setup_pre_neighbor() void FixRigidSmall::setup_pre_neighbor()
{ {
if (!staticflag || !infile) setup_bodies_static(); if (!staticflag || !infile) setup_bodies_static();
else pre_neighbor();
setup_bodies_dynamic(); setup_bodies_dynamic();
} }
@ -532,7 +538,6 @@ void FixRigidSmall::setup(int vflag)
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
imageint *image = atom->image;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double *xcm,*fcm,*tcm; double *xcm,*fcm,*tcm;
@ -555,7 +560,7 @@ void FixRigidSmall::setup(int vflag)
fcm[1] += f[i][1]; fcm[1] += f[i][1];
fcm[2] += f[i][2]; fcm[2] += f[i][2];
domain->unmap(x[i],image[i],unwrap); domain->unmap(x[i],xcmimage[i],unwrap);
xcm = b->xcm; xcm = b->xcm;
dx = unwrap[0] - xcm[0]; dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1]; dy = unwrap[1] - xcm[1];
@ -740,7 +745,6 @@ void FixRigidSmall::final_integrate()
// sum over atoms to get force and torque on rigid body // sum over atoms to get force and torque on rigid body
imageint *image = atom->image;
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
@ -765,7 +769,7 @@ void FixRigidSmall::final_integrate()
fcm[1] += f[i][1]; fcm[1] += f[i][1];
fcm[2] += f[i][2]; fcm[2] += f[i][2];
domain->unmap(x[i],image[i],unwrap); domain->unmap(x[i],xcmimage[i],unwrap);
xcm = b->xcm; xcm = b->xcm;
dx = unwrap[0] - xcm[0]; dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1]; dy = unwrap[1] - xcm[1];
@ -867,6 +871,26 @@ void FixRigidSmall::final_integrate_respa(int ilevel, int iloop)
final_integrate(); final_integrate();
} }
/* ----------------------------------------------------------------------
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 adjust rgdimage flags of all atoms in bodies for two effects
(1) change in true image flags due to pbc() call during exchange
(2) change in imagebody due to xcm remap
rgdimage flags are always -1,0,-1 so that body can be unwrapped
around in-box xcm and stay close to simulation box
if don't do this, then a body could end up very far away
when unwrapped by true image flags, and 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 probem?
------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
atom to processor assignments have changed, atom to processor assignments have changed,
so acquire ghost bodies and setup atom2body so acquire ghost bodies and setup atom2body
@ -883,31 +907,6 @@ void FixRigidSmall::final_integrate_respa(int ilevel, int iloop)
void FixRigidSmall::pre_neighbor() void FixRigidSmall::pre_neighbor()
{ {
// remap xcm and image flags of each body as needed
imageint original,oldimage,newimage;
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
original = b->image;
domain->remap(b->xcm,b->image);
if (original == b->image) b->remapflag[3] = 0;
else {
oldimage = original & IMGMASK;
newimage = b->image & IMGMASK;
b->remapflag[0] = newimage - oldimage;
oldimage = (original >> IMGBITS) & IMGMASK;
newimage = (b->image >> IMGBITS) & IMGMASK;
b->remapflag[1] = newimage - oldimage;
oldimage = original >> IMG2BITS;
newimage = b->image >> IMG2BITS;
b->remapflag[2] = newimage - oldimage;
b->remapflag[3] = 1;
}
}
// acquire ghost bodies via forward comm // acquire ghost bodies via forward comm
// also gets new remapflags needed for adjusting atom image flags // also gets new remapflags needed for adjusting atom image flags
// reset atom2body for owned atoms // reset atom2body for owned atoms
@ -919,40 +918,42 @@ void FixRigidSmall::pre_neighbor()
reset_atom2body(); reset_atom2body();
//check(4); //check(4);
// adjust image flags of any atom in a rigid body whose xcm was remapped for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
domain->remap(b->xcm,b->image);
}
image_shift();
}
/* ----------------------------------------------------------------------
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()
{
int ibody;
imageint tdim,bdim,xdim[3];
imageint *image = atom->image; imageint *image = atom->image;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
imageint idim,otherdims;
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue; if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]]; Body *b = &body[atom2body[i]];
if (b->remapflag[3] == 0) continue; 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;
if (b->remapflag[0]) { xcmimage[i] = (xdim[2] << IMG2BITS) | (xdim[1] << IMGBITS) | xdim[0];
idim = image[i] & IMGMASK;
otherdims = image[i] ^ idim;
idim -= b->remapflag[0];
idim &= IMGMASK;
image[i] = otherdims | idim;
}
if (b->remapflag[1]) {
idim = (image[i] >> IMGBITS) & IMGMASK;
otherdims = image[i] ^ (idim << IMGBITS);
idim -= b->remapflag[1];
idim &= IMGMASK;
image[i] = otherdims | (idim << IMGBITS);
}
if (b->remapflag[2]) {
idim = image[i] >> IMG2BITS;
otherdims = image[i] ^ (idim << IMG2BITS);
idim -= b->remapflag[2];
idim &= IMGMASK;
image[i] = otherdims | (idim << IMG2BITS);
}
} }
} }
@ -1091,7 +1092,6 @@ void FixRigidSmall::set_xv()
double xz = domain->xz; double xz = domain->xz;
double yz = domain->yz; double yz = domain->yz;
imageint *image = atom->image;
double **x = atom->x; double **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **f = atom->f; double **f = atom->f;
@ -1106,9 +1106,9 @@ void FixRigidSmall::set_xv()
if (atom2body[i] < 0) continue; if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]]; Body *b = &body[atom2body[i]];
xbox = (image[i] & IMGMASK) - IMGMAX; xbox = (xcmimage[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX; zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX;
// save old positions and velocities for virial // save old positions and velocities for virial
@ -1264,7 +1264,6 @@ void FixRigidSmall::set_v()
double *rmass = atom->rmass; double *rmass = atom->rmass;
double *mass = atom->mass; double *mass = atom->mass;
int *type = atom->type; int *type = atom->type;
imageint *image = atom->image;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
// set v of each atom // set v of each atom
@ -1300,9 +1299,9 @@ void FixRigidSmall::set_v()
fc1 = massone*(v[i][1] - v1)/dtf - f[i][1]; fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; fc2 = massone*(v[i][2] - v2)/dtf - f[i][2];
xbox = (image[i] & IMGMASK) - IMGMAX; xbox = (xcmimage[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX; zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX;
if (triclinic == 0) { if (triclinic == 0) {
x0 = x[i][0] + xbox*xprd; x0 = x[i][0] + xbox*xprd;
@ -1771,6 +1770,15 @@ void FixRigidSmall::setup_bodies_static()
} }
} }
// first-time setting of body xcmimage flags = true image flags
if (!staticflag) {
imageint *image = atom->image;
for (i = 0; i < nlocal; i++)
if (bodytag[i] >= 0) xcmimage[i] = image[i];
else xcmimage[i] = 0;
}
// acquire ghost bodies via forward comm // acquire ghost bodies via forward comm
// set atom2body for ghost atoms via forward comm // set atom2body for ghost atoms via forward comm
// set atom2body for other owned atoms via reset_atom2body() // set atom2body for other owned atoms via reset_atom2body()
@ -1783,7 +1791,6 @@ void FixRigidSmall::setup_bodies_static()
// compute mass & center-of-mass of each rigid body // compute mass & center-of-mass of each rigid body
double **x = atom->x; double **x = atom->x;
imageint *image = atom->image;
double *xcm; double *xcm;
@ -1803,7 +1810,7 @@ void FixRigidSmall::setup_bodies_static()
if (rmass) massone = rmass[i]; if (rmass) massone = rmass[i];
else massone = mass[type[i]]; else massone = mass[type[i]];
domain->unmap(x[i],image[i],unwrap); domain->unmap(x[i],xcmimage[i],unwrap);
xcm = b->xcm; xcm = b->xcm;
xcm[0] += unwrap[0] * massone; xcm[0] += unwrap[0] * massone;
xcm[1] += unwrap[1] * massone; xcm[1] += unwrap[1] * massone;
@ -1833,14 +1840,16 @@ void FixRigidSmall::setup_bodies_static()
readfile(0,NULL,inbody); readfile(0,NULL,inbody);
} }
// set image flags for each rigid body to default values // one-time set of rigid body image flags to default values
// then remap the xcm of each body back into simulation box if needed // staticflag insures this is only done once, not on successive runs
// staticflag check insures this in only done once, not on successive runs // then remap the xcm of each body back into simulation box
// and reset body xcmimage flags via pre_neighbor()
if (!staticflag) if (!staticflag) {
for (ibody = 0; ibody < nlocal_body; ibody++) for (ibody = 0; ibody < nlocal_body; ibody++)
body[ibody].image = ((imageint) IMGMAX << IMG2BITS) | body[ibody].image = ((imageint) IMGMAX << IMG2BITS) |
((imageint) IMGMAX << IMGBITS) | IMGMAX; ((imageint) IMGMAX << IMGBITS) | IMGMAX;
}
pre_neighbor(); pre_neighbor();
@ -1859,7 +1868,7 @@ void FixRigidSmall::setup_bodies_static()
if (atom2body[i] < 0) continue; if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]]; Body *b = &body[atom2body[i]];
domain->unmap(x[i],image[i],unwrap); domain->unmap(x[i],xcmimage[i],unwrap);
xcm = b->xcm; xcm = b->xcm;
dx = unwrap[0] - xcm[0]; dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1]; dy = unwrap[1] - xcm[1];
@ -2015,7 +2024,7 @@ void FixRigidSmall::setup_bodies_static()
Body *b = &body[atom2body[i]]; Body *b = &body[atom2body[i]];
domain->unmap(x[i],image[i],unwrap); domain->unmap(x[i],xcmimage[i],unwrap);
xcm = b->xcm; xcm = b->xcm;
delta[0] = unwrap[0] - xcm[0]; delta[0] = unwrap[0] - xcm[0];
delta[1] = unwrap[1] - xcm[1]; delta[1] = unwrap[1] - xcm[1];
@ -2203,7 +2212,6 @@ void FixRigidSmall::setup_bodies_dynamic()
double *rmass = atom->rmass; double *rmass = atom->rmass;
double *mass = atom->mass; double *mass = atom->mass;
int *type = atom->type; int *type = atom->type;
imageint *image = atom->image;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double *xcm,*vcm,*acm; double *xcm,*vcm,*acm;
@ -2229,7 +2237,7 @@ void FixRigidSmall::setup_bodies_dynamic()
vcm[1] += v[i][1] * massone; vcm[1] += v[i][1] * massone;
vcm[2] += v[i][2] * massone; vcm[2] += v[i][2] * massone;
domain->unmap(x[i],image[i],unwrap); domain->unmap(x[i],xcmimage[i],unwrap);
xcm = b->xcm; xcm = b->xcm;
dx = unwrap[0] - xcm[0]; dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1]; dy = unwrap[1] - xcm[1];
@ -2525,6 +2533,7 @@ void FixRigidSmall::grow_arrays(int nmax)
memory->grow(bodyown,nmax,"rigid/small:bodyown"); memory->grow(bodyown,nmax,"rigid/small:bodyown");
memory->grow(bodytag,nmax,"rigid/small:bodytag"); memory->grow(bodytag,nmax,"rigid/small:bodytag");
memory->grow(atom2body,nmax,"rigid/small:atom2body"); memory->grow(atom2body,nmax,"rigid/small:atom2body");
memory->grow(xcmimage,nmax,"rigid/small:xcmimage");
memory->grow(displace,nmax,3,"rigid/small:displace"); memory->grow(displace,nmax,3,"rigid/small:displace");
if (extended) { if (extended) {
memory->grow(eflags,nmax,"rigid/small:eflags"); memory->grow(eflags,nmax,"rigid/small:eflags");
@ -2540,6 +2549,7 @@ void FixRigidSmall::grow_arrays(int nmax)
void FixRigidSmall::copy_arrays(int i, int j, int delflag) void FixRigidSmall::copy_arrays(int i, int j, int delflag)
{ {
bodytag[j] = bodytag[i]; bodytag[j] = bodytag[i];
xcmimage[j] = xcmimage[i];
displace[j][0] = displace[i][0]; displace[j][0] = displace[i][0];
displace[j][1] = displace[i][1]; displace[j][1] = displace[i][1];
displace[j][2] = displace[i][2]; displace[j][2] = displace[i][2];
@ -2579,6 +2589,7 @@ void FixRigidSmall::set_arrays(int i)
bodyown[i] = -1; bodyown[i] = -1;
bodytag[i] = 0; bodytag[i] = 0;
atom2body[i] = -1; atom2body[i] = -1;
xcmimage[i] = 0;
displace[i][0] = 0.0; displace[i][0] = 0.0;
displace[i][1] = 0.0; displace[i][1] = 0.0;
displace[i][2] = 0.0; displace[i][2] = 0.0;
@ -2670,13 +2681,14 @@ void FixRigidSmall::set_molecule(int nlocalprev, tagint tagprev, int imol,
int FixRigidSmall::pack_exchange(int i, double *buf) int FixRigidSmall::pack_exchange(int i, double *buf)
{ {
buf[0] = ubuf(bodytag[i]).d; buf[0] = ubuf(bodytag[i]).d;
buf[1] = displace[i][0]; buf[1] = ubuf(xcmimage[i]).d;
buf[2] = displace[i][1]; buf[2] = displace[i][0];
buf[3] = displace[i][2]; buf[3] = displace[i][1];
buf[4] = displace[i][2];
// extended attribute info // extended attribute info
int m = 4; int m = 5;
if (extended) { if (extended) {
buf[m++] = eflags[i]; buf[m++] = eflags[i];
for (int j = 0; j < orientflag; j++) for (int j = 0; j < orientflag; j++)
@ -2714,13 +2726,14 @@ int FixRigidSmall::pack_exchange(int i, double *buf)
int FixRigidSmall::unpack_exchange(int nlocal, double *buf) int FixRigidSmall::unpack_exchange(int nlocal, double *buf)
{ {
bodytag[nlocal] = (tagint) ubuf(buf[0]).i; bodytag[nlocal] = (tagint) ubuf(buf[0]).i;
displace[nlocal][0] = buf[1]; xcmimage[nlocal] = (imageint) ubuf(buf[1]).i;
displace[nlocal][1] = buf[2]; displace[nlocal][0] = buf[2];
displace[nlocal][2] = buf[3]; displace[nlocal][1] = buf[3];
displace[nlocal][2] = buf[4];
// extended attribute info // extended attribute info
int m = 4; int m = 5;
if (extended) { if (extended) {
eflags[nlocal] = static_cast<int> (buf[m++]); eflags[nlocal] = static_cast<int> (buf[m++]);
for (int j = 0; j < orientflag; j++) for (int j = 0; j < orientflag; j++)
@ -3327,7 +3340,8 @@ double FixRigidSmall::compute_scalar()
double FixRigidSmall::memory_usage() double FixRigidSmall::memory_usage()
{ {
int nmax = atom->nmax; int nmax = atom->nmax;
double bytes = 2 * nmax * sizeof(int); double bytes = nmax*2 * sizeof(int);
bytes += nmax * sizeof(imageint);
bytes += nmax*3 * sizeof(double); bytes += nmax*3 * sizeof(double);
bytes += maxvatom*6 * sizeof(double); // vatom bytes += maxvatom*6 * sizeof(double); // vatom
if (extended) { if (extended) {

16
src/RIGID/fix_rigid_small.h
View File

@ -117,15 +117,16 @@ class FixRigidSmall : public Fix {
// ID = tag of atom that owns body // ID = tag of atom that owns body
int *atom2body; // index of owned/ghost body this atom is in, -1 if not int *atom2body; // index of owned/ghost body this atom is in, -1 if not
// can point to original or any image of the body // can point to original or any image of the body
imageint *xcmimage; // internal image flags for atoms in rigid bodies
// set relative to in-box xcm of each body
double **displace; // displacement of each atom in body coords double **displace; // displacement of each atom in body coords
int *eflags; // flags for extended particles
double **orient; // orientation vector of particle wrt rigid body
double **dorient; // orientation of dipole mu wrt rigid body
int *eflags; // flags for extended particles int extended; // 1 if any particles have extended attributes
double **orient; // orientation vector of particle wrt rigid body int orientflag; // 1 if particles store spatial orientation
double **dorient; // orientation of dipole mu wrt rigid body int dorientflag; // 1 if particles store dipole orientation
int extended; // 1 if any particles have extended attributes
int orientflag; // 1 if particles store spatial orientation
int dorientflag; // 1 if particles store dipole orientation
int POINT,SPHERE,ELLIPSOID,LINE,TRIANGLE,DIPOLE; // bitmasks for eflags int POINT,SPHERE,ELLIPSOID,LINE,TRIANGLE,DIPOLE; // bitmasks for eflags
int OMEGA,ANGMOM,TORQUE; int OMEGA,ANGMOM,TORQUE;
@ -182,6 +183,7 @@ class FixRigidSmall : public Fix {
double *rsqclose; double *rsqclose;
double rsqfar; double rsqfar;
void image_shift();
void set_xv(); void set_xv();
void set_v(); void set_v();
void create_bodies(); void create_bodies();