add gao-weber potentials (regular and with ZBL core) with SiC potential files

# DATE: 2016-05-06 CONTRIBUTOR: German Samolyuk, samolyuk@gmail.com CITATION: ???

# Gao-Weber parameters for various elements and mixtures

# multiple entries can be added to this file, LAMMPS reads the ones it needs

# these entries are in LAMMPS "metal" units:

# format of a single entry (one or more lines):

#   element 1, element 2, element 3, 

#        m,  gamma, lambda3, c,     d,       h,          n,        beta,   lambda2, X_ij*B,  R,   D, lambda1, A

#E1 E2 E3 m gamma lambda3 c d h n beta lambda2 B R D lambda1 A

Si Si Si 1 0.013318 0 14 2.1 -1 0.78000 1 1.80821400248640 632.658058300867 2.35 0.15 2.38684248328205 1708.79738703139

Si Si C  1 0.013318 0 14 2.1 -1 0.78000 1 1.80821400248640 632.658058300867 2.35 0.15 2.38684248328205 1708.79738703139

Si C  Si 1 0.013318 0 14 2.1 -1 0.78000 1 1.96859970919818 428.946015420752 2.35 0.15 3.03361215187440 1820.05673775234

C  Si Si 1 0.011304 0 19 2.5 -1 0.80468 1 1.96859970919818 428.946015420752 2.35 0.15 3.03361215187440 1820.05673775234

C  C  Si 1 0.011304 0 19 2.5 -1 0.80469 1 1.76776695296637 203.208547714849 2.35 0.15 2.54558441227157 458.510465798439

C  Si C  1 0.011304 0 19 2.5 -1 0.80469 1 1.96859970919818 428.946015420752 2.35 0.15 3.03361215187440 1820.05673775234

Si C  C  1 0.013318 0 14 2.1 -1 0.78000 1 1.96859970919818 428.946015420752 2.35 0.15 3.03361215187440 1820.05673775234

C  C  C  1 0.011304 0 19 2.5 -1 0.80469 1 1.76776695296637 203.208547714849 2.35 0.15 2.54558441227157 458.510465798439

# DATE: 2016-05-06 CONTRIBUTOR: German Samolyuk, samolyuk@gmail.com CITATION: ???

# Gao-Weber parameters for various elements and mixtures

# multiple entries can be added to this file, LAMMPS reads the ones it needs

# these entries are in LAMMPS "metal" units:

# format of a single entry (one or more lines):

#   element 1, element 2, element 3, 

#        m,  gamma, lambda3, c,     d,       h,          n,        beta,   lambda2, X_ij*B,  R,   D, lambda1, A

#E1 E2 E3 m gamma lambda3 c d h n beta lambda2 B R D lambda1 A Z_i, Z_j, ZBLcut, ZBLexpscale

Si Si Si 1 0.013318 0 14 2.1 -1 0.78000 1 1.80821400248640 632.658058300867 2.35 0.15 2.38684248328205 1708.79738703139 14 14 .95   14

Si Si C  1 0.013318 0 14 2.1 -1 0.78000 1 1.80821400248640 632.658058300867 2.35 0.15 2.38684248328205 1708.79738703139 14 14 .95   14

Si C  Si 1 0.013318 0 14 2.1 -1 0.78000 1 1.96859970919818 428.946015420752 2.35 0.15 3.03361215187440 1820.05673775234 14 6  .95   14

C  Si Si 1 0.011304 0 19 2.5 -1 0.80468 1 1.96859970919818 428.946015420752 2.35 0.15 3.03361215187440 1820.05673775234 6  14 .95   14

C  C  Si 1 0.011304 0 19 2.5 -1 0.80469 1 1.76776695296637 203.208547714849 2.35 0.15 2.54558441227157 458.510465798439 6  6  .95   14

C  Si C  1 0.011304 0 19 2.5 -1 0.80469 1 1.96859970919818 428.946015420752 2.35 0.15 3.03361215187440 1820.05673775234 6  14 .95   14

