bug fix for fix GCMC w/ fix shake, enhance of fix wall/gran/region with restarting

This commit is contained in:
Steve Plimpton 2016-10-13 16:55:53 -06:00
parent 63e71cd45b
commit 2c8a7a318a
18 changed files with 254 additions and 145 deletions

View File

@ -91,6 +91,3 @@ This will install virtualenv from the Python Package Index.
----------------
Installing prerequisites for PDF build

View File

@ -72,9 +72,10 @@ compute peratom all pe/atom
compute pe all reduce sum c_peratom
thermo_style custom step temp etotal press pe c_pe :pre
NOTE: The per-atom energy does not any Lennard-Jones tail corrections
invoked by the "pair_modify tail yes"_pair_modify.html command, since
those are global contributions to the system energy.
NOTE: The per-atom energy does not include any Lennard-Jones tail
corrections to the energy added by the "pair_modify tail
yes"_pair_modify.html command, since those are contributions to the
global system energy.
[Output info:]

View File

@ -133,6 +133,11 @@ compute p all reduce sum c_peratom\[1\] c_peratom\[2\] c_peratom\[3\]
variable press equal -(c_p\[1\]+c_p\[2\]+c_p\[3\])/(3*vol)
thermo_style custom step temp etotal press v_press :pre
NOTE: The per-atom stress does not include any Lennard-Jones tail
corrections to the pressure added by the "pair_modify tail
yes"_pair_modify.html command, since those are contributions to the
global system pressure.
[Output info:]
This compute calculates a per-atom array with 6 columns, which can be

View File

@ -163,16 +163,20 @@ info on how to re-specify a fix in an input script that reads a
restart file, so that the operation of the fix continues in an
uninterrupted fashion.
Note that info about region definitions is NOT included in restart
files. So you must re-define your region and if it is a moving
region, define its motion attributes in a way that is consistent with
the simulation that wrote the restart file. In particular, if you
want to change its motion attributes (e.g. its velocity), then you
should insure the postition/orientation of the region at the initial
restart timestep is the same as it was on the timestep the restart
file was written. If this is not possible, then you may need to
ignore info in the restart file by defining a new fix wall/gran/region
command in your restart script (e.g. with a different fix ID).
NOTE: Information about region definitions is NOT included in restart
files, as discussed on the "read_restart"_read_restart.html doc page.
So you must re-define your region and if it is a moving region, define
its motion attributes in a way that is consistent with the simulation
that wrote the restart file. In particular, if you want to change the
region motion attributes (e.g. its velocity), then you should ensure
the postition/orientation of the region at the initial restart
timestep is the same as it was on the timestep the restart file was
written. If this is not possible, you may need to ignore info in the
restart file by defining a new fix wall/gran/region command in your
restart script, e.g. with a different fix ID. Or if you want to keep
the shear history info but discard the region motion information, you
can use the same fix ID for fix wall/gran/region, but assign it a
region with a different region ID.
None of the "fix_modify"_fix_modify.html options are relevant to this
fix. No global or per-atom quantities are stored by this fix for

View File

@ -197,12 +197,13 @@ was used in the input script that wrote the restart file.
If a match is found, LAMMPS prints a message indicating that the fix
is being re-enabled. If no match is found before the first run or
minimization is performed by the new script, the "state" information
for the saved fix is discarded. LAMMPS will also print a list of
fixes for which the information is being discarded. See the doc pages
for individual fixes for info on which ones can be restarted in this
manner. Note that fixes which are created internally by other LAMMPS
commands (computes, fixes, etc) will have style names which are
all-capitalized, and IDs which are generated internally.
for the saved fix is discarded. At the time the discard occurs,
LAMMPS will also print a list of fixes for which the information is
being discarded. See the doc pages for individual fixes for info on
which ones can be restarted in this manner. Note that fixes which are
created internally by other LAMMPS commands (computes, fixes, etc)
will have style names which are all-capitalized, and IDs which are
generated internally.
Likewise, the "computes"_fix.html used for a simulation are not stored
in the restart file. This means the new input script should specify
@ -219,12 +220,16 @@ described in the previous paragraph, so that the compute can continue
its calculations in a consistent manner.
NOTE: There are a handful of commands which can be used before or
between runs which require a system initialization. Examples include
the "balance", "displace_atoms", and "delete_atoms" commands. This is
because they may migrate atoms to new processors. Thus they will also
discard unused "state" information from fixes. This means that, if
desired, you must re-specify the relevant fixes and computes (which
create fixes) before those commands are used.
between runs which may require a system initialization. Examples
include the "balance", "displace_atoms", "delete_atoms", "set" (some
options), and "velocity" (some options) commands. This is because
they can migrate atoms to new processors. Thus they will also discard
unused "state" information from fixes. You will know the discard has
occurred because a list of discarded fixes will be printed to the
screen and log file, as explained above. This means that if you wish
to retain that info in a restarted run, you must re-specify the
relevant fixes and computes (which create fixes) before those commands
are used.
Some pair styles, like the "granular pair styles"_pair_gran.html, also
use a fix to store "state" information that persists from timestep to

