additional fix external hooks for calling programs

dir = {x} or {y} or {z} :l

shockvel = shock velocity (strictly positive, distance/time units) :l

zero or more keyword value pairs may be appended :l

keyword = {q} or {mu} or {p0} or {v0} or {e0} or {tscale} :l

keyword = {q} or {mu} or {p0} or {v0} or {e0} or {tscale} or {beta} or {dftb} :l

  {q} value = cell mass-like parameter (mass^2/distance^4 units)

  {mu} value = artificial viscosity (mass/length/time units)

  {p0} value = initial pressure in the shock equations (pressure units)

  {v0} value = initial simulation cell volume in the shock equations (distance^3 units)

  {e0} value = initial total energy (energy units)

  {tscale} value = reduction in initial temperature (unitless fraction between 0.0 and 1.0) :pre

  {tscale} value = reduction in initial temperature (unitless fraction between 0.0 and 1.0) 

  {dftb} value = {yes} or {no} for whether using MSST in conjunction with DFTB+

  {beta} value = scale factor on energy contribution of DFTB+ :pre

:ule

[Examples:]

fix 1 all msst y 100.0 q 1.0e5 mu 1.0e5

fix 2 all msst z 50.0 q 1.0e4 mu 1.0e4  v0 4.3419e+03 p0 3.7797e+03 e0 -9.72360e+02 tscale 0.01 :pre

fix 2 all msst z 50.0 q 1.0e4 mu 1.0e4  v0 4.3419e+03 p0 3.7797e+03 e0 -9.72360e+02 tscale 0.01

fix 1 all msst y 100.0 q 1.0e5 mu 1.0e5 dftb yes beta 0.5 :pre

[Description:]

symmetry to equilibrate to the shock Hugoniot and Rayleigh line more

rapidly in such cases.

{tscale} is a factor between 0 and 1 that determines what fraction of

thermal kinetic energy is converted to compressive strain kinetic

energy at the start of the simulation.  Setting this parameter to a

non-zero value may assist in compression at the start of simulations

where it is slow to occur.

The keyword {tscale} is a factor between 0 and 1 that determines what

fraction of thermal kinetic energy is converted to compressive strain

kinetic energy at the start of the simulation.  Setting this parameter

to a non-zero value may assist in compression at the start of

simulations where it is slow to occur.

If keywords {e0}, {p0},or {v0} are not supplied, these quantities will

be calculated on the first step, after the energy specified by

Parrinello-Rahman boundary conditions (tilted box) are supported by

LAMMPS, but are not implemented for MSST.

This fix computes a temperature and pressure each timestep. To do

this, the fix creates its own computes of style "temp" and "pressure",

as if these commands had been issued:

This fix computes a temperature and pressure and potential energy each

timestep. To do this, the fix creates its own computes of style "temp"

"pressure", and "pe", as if these commands had been issued:

compute fix-ID_temp group-ID temp

compute fix-ID_press group-ID pressure fix-ID_temp :pre

compute fix-ID_MSST_temp all temp

compute fix-ID_MSST_press all pressure fix-ID_MSST_temp :pre

compute fix-ID_MSST_pe all pe :pre

See the "compute temp"_compute_temp.html and "compute

pressure"_compute_pressure.html commands for details.  Note that the

IDs of the new computes are the fix-ID + underscore + "temp" or fix_ID

+ underscore + "press".  The group for the new computes is "all".

IDs of the new computes are the fix-ID + "_MSST_temp" or "_MSST_press"

or "_MSST_pe".  The group for the new computes is "all".

:line

The {dftb} and {beta} keywords are to allow this fix to be used when

LAMMPS is being driven by DFTB+, a density-functional tight-binding

code.

If the keyword {dftb} is used with a value of {yes}, then the MSST

equations are altered to account for an energy contribution compute by

DFTB+.  In this case, you must define a "fix

external"_fix_external.html command in your input script, which is

used to callback to DFTB+ during the LAMMPS timestepping.  DFTB+ will

communicate its info to LAMMPS via that fix.

The keyword {beta} is a scale factor on the DFTB+ energy contribution.

The value of {beta} must be between 0.0 and 1.0 inclusive.  A value of

0.0 means no contribution, a value of 1.0 means a full contribution.

(July 2017) More information about these keywords and the use of

LAMMPS with DFTB+ will be added to the LAMMMPS documention soon.

:line

[Restart, fix_modify, output, run start/stop, minimize info:]

[Default:]

The keyword defaults are q = 10, mu = 0, tscale = 0.01. p0, v0, and e0

are calculated on the first step.

The keyword defaults are q = 10, mu = 0, tscale = 0.01, dftb = no,

beta = 0.0.  Note that p0, v0, and e0 are calculated on the first

timestep.

:line

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

   Contributing authors: Laurence Fried (LLNL), Evan Reed (LLNL, Stanford)

     implementation of the Multi-Scale Shock Method

   See Reed, Fried, Joannopoulos, Phys Rev Lett, 90, 235503 (2003)

   see Reed, Fried, Joannopoulos, Phys Rev Lett, 90, 235503 (2003)

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

#include <string.h>

#include "comm.h"

#include "output.h"

#include "modify.h"

#include "fix_external.h"

#include "compute.h"

#include "kspace.h"

#include "update.h"

#include "respa.h"

#include "domain.h"

#include "error.h"

#include "thermo.h"

#include "memory.h"

#include "error.h"

using namespace LAMMPS_NS;

using namespace FixConst;

FixMSST::FixMSST(LAMMPS *lmp, int narg, char **arg) :

  Fix(lmp, narg, arg), old_velocity(NULL), rfix(NULL), 

  id_temp(NULL), id_press(NULL), id_pe(NULL), temperature(NULL), 

  pressure(NULL), pe(NULL), atoms_allocated(0)

  pressure(NULL), pe(NULL)

