created vashishta/table variant

pair_style vashishta command :h3

pair_style vashishta/omp command :h3

pair_style vashishta/table command :h3

[Syntax:]

pair_style vashishta :pre

pair_style style args :pre

style = {vashishta} or {vashishta/table}

args = list of arguments for a particular style :ul

  {vashishta} args = none

  {vashishta/table} args = Ntable cutinner

    Ntable = # of tabulation points

    cutinner = tablulate from cutinner to cutoff :pre

[Examples:]

pair_style vashishta

pair_coeff * * SiC.vashishta Si C :pre

pair_style vashishta/table 100000 0.2

pair_coeff * * SiC.vashishta Si C :pre

[Description:]

The {vashishta} style computes the combined 2-body and 3-body 

family of potentials developed in the group of Vashishta and

co-workers. By combining repulsive, screened Coulombic, 

screened charge-dipole, and dispersion interactions with a 

bond-angle energy based on the Stillinger-Weber potential, 

this potential has been used to describe a variety of inorganic

compounds, including SiO2 "Vashishta1990"_#Vashishta1990, 

SiC "Vashishta2007"_#Vashishta2007,

The {vashishta} and {vashishta/table} styles compute the combined

2-body and 3-body family of potentials developed in the group of

Vashishta and co-workers. By combining repulsive, screened Coulombic,

screened charge-dipole, and dispersion interactions with a bond-angle

energy based on the Stillinger-Weber potential, this potential has

been used to describe a variety of inorganic compounds, including SiO2

"Vashishta1990"_#Vashishta1990, SiC "Vashishta2007"_#Vashishta2007,

and InP "Branicio2009"_#Branicio2009.

The potential for the energy U of a system of atoms is

is over all neighbors J and K within a cut-off distance = {r0},

where the exponential screening function becomes zero.

Only a single pair_coeff command is used with the {vashishta} style which

specifies a Vashishta potential file with parameters for all

needed elements.  These are mapped to LAMMPS atom types by specifying

N additional arguments after the filename in the pair_coeff command,

The {vashishta} style computes these formulas analytically.  The

{vashishta/table} style tabulates the analytic values for {Ntable}

points from cutinner to the cutoff of the potential.  The points are

equally spaced in R^2 space from cutinner^2 to cutoff^2.  For the

two-body term in the above equation, a linear interpolation for each

pairwise distance between adjacent points in the table.  In practice

the tabulated version can run 3-5x faster than the analytic version

with little loss of accuracy for Ntable = ???.

Only a single pair_coeff command is used with either style which

specifies a Vashishta potential file with parameters for all needed

elements.  These are mapped to LAMMPS atom types by specifying N

additional arguments after the filename in the pair_coeff command,

where N is the number of LAMMPS atom types:

filename

C

costheta0 :ul

The non-annotated parameters are unitless.

The Vashishta potential file must contain entries for all the

elements listed in the pair_coeff command.  It can also contain

entries for additional elements not being used in a particular

simulation; LAMMPS ignores those entries.

For a single-element simulation, only a single entry is required

(e.g. SiSiSi).  For a two-element simulation, the file must contain 8

entries (for SiSiSi, SiSiC, SiCSi, SiCC, CSiSi, CSiC, CCSi, CCC), that

specify parameters for all permutations of the two elements

interacting in three-body configurations.  Thus for 3 elements, 27

entries would be required, etc.

Depending on the particular version of the Vashishta potential,

the values of these parameters may be keyed to the identities of

zero, one, two, or three elements.

In order to make the input file format unambiguous, general,

and simple to code,

LAMMPS uses a slightly confusing method for specifying parameters.

All parameters are divided into two classes: two-body and three-body.

Two-body and three-body parameters are handled differently,

as described below.

The two-body parameters are H, eta, lambda1, D, lambda4, W, rc, gamma, and r0.

They appear in the above formulae with two subscripts.

The parameters Zi and Zj are also classified as two-body parameters,

even though they only have 1 subscript.

The three-body parameters are B, C, costheta0.

They appear in the above formulae with three subscripts.

Two-body and three-body parameters are handled differently,

as described below.

The first element in each entry is the center atom

in a three-body interaction, while the second and third elements

are two neighbor atoms. Three-body parameters for a central atom I