Si C  C  1 0.013318 0 14 2.1 -1 0.78000 1 1.96859970919818 428.946015420752 2.35 0.15 3.03361215187440 1820.05673775234 14 6  .95   14

C  C  C  1 0.011304 0 19 2.5 -1 0.80469 1 1.76776695296637 203.208547714849 2.35 0.15 2.54558441227157 458.510465798439 6  6  .95   14

/* -*- c++ -*- ----------------------------------------------------------

   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator

   http://lammps.sandia.gov, Sandia National Laboratories

   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract

   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains

   certain rights in this software.  This software is distributed under

   the GNU General Public License.

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

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

#ifdef PAIR_CLASS

PairStyle(gw,PairGW)

#else

#ifndef LMP_PAIR_GW_H

#define LMP_PAIR_GW_H

#include "pair.h"

namespace LAMMPS_NS {

class PairGW : public Pair {

 public:

  PairGW(class LAMMPS *);

  virtual ~PairGW();

  virtual void compute(int, int);

  void settings(int, char **);

  void coeff(int, char **);

  void init_style();

  double init_one(int, int);

 protected:

  struct Param {

    double lam1,lam2,lam3;

    double c,d,h;

    double gamma,powerm;

    double powern,beta;

    double biga,bigb,bigd,bigr;

    double cut,cutsq;

    double c1,c2,c3,c4;

    int ielement,jelement,kelement;

    int powermint;

    double Z_i,Z_j;

    double ZBLcut,ZBLexpscale;

  };

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

  char **elements;              // names of unique elements

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

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

  double cutmax;                // max cutoff for all elements

  int nelements;                // # of unique elements

  int nparams;                  // # of stored parameter sets

  int maxparam;                 // max # of parameter sets

  int **pages;                     // neighbor list pages

  int maxlocal;                    // size of numneigh, firstneigh arrays

  int maxpage;                     // # of pages currently allocated

  int pgsize;                      // size of neighbor page

  int oneatom;                     // max # of neighbors for one atom

  int *GW_numneigh;             // # of pair neighbors for each atom

  int **GW_firstneigh;          // ptr to 1st neighbor of each atom

  void GW_neigh();

  void add_pages(int howmany = 1);

  void allocate();

  virtual void read_file(char *);

  void setup_params();

  virtual void repulsive(Param *, double, double &, int, double &);

  double zeta(Param *, double, double, double *, double *);

  virtual void force_zeta(Param *, double, double, double &,

                          double &, int, double &);

  void attractive(Param *, double, double, double, double *, double *,

                  double *, double *, double *);

  double gw_fc(double, Param *);

  double gw_fc_d(double, Param *);

  virtual double gw_fa(double, Param *);

  virtual double gw_fa_d(double, Param *);

  double gw_bij(double, Param *);

  double gw_bij_d(double, Param *);

  void gw_zetaterm_d(double, double *, double, double *, double,

                               double *, double *, double *, Param *);

  void costheta_d(double *, double, double *, double,

                  double *, double *, double *);

  // inlined functions for efficiency

  inline double gw_gijk(const double costheta,

                          const Param * const param) const {

    const double gw_c = param->c * param->c;

    const double gw_d = param->d * param->d;

    const double hcth = param->h - costheta;

	  //printf("gw_gijk: gw_c=%f gw_d=%f hcth=%f=%f-%f\n", gw_c, gw_d, hcth, param->h, costheta);

    return param->gamma*(1.0 + gw_c/gw_d - gw_c / (gw_d + hcth*hcth));

  }

  inline double gw_gijk_d(const double costheta,

                            const Param * const param) const {

    const double gw_c = param->c * param->c;

    const double gw_d = param->d * param->d;

    const double hcth = param->h - costheta;

    const double numerator = -2.0 * gw_c * hcth;

    const double denominator = 1.0/(gw_d + hcth*hcth);

    return param->gamma*numerator*denominator*denominator;

  }

  inline double vec3_dot(const double x[3], const double y[3]) const {

    return x[0]*y[0] + x[1]*y[1] + x[2]*y[2];

  }

  inline void vec3_add(const double x[3], const double y[3],

                       double * const z) const {

    z[0] = x[0]+y[0];  z[1] = x[1]+y[1];  z[2] = x[2]+y[2];

  }

  inline void vec3_scale(const double k, const double x[3],

                         double y[3]) const {

    y[0] = k*x[0];  y[1] = k*x[1];  y[2] = k*x[2];

  }

  inline void vec3_scaleadd(const double k, const double x[3],

                            const double y[3], double * const z) const {

    z[0] = k*x[0]+y[0];

    z[1] = k*x[1]+y[1];

    z[2] = k*x[2]+y[2];

  }

};

}

#endif

#endif

/* ERROR/WARNING messages:

E: Illegal ... command

Self-explanatory.  Check the input script syntax and compare to the

documentation for the command.  You can use -echo screen as a

command-line option when running LAMMPS to see the offending line.

E: Incorrect args for pair coefficients

Self-explanatory.  Check the input script or data file.

E: Pair style GW requires atom IDs

This is a requirement to use the GW potential.

E: Pair style GW requires newton pair on

See the newton command.  This is a restriction to use the GW

potential.

E: All pair coeffs are not set

All pair coefficients must be set in the data file or by the

pair_coeff command before running a simulation.

E: Cannot open GW potential file %s

The specified GW potential file cannot be opened.  Check that the

path and name are correct.

E: Incorrect format in GW potential file

Incorrect number of words per line in the potential file.

E: Illegal GW parameter

One or more of the coefficients defined in the potential file is

invalid.

E: Potential file has duplicate entry

The potential file for a SW or GW potential has more than

one entry for the same 3 ordered elements.

E: Potential file is missing an entry

The potential file for a SW or GW potential does not have a

needed entry.

*/

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

   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator

   http://lammps.sandia.gov, Sandia National Laboratories

   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract

   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains

   certain rights in this software.  This software is distributed under

   the GNU General Public License.

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

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

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

   Contributing author: German Samolyuk (ORNL)

   Based on PairTersoffZBL by Aidan Thompson (SNL) and David Farrell (NWU)

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

