diff --git a/doc/src/fix_adapt.rst b/doc/src/fix_adapt.rst index 78677a0e2d..5745c0f69e 100644 --- a/doc/src/fix_adapt.rst +++ b/doc/src/fix_adapt.rst @@ -35,7 +35,7 @@ Syntax v_name = variable with name that calculates value of aparam * zero or more keyword/value pairs may be appended -* keyword = *scale* or *reset* +* keyword = *scale* or *reset* or *mass* .. parsed-literal:: @@ -45,6 +45,9 @@ Syntax *reset* value = *no* or *yes* *no* = values will remain altered at the end of a run *yes* = reset altered values to their original values at the end of a run + *mass* value = *no* or *yes* + *no* = mass is not altered by changes in diameter + *yes* = mass is altered by changes in diameter Examples """""""" @@ -53,7 +56,7 @@ Examples fix 1 all adapt 1 pair soft a 1 1 v_prefactor fix 1 all adapt 1 pair soft a 2* 3 v_prefactor - fix 1 all adapt 1 pair lj/cut epsilon * * v_scale1 coul/cut scale 3 3 v_scale2 scale yes reset yes + fix 1 all adapt 1 pair lj/cut epsilon * * v_scale1 pair coul/cut scale 3 3 v_scale2 scale yes reset yes fix 1 all adapt 10 atom diameter v_size variable ramp_up equal "ramp(0.01,0.5)" @@ -86,12 +89,13 @@ the end of a simulation. Even if *reset* is specified as *yes*\ , a restart file written during a simulation will contain the modified settings. -If the *scale* keyword is set to *no*\ , then the value of the altered -parameter will be whatever the variable generates. If the *scale* -keyword is set to *yes*\ , then the value of the altered parameter will -be the initial value of that parameter multiplied by whatever the -variable generates. I.e. the variable is now a "scale factor" applied -in (presumably) a time-varying fashion to the parameter. +If the *scale* keyword is set to *no*\ , which is the default, then +the value of the altered parameter will be whatever the variable +generates. If the *scale* keyword is set to *yes*\ , then the value +of the altered parameter will be the initial value of that parameter +multiplied by whatever the variable generates. I.e. the variable is +now a "scale factor" applied in (presumably) a time-varying fashion to +the parameter. Note that whether scale is *no* or *yes*\ , internally, the parameters themselves are actually altered by this fix. Make sure you use the @@ -107,16 +111,17 @@ style supports it. Note that the :doc:`pair_style ` and to specify these parameters initially; the fix adapt command simply overrides the parameters. -The *pstyle* argument is the name of the pair style. If :doc:`pair_style hybrid or hybrid/overlay ` is used, *pstyle* should be -a sub-style name. If there are multiple sub-styles using the same -pair style, then *pstyle* should be specified as "style:N" where N is -which instance of the pair style you wish to adapt, e.g. the first, -second, etc. For example, *pstyle* could be specified as "soft" or -"lubricate" or "lj/cut:1" or "lj/cut:2". The *pparam* argument is the -name of the parameter to change. This is the current list of pair -styles and parameters that can be varied by this fix. See the doc -pages for individual pair styles and their energy formulas for the -meaning of these parameters: +The *pstyle* argument is the name of the pair style. If +:doc:`pair_style hybrid or hybrid/overlay ` is used, +*pstyle* should be a sub-style name. If there are multiple +sub-styles using the same pair style, then *pstyle* should be specified +as "style:N" where N is which instance of the pair style you wish to +adapt, e.g. the first, second, etc. For example, *pstyle* could be +specified as "soft" or "lubricate" or "lj/cut:1" or "lj/cut:2". The +*pparam* argument is the name of the parameter to change. This is the +current list of pair styles and parameters that can be varied by this +fix. See the doc pages for individual pair styles and their energy +formulas for the meaning of these parameters: +---------------------------------------------------------------------+--------------------------------------------------+-------------+ | :doc:`born ` | a,b,c | type pairs | @@ -234,31 +239,32 @@ the coefficients for the symmetric J,I interaction to the same values. A wild-card asterisk can be used in place of or in conjunction with the I,J arguments to set the coefficients for multiple pairs of atom -types. This takes the form "\*" or "\*n" or "n\*" or "m\*n". If N = the -number of atom types, then an asterisk with no numeric values means -all types from 1 to N. A leading asterisk means all types from 1 to n -(inclusive). A trailing asterisk means all types from n to N +types. This takes the form "\*" or "\*n" or "n\*" or "m\*n". If N = +the number of atom types, then an asterisk with no numeric values +means all types from 1 to N. A leading asterisk means all types from +1 to n (inclusive). A trailing asterisk means all types from n to N (inclusive). A middle asterisk means all types from m to n (inclusive). Note that only type pairs with I <= J are considered; if asterisks imply type pairs where J < I, they are ignored. -IMPROTANT NOTE: If :doc:`pair_style hybrid or hybrid/overlay ` is being used, then the *pstyle* will -be a sub-style name. You must specify I,J arguments that correspond -to type pair values defined (via the :doc:`pair_coeff ` -command) for that sub-style. +IMPROTANT NOTE: If :doc:`pair_style hybrid or hybrid/overlay +` is being used, then the *pstyle* will be a sub-style +name. You must specify I,J arguments that correspond to type pair +values defined (via the :doc:`pair_coeff ` command) for +that sub-style. The *v_name* argument for keyword *pair* is the name of an -:doc:`equal-style variable ` which will be evaluated each time -this fix is invoked to set the parameter to a new value. It should be -specified as v_name, where name is the variable name. Equal-style -variables can specify formulas with various mathematical functions, -and include :doc:`thermo_style ` command keywords for the -simulation box parameters and timestep and elapsed time. Thus it is -easy to specify parameters that change as a function of time or span -consecutive runs in a continuous fashion. For the latter, see the -*start* and *stop* keywords of the :doc:`run ` command and the -*elaplong* keyword of :doc:`thermo_style custom ` for -details. +:doc:`equal-style variable ` which will be evaluated each +time this fix is invoked to set the parameter to a new value. It +should be specified as v_name, where name is the variable name. +Equal-style variables can specify formulas with various mathematical +functions, and include :doc:`thermo_style ` command +keywords for the simulation box parameters and timestep and elapsed +time. Thus it is easy to specify parameters that change as a function +of time or span consecutive runs in a continuous fashion. For the +latter, see the *start* and *stop* keywords of the :doc:`run ` +command and the *elaplong* keyword of :doc:`thermo_style custom +` for details. For example, these commands would change the prefactor coefficient of the :doc:`pair_style soft ` potential from 10.0 to 30.0 in a @@ -319,23 +325,28 @@ The *atom* keyword enables various atom properties to be changed. The current list of atom parameters that can be varied by this fix: * charge = charge on particle -* diameter, or, diameter/disc = diameter of particle +* diameter or diameter/disc = diameter of particle The *v_name* argument of the *atom* keyword is the name of an -:doc:`equal-style variable ` which will be evaluated each time -this fix is invoked to set, or scale the parameter to a new value. -It should be specified as v_name, where name is the variable name. See the -discussion above describing the formulas associated with equal-style -variables. The new value is assigned to the corresponding attribute -for all atoms in the fix group. +:doc:`equal-style variable ` which will be evaluated each +time this fix is invoked to set, or scale the parameter to a new +value. It should be specified as v_name, where name is the variable +name. See the discussion above describing the formulas associated +with equal-style variables. The new value is assigned to the +corresponding attribute for all atoms in the fix group. If the atom parameter is *diameter* and per-atom density and per-atom -mass are defined for particles (e.g. :doc:`atom_style granular `), then the mass of each particle is also -changed when the diameter changes. The mass is set from the particle volume -for 3d systems (density is assumed to stay constant). For 2d, the default is -for LAMMPS to model particles with a radius attribute as spheres. -However, if the atom parameter is *diameter/disc*, then the mass is -set from the particle area (the density is assumed to be in mass/distance^2 units). +mass are defined for particles (e.g. :doc:`atom_style granular +`), then the mass of each particle is, by default, also +changed when the diameter changes. The mass is set from the particle +volume for 3d systems (density is assumed to stay constant). For 2d, +the default is for LAMMPS to model particles with a radius attribute +as spheres. However, if the atom parameter is *diameter/disc*, then the +mass is set from the particle area (the density is assumed to be in +mass/distance^2 units). The mass of the particle may also be kept constant +if the *mass* keyword is set to *no*. This can be useful to account for +diameter changes that do not involve mass changes, e.g., thermal expansion. + For example, these commands would shrink the diameter of all granular particles in the "center" group from 1.0 to 0.1 in a linear fashion @@ -348,13 +359,63 @@ over the course of a 1000-step simulation: ---------- +This fix can be used in long simulations which are restarted one or +more times to continuously adapt simulation parameters, but it must be +done carefully. There are two issues to consider. The first is how +to adapt the parameters in a continuous manner from one simulation to +the next. The second is how, if desired, to reset the parameters to +their original values at the end of the last restarted run. + +Note that all the parameters changed by this fix are written into a +restart file in their current changed state. A new restarted +simulation does not know their original time=0 values, unless the +input script explicitly resets the parameters (after the restart file +is read), to their original values. + +Also note, that the time-dependent variable(s) used in the restart +script should typically be written as a function of time elapsed since +the original simulation began. + +With this in mind, if the *scale* keyword is set to *no* (the default) +in a restarted simulation, original parameters are not needed. The +adapted parameters should seamlessly continue their variation relative +to the preceding simulation. + +If the *scale* keyword is set to *yes*, then the input script should +typically reset the parameters being adapted to their original values, +so that the scaling formula specified by the variable will operate +correctly. An exception is if the *atom* keyword is being used with +*scale yes*. In this case, information is added to the restart file +so that per-atom properties in the new run will automatically be +scaled relative to their original values. This will only work if the +fix adapt command specified in the restart script has the same ID as +the one used in the original script. + +In a restarted run, if the *reset* keyword is set to *yes*, and the +run ends in this script (as opposed to just writing more restart +files, parameters will be restored to the values they were at the +beginning of the run command in the restart script. Which as +explained above, may or may not be the original values of the +parameters. Again, an exception is if the *atom* keyword is being +used with *reset yes* (in all the runs). In that case, the original +per-atom parameters are stored in the restart file, and will be +restored when the restarted run finally completes. + +---------- + **Restart, fix_modify, output, run start/stop, minimize info:** -No information about this fix is written to :doc:`binary restart files `. None of the :doc:`fix_modify ` options -are relevant to this fix. No global or per-atom quantities are stored -by this fix for access by various :doc:`output commands `. -No parameter of this fix can be used with the *start/stop* keywords of -the :doc:`run ` command. This fix is not invoked during :doc:`energy minimization `. +If the *atom* keyword is used and the *scale* or *reset* keyword is +set to *yes*, then this fix writes information to a restart file so +that in a restarted run scaling can continue in a seamless manner +and/or the per-atom values can be restored, as explained above. + +None of the :doc:`fix_modify ` options are relevant to +this fix. No global or per-atom quantities are stored by this fix for +access by various :doc:`output commands `. No parameter +of this fix can be used with the *start/stop* keywords of the +:doc:`run ` command. This fix is not invoked during :doc:`energy +minimization `. For :doc:`rRESPA time integration `, this fix changes parameters on the outermost rRESPA level. @@ -371,4 +432,4 @@ Related commands Default """"""" -The option defaults are scale = no, reset = no. +The option defaults are scale = no, reset = no, mass = yes. diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp index f241d4ea36..e7254e77d7 100644 --- a/src/fix_adapt.cpp +++ b/src/fix_adapt.cpp @@ -108,6 +108,7 @@ nadapt(0), id_fix_diam(NULL), id_fix_chg(NULL), adapt(NULL) } else error->all(FLERR,"Illegal fix adapt command"); nadapt++; iarg += 6; + } else if (strcmp(arg[iarg],"bond") == 0 ){ if (iarg+5 > narg) error->all(FLERR, "Illegal fix adapt command"); adapt[nadapt].which = BOND; @@ -127,6 +128,7 @@ nadapt(0), id_fix_diam(NULL), id_fix_chg(NULL), adapt(NULL) } else error->all(FLERR,"Illegal fix adapt command"); nadapt++; iarg += 5; + } else if (strcmp(arg[iarg],"kspace") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); adapt[nadapt].which = KSPACE; @@ -137,10 +139,12 @@ nadapt(0), id_fix_diam(NULL), id_fix_chg(NULL), adapt(NULL) } else error->all(FLERR,"Illegal fix adapt command"); nadapt++; iarg += 2; + } else if (strcmp(arg[iarg],"atom") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command"); adapt[nadapt].which = ATOM; - if (strcmp(arg[iarg+1],"diameter") == 0 || strcmp(arg[iarg+1],"diameter/disc") == 0) { + if (strcmp(arg[iarg+1],"diameter") == 0 || + strcmp(arg[iarg+1],"diameter/disc") == 0) { adapt[nadapt].aparam = DIAMETER; diamflag = 1; discflag = 0; @@ -163,6 +167,7 @@ nadapt(0), id_fix_diam(NULL), id_fix_chg(NULL), adapt(NULL) resetflag = 0; scaleflag = 0; + massflag = 1; while (iarg < narg) { if (strcmp(arg[iarg],"reset") == 0) { @@ -177,9 +182,22 @@ nadapt(0), id_fix_diam(NULL), id_fix_chg(NULL), adapt(NULL) else if (strcmp(arg[iarg+1],"yes") == 0) scaleflag = 1; else error->all(FLERR,"Illegal fix adapt command"); iarg += 2; + } else if (strcmp(arg[iarg],"mass") == 0) { + if (iarg+2 > narg)error->all(FLERR,"Illegal fix adapt command"); + if (strcmp(arg[iarg+1],"no") == 0) massflag = 0; + else if (strcmp(arg[iarg+1],"yes") == 0) massflag = 1; + else error->all(FLERR,"Illegal fix adapt command"); + iarg += 2; } else error->all(FLERR,"Illegal fix adapt command"); } + // if scaleflag set with diameter or charge adaptation, + // then previous step scale factors are written to restart file + // initialize them here in case one is used and other is never defined + + if (scaleflag && (diamflag || chgflag)) restart_global = 1; + previous_diam_scale = previous_chg_scale = 1.0; + // allocate pair style arrays int n = atom->ntypes; @@ -238,6 +256,7 @@ int FixAdapt::setmask() void FixAdapt::post_constructor() { + if (!resetflag) return; if (!diamflag && !chgflag) return; // new id = fix-ID + FIX_STORE_ATTRIBUTE @@ -432,15 +451,19 @@ void FixAdapt::init() error->all(FLERR,"Fix adapt requires atom attribute diameter"); if (!atom->rmass_flag) error->all(FLERR,"Fix adapt requires atom attribute mass"); - if (discflag && domain->dimension!=2) + if (discflag && domain->dimension != 2) error->all(FLERR,"Fix adapt requires 2d simulation"); + if (!restart_reset) previous_diam_scale = 1.0; } if (ad->aparam == CHARGE) { if (!atom->q_flag) error->all(FLERR,"Fix adapt requires atom attribute charge"); + if (!restart_reset) previous_chg_scale = 1.0; } } } + + if (restart_reset) restart_reset = 0; // make copy of original pair/bond array values @@ -570,48 +593,54 @@ void FixAdapt::change_settings() } else if (ad->which == ATOM) { - // reset radius from diameter - // also scale rmass to new value - + // reset radius to new value, for both owned and ghost atoms + // also reset rmass to new value assuming density remains constant + // for scaleflag, previous_diam_scale is the scale factor on previous step + if (ad->aparam == DIAMETER) { - double density; - - // Get initial diameter if `scale` keyword is used - - double *vec = fix_diam->vstore; + double density,scale; double *radius = atom->radius; double *rmass = atom->rmass; int *mask = atom->mask; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; - for (i = 0; i < nall; i++) + if (scaleflag) scale = value / previous_diam_scale; + + for (i = 0; i < nall; i++) { if (mask[i] & groupbit) { - if (discflag) density = rmass[i] / (MY_PI * radius[i]*radius[i]); - else density = rmass[i] / (4.0*MY_PI/3.0 * - radius[i]*radius[i]*radius[i]); - if (scaleflag) radius[i] = value * vec[i]; + if (massflag) { + if (!scaleflag) scale = 0.5*value / radius[i]; + if (discflag) rmass[i] *= scale*scale; + else rmass[i] *= scale*scale*scale; + } + if (scaleflag) radius[i] *= scale; else radius[i] = 0.5*value; - if (discflag) rmass[i] = MY_PI * radius[i]*radius[i] * density; - else rmass[i] = 4.0*MY_PI/3.0 * - radius[i]*radius[i]*radius[i] * density; } + } + + if (scaleflag) previous_diam_scale = value; + + // reset charge to new value, for both owned and ghost atoms + // for scaleflag, previous_chg_scale is the scale factor on previous step } else if (ad->aparam == CHARGE) { - - // Get initial charge if `scale` keyword is used - - double *vec = fix_chg->vstore; + double scale; double *q = atom->q; int *mask = atom->mask; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; - for (i = 0; i < nall; i++) + if (scaleflag) scale = value / previous_chg_scale; + + for (i = 0; i < nall; i++) { if (mask[i] & groupbit) { - if (scaleflag) q[i] = value * vec[i]; + if (scaleflag) q[i] *= scale; else q[i] = value; } + } + + if (scaleflag) previous_chg_scale = value; } } } @@ -672,23 +701,24 @@ void FixAdapt::restore_settings() } else if (ad->which == ATOM) { if (diamflag) { - double density; + double density,scale; double *vec = fix_diam->vstore; double *radius = atom->radius; double *rmass = atom->rmass; int *mask = atom->mask; int nlocal = atom->nlocal; + + if (scaleflag) scale = previous_diam_scale; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if(discflag) density = rmass[i] / (MY_PI * radius[i]*radius[i]); - else density = rmass[i] / (4.0*MY_PI/3.0 * - radius[i]*radius[i]*radius[i]); + if (massflag) { + if (!scaleflag) scale = vec[i] / radius[i]; + if (discflag) rmass[i] *= scale*scale; + else rmass[i] *= scale*scale*scale; + } radius[i] = vec[i]; - if(discflag) rmass[i] = MY_PI * radius[i]*radius[i] * density; - else rmass[i] = 4.0*MY_PI/3.0 * - radius[i]*radius[i]*radius[i] * density; } } if (chgflag) { @@ -717,3 +747,28 @@ void FixAdapt::set_arrays(int i) if (fix_diam) fix_diam->vstore[i] = atom->radius[i]; if (fix_chg) fix_chg->vstore[i] = atom->q[i]; } + +/* ---------------------------------------------------------------------- + write scale factors for diameter and charge to restart file +------------------------------------------------------------------------- */ + +void FixAdapt::write_restart(FILE *fp) +{ + int size = 2*sizeof(double); + + fwrite(&size,sizeof(int),1,fp); + fwrite(&previous_diam_scale,sizeof(double),1,fp); + fwrite(&previous_chg_scale,sizeof(double),1,fp); +} + +/* ---------------------------------------------------------------------- + use scale factors from restart file to restart the Fix +------------------------------------------------------------------------- */ + +void FixAdapt::restart(char *buf) +{ + double *dbuf = (double *) buf; + + previous_diam_scale = dbuf[0]; + previous_chg_scale = dbuf[1]; +} diff --git a/src/fix_adapt.h b/src/fix_adapt.h index 26def03a1d..1872e8b3a0 100644 --- a/src/fix_adapt.h +++ b/src/fix_adapt.h @@ -40,13 +40,16 @@ class FixAdapt : public Fix { void setup_pre_force_respa(int,int); void pre_force_respa(int,int,int); void set_arrays(int); + void write_restart(FILE *); + void restart(char *); private: - int nadapt,resetflag,scaleflag; + int nadapt,resetflag,scaleflag,massflag; int anypair, anybond; int nlevels_respa; char *id_fix_diam,*id_fix_chg; class FixStore *fix_diam,*fix_chg; + double previous_diam_scale,previous_chg_scale; int discflag; struct Adapt { diff --git a/src/modify.cpp b/src/modify.cpp index f367def0f4..f7fd18d55c 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -909,6 +909,7 @@ void Modify::add_fix(int narg, char **arg, int trysuffix) strcmp(style_restart_global[i],fix[ifix]->style) == 0) { fix[ifix]->restart(state_restart_global[i]); used_restart_global[i] = 1; + fix[ifix]->restart_reset = 1; if (comm->me == 0) { if (screen) fprintf(screen,"Resetting global fix info from restart file:\n");