and two neighbors J and K are taken from the IJK entry.

Note that even though three-body parameters do not depend on the order of

J and K, LAMMPS stores three-body parameters for both IJK and IKJ.

The user must ensure that these values are equal.     

Two-body parameters for an atom I interacting with atom J are taken from

the IJJ entry, where the 2nd and 3rd

elements are the same. Thus the two-body parameters 

for Si interacting with C come from the SiCC entry. Note that even

though two-body parameters (except possibly gamma and r0 in U3) 

do not depend on the order of the two elements, 

LAMMPS will get the Si-C value from the SiCC entry

and the C-Si value from the CSiSi entry. The user must ensure

that these values are equal. Two-body parameters appearing

in entries where the 2nd and 3rd elements are different are

stored but never used. It is good practice to enter zero for

these values. Note that the three-body function U3 above 

contains the two-body parameters gamma and r0. So U3 for a 

central C atom bonded to an Si atom and a second C atom

will take three-body parameters from the CSiC entry, but

The non-annotated parameters are unitless.  The Vashishta potential

file must contain entries for all the elements listed in the

pair_coeff command.  It can also contain entries for additional

elements not being used in a particular simulation; LAMMPS ignores

those entries.  For a single-element simulation, only a single entry

is required (e.g. SiSiSi).  For a two-element simulation, the file

must contain 8 entries (for SiSiSi, SiSiC, SiCSi, SiCC, CSiSi, CSiC,

CCSi, CCC), that specify parameters for all permutations of the two

elements interacting in three-body configurations.  Thus for 3

elements, 27 entries would be required, etc.

Depending on the particular version of the Vashishta potential, the

values of these parameters may be keyed to the identities of zero,

one, two, or three elements.  In order to make the input file format

unambiguous, general, and simple to code, LAMMPS uses a slightly

confusing method for specifying parameters.  All parameters are

divided into two classes: two-body and three-body.  Two-body and

three-body parameters are handled differently, as described below.

The two-body parameters are H, eta, lambda1, D, lambda4, W, rc, gamma,

and r0.  They appear in the above formulae with two subscripts.  The

parameters Zi and Zj are also classified as two-body parameters, even

though they only have 1 subscript.  The three-body parameters are B,

C, costheta0.  They appear in the above formulae with three

subscripts.  Two-body and three-body parameters are handled

differently, as described below.

The first element in each entry is the center atom in a three-body

interaction, while the second and third elements are two neighbor

atoms. Three-body parameters for a central atom I and two neighbors J

and K are taken from the IJK entry.  Note that even though three-body

parameters do not depend on the order of J and K, LAMMPS stores

three-body parameters for both IJK and IKJ.  The user must ensure that

these values are equal.  Two-body parameters for an atom I interacting

with atom J are taken from the IJJ entry, where the 2nd and 3rd

elements are the same. Thus the two-body parameters for Si interacting

with C come from the SiCC entry. Note that even though two-body

parameters (except possibly gamma and r0 in U3) do not depend on the

order of the two elements, LAMMPS will get the Si-C value from the

SiCC entry and the C-Si value from the CSiSi entry. The user must

ensure that these values are equal. Two-body parameters appearing in

entries where the 2nd and 3rd elements are different are stored but

never used. It is good practice to enter zero for these values. Note

that the three-body function U3 above contains the two-body parameters

gamma and r0. So U3 for a central C atom bonded to an Si atom and a

second C atom will take three-body parameters from the CSiC entry, but

two-body parameters from the CCC and CSiSi entries.

:line

described above from values in the potential file.

This pair style does not support the "pair_modify"_pair_modify.html

{shift} and {tail} options, but it supports the {table} option.

Turning on the {table} option will use a linear interpolation table

for force and energy of the two-body term with {2**N} entries.

Tabulation can reduce the computational cost by a factor of 2-3

with 20-bit to 12-bit table sizes, respectively, at a small to

moderate reduction in precision. It is not recommended to use

smaller than 12-bit tables.

shift, table, and tail options.

This pair style does not write its information to "binary restart

files"_restart.html, since it is stored in potential files.  Thus, you

This pair style requires the "newton"_newton.html setting to be "on"

for pair interactions.

The Vashishta potential files provided with LAMMPS (see the

potentials directory) are parameterized for metal "units"_units.html.

