forked from lijiext/lammps
enable fix rigid commands to add gravity to COM of rigid bodies
This commit is contained in:
Steve Plimpton
parent
b1668f28dd
commit
5859b18d2f
|
|
@ -52,8 +52,8 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
|||
torque(NULL), quat(NULL), imagebody(NULL), fflag(NULL),
|
||||
tflag(NULL), langextra(NULL), sum(NULL), all(NULL),
|
||||
remapflag(NULL), xcmimage(NULL), eflags(NULL), orient(NULL),
|
||||
dorient(NULL), id_dilate(NULL), random(NULL), avec_ellipsoid(NULL),
|
||||
avec_line(NULL), avec_tri(NULL)
|
||||
dorient(NULL), id_dilate(NULL), id_gravity(NULL), random(NULL),
|
||||
avec_ellipsoid(NULL), avec_line(NULL), avec_tri(NULL)
|
||||
{
|
||||
int i,ibody;
|
||||
|
||||
|
|
@ -124,14 +124,17 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
|||
if (custom_flag) {
|
||||
if (narg < 5) error->all(FLERR,"Illegal fix rigid command");
|
||||
|
||||
// determine whether atom-style variable or atom property is used.
|
||||
// determine whether atom-style variable or atom property is used
|
||||
|
||||
if (strstr(arg[4],"i_") == arg[4]) {
|
||||
int is_double=0;
|
||||
int custom_index = atom->find_custom(arg[4]+2,is_double);
|
||||
if (custom_index == -1)
|
||||
error->all(FLERR,"Fix rigid custom requires previously defined property/atom");
|
||||
error->all(FLERR,"Fix rigid custom requires "
|
||||
"previously defined property/atom");
|
||||
else if (is_double)
|
||||
error->all(FLERR,"Fix rigid custom requires integer-valued property/atom");
|
||||
error->all(FLERR,"Fix rigid custom requires "
|
||||
"integer-valued property/atom");
|
||||
int minval = INT_MAX;
|
||||
int *value = atom->ivector[custom_index];
|
||||
for (i = 0; i < nlocal; i++)
|
||||
|
|
@ -163,12 +166,15 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
|||
if (mask[i] & groupbit)
|
||||
molecule[i] = (tagint)((tagint)value[i] - minval + 1);
|
||||
delete[] value;
|
||||
|
||||
} else error->all(FLERR,"Unsupported fix rigid custom property");
|
||||
|
||||
} else {
|
||||
if (atom->molecule_flag == 0)
|
||||
error->all(FLERR,"Fix rigid molecule requires atom attribute molecule");
|
||||
molecule = atom->molecule;
|
||||
}
|
||||
|
||||
iarg = 4 + custom_flag;
|
||||
|
||||
tagint maxmol_tag = -1;
|
||||
|
|
@ -297,6 +303,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
|||
}
|
||||
|
||||
// number of linear rigid bodies is counted later
|
||||
|
||||
nlinear = 0;
|
||||
|
||||
// parse optional args
|
||||
|
|
@ -308,12 +315,14 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
|||
tstat_flag = 0;
|
||||
pstat_flag = 0;
|
||||
allremap = 1;
|
||||
id_dilate = NULL;
|
||||
t_chain = 10;
|
||||
t_iter = 1;
|
||||
t_order = 3;
|
||||
p_chain = 10;
|
||||
|
||||
inpfile = NULL;
|
||||
id_gravity = NULL;
|
||||
id_dilate = NULL;
|
||||
|
||||
pcouple = NONE;
|
||||
pstyle = ANISO;
|
||||
|
|
@ -543,12 +552,20 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
|||
iarg += 2;
|
||||
|
||||
} else if (strcmp(arg[iarg],"reinit") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command");
|
||||
if (strcmp("yes",arg[iarg+1]) == 0) reinitflag = 1;
|
||||
else if (strcmp("no",arg[iarg+1]) == 0) reinitflag = 0;
|
||||
else error->all(FLERR,"Illegal fix rigid command");
|
||||
iarg += 2;
|
||||
|
||||
} else if (strcmp(arg[iarg],"gravity") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command");
|
||||
delete [] id_gravity;
|
||||
int n = strlen(arg[iarg+1]) + 1;
|
||||
id_gravity = new char[n];
|
||||
strcpy(id_gravity,arg[iarg+1]);
|
||||
iarg += 2;
|
||||
|
||||
} else error->all(FLERR,"Illegal fix rigid command");
|
||||
}
|
||||
|
||||
|
|
@ -621,6 +638,9 @@ FixRigid::~FixRigid()
|
|||
|
||||
delete random;
|
||||
delete [] inpfile;
|
||||
delete [] id_dilate;
|
||||
delete [] id_gravity;
|
||||
|
||||
memory->destroy(mol2body);
|
||||
memory->destroy(body2mol);
|
||||
|
||||
|
|
@ -687,7 +707,7 @@ void FixRigid::init()
|
|||
avec_tri = (AtomVecTri *) atom->style_match("tri");
|
||||
|
||||
// warn if more than one rigid fix
|
||||
// if earlyflag, warn if any post-force fixes come after POEMS fix
|
||||
// if earlyflag, warn if any post-force fixes come after a rigid fix
|
||||
|
||||
int count = 0;
|
||||
for (i = 0; i < modify->nfix; i++)
|
||||
|
|
@ -701,7 +721,8 @@ void FixRigid::init()
|
|||
if (rflag && (modify->fmask[i] & POST_FORCE) &&
|
||||
!modify->fix[i]->rigid_flag) {
|
||||
char str[128];
|
||||
snprintf(str,128,"Fix %s alters forces after fix rigid",modify->fix[i]->id);
|
||||
snprintf(str,128,"Fix %s alters forces after fix rigid",
|
||||
modify->fix[i]->id);
|
||||
error->warning(FLERR,str);
|
||||
}
|
||||
}
|
||||
|
|
@ -713,12 +734,24 @@ void FixRigid::init()
|
|||
if (strcmp(modify->fix[i]->style,"npt") == 0) break;
|
||||
if (strcmp(modify->fix[i]->style,"nph") == 0) break;
|
||||
}
|
||||
|
||||
if (i < modify->nfix) {
|
||||
for (int j = i; j < modify->nfix; j++)
|
||||
if (strcmp(modify->fix[j]->style,"rigid") == 0)
|
||||
error->all(FLERR,"Rigid fix must come before NPT/NPH fix");
|
||||
}
|
||||
|
||||
// add gravity forces based on gravity vector from fix
|
||||
|
||||
if (id_gravity) {
|
||||
int ifix = modify->find_fix(id_gravity);
|
||||
if (ifix < 0) error->all(FLERR,"Fix rigid cannot find fix gravity ID");
|
||||
if (strcmp(modify->fix[ifix]->style,"gravity") != 0)
|
||||
error->all(FLERR,"Fix rigid gravity fix is invalid");
|
||||
int tmp;
|
||||
gvec = (double *) modify->fix[ifix]->extract("gvec",tmp);
|
||||
}
|
||||
|
||||
// timestep info
|
||||
|
||||
dtv = update->dt;
|
||||
|
|
@ -1061,8 +1094,17 @@ void FixRigid::compute_forces_and_torques()
|
|||
torque[ibody][1] = all[ibody][4] + langextra[ibody][4];
|
||||
torque[ibody][2] = all[ibody][5] + langextra[ibody][5];
|
||||
}
|
||||
}
|
||||
|
||||
// add gravity force to COM of each body
|
||||
|
||||
if (id_gravity) {
|
||||
for (ibody = 0; ibody < nbody; ibody++) {
|
||||
fcm[ibody][0] += gvec[0]*masstotal[ibody];
|
||||
fcm[ibody][1] += gvec[1]*masstotal[ibody];
|
||||
fcm[ibody][2] += gvec[2]*masstotal[ibody];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
|
|
|
|||
|
|
@ -66,7 +66,8 @@ class FixRigid : public Fix {
|
|||
double *step_respa;
|
||||
int triclinic;
|
||||
|
||||
char *inpfile; // file to read rigid body attributes from
|
||||
char *inpfile; // file to read rigid body attributes from
|
||||
|
||||
int rstyle; // SINGLE,MOLECULE,GROUP
|
||||
int setupflag; // 1 if body properties are setup, else 0
|
||||
int earlyflag; // 1 if forces/torques computed at post_force()
|
||||
|
|
@ -131,6 +132,9 @@ class FixRigid : public Fix {
|
|||
int dilate_group_bit; // mask for dilation group
|
||||
char *id_dilate; // group name to dilate
|
||||
|
||||
char *id_gravity; // ID of fix gravity command to add gravity forces
|
||||
double *gvec; // ptr to gravity vector inside the fix
|
||||
|
||||
class RanMars *random;
|
||||
class AtomVecEllipsoid *avec_ellipsoid;
|
||||
class AtomVecLine *avec_line;
|
||||
|
|
|
|||
|
|
@ -57,7 +57,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
|
|||
xcmimage(NULL), displace(NULL), eflags(NULL), orient(NULL), dorient(NULL),
|
||||
avec_ellipsoid(NULL), avec_line(NULL), avec_tri(NULL), counts(NULL),
|
||||
itensor(NULL), mass_body(NULL), langextra(NULL), random(NULL),
|
||||
id_dilate(NULL), onemols(NULL)
|
||||
id_dilate(NULL), id_gravity(NULL), onemols(NULL)
|
||||
{
|
||||
int i;
|
||||
|
||||
|
|
@ -107,7 +107,8 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
|
|||
bodyID = new tagint[nlocal];
|
||||
customflag = 1;
|
||||
|
||||
// determine whether atom-style variable or atom property is used.
|
||||
// determine whether atom-style variable or atom property is used
|
||||
|
||||
if (strstr(arg[4],"i_") == arg[4]) {
|
||||
int is_double=0;
|
||||
int custom_index = atom->find_custom(arg[4]+2,is_double);
|
||||
|
|
@ -356,6 +357,13 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
|
|||
p_chain = force->numeric(FLERR,arg[iarg+1]);
|
||||
iarg += 2;
|
||||
|
||||
} else if (strcmp(arg[iarg],"gravity") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
|
||||
delete [] id_gravity;
|
||||
int n = strlen(arg[iarg+1]) + 1;
|
||||
id_gravity = new char[n];
|
||||
strcpy(id_gravity,arg[iarg+1]);
|
||||
iarg += 2;
|
||||
|
||||
} else error->all(FLERR,"Illegal fix rigid/small command");
|
||||
}
|
||||
|
|
@ -515,6 +523,8 @@ FixRigidSmall::~FixRigidSmall()
|
|||
|
||||
delete random;
|
||||
delete [] inpfile;
|
||||
delete [] id_dilate;
|
||||
delete [] id_gravity;
|
||||
|
||||
memory->destroy(langextra);
|
||||
memory->destroy(mass_body);
|
||||
|
|
@ -543,6 +553,7 @@ void FixRigidSmall::init()
|
|||
triclinic = domain->triclinic;
|
||||
|
||||
// warn if more than one rigid fix
|
||||
// if earlyflag, warn if any post-force fixes come after a rigid fix
|
||||
|
||||
int count = 0;
|
||||
for (i = 0; i < modify->nfix; i++)
|
||||
|
|
@ -575,6 +586,17 @@ void FixRigidSmall::init()
|
|||
error->all(FLERR,"Rigid fix must come before NPT/NPH fix");
|
||||
}
|
||||
|
||||
// add gravity forces based on gravity vector from fix
|
||||
|
||||
if (id_gravity) {
|
||||
int ifix = modify->find_fix(id_gravity);
|
||||
if (ifix < 0) error->all(FLERR,"Fix rigid/small cannot find fix gravity ID");
|
||||
if (strcmp(modify->fix[ifix]->style,"gravity") != 0)
|
||||
error->all(FLERR,"Fix rigid/small gravity fix is invalid");
|
||||
int tmp;
|
||||
gvec = (double *) modify->fix[ifix]->extract("gvec",tmp);
|
||||
}
|
||||
|
||||
// timestep info
|
||||
|
||||
dtv = update->dt;
|
||||
|
|
@ -954,8 +976,20 @@ void FixRigidSmall::compute_forces_and_torques()
|
|||
tcm[2] += langextra[ibody][5];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// add gravity force to COM of each body
|
||||
|
||||
if (id_gravity) {
|
||||
double mass;
|
||||
for (ibody = 0; ibody < nlocal_body; ibody++) {
|
||||
mass = body[ibody].mass;
|
||||
fcm = body[ibody].fcm;
|
||||
fcm[0] += gvec[0]*mass;
|
||||
fcm[1] += gvec[1]*mass;
|
||||
fcm[2] += gvec[2]*mass;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
|
|
|
|||
|
|
@ -164,6 +164,9 @@ class FixRigidSmall : public Fix {
|
|||
int dilate_group_bit; // mask for dilation group
|
||||
char *id_dilate; // group name to dilate
|
||||
|
||||
char *id_gravity; // ID of fix gravity command to add gravity forces
|
||||
double *gvec; // ptr to gravity vector inside the fix
|
||||
|
||||
double p_current[3],p_target[3];
|
||||
|
||||
// molecules added on-the-fly as rigid bodies
|
||||
|
|
|
|||
|
|
@ -37,7 +37,8 @@ enum{CONSTANT,EQUAL};
|
|||
|
||||
FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg),
|
||||
mstr(NULL), vstr(NULL), pstr(NULL), tstr(NULL), xstr(NULL), ystr(NULL), zstr(NULL)
|
||||
mstr(NULL), vstr(NULL), pstr(NULL), tstr(NULL),
|
||||
xstr(NULL), ystr(NULL), zstr(NULL)
|
||||
{
|
||||
if (narg < 5) error->all(FLERR,"Illegal fix gravity command");
|
||||
|
||||
|
|
@ -61,8 +62,10 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
|
|||
mstyle = CONSTANT;
|
||||
}
|
||||
|
||||
int iarg;
|
||||
|
||||
if (strcmp(arg[4],"chute") == 0) {
|
||||
if (narg != 6) error->all(FLERR,"Illegal fix gravity command");
|
||||
if (narg < 6) error->all(FLERR,"Illegal fix gravity command");
|
||||
style = CHUTE;
|
||||
if (strstr(arg[5],"v_") == arg[5]) {
|
||||
int n = strlen(&arg[5][2]) + 1;
|
||||
|
|
@ -73,9 +76,10 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
|
|||
vert = force->numeric(FLERR,arg[5]);
|
||||
vstyle = CONSTANT;
|
||||
}
|
||||
iarg = 6;
|
||||
|
||||
} else if (strcmp(arg[4],"spherical") == 0) {
|
||||
if (narg != 7) error->all(FLERR,"Illegal fix gravity command");
|
||||
if (narg < 7) error->all(FLERR,"Illegal fix gravity command");
|
||||
style = SPHERICAL;
|
||||
if (strstr(arg[5],"v_") == arg[5]) {
|
||||
int n = strlen(&arg[5][2]) + 1;
|
||||
|
|
@ -95,9 +99,10 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
|
|||
theta = force->numeric(FLERR,arg[6]);
|
||||
tstyle = CONSTANT;
|
||||
}
|
||||
iarg = 7;
|
||||
|
||||
} else if (strcmp(arg[4],"vector") == 0) {
|
||||
if (narg != 8) error->all(FLERR,"Illegal fix gravity command");
|
||||
if (narg < 8) error->all(FLERR,"Illegal fix gravity command");
|
||||
style = VECTOR;
|
||||
if (strstr(arg[5],"v_") == arg[5]) {
|
||||
int n = strlen(&arg[5][2]) + 1;
|
||||
|
|
@ -126,9 +131,23 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
|
|||
zdir = force->numeric(FLERR,arg[7]);
|
||||
zstyle = CONSTANT;
|
||||
}
|
||||
iarg = 8;
|
||||
|
||||
} else error->all(FLERR,"Illegal fix gravity command");
|
||||
|
||||
// optional keywords
|
||||
|
||||
disable = 0;
|
||||
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"disable") == 0) {
|
||||
disable = 1;
|
||||
iarg++;
|
||||
} else error->all(FLERR,"Illegal fix gravity command");
|
||||
}
|
||||
|
||||
// initializations
|
||||
|
||||
degree2rad = MY_PI/180.0;
|
||||
time_origin = update->ntimestep;
|
||||
|
||||
|
|
@ -266,6 +285,12 @@ void FixGravity::post_force(int /*vflag*/)
|
|||
set_acceleration();
|
||||
}
|
||||
|
||||
// just exit if application of force is disabled
|
||||
|
||||
if (disable) return;
|
||||
|
||||
// apply gravity force to each particle
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *rmass = atom->rmass;
|
||||
|
|
@ -338,9 +363,9 @@ void FixGravity::set_acceleration()
|
|||
}
|
||||
}
|
||||
|
||||
xacc = magnitude*xgrav;
|
||||
yacc = magnitude*ygrav;
|
||||
zacc = magnitude*zgrav;
|
||||
gvec[0] = xacc = magnitude*xgrav;
|
||||
gvec[1] = yacc = magnitude*ygrav;
|
||||
gvec[2] = zacc = magnitude*zgrav;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
|
@ -357,3 +382,16 @@ double FixGravity::compute_scalar()
|
|||
}
|
||||
return egrav_all;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
extract current gravity direction vector
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void *FixGravity::extract(const char *name, int &dim)
|
||||
{
|
||||
if (strcmp(name,"gvec") == 0) {
|
||||
dim = 1;
|
||||
return (void *) gvec;
|
||||
}
|
||||
return NULL;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -36,9 +36,10 @@ class FixGravity : public Fix {
|
|||
virtual void post_force(int);
|
||||
virtual void post_force_respa(int, int, int);
|
||||
double compute_scalar();
|
||||
void *extract(const char *, int &);
|
||||
|
||||
protected:
|
||||
int style;
|
||||
int style,disable;
|
||||
double magnitude;
|
||||
double vert,phi,theta;
|
||||
double xdir,ydir,zdir;
|
||||
|
|
@ -46,6 +47,8 @@ class FixGravity : public Fix {
|
|||
double degree2rad;
|
||||
int ilevel_respa;
|
||||
int time_origin;
|
||||
double gvec[3];
|
||||
|
||||
int eflag;
|
||||
double egrav,egrav_all;
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue