Finished fix nphug

   See the README file in the top-level LAMMPS directory.

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

/*

#include "string.h"

#include "stdlib.h"

#include "fix_nphug.h"

#include "modify.h"

#include "error.h"

#include "update.h"

#include "compute.h"

#include "force.h"

#include "domain.h"

#include "group.h"

#include "math.h"

#include "memory.h"

#include "comm.h"

#include "math.h"

  This fix applies the NPHug Hugoniostat method of Ravelo et al.

using namespace LAMMPS_NS;

(Ravelo, Holian, Germann, and Lomdahl, PRB 70 014103 (2004))

enum{ISO,ANISO,TRICLINIC}; // same as fix_nh.cpp

It uses the Nose-Hoover thermostat and barostat (fix_nh.html).

The Nose-Hoover barostat is used to compress the system 

to a specified final stress state. This is done either

hydrostatically (using keyword iso, aniso, or tri) or uniaxially

(using keywords x, y, or z).  In the hydrostatic case,

the cell dimensions change dynamically so that the average axial stress

in all three directions converges towards the specified target value. 

In the uniaxial case, the chosen cell dimension changes dynamically 

so that the average

axial stress in that direction converges towards the target value. The

other two cell dimensions are kept fixed (zero lateral strain).

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

This leads to the following restrictions on the keywords:

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

  FixNH(lmp, narg, arg)

{

- The specified initial and

final target pressures must be the same.

  // extend vector of base-class computes

- Only one of the following keywords may be used:

iso, aniso, tri, x, y, z, 

  size_vector += 3;

- The keywords xy, xz, yz may not be used.

  // turn off deviatoric flag and remove strain energy from vector

- The only admissible value for the couple keyword is xyz,

which has the same effect as keyword iso

  deviatoric_flag = 0;

  size_vector -= 1;

- The drag parameter is proportional to the beta_H and beta_p

damping coefficients in the Ravelo paper.

  // use initial state as reference state

- The temp keyword serves only to set the value of tdamp. The initial

and final target temperatures are ignored. 

  v0_set = p0_set = e0_set = 0;

- The values of tdamp and pdamp are inversely proportional to the

coupling rate nu_H and nu_p in the Ravelo paper

  // check pressure settings

- All other keywords function in the same way. 

  if (p_start[0] != p_stop[0] ||

      p_start[1] != p_stop[1] ||  

      p_start[2] != p_stop[2])

    error->all("Invalid argument for fix nphug");

*/

  // uniaxial = 0 means hydrostatic compression

  // uniaxial = 1 means uniaxial compression

  //      in x, y, or z (idir = 0, 1, or 2)

  // isotropic hydrostatic compression

  if (pstyle == ISO) {

    uniaxial = 0;

    // anisotropic compression

  } else if (pstyle == ANISO) {

#include "string.h"

#include "fix_nphug.h"

#include "modify.h"

#include "error.h"

#include "update.h"

#include "compute.h"

#include "force.h"

#include "domain.h"

    // anisotropic hydrostatic compression

using namespace LAMMPS_NS;

    if (p_start[0] == p_start[1] &&  

	p_start[0] == p_start[2] )

      uniaxial = 0;

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

    // uniaxial compression

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

  FixNH(lmp, narg, arg)

{

  // Hard-code to use initial state of system

    else if (p_flag[0] == 1 && p_flag[1] == 0 

	&& p_flag[2] == 0) {

      uniaxial = 1;

      idir = 0;

    } else if (p_flag[0] == 0 && p_flag[1] == 1 

	   && p_flag[2] == 0) {

      uniaxial = 1;

      idir = 1;

    } else if (p_flag[0] == 0 && p_flag[1] == 0 

	       && p_flag[2] == 1) {

      uniaxial = 1;

      idir = 2;

  v0_set = p0_set = e0_set = 0;

    } else error->all("Invalid argument for fix nphug");

  // Hard-code to use z direction

    // triclinic hydrostatic compression

  direction = 2;

  if (p_flag[0] == 1 || p_flag[1] == 1 ||

      p_flag[3] == 1 || p_flag[4] == 1 || p_flag[5] == 1)

    error->all("Only pressure control in z direction to be used with fix nphug");

  if (p_flag[2] == 0)

    error->all("Pressure control in z direction must be used with fix nphug");

  } else if (pstyle == TRICLINIC) {

    if (p_start[0] == p_start[1] &&  

	p_start[0] == p_start[2] &&

	p_start[3] == 0.0 &&  

	p_start[4] == 0.0 &&  

	p_start[5] == 0.0 )

      uniaxial = 0;

    else error->all("Invalid argument for fix nphug");

  }

  if (!tstat_flag)

    error->all("Temperature control must be used with fix nphug");

  } 

  if ( p0_set == 0 ) {

    p0 = p_current[direction];

    p0_set = 1;

    if (uniaxial == 1)

      p0 = p_current[idir];

    else

      p0 = (p_current[0]+p_current[1]+p_current[2])/3.0;

  }

  if ( e0_set == 0 ) {

    e0_set = 1;

  }

  double masstot = group->mass(igroup);

  rho0 = nktv2p*force->mvv2e*masstot/v0;

  t_target = 0.01;

  pe->addstep(update->ntimestep+1);

}

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

{

  t_target = t_current + compute_hugoniot();

  ke_target = tdof * boltz * t_target;

  if (ke_target < 0.0) ke_target = 0.0;

  // If t_target is very small, need to choose 

  // more reasonable value for use by barostat and 

  // thermostat masses. ke_target is left as is.

  if (t_target <= 1.0e-6) {

    if (strcmp(update->unit_style,"lj") == 0) t0 = 1.0;

    else t0 = 300.0;

  }

  pressure->addstep(update->ntimestep+1);

  pe->addstep(update->ntimestep+1);

}

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

double FixNPHug::compute_etotal()

  epot = pe->compute_scalar();

  if (thermo_energy) epot -= compute_scalar();

  ekin = temperature->compute_scalar();

  ekin *= 0.5 * temperature->dof * force->boltz;

  ekin *= 0.5 * tdof * force->boltz;

  etot = epot+ekin;

  return etot;

}

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

   Computes the deviation of the current point 

   from the Hugoniot in energy units.

   from the Hugoniot in temperature units.

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

double FixNPHug::compute_hugoniot()

  e = compute_etotal();

  temperature->compute_vector();

  if (uniaxial == 1) {

    pressure->compute_vector();

  p = pressure->vector[direction];

    p = pressure->vector[idir];

  } else 

    p = pressure->compute_scalar();

  v = compute_vol();

  dhugo = (0.5 * (p + p0 ) * ( v0 - v)) / 

    force->nktv2p + e0 - e;

  dhugo /= tdof * boltz;

  return dhugo;

}

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

   Compute shock velocity is distance/time units

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

double FixNPHug::compute_us()

{

  double v,p;

  double eps,us;

  temperature->compute_vector();

  if (uniaxial == 1) {

    pressure->compute_vector();

    p = pressure->vector[idir];

  } else 

    p = pressure->compute_scalar();

  v = compute_vol();

  // Us^2 = (p-p0)/(rho0*eps)

  eps = 1.0 - v/v0;

  if (eps < 1.0e-10) us = 0.0;

  else if (p < p0) us = 0.0; 

  else us = sqrt((p-p0)/(rho0*eps));

  return us;

}

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

   Compute particle velocity is distance/time units

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

double FixNPHug::compute_up()

{

  double v;

  double eps,us,up;

  v = compute_vol();

  us = compute_us();

  // u = eps*Us

  eps = 1.0 - v/v0;

  up = us*eps;

  return up;

}

// look for index in local class

// if index not found, look in base class

double FixNPHug::compute_vector(int n)

{

  int ilen;

  // n = 0: Hugoniot energy difference (temperature units)

  ilen = 1;

  if (n < ilen) return compute_hugoniot();

  n -= ilen;

  // n = 1: Shock velocity

  ilen = 1;

  if (n < ilen) return compute_us();

  n -= ilen;

  // n = 2: Particle velocity

  ilen = 1;

  if (n < ilen) return compute_up();

  n -= ilen;

  // index not found, look in base class

  FixNH::compute_vector(n);

}

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

   pack entire state of Fix into one write 

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

void FixNPHug::write_restart(FILE *fp)

{

  int nsize = 3;

  // Add on the base class size

  int nsizetot = nsize + FixNH::size_restart();

  double *list;

  memory->create(list,nsize,"nphug:list");

  int n = 0;

  list[n++] = e0;

  list[n++] = v0;

  list[n++] = p0;

  if (comm->me == 0) {

    int sizetot = nsizetot * sizeof(double);

    fwrite(&sizetot,sizeof(int),1,fp);

    fwrite(list,sizeof(double),nsize,fp);

  }

  memory->destroy(list);

  // Write the base class data without size

  FixNH::write_restart_nosize(fp);

}

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

   use state info from restart file to restart the Fix 

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

void FixNPHug::restart(char *buf)

{

  int n = 0;

  double *list = (double *) buf;

  e0 = list[n++];

  v0 = list[n++];

  p0 = list[n++];

  e0_set = 1;

  v0_set = 1;

  p0_set = 1;

  // Call the base class function

  buf += n*sizeof(double);

  FixNH::restart(buf);

}

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

int FixNPHug::modify_param(int narg, char **arg)

{

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

    if (narg < 2) error->all("Illegal fix nphug command");

    e0 = atof(arg[1]);

    e0_set = 1;

    return 2;

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

    if (narg < 2) error->all("Illegal fix nphug command");

    v0 = atof(arg[1]);

    v0_set = 1;

    return 2;

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

    if (narg < 2) error->all("Illegal fix nphug command");

    p0 = atof(arg[1]);

    p0_set = 1;

    return 2;

  }

  return 0;

}

  ~FixNPHug();

  void init();

  void setup(int);

  int modify_param(int, char **);

  void write_restart(FILE *);

  void restart(char *);

 private:

  class Compute *pe;               // PE compute pointer

  void compute_temp_target();

  double compute_vector(int);

  double compute_etotal();

  double compute_vol();

  double compute_hugoniot();

  double compute_us();

  double compute_up();

  char *id_pe;

  int peflag;

  int v0_set,p0_set,e0_set;

  double v0,p0,e0;

  int direction;

  double v0,p0,e0,rho0;

  int idir;

  int uniaxial;

};

}