diff --git a/src/MANYBODY/pair_gw_zbl.cpp b/src/MANYBODY/pair_gw_zbl.cpp index 27ab1a3bde..514892d527 100644 --- a/src/MANYBODY/pair_gw_zbl.cpp +++ b/src/MANYBODY/pair_gw_zbl.cpp @@ -28,6 +28,8 @@ #include "memory.h" #include "error.h" #include "utils.h" +#include "tokenizer.h" +#include "potential_file_reader.h" #include "math_const.h" using namespace LAMMPS_NS; @@ -61,146 +63,101 @@ PairGWZBL::PairGWZBL(LAMMPS *lmp) : PairGW(lmp) 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]; - snprintf(str,128,"Cannot open GW potential file %s",file); - error->one(FLERR,str); - } - } + PotentialFileReader reader(lmp, file, "GW/ZBL"); + char * line; - // read each line out of file, skipping blank lines or leading '#' - // store line of params if all 3 element tags are in element list + while(line = reader.next_line(NPARAMS_PER_LINE)) { + try { + ValueTokenizer values(line, " \t\n\r\f"); - int n,nwords,ielement,jelement,kelement; - char line[MAXLINE],*ptr; - int eof = 0; + std::string iname = values.next_string(); + std::string jname = values.next_string(); + std::string kname = values.next_string(); - 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); + // ielement,jelement,kelement = 1st args + // if all 3 args are in element list, then parse this line + // else skip to next entry in file + int ielement, jelement, kelement; - // strip comment, skip line if blank + for (ielement = 0; ielement < nelements; ielement++) + if (iname == elements[ielement]) break; + if (ielement == nelements) continue; + for (jelement = 0; jelement < nelements; jelement++) + if (jname == elements[jelement]) break; + if (jelement == nelements) continue; + for (kelement = 0; kelement < nelements; kelement++) + if (kname == elements[kelement]) break; + if (kelement == nelements) continue; - if ((ptr = strchr(line,'#'))) *ptr = '\0'; - nwords = utils::count_words(line); - if (nwords == 0) continue; + // load up parameter settings and error check their values - // concatenate additional lines until have params_per_line words + if (nparams == maxparam) { + maxparam += DELTA; + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), + "pair:params"); + } - 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; + params[nparams].ielement = ielement; + params[nparams].jelement = jelement; + params[nparams].kelement = kelement; + params[nparams].powerm = values.next_double(); + params[nparams].gamma = values.next_double(); + params[nparams].lam3 = values.next_double(); + params[nparams].c = values.next_double(); + params[nparams].d = 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].Z_i = values.next_double(); + params[nparams].Z_j = values.next_double(); + params[nparams].ZBLcut = values.next_double(); + params[nparams].ZBLexpscale = values.next_double(); + params[nparams].powermint = int(params[nparams].powerm); + } catch (TokenizerException & e) { + error->one(FLERR, e.what()); } - 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); + + // currently only allow m exponent of 1 or 3 + 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->one(FLERR,"Illegal GW parameter"); + + nparams++; } - - 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; + 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"); + } + + MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world); } /* ---------------------------------------------------------------------- */ diff --git a/src/MANYBODY/pair_gw_zbl.h b/src/MANYBODY/pair_gw_zbl.h index 62b3f9a43f..95129c4eb8 100644 --- a/src/MANYBODY/pair_gw_zbl.h +++ b/src/MANYBODY/pair_gw_zbl.h @@ -29,6 +29,8 @@ class PairGWZBL : public PairGW { PairGWZBL(class LAMMPS *); ~PairGWZBL() {} + static const int NPARAMS_PER_LINE = 21; + private: double global_a_0; // Bohr radius for Coulomb repulsion double global_epsilon_0; // permittivity of vacuum for Coulomb repulsion diff --git a/unittest/formats/test_potential_file_reader.cpp b/unittest/formats/test_potential_file_reader.cpp index e948b23222..d5d9c61bbe 100644 --- a/unittest/formats/test_potential_file_reader.cpp +++ b/unittest/formats/test_potential_file_reader.cpp @@ -11,6 +11,7 @@ #include "MANYBODY/pair_tersoff_mod_c.h" #include "MANYBODY/pair_tersoff_zbl.h" #include "MANYBODY/pair_gw.h" +#include "MANYBODY/pair_gw_zbl.h" #include @@ -24,6 +25,7 @@ const int LAMMPS_NS::PairTersoffMOD::NPARAMS_PER_LINE; const int LAMMPS_NS::PairTersoffMODC::NPARAMS_PER_LINE; const int LAMMPS_NS::PairTersoffZBL::NPARAMS_PER_LINE; const int LAMMPS_NS::PairGW::NPARAMS_PER_LINE; +const int LAMMPS_NS::PairGWZBL::NPARAMS_PER_LINE; class PotenialFileReaderTest : public ::testing::Test { protected: @@ -117,6 +119,15 @@ TEST_F(PotenialFileReaderTest, GW) { ASSERT_EQ(utils::count_words(line), PairGW::NPARAMS_PER_LINE); } +TEST_F(PotenialFileReaderTest, GWZBL) { + ::testing::internal::CaptureStdout(); + PotentialFileReader reader(lmp, "SiC.gw.zbl", "GWZBL"); + ::testing::internal::GetCapturedStdout(); + + auto line = reader.next_line(PairGWZBL::NPARAMS_PER_LINE); + ASSERT_EQ(utils::count_words(line), PairGWZBL::NPARAMS_PER_LINE); +} + int main(int argc, char **argv) { MPI_Init(&argc, &argv);