You can use the Vashishta potential with any LAMMPS units, but you would need

to create your own Vashishta potential file with coefficients listed in the

appropriate units if your simulation doesn't use "metal" units.

The Vashishta potential files provided with LAMMPS (see the potentials

directory) are parameterized for metal "units"_units.html.  You can

use the Vashishta potential with any LAMMPS units, but you would need

to create your own Vashishta potential file with coefficients listed

in the appropriate units if your simulation doesn't use "metal" units.

[Related commands:]

:line

:link(Vashishta1990)

[(Vashishta1990)] P. Vashishta, R. K. Kalia, J. P. Rino, Phys. Rev. B 41, 12197 (1990).

[(Vashishta1990)] P. Vashishta, R. K. Kalia, J. P. Rino, Phys. Rev. B

41, 12197 (1990).

:link(Vashishta2007)

[(Vashishta2007)] P. Vashishta, R. K. Kalia, A. Nakano, J. P. Rino. J. Appl. Phys. 101, 103515 (2007).

[(Vashishta2007)] P. Vashishta, R. K. Kalia, A. Nakano,

J. P. Rino. J. Appl. Phys. 101, 103515 (2007).

:link(Branicio2009)

[(Branicio2009)] Branicio, Rino, Gan and Tsuzuki, J. Phys Condensed Matter 21 (2009) 095002

[(Branicio2009)] Branicio, Rino, Gan and Tsuzuki, J. Phys Condensed

Matter 21 (2009) 095002

# choose potential

pair_style	vashishta

pair_style 	vashishta/table 100000 0.2

pair_coeff 	* * InP.vashishta In P

pair_modify     table 16

# setup neighbor style

velocity	all create 2000.0 277387 mom yes

displace_atoms	all move 0.05 0.9 0.4 units box

pair_style 	vashishta

pair_style 	vashishta/table 100000 0.2

pair_coeff	* *  SiO.1990.vashishta Si O

pair_modify	table 16

neighbor	0.3 bin

neigh_modify	delay 10

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

   Contributing author:  Yongnan Xiong (HNU), xyn@hnu.edu.cn

                         Aidan Thompson (SNL)

                         Anders Hafreager (UiO), andershaf@gmail.com

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

#include <math.h>

#include "neigh_list.h"

#include "memory.h"

#include "error.h"

using namespace LAMMPS_NS;

#define MAXLINE 1024

  params = NULL;

  elem2param = NULL;

  map = NULL;

  neigh3BodyMax = 0;

  neigh3BodyCount = NULL; 

  neigh3Body = NULL;

  useTable = false;

  nTablebits = 12;

  forceTable = NULL;

  potentialTable = NULL;

}

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

  memory->destroy(params);

  memory->destroy(elem2param);

  memory->destroy(forceTable);

  memory->destroy(potentialTable);

  memory->destroy(neigh3BodyCount);

  memory->destroy(neigh3Body);

  if (allocated) {

    memory->destroy(setflag);

    memory->destroy(cutsq);

  }

}

void PairVashishta::modify_params(int narg, char **arg)

{

  if (narg == 0) error->all(FLERR,"Illegal pair_modify command");

  int iarg = 0;

  while (iarg < narg) {

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

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

      nTablebits = force->inumeric(FLERR,arg[iarg+1]);

      if ((nTablebits < 0) || (nTablebits > sizeof(int)*CHAR_BIT-1))

        error->all(FLERR,"Illegal interpolation table size");

      if (nTablebits == 0) {

        useTable = false;

      } else {

        useTable = true;

      }

      iarg += 2;

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

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

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

      if (tabinner <= 0.0)

        error->all(FLERR,"Illegal inner cutoff for tabulation");

      iarg += 2;

    }

  }

  updateTables();

}

void PairVashishta::validateNeigh3Body() {

  if (atom->nlocal > neigh3BodyMax) {

        neigh3BodyMax = atom->nmax;

        memory->destroy(neigh3BodyCount);

        memory->destroy(neigh3Body);

        memory->create(neigh3BodyCount,neigh3BodyMax,  "pair:vashishta:neigh3BodyCount");

        memory->create(neigh3Body,neigh3BodyMax, 1000, "pair:vashishta:neigh3Body"); // TODO: pick this number more wisely? 

    }

}