#include <math.h>

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include "pair_gw_zbl.h"

#include "atom.h"

#include "update.h"

#include "neighbor.h"

#include "neigh_list.h"

#include "neigh_request.h"

#include "force.h"

#include "comm.h"

#include "memory.h"

#include "error.h"

#include "math_const.h"

using namespace LAMMPS_NS;

using namespace MathConst;

#define MAXLINE 1024

#define DELTA 4

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

PairGWZBL::PairGWZBL(LAMMPS *lmp) : PairGW(lmp)

{

  // hard-wired constants in metal or real units

  // a0 = Bohr radius

  // epsilon0 = permittivity of vacuum = q / energy-distance units

  // e = unit charge

  // 1 Kcal/mole = 0.043365121 eV

  if (strcmp(update->unit_style,"metal") == 0) {

    global_a_0 = 0.529;

    global_epsilon_0 = 0.00552635;

    global_e = 1.0;

  } else if (strcmp(update->unit_style,"real") == 0) {

    global_a_0 = 0.529;

    global_epsilon_0 = 0.00552635 * 0.043365121;

    global_e = 1.0;

  } else error->all(FLERR,"Pair gw/zbl requires metal or real units");

}

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

void PairGWZBL::read_file(char *file)

{

  int params_per_line = 21;

  char **words = new char*[params_per_line+1];

  memory->sfree(params);

  params = NULL;

  nparams = maxparam = 0;

  // open file on proc 0

  FILE *fp;

  if (comm->me == 0) {

    fp = force->open_potential(file);

    if (fp == NULL) {

      char str[128];

      sprintf(str,"Cannot open GW potential file %s",file);

      error->one(FLERR,str);

    }

  }

  // read each line out of file, skipping blank lines or leading '#'

  // store line of params if all 3 element tags are in element list

  int n,nwords,ielement,jelement,kelement;

  char line[MAXLINE],*ptr;

  int eof = 0;

  while (1) {

    if (comm->me == 0) {

      ptr = fgets(line,MAXLINE,fp);

      if (ptr == NULL) {

        eof = 1;

        fclose(fp);

      } else n = strlen(line) + 1;

    }

    MPI_Bcast(&eof,1,MPI_INT,0,world);

    if (eof) break;

    MPI_Bcast(&n,1,MPI_INT,0,world);

    MPI_Bcast(line,n,MPI_CHAR,0,world);

    // strip comment, skip line if blank

    if ((ptr = strchr(line,'#'))) *ptr = '\0';

    nwords = atom->count_words(line);

    if (nwords == 0) continue;

    // concatenate additional lines until have params_per_line words

    while (nwords < params_per_line) {

      n = strlen(line);

      if (comm->me == 0) {

        ptr = fgets(&line[n],MAXLINE-n,fp);

        if (ptr == NULL) {

          eof = 1;

          fclose(fp);

        } else n = strlen(line) + 1;

      }

      MPI_Bcast(&eof,1,MPI_INT,0,world);

      if (eof) break;

      MPI_Bcast(&n,1,MPI_INT,0,world);

      MPI_Bcast(line,n,MPI_CHAR,0,world);

      if ((ptr = strchr(line,'#'))) *ptr = '\0';

      nwords = atom->count_words(line);

    }

    if (nwords != params_per_line)

      error->all(FLERR,"Incorrect format in GW potential file");

    // words = ptrs to all words in line

    nwords = 0;

    words[nwords++] = strtok(line," \t\n\r\f");

    while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;

    // ielement,jelement,kelement = 1st args

    // if all 3 args are in element list, then parse this line

    // else skip to next line

    for (ielement = 0; ielement < nelements; ielement++)

      if (strcmp(words[0],elements[ielement]) == 0) break;

    if (ielement == nelements) continue;

    for (jelement = 0; jelement < nelements; jelement++)

      if (strcmp(words[1],elements[jelement]) == 0) break;

    if (jelement == nelements) continue;

    for (kelement = 0; kelement < nelements; kelement++)

      if (strcmp(words[2],elements[kelement]) == 0) break;

    if (kelement == nelements) continue;

    // load up parameter settings and error check their values

    if (nparams == maxparam) {

      maxparam += DELTA;

      params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),

                                          "pair:params");

    }

    params[nparams].ielement = ielement;

    params[nparams].jelement = jelement;

    params[nparams].kelement = kelement;

    params[nparams].powerm = atof(words[3]);

    params[nparams].gamma = atof(words[4]);

    params[nparams].lam3 = atof(words[5]);

    params[nparams].c = atof(words[6]);

    params[nparams].d = atof(words[7]);

    params[nparams].h = atof(words[8]);

    params[nparams].powern = atof(words[9]);

    params[nparams].beta = atof(words[10]);

    params[nparams].lam2 = atof(words[11]);

    params[nparams].bigb = atof(words[12]);

    params[nparams].bigr = atof(words[13]);

    params[nparams].bigd = atof(words[14]);

    params[nparams].lam1 = atof(words[15]);

    params[nparams].biga = atof(words[16]);

    params[nparams].Z_i = atof(words[17]);

    params[nparams].Z_j = atof(words[18]);

    params[nparams].ZBLcut = atof(words[19]);

    params[nparams].ZBLexpscale = atof(words[20]);

    // currently only allow m exponent of 1 or 3

    params[nparams].powermint = int(params[nparams].powerm);

    if (

        params[nparams].lam3 < 0.0 || params[nparams].c < 0.0 ||

        params[nparams].d < 0.0 || params[nparams].powern < 0.0 ||

        params[nparams].beta < 0.0 || params[nparams].lam2 < 0.0 ||

        params[nparams].bigb < 0.0 || params[nparams].bigr < 0.0 ||

        params[nparams].bigd < 0.0 ||

        params[nparams].bigd > params[nparams].bigr ||

        params[nparams].lam3 < 0.0 || params[nparams].biga < 0.0 ||

        params[nparams].powerm - params[nparams].powermint != 0.0 ||

        (params[nparams].powermint != 3 && params[nparams].powermint != 1) ||

        params[nparams].gamma < 0.0 ||

        params[nparams].Z_i < 1.0 || params[nparams].Z_j < 1.0 ||

        params[nparams].ZBLcut < 0.0 || params[nparams].ZBLexpscale < 0.0)

      error->all(FLERR,"Illegal GW parameter");

    nparams++;

  }

  delete [] words;

}

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

