Unverified Commit a4db35e1 authored by Richard Berger's avatar Richard Berger
Browse files

Update pair tersoff mod

parent 8329a549
Loading
Loading
Loading
Loading
+84 −127
Original line number Diff line number Diff line
@@ -29,12 +29,13 @@
#include "memory.h"
#include "error.h"
#include "utils.h"
#include "tokenizer.h"
#include "potential_file_reader.h"

using namespace LAMMPS_NS;
using namespace MathConst;
using namespace MathSpecial;

#define MAXLINE 1024
#define DELTA 4

/* ---------------------------------------------------------------------- */
@@ -45,93 +46,40 @@ PairTersoffMOD::PairTersoffMOD(LAMMPS *lmp) : PairTersoff(lmp) {}

void PairTersoffMOD::read_file(char *file)
{
  int params_per_line = NPARAMS_PER_LINE;
  char **words = new char*[params_per_line+1];

  memory->sfree(params);
  params = NULL;
  params = nullptr;
  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];
      snprintf(str,128,"Cannot open Tersoff 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 = utils::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 = utils::count_words(line);
    }
    PotentialFileReader reader(lmp, file, "Tersoff");
    char * line;

    if (nwords != params_per_line)
      error->all(FLERR,"Incorrect format in Tersoff potential file");
    while(line = reader.next_line(NPARAMS_PER_LINE)) {
      try {
        ValueTokenizer values(line, " \t\n\r\f");

    // 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;
        std::string iname = values.next_string();
        std::string jname = values.next_string();
        std::string kname = values.next_string();

        // ielement,jelement,kelement = 1st args
        // if all 3 args are in element list, then parse this line
    // else skip to next line
        // else skip to next entry in file
        int ielement, jelement, kelement;

        for (ielement = 0; ielement < nelements; ielement++)
      if (strcmp(words[0],elements[ielement]) == 0) break;
          if (iname == elements[ielement]) break;
        if (ielement == nelements) continue;
        for (jelement = 0; jelement < nelements; jelement++)
      if (strcmp(words[1],elements[jelement]) == 0) break;
          if (jname == elements[jelement]) break;
        if (jelement == nelements) continue;
        for (kelement = 0; kelement < nelements; kelement++)
      if (strcmp(words[2],elements[kelement]) == 0) break;
          if (kname == elements[kelement]) break;
        if (kelement == nelements) continue;


        // load up parameter settings and error check their values

        if (nparams == maxparam) {
@@ -143,28 +91,29 @@ void PairTersoffMOD::read_file(char *file)
        params[nparams].ielement = ielement;
        params[nparams].jelement = jelement;
        params[nparams].kelement = kelement;
    params[nparams].powerm = atof(words[3]);
    params[nparams].lam3 = atof(words[4]);
    params[nparams].h = atof(words[5]);
    params[nparams].powern = atof(words[6]);
    params[nparams].beta = atof(words[7]);
    params[nparams].lam2 = atof(words[8]);
    params[nparams].bigb = atof(words[9]);
    params[nparams].bigr = atof(words[10]);
    params[nparams].bigd = atof(words[11]);
    params[nparams].lam1 = atof(words[12]);
    params[nparams].biga = atof(words[13]);
    params[nparams].powern_del = atof(words[14]);
    params[nparams].c1 = atof(words[15]);
    params[nparams].c2 = atof(words[16]);
    params[nparams].c3 = atof(words[17]);
    params[nparams].c4 = atof(words[18]);
    params[nparams].c5 = atof(words[19]);

    // currently only allow m exponent of 1 or 3

        params[nparams].powerm     = values.next_double();
        params[nparams].lam3       = values.next_double();
        params[nparams].h          = values.next_double();
        params[nparams].powern     = values.next_double();
        params[nparams].beta       = values.next_double();
        params[nparams].lam2       = values.next_double();
        params[nparams].bigb       = values.next_double();
        params[nparams].bigr       = values.next_double();
        params[nparams].bigd       = values.next_double();
        params[nparams].lam1       = values.next_double();
        params[nparams].biga       = values.next_double();
        params[nparams].powern_del = values.next_double();
        params[nparams].c1         = values.next_double(); 
        params[nparams].c2         = values.next_double(); 
        params[nparams].c3         = values.next_double(); 
        params[nparams].c4         = values.next_double(); 
        params[nparams].c5         = values.next_double(); 
        params[nparams].powermint = int(params[nparams].powerm);
      } catch (TokenizerException & e) {
        error->one(FLERR, e.what());
      }

      // currently only allow m exponent of 1 or 3
      if (params[nparams].powern < 0.0 ||
          params[nparams].beta < 0.0 ||
          params[nparams].lam2 < 0.0 ||
@@ -182,8 +131,16 @@ void PairTersoffMOD::read_file(char *file)

      nparams++;
    }
  }

  MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
  MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);

  if(comm->me != 0) {
    params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params");
  }

  delete [] words;
  MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world);
}

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