void PairVashishta::updateTables()

{

  memory->destroy(forceTable);

  memory->destroy(potentialTable);

  forceTable = NULL;

  potentialTable = NULL;

  if (!useTable) return;

  int ntable = 1<<nTablebits;

  tabinnersq = tabinner*tabinner;

  deltaR2 = (cutmax*cutmax - tabinnersq) / (ntable-1);

  oneOverDeltaR2 = 1.0/deltaR2;

  memory->create(forceTable,nelements,nelements,ntable+1,"pair:vashishta:forceTable");

  memory->create(potentialTable,nelements,nelements,ntable+1,"pair:vashishta:potentialTable");

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

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

      int ijparam = elem2param[i][j][j];

      for(int idx=0; idx <= ntable; idx++) {

        double rsq = tabinnersq + idx*deltaR2;

        double fpair, eng;

        // call analytical version of two-body term function

        twobody(params[ijparam], rsq, fpair, 1, eng, false);

        forceTable[i][j][idx] = fpair;

        potentialTable[i][j][idx] = eng;

      }

    }

  }

}

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

void PairVashishta::compute(int eflag, int vflag)

{

  validateNeigh3Body();

  int i,j,k,ii,jj,kk,inum,jnum,jnumm1;

  int itype,jtype,ktype,ijparam,ikparam,ijkparam;

  tagint itag,jtag;

    ytmp = x[i][1];

    ztmp = x[i][2];

    neigh3BodyCount[i] = 0; // Reset the 3-body neighbor list

    // two-body interactions, skip half of them

    jlist = firstneigh[i];

      j &= NEIGHMASK;

      jtag = tag[j];

      jtype = map[type[j]];

      delx = xtmp - x[j][0];

      dely = ytmp - x[j][1];

      delz = ztmp - x[j][2];

      rsq = delx*delx + dely*dely + delz*delz;

      ijparam = elem2param[itype][jtype][jtype];

      if (rsq <= params[ijparam].cutsq2) {

          neigh3Body[i][neigh3BodyCount[i]] = j;

          neigh3BodyCount[i]++;

      }

      if (rsq > params[ijparam].cutsq) continue;

      if (itag > jtag) {

        if ((itag+jtag) % 2 == 0) continue;

      } else if (itag < jtag) {

        if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;

      }

      twobody(params[ijparam],rsq,fpair,eflag,evdwl, useTable);

      jtype = map[type[j]];

      delx = xtmp - x[j][0];

      dely = ytmp - x[j][1];

      delz = ztmp - x[j][2];

      rsq = delx*delx + dely*dely + delz*delz;

      ijparam = elem2param[itype][jtype][jtype];

      if (rsq > params[ijparam].cutsq) continue;

      twobody(&params[ijparam],rsq,fpair,eflag,evdwl);

      f[i][0] += delx*fpair;

      f[i][1] += dely*fpair;

      			   evdwl,0.0,fpair,delx,dely,delz);

    }

    jlist = neigh3Body[i];

    jnum = neigh3BodyCount[i];

    jnumm1 = jnum - 1;

    for (jj = 0; jj < jnumm1; jj++) {

    if (params[m].cut > cutmax) cutmax = params[m].cut;

    if (params[m].r0 > cutmax) cutmax = params[m].r0;

  }

  updateTables();

}

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

void PairVashishta::twobody(const Param &param, double rsq, double &fforce,

                     int eflag, double &eng, bool tabulated)

void PairVashishta::twobody(Param *param, double rsq, double &fforce,

                            int eflag, double &eng)