View File

@ -53,7 +53,6 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
if (!atom->sphere_flag)
error->all(FLERR,"Fix wall/gran requires atom style sphere");
restart_peratom = 1;
create_attribute = 1;
// set interaction style
@ -65,8 +64,8 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
//else if (strcmp(arg[3],"bonded/history") == 0) pairstyle = BONDED_HISTORY;
else error->all(FLERR,"Invalid fix wall/gran interaction style");
history = 1;
if (pairstyle == HOOKE) history = 0;
history = restart_peratom = 1;
if (pairstyle == HOOKE) history = restart_peratom = 0;
// wall/particle coefficients

View File

@ -101,20 +101,27 @@ void FixWallGranRegion::init()
error->all(FLERR,"Region ID for fix wall/gran/region does not exist");
region = domain->regions[iregion];
// region displacement and orientation theta at previous step
// check if region properties changed between runs
// reset if restart info was inconsistent
if (motion_resetflag) {
if (comm->me == 0) {
char str[128];
sprintf(str,"Properties for region %s do not match restart file, "
"resetting its motion",idregion);
error->warning(FLERR,str);
}
if (strcmp(idregion,region->id) != 0 ||
strcmp(region_style,region->style) != 0 ||
nregion != region->nregion) {
char str[256];
sprintf(str,"Region properties for region %s changed between runs, "
"resetting its motion",idregion);
error->warning(FLERR,str);
region->reset_vel();
}
}
}
if (motion_resetflag){
char str[256];
sprintf(str,"Region properties for region %s are inconsistent "
"with restart file, resetting its motion",idregion);
error->warning(FLERR,str);
region->reset_vel();
}
}
/* ---------------------------------------------------------------------- */
@ -483,19 +490,9 @@ int FixWallGranRegion::size_restart(int nlocal)
void FixWallGranRegion::write_restart(FILE *fp)
{
if (comm->me) return;
int size_id_str = (strlen(region->id) + 1) * sizeof(char);
int size_style_str = (strlen(region->style) + 1) * sizeof(char);
int size_tot = sizeof(int) + size_id_str +
sizeof(int) + size_style_str + sizeof(int) +
region->size_restart*sizeof(double);
fwrite(&size_tot,sizeof(int),1,fp);
fwrite(&size_id_str,sizeof(int),1,fp);
fwrite(region->id,sizeof(char),size_id_str,fp);
fwrite(&size_style_str,sizeof(int),1,fp);
fwrite(region->style,sizeof(char),size_style_str,fp);
fwrite(&region->nregion,sizeof(int),1,fp);
int len = 0;
region->length_restart_string(len);
fwrite(&len, sizeof(int),1,fp);
region->write_restart(fp);
}
@ -505,44 +502,6 @@ void FixWallGranRegion::write_restart(FILE *fp)
void FixWallGranRegion::restart(char *buf)
{
int n = 0;
int size_id_str = buf[n];
n += sizeof(int);
char *region_id_restart = new char[size_id_str];
for (int i = 0; i < size_id_str; i++){
region_id_restart[i] = buf[n++];
}
int size_style_str = buf[n];
n += sizeof(int);
char *region_style_restart = new char[size_style_str];
for (int i = 0; i < size_style_str; i++)
region_style_restart[i] = buf[n++];
int nregion_restart = buf[n];
n += sizeof(int);
if (check_consistent_region(region,region_id_restart,
region_style_restart,nregion_restart))
region->restart(buf,n);
else motion_resetflag = 1;
delete [] region_id_restart;
delete [] region_style_restart;
}
/* ----------------------------------------------------------------------
check that region id/style/number of sub-regions are consistent
------------------------------------------------------------------------- */
int FixWallGranRegion::check_consistent_region(Region *region,
char* region_id,
char* region_style, int nregion)
{
if (strcmp(region_id, region->id) != 0 ||
strcmp(region_style, region->style) != 0 ||
nregion != region->nregion)
return 0;
return 1;
int n = 0;
if (!region->restart(buf,n)) motion_resetflag = 1;
}

