creating clean fix npt/cauchy package and docs

fix nvt/kk command :h3

fix nvt/omp command :h3

fix npt command :h3

fix npt/cauchy command :h3

fix npt/intel command :h3

fix npt/kk command :h3

fix npt/omp command :h3

ID, group-ID are documented in "fix"_fix.html command :ulb,l

style_name = {nvt} or {npt} or {nph} :l

one or more keyword/value pairs may be appended :l

keyword = {temp} or {iso} or {aniso} or {tri} or {x} or {y} or {z} or {xy} or {yz} or {xz} or {couple} or {tchain} or {pchain} or {mtk} or {tloop} or {ploop} or {nreset} or {drag} or {dilate} or {scalexy} or {scaleyz} or {scalexz} or {flip} or {fixedpoint} or {update}

keyword = {temp} or {iso} or {aniso} or {tri} or {x} or {y} or {z} or {xy} or {yz} or {xz} or {couple} or {tchain} or {pchain} or {mtk} or {tloop} or {ploop} or {nreset} or {drag} or {dilate} or {scalexy} or {scaleyz} or {scalexz} or {flip} or {fixedpoint} or {update} or {alpha} or {continue}

  {temp} values = Tstart Tstop Tdamp

    Tstart,Tstop = external temperature at start/end of run

    Tdamp = temperature damping parameter (time units)

    x,y,z = perform barostat dilation/contraction around this point (distance units)

  {update} value = {dipole} or {dipole/dlm}

    dipole = update dipole orientation (only for sphere variants)

    dipole/dlm = use DLM integrator to update dipole orientation (only for sphere variants) :pre

    dipole/dlm = use DLM integrator to update dipole orientation (only for sphere variants)

   {alpha} value = alpha

    alpha = strength of Cauchystat control parameter.  This keyword is only recognized by npt/cauchy and requires the USER-CAUCHY package.

   {continue} value = {yes} or {no} = continue Cauchystat run using previously converged parameters.  This keyword is only recognized by npt/cauchy and requires the USER-CAUCHY package.

:pre

:ule

[Examples:]

:line

The thermostat parameters for fix styles {nvt} and {npt} is specified

The thermostat parameters for fix styles {nvt} and {npt} are specified

using the {temp} keyword.  Other thermostat-related keywords are

{tchain}, {tloop} and {drag}, which are discussed below.

{tchain}, {tloop}, and {drag}, which are discussed below.

The thermostat is applied to only the translational degrees of freedom

for the particles.  The translational degrees of freedom can also have

simulation.

Other barostat-related keywords are {pchain}, {mtk}, {ploop},

{nreset}, {drag}, and {dilate}, which are discussed below.

{nreset}, {drag}, {dilate}, {alpha} and {continue}, which are discussed below.

Orthogonal simulation boxes have 3 adjustable dimensions (x,y,z).

Triclinic (non-orthogonal) simulation boxes have 6 adjustable

Suzuki-Yoshida scheme "(Tuckerman)"_#nh-Tuckerman.  The keyword {ploop}

does the same thing for the barostat thermostat.

:line

By default, the barostat algorithm controls the Second-Piola Kirchhoff

stress, which is a stress measure referred to the undeformed (initial)

simulation box.  If the box deforms substantially during the

equilibration, the difference between the set values and the final

true (Cauchy) stresses can be considerable.  There are two ways to

deal with this: the {nreset} keyword and the {npt/cauchy} fix style

provided by the USER-CAUCHY package.

The keyword {nreset} controls how often the reference dimensions used

to define the strain energy are reset.  If this keyword is not used,

or is given a value of zero, then the reference dimensions are set to

means that every {nstep} timesteps, the reference dimensions are set

to those of the current simulation domain.

Installing the USER-CAUCHY package enables the {npt/cauchy} fix, which

implements the Cauchystat as per Miller et al. (Miller)_"#nh-Miller"

to directly control the Cauchy stress.  For this fix style, {alpha} is a

non-dimensional parameter, typically set to 0.001 or 0.01, that

determines how aggresively the algorithm drives the system towards the

set Cauchy stresses.  Larger values of {alpha} will modify the system

more quickly, but can lead to instabilities.  Smaller values will lead

to longer convergence time.  Since {alpha} also influences how much

the stress fluctuations deviate from the equilibrium fluctuations, it

should be set as small as possible.

Setting {alpha} to zero is not permitted.  To "turn off" the

Cauchystat control and thus restore the equilibrium stress

fluctuations, two subsequent fixes should be used.  In the first, the

{npt/cauchy} fix is used and the simulation box equilibrates to the correct

shape for the desired stresses.  In the second, the regular {npt} fix is

used, and it will restore the original

Parrinello-Rahman algorithm, but now with the correct simulation box

shape determined in the first fix.

A {continue} value of {yes} is used for the {npt/cauchy} fix to signal that