void PairGWZBL::repulsive(Param *param, double rsq, double &fforce,

                               int eflag, double &eng)

{

  double r,tmp_fc,tmp_fc_d,tmp_exp;

  // GW repulsive portion

  r = sqrt(rsq);

  tmp_fc = gw_fc(r,param);

  tmp_fc_d = gw_fc_d(r,param);

  tmp_exp = exp(-param->lam1 * r);

  double fforce_gw = param->biga * tmp_exp * (tmp_fc_d - tmp_fc*param->lam1);

  double eng_gw = tmp_fc * param->biga * tmp_exp;

  // ZBL repulsive portion

  double esq = pow(global_e,2.0);

  double a_ij = (0.8854*global_a_0) /

    (pow(param->Z_i,0.23) + pow(param->Z_j,0.23));

  double premult = (param->Z_i * param->Z_j * esq)/(4.0*MY_PI*global_epsilon_0);

  double r_ov_a = r/a_ij;

  double phi = 0.1818*exp(-3.2*r_ov_a) + 0.5099*exp(-0.9423*r_ov_a) +

    0.2802*exp(-0.4029*r_ov_a) + 0.02817*exp(-0.2016*r_ov_a);

  double dphi = (1.0/a_ij) * (-3.2*0.1818*exp(-3.2*r_ov_a) -

                              0.9423*0.5099*exp(-0.9423*r_ov_a) -

                              0.4029*0.2802*exp(-0.4029*r_ov_a) -

                              0.2016*0.02817*exp(-0.2016*r_ov_a));

  double fforce_ZBL = premult*-phi/rsq + premult*dphi/r;

  double eng_ZBL = premult*(1.0/r)*phi;

  // combine two parts with smoothing by Fermi-like function

  fforce = -(-F_fermi_d(r,param) * eng_ZBL +

             (1.0 - F_fermi(r,param))*fforce_ZBL +

             F_fermi_d(r,param)*eng_gw + F_fermi(r,param)*fforce_gw) / r;

  if (eflag)

    eng = (1.0 - F_fermi(r,param))*eng_ZBL + F_fermi(r,param)*eng_gw;

}

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

double PairGWZBL::gw_fa(double r, Param *param)

{

  if (r > param->bigr + param->bigd) return 0.0;

  return -param->bigb * exp(-param->lam2 * r) * gw_fc(r,param) *

    F_fermi(r,param);

}

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

double PairGWZBL::gw_fa_d(double r, Param *param)

{

  if (r > param->bigr + param->bigd) return 0.0;

  return param->bigb * exp(-param->lam2 * r) *

    (param->lam2 * gw_fc(r,param) * F_fermi(r,param) -

     gw_fc_d(r,param) * F_fermi(r,param) - gw_fc(r,param) *

     F_fermi_d(r,param));

}

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

   Fermi-like smoothing function

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

double PairGWZBL::F_fermi(double r, Param *param)

{

  return 1.0 / (1.0 + exp(-param->ZBLexpscale*(r-param->ZBLcut)));

}

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

   Fermi-like smoothing function derivative with respect to r

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

double PairGWZBL::F_fermi_d(double r, Param *param)

{

  return param->ZBLexpscale*exp(-param->ZBLexpscale*(r-param->ZBLcut)) /

    pow(1.0 + exp(-param->ZBLexpscale*(r-param->ZBLcut)),2.0);

}