diff --git a/doc/src/fix_adapt.rst b/doc/src/fix_adapt.rst index fbc4aa4661..5ba9ed0d26 100644 --- a/doc/src/fix_adapt.rst +++ b/doc/src/fix_adapt.rst @@ -53,7 +53,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 +86,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 +108,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-stylesusing 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 +236,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 +322,25 @@ 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 +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). +set from the particle area (the density is assumed to be in +mass/distance^2 units). 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 +353,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 explicity 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 preceeding 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. diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp index a29e1bfebb..ef2cb3e382 100644 --- a/src/fix_adapt.cpp +++ b/src/fix_adapt.cpp @@ -50,7 +50,6 @@ nadapt(0), id_fix_diam(NULL), id_fix_chg(NULL), adapt(NULL) dynamic_group_allow = 1; create_attribute = 1; - restart_global = 1; // count # of adaptations @@ -109,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; @@ -128,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; @@ -138,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; @@ -181,6 +184,13 @@ nadapt(0), id_fix_diam(NULL), id_fix_chg(NULL), adapt(NULL) } 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; @@ -434,17 +444,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) diam_scale = 1.0; + 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) chg_scale = 1.0; + if (!restart_reset) previous_chg_scale = 1.0; } } } + + if (restart_reset) restart_reset = 0; // make copy of original pair/bond array values @@ -479,8 +491,6 @@ void FixAdapt::init() if (strstr(update->integrate_style,"respa")) nlevels_respa = ((Respa *) update->integrate)->nlevels; - - if(restart_reset) restart_reset = 0; } /* ---------------------------------------------------------------------- */ @@ -576,43 +586,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; - + 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 (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 * radius[i] / diam_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) diam_scale = value; + 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] *= 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) { + 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 * q[i] / chg_scale; + if (scaleflag) q[i] *= scale; else q[i] = value; } - if (scaleflag) chg_scale = value; + } + + if (scaleflag) previous_chg_scale = value; } } } @@ -683,11 +704,11 @@ void FixAdapt::restore_settings() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if(discflag) density = rmass[i] / (MY_PI * radius[i]*radius[i]); + 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]); radius[i] = vec[i]; - if(discflag) rmass[i] = MY_PI * radius[i]*radius[i] * density; + 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; } @@ -728,11 +749,10 @@ void FixAdapt::write_restart(FILE *fp) int size = 2*sizeof(double); fwrite(&size,sizeof(int),1,fp); - fwrite(&diam_scale,sizeof(double),1,fp); - fwrite(&chg_scale,sizeof(double),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 ------------------------------------------------------------------------- */ @@ -741,6 +761,6 @@ void FixAdapt::restart(char *buf) { double *dbuf = (double *) buf; - diam_scale = dbuf[0]; - chg_scale = dbuf[1]; -} \ No newline at end of file + previous_diam_scale = dbuf[0]; + previous_chg_scale = dbuf[1]; +} diff --git a/src/fix_adapt.h b/src/fix_adapt.h index 1664a09b11..390a8f47a8 100644 --- a/src/fix_adapt.h +++ b/src/fix_adapt.h @@ -49,7 +49,7 @@ class FixAdapt : public Fix { int nlevels_respa; char *id_fix_diam,*id_fix_chg; class FixStore *fix_diam,*fix_chg; - double diam_scale,chg_scale; + double previous_diam_scale,previous_chg_scale; int discflag; struct Adapt {