From 00332d299b2bbd60129c2ad70fd0913d719e15b7 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 24 Jun 2020 21:05:00 -0400 Subject: [PATCH] add unit conversion support to pair style tersoff/table --- src/USER-MISC/pair_tersoff_table.cpp | 218 ++++++++---------- src/USER-MISC/pair_tersoff_table.h | 4 +- unittest/formats/test_pair_unit_convert.cpp | 43 ++++ .../formats/test_potential_file_reader.cpp | 13 ++ 4 files changed, 154 insertions(+), 124 deletions(-) diff --git a/src/USER-MISC/pair_tersoff_table.cpp b/src/USER-MISC/pair_tersoff_table.cpp index 8ab64b1d5a..be8e4d701a 100644 --- a/src/USER-MISC/pair_tersoff_table.cpp +++ b/src/USER-MISC/pair_tersoff_table.cpp @@ -33,6 +33,8 @@ #include "comm.h" #include "memory.h" #include "utils.h" +#include "tokenizer.h" +#include "potential_file_reader.h" #include "error.h" @@ -59,6 +61,7 @@ PairTersoffTable::PairTersoffTable(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); nelements = 0; elements = nullptr; @@ -840,146 +843,115 @@ double PairTersoffTable::init_one(int i, int j) void PairTersoffTable::read_file(char *file) { - int params_per_line = 17; - 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); - } - } + PotentialFileReader reader(lmp, file, "TersoffTable", unit_convert_flag); + char *line; + + // transparently convert units for supported conversions - // read each set of params from potential file - // one set of params can span multiple lines - // store params if all 3 element tags are in element list + int unit_convert = reader.get_unit_convert(); + double conversion_factor = utils::get_conversion_factor(utils::ENERGY, + unit_convert); - int n,nwords,ielement,jelement,kelement; - char line[MAXLINE],*ptr; - int eof = 0; + while((line = reader.next_line(NPARAMS_PER_LINE))) { - 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); + try { + ValueTokenizer values(line); - // strip comment, skip line if blank + std::string iname = values.next_string(); + std::string jname = values.next_string(); + std::string kname = values.next_string(); - if ((ptr = strchr(line,'#'))) *ptr = '\0'; - nwords = utils::count_words(line); - if (nwords == 0) continue; + // ielement,jelement,kelement = 1st args + // if all 3 args are in element list, then parse this line + // else skip to next entry in file - // concatenate additional lines until have params_per_line words + int ielement, jelement, kelement; - 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; + 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; + + // load up parameter settings and error check their values + + if (nparams == maxparam) { + maxparam += DELTA; + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), + "pair:params"); + } + + // some parameters are not used since only Tersoff_2 is implemented + + params[nparams].ielement = ielement; + params[nparams].jelement = jelement; + params[nparams].kelement = kelement; + params[nparams].powerm = values.next_double(); // not used + params[nparams].gamma = values.next_double(); // not used + params[nparams].lam3 = values.next_double(); // not used + 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(); + double bigr = values.next_double(); + double bigd = values.next_double(); + params[nparams].cutoffR = bigr - bigd; + params[nparams].cutoffS = bigr + bigd; + params[nparams].lam1 = values.next_double(); + params[nparams].biga = values.next_double(); + + if (unit_convert) { + params[nparams].biga *= conversion_factor; + params[nparams].bigb *= conversion_factor; + } + } 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); + + if (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].cutoffR < 0.0 || + params[nparams].cutoffS < 0.0 || + params[nparams].cutoffR > params[nparams].cutoffS || + params[nparams].lam1 < 0.0 || + params[nparams].biga < 0.0 + ) error->one(FLERR,"Illegal Tersoff parameter"); + + // only tersoff_2 parametrization is implemented + + if (params[nparams].gamma != 1.0 || params[nparams].lam3 != 0.0) + error->one(FLERR,"Currently the tersoff/table pair_style only " + "implements the Tersoff_2 parametrization"); + nparams++; } - - if (nwords != params_per_line) - error->all(FLERR,"Incorrect format in Tersoff 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 entry in file - - 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]); // not used (only tersoff_2 is implemented) - params[nparams].gamma = atof(words[4]); // not used (only tersoff_2 is implemented) - params[nparams].lam3 = atof(words[5]); // not used (only tersoff_2 is implemented) - 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]); - - // current implementation is based on functional form - // of tersoff_2 as reported in the reference paper - - double bigr = atof(words[13]); - double bigd = atof(words[14]); - params[nparams].cutoffR = bigr - bigd; - params[nparams].cutoffS = bigr + bigd; - params[nparams].lam1 = atof(words[15]); - params[nparams].biga = atof(words[16]); - - if (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].cutoffR < 0.0 || - params[nparams].cutoffS < 0.0 || - params[nparams].cutoffR > params[nparams].cutoffS || - params[nparams].lam1 < 0.0 || - params[nparams].biga < 0.0 - ) error->all(FLERR,"Illegal Tersoff parameter"); - - // only tersoff_2 parametrization is implemented - if (params[nparams].gamma != 1.0 || params[nparams].lam3 != 0.0) - error->all(FLERR,"Current tersoff/table pair_style implements only tersoff_2 parametrization"); - 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/USER-MISC/pair_tersoff_table.h b/src/USER-MISC/pair_tersoff_table.h index cd5549c345..faa704191c 100644 --- a/src/USER-MISC/pair_tersoff_table.h +++ b/src/USER-MISC/pair_tersoff_table.h @@ -40,8 +40,10 @@ class PairTersoffTable : public Pair { virtual void compute(int, int); void settings(int, char **); void coeff(int, char **); - double init_one(int, int); void init_style(); + double init_one(int, int); + + static const int NPARAMS_PER_LINE = 17; protected: struct Param { diff --git a/unittest/formats/test_pair_unit_convert.cpp b/unittest/formats/test_pair_unit_convert.cpp index 2583ee474b..e77d7e535b 100644 --- a/unittest/formats/test_pair_unit_convert.cpp +++ b/unittest/formats/test_pair_unit_convert.cpp @@ -214,6 +214,49 @@ TEST_F(PairUnitConvertTest, tersoff) EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error)); } +TEST_F(PairUnitConvertTest, tersoff_table) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "tersoff/table")) GTEST_SKIP(); + + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("units metal"); + lmp->input->one("read_data test_pair_unit_convert.data"); + lmp->input->one("pair_style tersoff/table"); + lmp->input->one("pair_coeff * * SiC.tersoff Si C"); + lmp->input->one("run 0 post no"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + // copy pressure, energy, and force from first step + double pold; + lmp->output->thermo->evaluate_keyword("press", &pold); + double eold = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul; + double **f = lmp->atom->f; + for (int i = 0; i < 4; ++i) + for (int j = 0; j < 3; ++j) + fold[i][j] = f[i][j]; + + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("clear"); + lmp->input->one("units real"); + lmp->input->one("read_data test_pair_unit_convert.data"); + lmp->input->one("pair_style tersoff/table"); + lmp->input->one("pair_coeff * * SiC.tersoff Si C"); + lmp->input->one("run 0 post no"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + double pnew; + lmp->output->thermo->evaluate_keyword("press", &pnew); + EXPECT_NEAR(pold, p_convert * pnew, fabs(pnew * rel_error)); + double enew = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul; + EXPECT_NEAR(ev_convert * eold, enew, fabs(enew * rel_error)); + + f = lmp->atom->f; + for (int i = 0; i < 4; ++i) + for (int j = 0; j < 3; ++j) + EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error)); +} + TEST_F(PairUnitConvertTest, sw) { // check if the prerequisite pair style is available diff --git a/unittest/formats/test_potential_file_reader.cpp b/unittest/formats/test_potential_file_reader.cpp index 54dc26e4dc..933fa66b96 100644 --- a/unittest/formats/test_potential_file_reader.cpp +++ b/unittest/formats/test_potential_file_reader.cpp @@ -23,6 +23,7 @@ #include "MANYBODY/pair_tersoff_mod_c.h" #include "MANYBODY/pair_tersoff_zbl.h" #include "MANYBODY/pair_vashishta.h" +#include "USER-MISC/pair_tersoff_table.h" #include "input.h" #include "lammps.h" #include "potential_file_reader.h" @@ -50,6 +51,7 @@ const int LAMMPS_NS::PairGW::NPARAMS_PER_LINE; const int LAMMPS_NS::PairGWZBL::NPARAMS_PER_LINE; const int LAMMPS_NS::PairNb3bHarmonic::NPARAMS_PER_LINE; const int LAMMPS_NS::PairVashishta::NPARAMS_PER_LINE; +const int LAMMPS_NS::PairTersoffTable::NPARAMS_PER_LINE; class PotentialFileReaderTest : public ::testing::Test { protected: @@ -140,6 +142,17 @@ TEST_F(PotentialFileReaderTest, TersoffModC) ASSERT_EQ(utils::count_words(line), PairTersoffMODC::NPARAMS_PER_LINE); } +TEST_F(PotentialFileReaderTest, TersoffTable) +{ + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("units metal"); + PotentialFileReader reader(lmp, "Si.tersoff", "TersoffTable"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + auto line = reader.next_line(PairTersoffTable::NPARAMS_PER_LINE); + ASSERT_EQ(utils::count_words(line), PairTersoffTable::NPARAMS_PER_LINE); +} + TEST_F(PotentialFileReaderTest, TersoffZBL) { if (!verbose) ::testing::internal::CaptureStdout();