{

  if (tabulated) {

    if (rsq < tabinnersq) {

      sprintf(estr,"Pair distance < table inner cutoff: " 

              "ijtype %d %d dist %g",param.ielement+1,param.jelement+1,sqrt(rsq));

      error->one(FLERR,estr);

    }

    const int tableIndex = (rsq - tabinnersq)*oneOverDeltaR2;

    // double - int will only keep the 0.xxxx part

    const double fraction = (rsq - tabinnersq)*oneOverDeltaR2 - tableIndex;

    double force0 = forceTable[param.ielement][param.jelement][tableIndex];

    double force1 = forceTable[param.ielement][param.jelement][tableIndex+1];

    fforce = (1.0 - fraction)*force0 + fraction*force1; // force is linearly interpolated between the two values

    if (evflag) {

      double energy0 = potentialTable[param.ielement][param.jelement][tableIndex];

      double energy1 = potentialTable[param.ielement][param.jelement][tableIndex+1];

      eng = (1.0 - fraction)*energy0 + fraction*energy1;

    }

  } else {

  double r,rinvsq,r4inv,r6inv,reta,lam1r,lam4r,vc2,vc3;

  r = sqrt(rsq);

  rinvsq = 1.0/rsq;

  r4inv = rinvsq*rinvsq;

  r6inv = rinvsq*r4inv;

    reta = pow(r,-param.eta);

    lam1r = r*param.lam1inv;

    lam4r = r*param.lam4inv;

    vc2 = param.zizj * exp(-lam1r)/r;

    vc3 = param.mbigd * r4inv*exp(-lam4r);

    fforce = (param.dvrc*r

  	    - (4.0*vc3 + lam4r*vc3+param.big6w*r6inv

  	       - param.heta*reta - vc2 - lam1r*vc2)

  reta = pow(r,-param->eta);

  lam1r = r*param->lam1inv;

  lam4r = r*param->lam4inv;

  vc2 = param->zizj * exp(-lam1r)/r;

  vc3 = param->mbigd * r4inv*exp(-lam4r);

  fforce = (param->dvrc*r

	    - (4.0*vc3 + lam4r*vc3+param->big6w*r6inv

	       - param->heta*reta - vc2 - lam1r*vc2)

	    ) * rinvsq;

    if (eflag) eng = param.bigh*reta

  	       + vc2 - vc3 - param.bigw*r6inv

  	       - r*param.dvrc + param.c0;

  }

  if (eflag) eng = param->bigh*reta

	       + vc2 - vc3 - param->bigw*r6inv

	       - r*param->dvrc + param->c0;

}

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

  if (eflag) eng = facrad;

}

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

 *  add memory usage for tabulation

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

double PairVashishta::memory_usage()

{

  double bytes = Pair::memory_usage();

  if (useTable)

    bytes += 2*nelements*nelements*sizeof(double)*(1<<nTablebits);

  return bytes;

}

#else

#ifndef LMP_PAIR_Vashishta_H

#define LMP_PAIR_Vashishta_H

#ifndef LMP_PAIR_VASHISHITA_H

#define LMP_PAIR_VASHISHITA_H

#include "pair.h"

 public:

  PairVashishta(class LAMMPS *);

  virtual ~PairVashishta();

  virtual void modify_params(int, char **);

  virtual void compute(int, int);

  void settings(int, char **);

  virtual void settings(int, char **);

  void coeff(int, char **);

  virtual double init_one(int, int);

  virtual void init_style();

  virtual double memory_usage();

  double init_one(int, int);

  void init_style();

 protected:

  struct Param {

    int ielement,jelement,kelement;

  };

  bool useTable;

  int nTablebits;

  double deltaR2;

  double oneOverDeltaR2;

  double ***forceTable;          // Table containing the forces

  double ***potentialTable;      // Table containing the potential energies

  int neigh3BodyMax;            // Max size of short neighborlist

  int *neigh3BodyCount;         // Number of neighbors in short range 3 particle forces neighbor list

  int **neigh3Body;             // Neighborlist for short range 3 particle forces

  double cutmax;                // max cutoff for all elements

  int nelements;                // # of unique elements

  char **elements;              // names of unique elements

  int ***elem2param;            // mapping from element triplets to parameters

  int *map;                     // mapping from atom types to elements

  char estr[128];               // Char array to store error message with sprintf (i.e. twobody table)

  int nparams;                  // # of stored parameter sets

  int maxparam;                 // max # of parameter sets

  Param *params;                // parameter set for an I-J-K interaction

  virtual void allocate();

  void allocate();

  void read_file(char *);

  void setup_params();

  void updateTables();

  void validateNeigh3Body();

  void twobody(const Param &, double, double &, int, double &, bool tabulated=false);

  virtual void setup_params();

  void twobody(Param *, double, double &, int, double &);

  void threebody(Param *, Param *, Param *, double, double, double *, double *,

                 double *, double *, int, double &);

};