From 2acf71c3e2cb70cc444de6e75b168b13e43a2e6d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 25 Jun 2020 04:31:34 -0400 Subject: [PATCH] add unit conversion to table pair style --- src/pair_table.cpp | 46 +++++---- src/table_file_reader.cpp | 11 +- src/table_file_reader.h | 5 +- src/utils.cpp | 2 +- unittest/formats/test_pair_unit_convert.cpp | 105 ++++++++++++++++++++ 5 files changed, 140 insertions(+), 29 deletions(-) diff --git a/src/pair_table.cpp b/src/pair_table.cpp index 312420c354..f15d5a2cea 100644 --- a/src/pair_table.cpp +++ b/src/pair_table.cpp @@ -44,6 +44,7 @@ PairTable::PairTable(LAMMPS *lmp) : Pair(lmp) { ntables = 0; tables = NULL; + unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); } /* ---------------------------------------------------------------------- */ @@ -360,7 +361,13 @@ void PairTable::read_table(Table *tb, char *file, char *keyword) { TableFileReader reader(lmp, file, "pair"); - char * line = reader.find_section_start(keyword); + // transparently convert units for supported conversions + + int unit_convert = reader.get_unit_convert(); + double conversion_factor = utils::get_conversion_factor(utils::ENERGY, + unit_convert); + + char *line = reader.find_section_start(keyword); if (!line) { error->one(FLERR,"Did not find keyword in table file"); @@ -370,7 +377,7 @@ void PairTable::read_table(Table *tb, char *file, char *keyword) // allocate table arrays for file values line = reader.next_line(); - param_extract(tb,line); + param_extract(tb, line); memory->create(tb->rfile,tb->ninput,"pair:rfile"); memory->create(tb->efile,tb->ninput,"pair:efile"); memory->create(tb->ffile,tb->ninput,"pair:ffile"); @@ -404,8 +411,8 @@ void PairTable::read_table(Table *tb, char *file, char *keyword) ValueTokenizer values(line); values.next_int(); rfile = values.next_double(); - tb->efile[i] = values.next_double(); - tb->ffile[i] = values.next_double(); + tb->efile[i] = conversion_factor * values.next_double(); + tb->ffile[i] = conversion_factor * values.next_double(); } catch (TokenizerException & e) { ++cerror; } @@ -461,29 +468,26 @@ void PairTable::read_table(Table *tb, char *file, char *keyword) } } - if (ferror) { - std::string str = fmt::format("{} of {} force values in table are inconsistent with -dE/dr.\n" - " Should only be flagged at inflection points",ferror,tb->ninput); - error->warning(FLERR,str.c_str()); - } + if (ferror) + error->warning(FLERR,fmt::format("{} of {} force values in table are " + "inconsistent with -dE/dr.\n Should " + "only be flagged at inflection points", + ferror,tb->ninput)); // warn if re-computed distance values differ from file values - if (rerror) { - char str[128]; - sprintf(str,"%d of %d distance values in table with relative error\n" - " over %g to re-computed values",rerror,tb->ninput,EPSILONR); - error->warning(FLERR,str); - } + if (rerror) + error->warning(FLERR,fmt::format("{} of {} distance values in table with " + "relative error\n over {} to " + "re-computed values", + rerror,tb->ninput,EPSILONR)); // warn if data was read incompletely, e.g. columns were missing - if (cerror) { - char str[128]; - sprintf(str,"%d of %d lines in table were incomplete\n" - " or could not be parsed completely",cerror,tb->ninput); - error->warning(FLERR,str); - } + if (cerror) + error->warning(FLERR,fmt::format("{} of {} lines in table were " + "incomplete\n or could not be parsed " + "completely",cerror,tb->ninput)); } /* ---------------------------------------------------------------------- diff --git a/src/table_file_reader.cpp b/src/table_file_reader.cpp index 2f3317cb2a..525db24961 100644 --- a/src/table_file_reader.cpp +++ b/src/table_file_reader.cpp @@ -27,17 +27,18 @@ using namespace LAMMPS_NS; TableFileReader::TableFileReader(LAMMPS *lmp, - const std::string &filename, - const std::string &type) : - PotentialFileReader(lmp, filename, type + " table") + const std::string &filename, + const std::string &type, + const int auto_convert) : + PotentialFileReader(lmp, filename, type + " table", auto_convert) { } TableFileReader::~TableFileReader() { } -char * TableFileReader::find_section_start(const std::string & keyword) { - char * line = nullptr; +char *TableFileReader::find_section_start(const std::string & keyword) { + char *line = nullptr; while ((line = reader->next_line())) { ValueTokenizer values(line); std::string word = values.next_string(); diff --git a/src/table_file_reader.h b/src/table_file_reader.h index 209c649346..962ddf4209 100644 --- a/src/table_file_reader.h +++ b/src/table_file_reader.h @@ -24,10 +24,11 @@ namespace LAMMPS_NS { class TableFileReader : public PotentialFileReader { public: - TableFileReader(class LAMMPS *lmp, const std::string &filename, const std::string & type); + TableFileReader(class LAMMPS *lmp, const std::string &filename, + const std::string &type, const int auto_convert = 0); virtual ~TableFileReader(); - char * find_section_start(const std::string & keyword); + char *find_section_start(const std::string &keyword); }; } // namespace LAMMPS_NS diff --git a/src/utils.cpp b/src/utils.cpp index 3189b26045..df858dbb36 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -687,7 +687,7 @@ double utils::get_conversion_factor(const int property, const int conversion) } else if (conversion == METAL2REAL) { return 23.060549; } else if (conversion == REAL2METAL) { - return 0.04336410204284381954; + return 1.0/23.060549; } } return 0.0; diff --git a/unittest/formats/test_pair_unit_convert.cpp b/unittest/formats/test_pair_unit_convert.cpp index 8b96e6d46e..abc02fedd4 100644 --- a/unittest/formats/test_pair_unit_convert.cpp +++ b/unittest/formats/test_pair_unit_convert.cpp @@ -140,6 +140,10 @@ TEST_F(PairUnitConvertTest, lj_cut) lmp->input->one("read_data test_pair_unit_convert.data"); lmp->input->one("pair_style lj/cut 6.0"); lmp->input->one("pair_coeff * * 0.01014286346782117 2.0"); + remove("test.table.metal"); + lmp->input->one("pair_write 1 1 1000 r 0.1 6.0 test.table.metal lj_1_1"); + lmp->input->one("pair_write 1 2 1000 r 0.1 6.0 test.table.metal lj_1_2"); + lmp->input->one("pair_write 2 2 1000 r 0.1 6.0 test.table.metal lj_2_2"); lmp->input->one("run 0 post no"); if (!verbose) ::testing::internal::GetCapturedStdout(); @@ -158,6 +162,10 @@ TEST_F(PairUnitConvertTest, lj_cut) lmp->input->one("read_data test_pair_unit_convert.data"); lmp->input->one("pair_style lj/cut 6.0"); lmp->input->one("pair_coeff * * 0.2339 2.0"); + remove("test.table.real"); + lmp->input->one("pair_write 1 1 1000 r 0.1 6.0 test.table.real lj_1_1"); + lmp->input->one("pair_write 1 2 1000 r 0.1 6.0 test.table.real lj_1_2"); + lmp->input->one("pair_write 2 2 1000 r 0.1 6.0 test.table.real lj_2_2"); lmp->input->one("run 0 post no"); if (!verbose) ::testing::internal::GetCapturedStdout(); @@ -216,6 +224,101 @@ TEST_F(PairUnitConvertTest, sw) EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error)); } +TEST_F(PairUnitConvertTest, table_metal2real) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "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 table linear 1000"); + lmp->input->one("pair_coeff 1 1 test.table.metal lj_1_1"); + lmp->input->one("pair_coeff 1 2 test.table.metal lj_1_2"); + lmp->input->one("pair_coeff 2 2 test.table.metal lj_2_2"); + 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 table linear 1000"); + lmp->input->one("pair_coeff 1 1 test.table.metal lj_1_1"); + lmp->input->one("pair_coeff 1 2 test.table.metal lj_1_2"); + lmp->input->one("pair_coeff 2 2 test.table.metal lj_2_2"); + 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, table_real2metal) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "table")) GTEST_SKIP(); + + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("units real"); + lmp->input->one("read_data test_pair_unit_convert.data"); + lmp->input->one("pair_style table linear 1000"); + lmp->input->one("pair_coeff 1 1 test.table.real lj_1_1"); + lmp->input->one("pair_coeff 1 2 test.table.real lj_1_2"); + lmp->input->one("pair_coeff 2 2 test.table.real lj_2_2"); + 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 metal"); + lmp->input->one("read_data test_pair_unit_convert.data"); + lmp->input->one("pair_style table linear 1000"); + lmp->input->one("pair_coeff 1 1 test.table.real lj_1_1"); + lmp->input->one("pair_coeff 1 2 test.table.real lj_1_2"); + lmp->input->one("pair_coeff 2 2 test.table.real lj_2_2"); + lmp->input->one("run 0 post no"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + double pnew; + lmp->output->thermo->evaluate_keyword("press", &pnew); + EXPECT_NEAR(pold, 1.0/p_convert * pnew, fabs(pnew * rel_error)); + double enew = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul; + EXPECT_NEAR(1.0/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(1.0/ev_convert * fold[i][j], f[i][j], + fabs(f[i][j] * rel_error)); +} + TEST_F(PairUnitConvertTest, tersoff) { // check if the prerequisite pair style is available @@ -539,5 +642,7 @@ int main(int argc, char **argv) int rv = RUN_ALL_TESTS(); MPI_Finalize(); + remove("test.table.metal"); + remove("test.table.real"); return rv; }