{

  if (narg < 4) error->all(FLERR,"Illegal fix msst command");

  p0 = 0.0;

  v0 = 1.0;

  e0 = 0.0;

  TS_int = 0;

  T0S0 = 0.0;

  S_elec = 0.0;

  S_elec_1 = 0.0;

  S_elec_2 = 0.0;

  qmass = 1.0e1;

  mu = 0.0;

  v0_set = 0;

  e0_set = 0;

  tscale = 0.01;

  dftb = 0;

  beta = 0.0;

  if ( strcmp(arg[3],"x") == 0 )

    direction = 0;

  else if ( strcmp(arg[3],"y") == 0 )

    direction = 1;

  else if ( strcmp(arg[3],"z") == 0 )

    direction = 2;

  else {

    error->all(FLERR,"Illegal fix msst command");

  }

  if (strcmp(arg[3],"x") == 0) direction = 0;

  else if (strcmp(arg[3],"y") == 0) direction = 1;

  else if (strcmp(arg[3],"z") == 0) direction = 2;

  else error->all(FLERR,"Illegal fix msst command");

  velocity = force->numeric(FLERR,arg[4]);

  if ( velocity < 0 )

    error->all(FLERR,"Illegal fix msst command");

  if (velocity < 0) error->all(FLERR,"Illegal fix msst command");

  for ( int iarg = 5; iarg < narg; iarg++ ) {

  // optional args

  int iarg = 5;

  while (iarg < narg) {

    if (strcmp(arg[iarg],"q") == 0) {

      if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command");

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

      iarg++;

      iarg += 2;

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

      if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command");

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

      iarg++;

      iarg += 2;

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

      if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command");

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

      iarg++;

      p0_set = 1;

      iarg += 2;

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

      if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command");

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

      v0_set = 1;

      iarg++;

      iarg += 2;

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

      if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command");

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

      e0_set = 1;

      iarg++;

      iarg += 2;

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

      if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command");

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

      if (tscale < 0.0 || tscale > 1.0)

        error->all(FLERR,"Fix msst tscale must satisfy 0 <= tscale < 1");

      iarg++;

      iarg += 2;

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

      if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command");

      if (strcmp(arg[iarg+1],"yes") == 0) dftb = 1;

      else if (strcmp(arg[iarg+1],"yes") == 0) dftb = 0;

      else error->all(FLERR,"Illegal fix msst command");

      iarg += 2;

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

      if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command");

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

      if (beta < 0.0 || beta > 1.0) 

        error->all(FLERR,"Illegal fix msst command");

      iarg += 2;

    } else error->all(FLERR,"Illegal fix msst command");

  }

  // output MSST info

  if (comm->me == 0) {

    if (screen) {

      fprintf(screen,"MSST parameters:\n");

  if (domain->nonperiodic) error->all(FLERR,"Fix msst requires a periodic box");

  // create a new compute temp style

  // id = fix-ID + temp

  // create a new temperature compute

  // id = fix-ID + "MSST_temp"

  // compute group = all since pressure is always global (group all)

  //   and thus its KE/temperature contribution should use group all

  int n = strlen(id) + 6;

  int n = strlen(id) + 10;

  id_temp = new char[n];

  strcpy(id_temp,id);

  strcat(id_temp,"_temp");

  strcat(id_temp,"MSST_temp");

  char **newarg = new char*[3];

  newarg[0] = id_temp;

  delete [] newarg;

  tflag = 1;

  // create a new compute pressure style

  // id = fix-ID + press, compute group = all

  // create a new pressure compute

  // id = fix-ID + "MSST_press", compute group = all

  // pass id_temp as 4th arg to pressure constructor

  n = strlen(id) + 7;

  n = strlen(id) + 11;

  id_press = new char[n];

  strcpy(id_press,id);

  strcat(id_press,"_press");

  strcat(id_press,"MSST_press");

  newarg = new char*[4];

  newarg[0] = id_press;

  delete [] newarg;

  pflag = 1;

  // create a new compute potential energy compute

  // create a new potential energy compute

  // id = fix-ID + "MSST_pe", compute group = all

  n = strlen(id) + 4;

  n = strlen(id) + 8;

  id_pe = new char[n];

  strcpy(id_pe,id);

  strcat(id_pe,"_pe");

  strcat(id_pe,"MSST_pe");

  newarg = new char*[3];

  newarg[0] = id_pe;

  delete [] newarg;

  peflag = 1;

  // initialize the time derivative of the volume.

  omega[0] = omega[1] = omega[2] = 0.0;

  // initialize the time derivative of the volume

  omega[0] = omega[1] = omega[2] = 0.0;

  nrigid = 0;

  rfix = NULL;

  old_velocity = new double* [atom->nlocal];

  for ( int j = 0; j < atom->nlocal; j++ ) {

    old_velocity[j] = new double [3];

  }

  atoms_allocated = atom->nlocal;

  maxold = 0;

  old_velocity = NULL;

}

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

  delete [] id_press;

  delete [] id_pe;

  for ( int j = 0; j < atoms_allocated; j++ ) {

    delete [] old_velocity[j];

  }

  delete [] old_velocity;

  memory->destroy(old_velocity);

}

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

          strcmp(modify->fix[i]->style,"poems") == 0) rfix[nrigid++] = i;

  }

  // find fix external being used to drive LAMMPS from DFTB+

  if (dftb) {

    for (int i = 0; i < modify->nfix; i++)

      if (strcmp(modify->fix[i]->style,"external") == 0)

        fix_external = (FixExternal *) modify->fix[i];

    if (fix_external == NULL) 

      error->all(FLERR,"Fix msst dftb cannot be used w/out fix external");

  }

}

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

void FixMSST::initial_integrate(int vflag)

