Commit 06cc38e1 authored by Anders Hafreager's avatar Anders Hafreager Committed by Axel Kohlmeyer
Browse files

Fixed so tabulated pair_vashishta uses same pair_modify command style as other pair styles

parent 10ec14f0

src/MANYBODY/pair_vashishta.cpp

+58 −46
Original line number Diff line number Diff line
@@ -32,7 +32,6 @@ 
#include "neigh_list.h"

#include "memory.h"

#include "error.h"



using namespace LAMMPS_NS;



#define MAXLINE 1024
@@ -58,8 +57,8 @@ PairVashishta::PairVashishta(LAMMPS *lmp) : Pair(lmp) 
  neigh3BodyCount = NULL; 

  neigh3Body = NULL;

  

  useTable = true;

  tableSize = 65536;

  useTable = false;

  nTablebits = 12;

  forceTable = NULL;

  potentialTable = NULL;

}
@@ -90,25 +89,27 @@ PairVashishta::~PairVashishta() 


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

{

  if(narg < 2 || narg > 3) // We accept table yes/no [tableSize]

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

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

    error->all(FLERR,"Illegal pair_modify command. Expected 'table yes/no [tableSize]'");

  if(strcmp(arg[1], "yes") == 0) 

  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 > sizeof(float)*CHAR_BIT)

        error->all(FLERR,"Too many total bits for bitmapped lookup table");

      if(nTablebits != 0) 

        useTable = true;

  else if(strcmp(arg[1], "no") == 0) 

    useTable = false;

  else 

    error->all(FLERR,"Illegal pair_modify command. Expected 'table yes/no [tableSize]'");



  if(narg == 3) {

    tableSize = atoi(arg[2]);

    if(tableSize <= 0) 

      error->all(FLERR,"Illegal pair_modify command. Expected 'table yes/no [tableSize]', but tableSize has value 0.");

      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]);

      iarg += 2;

    }

  }



  createTable();

}

}



void PairVashishta::validateNeigh3Body() {

  if (atom->nlocal > neigh3BodyMax) {
@@ -123,27 +124,32 @@ void PairVashishta::validateNeigh3Body() { 
void PairVashishta::createTable()

{

  int ntypes = atom->ntypes+1;

  deltaR2 = cutmax*cutmax / (tableSize-1);

  tabinnersq = tabinner*tabinner;

  

  int ntable = 1;

  for (int i = 0; i < nTablebits; i++) ntable *= 2;



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

  oneOverDeltaR2 = 1.0/deltaR2;



  memory->destroy(forceTable);

  memory->destroy(potentialTable);

  memory->create(forceTable,ntypes,ntypes,tableSize+1,"pair:vashishta:forceTable");

  memory->create(potentialTable,ntypes,ntypes,tableSize+1,"pair:vashishta:potentialTable");

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

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



  for (int ii = 0; ii < ntypes; ii++) {

    int i = map[ii];

    for (int jj = 0; jj < ntypes; jj++) {

      int j = map[jj];

      if (i < 0 || j < 0 || ii == 0 || jj == 0) {

            for(int tableIndex=0; tableIndex<=tableSize; tableIndex++) {

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

            forceTable[ii][jj][tableIndex] = 0;

            potentialTable[ii][jj][tableIndex] = 0;

        }

      } else {

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

            for(int tableIndex=0; tableIndex<=tableSize; tableIndex++) {

                double rsq = tableIndex*deltaR2;

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

            double rsq = tabinnersq + tableIndex*deltaR2;

            double fpair, eng;

            twobody(&params[ijparam], rsq, fpair, 1, eng, false /* Don't ask for tabulated since we are generating it now*/ );

            forceTable[ii][jj][tableIndex] = fpair;
@@ -637,16 +643,22 @@ void PairVashishta::twobody(Param *param, double rsq, double &fforce, 
                     int eflag, double &eng, bool tabulated)

{

  if(tabulated) {

    int potentialTableIndex = rsq*oneOverDeltaR2;

    double fractionalRest = rsq*oneOverDeltaR2 - potentialTableIndex; // double - int will only keep the 0.xxxx part

    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);

    }



    int tableIndex = (rsq - tabinnersq)*oneOverDeltaR2;

    double fraction = (rsq - tabinnersq)*oneOverDeltaR2 - tableIndex; // double - int will only keep the 0.xxxx part



    double force0 = forceTable[param->ielement+1][param->jelement+1][potentialTableIndex];

    double force1 = forceTable[param->ielement+1][param->jelement+1][potentialTableIndex+1];

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

    double force0 = forceTable[param->ielement+1][param->jelement+1][tableIndex];

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

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

    if(evflag) {

        double energy0 = potentialTable[param->ielement+1][param->jelement+1][potentialTableIndex];

        double energy1 = potentialTable[param->ielement+1][param->jelement+1][potentialTableIndex+1];

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

        double energy0 = potentialTable[param->ielement+1][param->jelement+1][tableIndex];

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

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

    }

  } else {

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

src/MANYBODY/pair_vashishta.h

+4 −2
Original line number Diff line number Diff line
@@ -46,8 +46,9 @@ class PairVashishta : public Pair { 
    double lam1rc,lam4rc,vrcc2,vrcc3,vrc,dvrc,c0;

    int ielement,jelement,kelement;

  };



  bool useTable;

  int tableSize;

  int nTablebits;

  double deltaR2;

  double oneOverDeltaR2;

  double ***forceTable;          // Table containing the forces
@@ -62,6 +63,7 @@ class PairVashishta : public Pair { 
  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
@@ -71,7 +73,7 @@ class PairVashishta : public Pair { 
  void setup_params();

  void createTable();

  void validateNeigh3Body();

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

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

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

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

};

CRA Git | Maintained and supported by SUSTech CRA and CCSE