View File

@ -62,7 +62,6 @@ class FixWallGranRegion : public FixWallGran {
// vel info is to be reset
void update_contacts(int, int);
int check_consistent_region(Region *, char*, char*, int);
};
}

View File

@ -40,7 +40,8 @@ using namespace FixConst;
FixBondCreate::FixBondCreate(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
bondcount(NULL), partner(NULL), finalpartner(NULL), distsq(NULL), created(NULL), copy(NULL), random(NULL)
bondcount(NULL), partner(NULL), finalpartner(NULL), distsq(NULL),
created(NULL), copy(NULL), random(NULL)
{
if (narg < 8) error->all(FLERR,"Illegal fix bond/create command");

View File

@ -2108,6 +2108,12 @@ double FixGCMC::energy_full()
int eflag = 1;
int vflag = 0;
// clear forces so they don't accumulate over multiple
// calls within fix gcmc timestep, e.g. for fix shake
size_t nbytes = sizeof(double) * (atom->nlocal + atom->nghost);
if (nbytes) memset(&atom->f[0][0],0,3*nbytes);
if (modify->n_pre_force) modify->pre_force(vflag);
if (force->pair) force->pair->compute(eflag,vflag);
@ -2121,9 +2127,18 @@ double FixGCMC::energy_full()
if (force->kspace) force->kspace->compute(eflag,vflag);
// unlike Verlet, not performing a reverse_comm() or forces here
// b/c GCMC does not care about forces
// don't think it will mess up energy due to any post_force() fixes
if (modify->n_post_force) modify->post_force(vflag);
if (modify->n_end_of_step) modify->end_of_step();
// NOTE: all fixes with THERMO_ENERGY mask set and which
// operate at pre_force() or post_force() or end_of_step()
// and which user has enable via fix_modify thermo yes,
// will contribute to total MC energy via pe->compute_scalar()
update->eflag_global = update->ntimestep;
double total_energy = c_pe->compute_scalar();

View File

@ -54,8 +54,6 @@ FixRecenter::FixRecenter(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
dynamic_group_allow = 1;
/* ---------------------------------------------------------------------- */
if (strcmp(arg[3],"NULL") == 0) xflag = 0;
else if (strcmp(arg[3],"INIT") == 0) xinitflag = 1;
else xcom = force->numeric(FLERR,arg[3]);

View File

@ -707,7 +707,7 @@ void LAMMPS::init()
atom->init(); // atom must come after force and domain
// atom deletes extra array
// used by fix shear_history::unpack_restart()
// when force->pair->gran_history creates fix ??
// when force->pair->gran_history creates fix
// atom_vec init uses deform_vremap
modify->init(); // modify must come after update, force, atom, domain
neighbor->init(); // neighbor must come after force, modify

View File

@ -42,10 +42,6 @@ Region::Region(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
xstr = ystr = zstr = tstr = NULL;
dx = dy = dz = 0.0;
// used by set_velocity() even if no rotation specified
point[0] = point[1] = point[2] = 0.0;
size_restart = 5;
reset_vel();
copymode = 0;
@ -398,7 +394,6 @@ void Region::options(int narg, char **arg)
else error->all(FLERR,"Illegal region command");
}
// error check
if ((moveflag || rotateflag) &&
@ -524,42 +519,95 @@ void Region::velocity_contact(double *vwall, double *x, int ic)
{
Contact c = contact[ic];
double xc[3];
xc[0] = x[0] - contact[ic].delx;
xc[1] = x[1] - contact[ic].dely;
xc[2] = x[2] - contact[ic].delz;
vwall[0] = v[0] + omega[1]*(xc[2] - rpoint[2]) - omega[2]*(xc[1] - rpoint[1]);
vwall[1] = v[1] + omega[2]*(xc[0] - rpoint[0]) - omega[0]*(xc[2] - rpoint[2]);
vwall[2] = v[2] + omega[0]*(xc[1] - rpoint[1]) - omega[1]*(xc[0] - rpoint[0]);
if (varshape && contact[ic].varflag) velocity_contact_shape(vwall,xc);
vwall[0] = vwall[1] = vwall[2] = 0.0;
if (moveflag){
vwall[0] = v[0];
vwall[1] = v[1];
vwall[2] = v[2];
}
if (rotateflag){
xc[0] = x[0] - contact[ic].delx;
xc[1] = x[1] - contact[ic].dely;
xc[2] = x[2] - contact[ic].delz;
vwall[0] += omega[1]*(xc[2] - rpoint[2]) - omega[2]*(xc[1] - rpoint[1]);
vwall[1] += omega[2]*(xc[0] - rpoint[0]) - omega[0]*(xc[2] - rpoint[2]);
vwall[2] += omega[0]*(xc[1] - rpoint[1]) - omega[1]*(xc[0] - rpoint[0]);
}
if (varshape && contact[ic].varflag) velocity_contact_shape(vwall, xc);
}
/* ----------------------------------------------------------------------
increment length of restart buffer based on region info
used by restart of fix/wall/gran/region
------------------------------------------------------------------------- */
void Region::length_restart_string(int &n)
{
n += sizeof(int) + strlen(id)+1 +
sizeof(int) + strlen(style)+1 + sizeof(int) +
size_restart*sizeof(double);
}
/* ----------------------------------------------------------------------
region writes its current position/angle
region writes its current style, id, number of sub-regions
and position/angle
needed by fix/wall/gran/region to compute velocity by differencing scheme
------------------------------------------------------------------------- */
void Region::write_restart(FILE *fp)
{
int sizeid = (strlen(id)+1);
int sizestyle = (strlen(style)+1);
fwrite(&sizeid, sizeof(int), 1, fp);
fwrite(id, 1, sizeid, fp);
fwrite(&sizestyle, sizeof(int), 1, fp);
fwrite(style, 1, sizestyle, fp);
fwrite(&nregion,sizeof(int),1,fp);
fwrite(prev, sizeof(double), size_restart, fp);
}
/* ----------------------------------------------------------------------
region reads its previous position/angle
region reads style, id, number of sub-regions from restart file
if they match current region, also read previous position/angle
needed by fix/wall/gran/region to compute velocity by differencing scheme
------------------------------------------------------------------------- */
int Region::restart(char *buf, int n)
int Region::restart(char *buf, int &n)
{
int sizeid = buf[n];
n += sizeof(int);
char *restart_id = new char[sizeid];
for (int i = 0; i < sizeid; i++)
restart_id[i] = buf[n++];
if (strcmp(restart_id,id) != 0) return 0;
int sizestyle = buf[n];
n += sizeof(int);
char *restart_style = new char[sizestyle];
for (int i = 0; i < sizestyle; i++)
restart_style[i] = buf[n++];
if (strcmp(restart_style,style) != 0) return 0;
int restart_nreg = buf[n];
n += sizeof(int);
if (restart_nreg != nregion) return 0;
char *rlist = new char[size_restart*sizeof(double)];
for (int i = 0; i < size_restart*sizeof(double); i++)
rlist[i] = buf[n++];
for (int i = 0; i < size_restart; i++) {
for (int i = 0; i < size_restart; i++){
prev[i] = ((double *)rlist)[i];
}
delete [] rlist;
return n;
delete [] restart_id;
delete [] restart_style;
return 1;
}
/* ----------------------------------------------------------------------
@ -568,6 +616,5 @@ int Region::restart(char *buf, int n)
void Region::reset_vel()
{
for (int i = 0; i < size_restart; i++)
prev[i] = 0;
for (int i = 0; i < size_restart; i++) prev[i] = 0;
}

View File

@ -59,8 +59,7 @@ class Region : protected Pointers {
double rpoint[3]; // current origin of rotation axis
double omega[3]; // angular velocity
double rprev; // speed of time-dependent radius, if applicable
double xcenter[3]; // translated/rotated center of cylinder/sphere
// only used with varshape
double xcenter[3]; // translated/rotated center of cylinder/sphere (only used if varshape)
double prev[5]; // stores displacement (X3), angle and if
// necessary, region variable size (e.g. radius)
// at previous time step
@ -84,7 +83,8 @@ class Region : protected Pointers {
virtual void set_velocity();
void velocity_contact(double *, double *, int);
virtual void write_restart(FILE *);
virtual int restart(char *, int);
virtual int restart(char *, int&);
virtual void length_restart_string(int&);
virtual void reset_vel();
// implemented by each region, not called by other classes
@ -106,12 +106,11 @@ class Region : protected Pointers {
void options(int, char **);
void point_on_line_segment(double *, double *, double *, double *);
void forward_transform(double &, double &, double &);
double point[3],runit[3];
private:
char *xstr,*ystr,*zstr,*tstr;
int xvar,yvar,zvar,tvar;
double axis[3];
double axis[3],point[3],runit[3];
void inverse_transform(double &, double &, double &);
void rotate(double &, double &, double &, double);

View File

@ -91,11 +91,8 @@ RegIntersect::RegIntersect(LAMMPS *lmp, int narg, char **arg) :
// for near contacts and touching contacts
cmax = 0;
size_restart = 0;
for (int ilist = 0; ilist < nregion; ilist++){
for (int ilist = 0; ilist < nregion; ilist++)
cmax += regions[list[ilist]]->cmax;
size_restart += regions[list[ilist]]->size_restart;
}
contact = new Contact[cmax];
tmax = 0;
@ -274,6 +271,7 @@ void RegIntersect::pretransform()
regions[list[ilist]]->pretransform();
}
/* ----------------------------------------------------------------------
get translational/angular velocities of all subregions
------------------------------------------------------------------------- */
@ -285,6 +283,19 @@ void RegIntersect::set_velocity()
regions[list[ilist]]->set_velocity();
}
/* ----------------------------------------------------------------------
increment length of restart buffer based on region info
used by restart of fix/wall/gran/region
------------------------------------------------------------------------- */
void RegIntersect::length_restart_string(int& n)
{
n += sizeof(int) + strlen(id)+1 +
sizeof(int) + strlen(style)+1 + sizeof(int);
for (int ilist = 0; ilist < nregion; ilist++)
domain->regions[list[ilist]]->length_restart_string(n);
}
/* ----------------------------------------------------------------------
region writes its current position/angle
needed by fix/wall/gran/region to compute velocity by differencing scheme
@ -292,8 +303,17 @@ void RegIntersect::set_velocity()
void RegIntersect::write_restart(FILE *fp)
{
for (int ilist = 0; ilist < nregion; ilist++)
int sizeid = (strlen(id)+1);
int sizestyle = (strlen(style)+1);
fwrite(&sizeid, sizeof(int), 1, fp);
fwrite(id, 1, sizeid, fp);
fwrite(&sizestyle, sizeof(int), 1, fp);
fwrite(style, 1, sizestyle, fp);
fwrite(&nregion,sizeof(int),1,fp);
for (int ilist = 0; ilist < nregion; ilist++){
domain->regions[list[ilist]]->write_restart(fp);
}
}
/* ----------------------------------------------------------------------
@ -301,11 +321,33 @@ void RegIntersect::write_restart(FILE *fp)
needed by fix/wall/gran/region to compute velocity by differencing scheme
------------------------------------------------------------------------- */
int RegIntersect::restart(char *buf, int n)
int RegIntersect::restart(char *buf, int& n)
{
for (int ilist = 0; ilist < nregion; ilist++)
n = domain->regions[list[ilist]]->restart(buf, n);
return n;
int sizeid = buf[n];
n += sizeof(int);
char *restart_id = new char[sizeid];
for (int i = 0; i < sizeid; i++)
restart_id[i] = buf[n++];
if (strcmp(restart_id,id) != 0) return 0;
int sizestyle = buf[n];
n += sizeof(int);
char *restart_style = new char[sizestyle];
for (int i = 0; i < sizestyle; i++)
restart_style[i] = buf[n++];
if (strcmp(restart_style,style) != 0) return 0;
int restart_nreg = buf[n];
n += sizeof(int);
if (restart_nreg != nregion) return 0;
for (int ilist = 0; ilist < nregion; ilist++){
if (!domain->regions[list[ilist]]->restart(buf, n)){
return 0;
}
}
return 1;
}
/* ----------------------------------------------------------------------

View File

@ -35,8 +35,9 @@ class RegIntersect : public Region {
void shape_update();
void pretransform();
void set_velocity();
void length_restart_string(int&);
void write_restart(FILE *);
int restart(char *, int);
int restart(char *, int&);
void reset_vel();
private:

View File

@ -85,11 +85,8 @@ RegUnion::RegUnion(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
// for near contacts and touching contacts
cmax = 0;
size_restart = 0;
for (int ilist = 0; ilist < nregion; ilist++){
for (int ilist = 0; ilist < nregion; ilist++)
cmax += regions[list[ilist]]->cmax;
size_restart += regions[list[ilist]]->size_restart;
}
contact = new Contact[cmax];
tmax = 0;
@ -279,6 +276,20 @@ void RegUnion::set_velocity()
for (int ilist = 0; ilist < nregion; ilist++)
regions[list[ilist]]->set_velocity();
}
/* ----------------------------------------------------------------------
increment length of restart buffer based on region info
used by restart of fix/wall/gran/region
------------------------------------------------------------------------- */
void RegUnion::length_restart_string(int& n)
{
n += sizeof(int) + strlen(id)+1 +
sizeof(int) + strlen(style)+1 + sizeof(int);
for (int ilist = 0; ilist < nregion; ilist++)
domain->regions[list[ilist]]->length_restart_string(n);
}
/* ----------------------------------------------------------------------
region writes its current position/angle
needed by fix/wall/gran/region to compute velocity by differencing scheme
@ -286,8 +297,15 @@ void RegUnion::set_velocity()
void RegUnion::write_restart(FILE *fp)
{
int sizeid = (strlen(id)+1);
int sizestyle = (strlen(style)+1);
fwrite(&sizeid, sizeof(int), 1, fp);
fwrite(id, 1, sizeid, fp);
fwrite(&sizestyle, sizeof(int), 1, fp);
fwrite(style, 1, sizestyle, fp);
fwrite(&nregion,sizeof(int),1,fp);
for (int ilist = 0; ilist < nregion; ilist++)
domain->regions[list[ilist]]->write_restart(fp);
domain->regions[list[ilist]]->write_restart(fp);
}
/* ----------------------------------------------------------------------
@ -295,11 +313,30 @@ void RegUnion::write_restart(FILE *fp)
needed by fix/wall/gran/region to compute velocity by differencing scheme
------------------------------------------------------------------------- */
int RegUnion::restart(char *buf, int n)
int RegUnion::restart(char *buf, int &n)
{
for (int ilist = 0; ilist < nregion; ilist++)
n = domain->regions[list[ilist]]->restart(buf, n);
return n;
int sizeid = buf[n];
n += sizeof(int);
char *restart_id = new char[sizeid];
for (int i = 0; i < sizeid; i++)
restart_id[i] = buf[n++];
if (strcmp(restart_id,id) != 0) return 0;
int sizestyle = buf[n];
n += sizeof(int);
char *restart_style = new char[sizestyle];
for (int i = 0; i < sizestyle; i++)
restart_style[i] = buf[n++];
if (strcmp(restart_style,style) != 0) return 0;
int restart_nreg = buf[n];
n += sizeof(int);
if (restart_nreg != nregion) return 0;
for (int ilist = 0; ilist < nregion; ilist++){
if (!domain->regions[list[ilist]]->restart(buf, n)) return 0;
}
return 1;
}
/* ----------------------------------------------------------------------

View File

@ -35,10 +35,10 @@ class RegUnion : public Region {
void shape_update();
void pretransform();
void set_velocity();
void length_restart_string(int&);
void write_restart(FILE *);
int restart(char *, int);
int restart(char *, int&);
void reset_vel();
private:
char **idsub;
};