{

  int sd;

  double p_msst;                // MSST driving pressure.

  int i,k;

  double vol;

  double p_msst;                       // MSST driving pressure

  double vol,TS,TS_term,escale_term;

  int nlocal = atom->nlocal;

  int *mask = atom->mask;

  double **v = atom->v;

  double *mass = atom->mass;

  int *type = atom->type;

  double **x = atom->x;

  int sd = direction;

  // check to see if old_velocity is correctly allocated

  // realloc old_velocity if necessary

  check_alloc(nlocal);

  if (nlocal > maxold) {

    memory->destroy(old_velocity);

    maxold = atom->nmax;

    memory->create(old_velocity,maxold,3,"msst:old_velocity");

  }

  sd = direction;

  // for DFTB, extract TS_dftb from fix external

  // must convert energy to mv^2 units

  if (dftb) {

    TS_dftb = fix_external->compute_vector(1);

    TS = force->ftm2v*TS_dftb;

    if (update->ntimestep == 1) T0S0 = TS;

  } else {

    TS = 0.0;

    T0S0 = 0.0;

  }

  // compute new pressure and volume

  // compute new pressure and volume.

  temperature->compute_vector();

  pressure->compute_vector();

  couple();

  vol = compute_vol();

  // update S_elec terms and compute TS_dot via finite differences

  S_elec_2 = S_elec_1;

  S_elec_1 = S_elec;

  double Temp = temperature->compute_scalar();

  S_elec = TS/Temp;

  TS_dot = Temp*(3.0*S_elec-4.0*S_elec_1+S_elec_2)/(2.0*update->dt);

  TS_int += (update->dt*TS_dot);

  //TS_int += (update->dt*TS_dot)/total_mass;

  // compute etot + extra terms for conserved quantity 

  double e_scale = compute_etotal() + compute_scalar();

  // propagate the time derivative of

  // the volume 1/2 step at fixed vol, r, rdot.

  // the volume 1/2 step at fixed vol, r, rdot

  p_msst = nktv2p * mvv2e * velocity * velocity * total_mass *

    ( v0 - vol)/( v0 * v0);

    (qmass * nktv2p * mvv2e);

  double B = total_mass * mu / ( qmass * vol );

  // prevent blow-up of the volume.

  // prevent blow-up of the volume

  if ( vol > v0 && A > 0.0 ) {

    A = -A;

  }

  if (vol > v0 && A > 0.0) A = -A;

  // use taylor expansion to avoid singularity at B == 0.

  // use Taylor expansion to avoid singularity at B = 0

  if ( B * dthalf > 1.0e-06 ) {

    omega[sd] = ( omega[sd] + A * ( exp(B * dthalf) - 1.0 ) / B )

  }

  // propagate velocity sum 1/2 step by

  // temporarily propagating the velocities.

  // temporarily propagating the velocities

  velocity_sum = compute_vsum();

  if (dftb) {

    for (i = 0; i < nlocal; i++) {

      if (mask[i] & groupbit) {

        for ( k = 0; k < 3; k++ ) {

          double C = f[i][k] * force->ftm2v / mass[type[i]];

          TS_term = TS_dot/(mass[type[i]]*velocity_sum);

          escale_term = force->ftm2v*beta*(e0-e_scale) / 

            (mass[type[i]]*velocity_sum);

          double D = mu * omega[sd] * omega[sd] /

            (velocity_sum * mass[type[i]] * vol );

          D += escale_term - TS_term;

          old_velocity[i][k] = v[i][k];

          if ( k == direction ) D -= 2.0 * omega[sd] / vol;

          if ( fabs(dthalf * D) > 1.0e-06 ) {

            double expd = exp(D * dthalf);

            v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;

          } else {

            v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +

              0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf;

          }

        }

      }

    }

  } else {

    for (i = 0; i < nlocal; i++) {

      if (mask[i] & groupbit) {

        for ( k = 0; k < 3; k++ ) {

        }

      }

    }

  }

  velocity_sum = compute_vsum();

  // reset the velocities.

  // reset the velocities

  for (i = 0; i < nlocal; i++) {

    if (mask[i] & groupbit) {

      for ( k = 0; k < 3; k++ ) {

        v[i][k] = old_velocity[i][k];

      }

      v[i][0] = old_velocity[i][0];

      v[i][1] = old_velocity[i][1];

      v[i][2] = old_velocity[i][2];

    }

  }

  // propagate velocities 1/2 step using the new velocity sum.

  // propagate velocities 1/2 step using the new velocity sum

  if (dftb) {

    for (i = 0; i < nlocal; i++) {

      if (mask[i] & groupbit) {

        for ( k = 0; k < 3; k++ ) {

          double C = f[i][k] * force->ftm2v / mass[type[i]];

          TS_term = TS_dot/(mass[type[i]]*velocity_sum);

          escale_term = force->ftm2v*beta*(e0-e_scale) / 

            (mass[type[i]]*velocity_sum);

          double D = mu * omega[sd] * omega[sd] /

            (velocity_sum * mass[type[i]] * vol );

          D += escale_term - TS_term;

          if ( k == direction ) D -= 2.0 * omega[sd] / vol;

          if ( fabs(dthalf * D) > 1.0e-06 ) {

            double expd = exp(D * dthalf);

            v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;

          } else {

            v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +

              0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf;

          }

        }

      }

    }

  } else {

    for (i = 0; i < nlocal; i++) {

      if (mask[i] & groupbit) {

        for ( k = 0; k < 3; k++ ) {

        }

      }

    }

  }

  // propagate the volume 1/2 step.

  // propagate the volume 1/2 step

  double vol1 = vol + omega[sd] * dthalf;

  // rescale positions and change box size.

  // rescale positions and change box size

  dilation[sd] = vol1/vol;

  remap(0);

  // propagate particle positions 1 time step.

  // propagate particle positions 1 time step

  for (i = 0; i < nlocal; i++) {

    if (mask[i] & groupbit) {

    }

  }

  // propagate the volume 1/2 step.

  // propagate the volume 1/2 step

  double vol2 = vol1 + omega[sd] * dthalf;

  // rescale positions and change box size.

  // rescale positions and change box size

  dilation[sd] = vol2/vol1;

  remap(0);

void FixMSST::final_integrate()

{

  int i;

  double p_msst;                  // MSST driving pressure

  double TS,TS_term,escale_term;

  // v update only for atoms in MSST group

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  double vol = compute_vol();

  double p_msst;

  int sd = direction;

  // propagate particle velocities 1/2 step.

  // for DFTB, extract TS_dftb from fix external

  // must convert energy to mv^2 units

  if (dftb) {

    TS_dftb = fix_external->compute_vector(1);

    TS = force->ftm2v*TS_dftb;

  } else TS = 0.0;

  // compute etot + extra terms for conserved quantity 

  double e_scale = compute_etotal() + compute_scalar();

  // propagate particle velocities 1/2 step

  if (dftb) {

    for (i = 0; i < nlocal; i++) {

      if (mask[i] & groupbit) {

        for ( int k = 0; k < 3; k++ ) {

          double C = f[i][k] * force->ftm2v / mass[type[i]];

          TS_term = TS_dot/(mass[type[i]]*velocity_sum);

          escale_term = force->ftm2v*beta*(e0-e_scale) / 

            (mass[type[i]]*velocity_sum);

          double D = mu * omega[sd] * omega[sd] /

            (velocity_sum * mass[type[i]] * vol );

          D += escale_term - TS_term;

          if ( k == direction ) D -= 2.0 * omega[sd] / vol;

          if ( fabs(dthalf * D) > 1.0e-06 ) {

            double expd = exp(D * dthalf);

            v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;

          } else {

            v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +

              0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf;

          }

        }

      }

    }

  } else {

    for (i = 0; i < nlocal; i++) {

      if (mask[i] & groupbit) {

        for ( int k = 0; k < 3; k++ ) {

        }

      }

    }

  }

  // compute new pressure and volume.

  // compute new pressure and volume