diff --git a/doc/src/pair_tersoff_mod.rst b/doc/src/pair_tersoff_mod.rst index 7c969ad490..333f4a21bc 100644 --- a/doc/src/pair_tersoff_mod.rst +++ b/doc/src/pair_tersoff_mod.rst @@ -74,7 +74,7 @@ formulation of the V_ij term, where it contains an additional c0 term. .. math:: - V_{ij} & = f_C(r_{ij}) \left[ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) + c_0 \right] + V_{ij} = f_C(r_{ij}) \left[ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) + c_0 \right] The modified cutoff function :math:`f_C` proposed by :ref:`(Murty) ` and having a continuous second-order differential is employed. The diff --git a/src/MANYBODY/pair_tersoff_mod.cpp b/src/MANYBODY/pair_tersoff_mod.cpp index 441ba6f21e..6a60b0855c 100644 --- a/src/MANYBODY/pair_tersoff_mod.cpp +++ b/src/MANYBODY/pair_tersoff_mod.cpp @@ -53,9 +53,14 @@ void PairTersoffMOD::read_file(char *file) // open file on proc 0 if (comm->me == 0) { - PotentialFileReader reader(lmp, file, "Tersoff"); + PotentialFileReader reader(lmp, file, "TersoffMod", unit_convert_flag); char * line; + // transparently convert units for supported conversions + + int unit_convert = reader.get_unit_convert(); + double conversion_factor = utils::get_conversion_factor(utils::ENERGY, + unit_convert); while((line = reader.next_line(NPARAMS_PER_LINE))) { try { ValueTokenizer values(line); @@ -109,6 +114,11 @@ void PairTersoffMOD::read_file(char *file) params[nparams].c4 = values.next_double(); params[nparams].c5 = values.next_double(); params[nparams].powermint = int(params[nparams].powerm); + + if (unit_convert) { + params[nparams].biga *= conversion_factor; + params[nparams].bigb *= conversion_factor; + } } catch (TokenizerException & e) { error->one(FLERR, e.what()); } diff --git a/src/MANYBODY/pair_tersoff_mod_c.cpp b/src/MANYBODY/pair_tersoff_mod_c.cpp index 33e88c70b8..da14463112 100644 --- a/src/MANYBODY/pair_tersoff_mod_c.cpp +++ b/src/MANYBODY/pair_tersoff_mod_c.cpp @@ -44,9 +44,14 @@ void PairTersoffMODC::read_file(char *file) // open file on proc 0 if (comm->me == 0) { - PotentialFileReader reader(lmp, file, "Tersoff"); + PotentialFileReader reader(lmp, file, "TersoffModC", unit_convert_flag); char * line; + // transparently convert units for supported conversions + + int unit_convert = reader.get_unit_convert(); + double conversion_factor = utils::get_conversion_factor(utils::ENERGY, + unit_convert); while((line = reader.next_line(NPARAMS_PER_LINE))) { try { ValueTokenizer values(line); @@ -101,6 +106,12 @@ void PairTersoffMODC::read_file(char *file) params[nparams].c5 = values.next_double(); params[nparams].c0 = values.next_double(); params[nparams].powermint = int(params[nparams].powerm); + + if (unit_convert) { + params[nparams].biga *= conversion_factor; + params[nparams].bigb *= conversion_factor; + params[nparams].c0 *= conversion_factor; + } } catch (TokenizerException & e) { error->one(FLERR, e.what()); } diff --git a/unittest/formats/test_pair_unit_convert.cpp b/unittest/formats/test_pair_unit_convert.cpp index e77d7e535b..4cf0166f6c 100644 --- a/unittest/formats/test_pair_unit_convert.cpp +++ b/unittest/formats/test_pair_unit_convert.cpp @@ -214,6 +214,91 @@ TEST_F(PairUnitConvertTest, tersoff) EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error)); } +TEST_F(PairUnitConvertTest, tersoff_mod) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "tersoff/mod")) 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/mod"); + lmp->input->one("pair_coeff * * Si.tersoff.mod Si Si"); + 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/mod"); + lmp->input->one("pair_coeff * * Si.tersoff.mod Si Si"); + 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, tersoff_mod_c) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "tersoff/mod/c")) 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/mod/c"); + lmp->input->one("pair_coeff * * Si.tersoff.modc Si Si"); + 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/mod/c"); + lmp->input->one("pair_coeff * * Si.tersoff.modc Si Si"); + 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, tersoff_table) { // check if the prerequisite pair style is available