the fix is subsequent to a previous Cauchystat run, and

the intention is to continue from the converged stress state at the

end of the previous run.  This may be required, for example, when

implementing a multi-step loading/unloading sequence over several

fixes.

The {scaleyz}, {scalexz}, and {scalexy} keywords control whether or

not the corresponding tilt factors are scaled with the associated box

dimensions when barostatting triclinic periodic cells.  The default

chosen value, except for a time-dependent constant translation of

positions.

:line

If the {update} keyword is used with the {dipole} value, then the

orientation of the dipole moment of each particle is also updated

during the time integration.  This option should be used for models

The keyword defaults are tchain = 3, pchain = 3, mtk = yes, tloop =

ploop = 1, nreset = 0, drag = 0.0, dilate = all, couple = none,

alpha = 0.001, continue = no,

scaleyz = scalexz = scalexy = yes if periodic in 2nd dimension and

not coupled to barostat, otherwise no.

:link(nh-Dullweber)

[(Dullweber)] Dullweber, Leimkuhler and McLachlan, J Chem Phys, 107,

5840 (1997).

:link(nh-Miller)

[(Miller)] Miller, Tadmor, Gibson, Bernstein and Pavia, J Chem Phys,

144, 184107 (2016).

                     y 0.0 0.0 0.1 &

                     z 0.0 0.0 0.1 &

                     xy -10000.0 -10000.0 0.1 &

                     couple none cauchystat 0.001 yes

                     couple none alpha 0.001 continue no

info fixes

                     y 0.0 0.0 0.1 &

                     z 0.0 0.0 0.1 &

                     xy -10000.0 -10000.0 0.1 &

                     couple none

                     couple none alpha 0.001 continue yes

info fixes

run 10000

  // for CauchyStat

  usePK = 1;

  alpha=0.001;

  initRUN = 0;

  restartPK = 0;

  restart_global = 1;

        dlm_flag = 1;

      } else error->all(FLERR,"Illegal fix nvt/npt/nph command");

      iarg += 2;

    } else if (strcmp(arg[iarg],"cauchystat") == 0) {

      usePK = 0;

      if (iarg+3 > narg) error->all(FLERR,

				    "Illegal fix npt cauchystat command:"

				    " wrong number of arguments");

      if (strcmp(arg[iarg+2],"yes") != 0 && strcmp(arg[iarg+2],"no") != 0)

    } else if (strcmp(arg[iarg],"alpha") == 0) {

      alpha = force->numeric(FLERR,arg[iarg+1]);

      iarg += 2;

    } else if (strcmp(arg[iarg],"continue") == 0) {

      if (strcmp(arg[iarg+1],"yes") != 0 && strcmp(arg[iarg+1],"no") != 0)

        error->all(FLERR,"Illegal cauchystat continue value.  "

                   "Must be 'yes' or 'no'");

      alpha = force->numeric(FLERR,arg[iarg+1]);

      restartPK = !strcmp(arg[iarg+2],"yes");

      iarg += 3;

      restartPK = !strcmp(arg[iarg+1],"yes");

      iarg += 2;

    } else if (strcmp(arg[iarg],"fixedpoint") == 0) {

      if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");

      fixedpoint[0] = force->numeric(FLERR,arg[iarg+1]);

void FixCauchy::post_constructor()

{

  if (!usePK) CauchyStat_init();

  else CauchyStat_cleanup();

  CauchyStat_init();

}

/* ---------------------------------------------------------------------- */

  // CauchyStat: call CauchyStat to modify p_target[i] and p_hydro,

  // if CauchyStat enabled and pressure->vector computation has been initiated

  if ((usePK == 0) && (initRUN == 1)) CauchyStat();

  if (initRUN == 1) CauchyStat();

  if (initRUN == 0) {

    for (int i=0; i < 6; i++) {

      h_old[i]=domain->h[i];

#undef invH0

}

/* ----------------------------------------------------------------------

   cleanup when not using Cauchystat

------------------------------------------------------------------------- */

void FixCauchy::CauchyStat_cleanup()

{

  if (id_store) {

    modify->delete_fix(id_store);

    delete[] id_store;

  } else {

    int n = strlen(id) + 14;

    id_store = new char[n];

    strcpy(id_store,id);

    strcat(id_store,"_FIX_NH_STORE");

    modify->delete_fix(id_store);

    delete[] id_store;

  }

  id_store = NULL;

  restart_stored = -1;

}

/* ----------------------------------------------------------------------

   Cauchystat

------------------------------------------------------------------------- */

  double alpha;         // integration parameter for the cauchystat

  int initPK;           // 1 if setPK needs to be initialized either

                        // from cauchy or restart, else 0

  int usePK;            // 0 if use CauchyStat else 1

  int restartPK;        // Read PK stress from the previous run

  int restart_stored;   // values of PK stress from the previous step stored

  int initRUN;          // 0 if run not initialized