diff --git a/doc/src/pair_eim.rst b/doc/src/pair_eim.rst index 1cea2f0993..107edda78d 100644 --- a/doc/src/pair_eim.rst +++ b/doc/src/pair_eim.rst @@ -48,7 +48,7 @@ and :math:`sigma_i` are calculated as \sigma_i = & \sum_{j=i_1}^{i_N} q_j \cdot \psi_{ij} \left(r_{ij}\right) \\ E_i\left(q_i,\sigma_i\right) = & \frac{1}{2} \cdot q_i \cdot \sigma_i -where :math:`\eta_{ji} is a pairwise function describing electron flow from atom +where :math:`\eta_{ji}` is a pairwise function describing electron flow from atom I to atom J, and :math:`\psi_{ij}` is another pairwise function. The multi-body nature of the EIM potential is a result of the embedding energy term. A complete list of all the pair functions used in EIM is summarized @@ -63,7 +63,7 @@ below \right.\\ \eta_{ji} = & A_{\eta,ij}\left(\chi_j-\chi_i\right)f_c\left(r,r_{s,\eta,ij},r_{c,\eta,ij}\right) \\ \psi_{ij}\left(r\right) = & A_{\psi,ij}\exp\left(-\zeta_{ij}r\right)f_c\left(r,r_{s,\psi,ij},r_{c,\psi,ij}\right) \\ - f_{c}\left(r,r_p,r_c\right) = & 0.510204 \mathrm{erfc}\left[\frac{1.64498\left(2r-r_p-r_c\right)}{r_c-r_p}\right] - 0.010204 + f_{c}\left(r,r_p,r_c\right) = & 0.510204 \cdot \mathrm{erfc}\left[\frac{1.64498\left(2r-r_p-r_c\right)}{r_c-r_p}\right] - 0.010204 Here :math:`E_b, r_e, r_(c,\phi), \alpha, \beta, A_(\psi), \zeta, r_(s,\psi), r_(c,\psi), A_(\eta), r_(s,\eta), r_(c,\eta), \chi,` and pair function type diff --git a/potentials/MOH.nb3b.harmonic b/potentials/MOH.nb3b.harmonic index 6bdc263b07..acf8f79c1c 100644 --- a/potentials/MOH.nb3b.harmonic +++ b/potentials/MOH.nb3b.harmonic @@ -1,4 +1,4 @@ -# DATE: 2016-07-29 UNITS: metal CONTRIBUTOR: Todd Zeitler, tzeitle@sandia.gov CITATION: none +# DATE: 2016-07-29 UNITS: real CONTRIBUTOR: Todd Zeitler, tzeitle@sandia.gov CITATION: none # nb3b/harmonic (nonbonded 3-body harmonic) parameters for various elements # # multiple entries can be added to this file, LAMMPS reads the ones it needs diff --git a/src/MANYBODY/pair_eim.cpp b/src/MANYBODY/pair_eim.cpp index bdc08662f3..7b84de9de9 100644 --- a/src/MANYBODY/pair_eim.cpp +++ b/src/MANYBODY/pair_eim.cpp @@ -41,6 +41,7 @@ PairEIM::PairEIM(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); setfl = NULL; nmax = 0; @@ -477,7 +478,7 @@ void PairEIM::read_file(char *filename) // read potential file if( comm->me == 0) { - EIMPotentialFileReader reader(lmp, filename); + EIMPotentialFileReader reader(lmp, filename, unit_convert_flag); reader.get_global(setfl); @@ -1050,14 +1051,18 @@ double PairEIM::memory_usage() return bytes; } -EIMPotentialFileReader::EIMPotentialFileReader(LAMMPS * lmp, const std::string & filename) : +EIMPotentialFileReader::EIMPotentialFileReader(LAMMPS *lmp, + const std::string &filename, + const int auto_convert) : Pointers(lmp), filename(filename) { if (comm->me != 0) { error->one(FLERR, "EIMPotentialFileReader should only be called by proc 0!"); } - FILE * fp = force->open_potential(filename.c_str()); + int unit_convert = auto_convert; + FILE *fp = force->open_potential(filename.c_str(), &unit_convert); + conversion_factor = utils::get_conversion_factor(utils::ENERGY,unit_convert); if (fp == NULL) { error->one(FLERR, fmt::format("cannot open EIM potential file {}", filename)); @@ -1186,7 +1191,7 @@ void EIMPotentialFileReader::parse(FILE * fp) PairData data; data.rcutphiA = values.next_double(); data.rcutphiR = values.next_double(); - data.Eb = values.next_double(); + data.Eb = values.next_double() * conversion_factor; data.r0 = values.next_double(); data.alpha = values.next_double(); data.beta = values.next_double(); @@ -1194,7 +1199,7 @@ void EIMPotentialFileReader::parse(FILE * fp) data.Asigma = values.next_double(); data.rq = values.next_double(); data.rcutsigma = values.next_double(); - data.Ac = values.next_double(); + data.Ac = values.next_double() * conversion_factor; data.zeta = values.next_double(); data.rs = values.next_double(); @@ -1217,20 +1222,18 @@ void EIMPotentialFileReader::parse(FILE * fp) } } -void EIMPotentialFileReader::get_global(PairEIM::Setfl * setfl) { +void EIMPotentialFileReader::get_global(PairEIM::Setfl *setfl) { setfl->division = division; setfl->rbig = rbig; setfl->rsmall = rsmall; } -void EIMPotentialFileReader::get_element(PairEIM::Setfl * setfl, int i, const std::string & name) { - if (elements.find(name) == elements.end()) { - char str[128]; - snprintf(str, 128, "Element %s not defined in EIM potential file", name.c_str()); - error->one(FLERR, str); - } +void EIMPotentialFileReader::get_element(PairEIM::Setfl *setfl, int i, + const std::string &name) { + if (elements.find(name) == elements.end()) + error->one(FLERR,"Element " + name + " not defined in EIM potential file"); - ElementData & data = elements[name]; + ElementData &data = elements[name]; setfl->ielement[i] = data.ielement; setfl->mass[i] = data.mass; setfl->negativity[i] = data.negativity; @@ -1240,16 +1243,16 @@ void EIMPotentialFileReader::get_element(PairEIM::Setfl * setfl, int i, const st setfl->q0[i] = data.q0; } -void EIMPotentialFileReader::get_pair(PairEIM::Setfl * setfl, int ij, const std::string & elemA, const std::string & elemB) { +void EIMPotentialFileReader::get_pair(PairEIM::Setfl *setfl, int ij, + const std::string &elemA, + const std::string &elemB) { auto p = get_pair(elemA, elemB); - if (pairs.find(p) == pairs.end()) { - char str[128]; - snprintf(str, 128, "Pair (%s, %s) not defined in EIM potential file", elemA.c_str(), elemB.c_str()); - error->one(FLERR, str); - } + if (pairs.find(p) == pairs.end()) + error->one(FLERR,"Element pair (" + elemA + ", " + elemB + + ") is not defined in EIM potential file"); - PairData & data = pairs[p]; + PairData &data = pairs[p]; setfl->rcutphiA[ij] = data.rcutphiA; setfl->rcutphiR[ij] = data.rcutphiR; setfl->Eb[ij] = data.Eb; diff --git a/src/MANYBODY/pair_eim.h b/src/MANYBODY/pair_eim.h index 986db52453..1266ffd242 100644 --- a/src/MANYBODY/pair_eim.h +++ b/src/MANYBODY/pair_eim.h @@ -98,17 +98,21 @@ class EIMPotentialFileReader : protected Pointers { std::string filename; static const int MAXLINE = 1024; char line[MAXLINE]; + double conversion_factor; - void parse(FILE * fp); - char * next_line(FILE * fp); - std::pair get_pair(const std::string & a, const std::string & b); + void parse(FILE *fp); + char *next_line(FILE *fp); + std::pair get_pair(const std::string &a, + const std::string &b); public: - EIMPotentialFileReader(class LAMMPS* lmp, const std::string & filename); + EIMPotentialFileReader(class LAMMPS* lmp, const std::string &filename, + const int auto_convert=0); - void get_global(PairEIM::Setfl * setfl); - void get_element(PairEIM::Setfl * setfl, int i, const std::string & name); - void get_pair(PairEIM::Setfl * setfl, int ij, const std::string & elemA, const std::string & elemB); + void get_global(PairEIM::Setfl *setfl); + void get_element(PairEIM::Setfl *setfl, int i, const std::string &name); + void get_pair(PairEIM::Setfl *setfl, int ij, + const std::string &elemA, const std::string &elemB); private: // potential parameters diff --git a/src/MANYBODY/pair_gw.cpp b/src/MANYBODY/pair_gw.cpp index e9604523b5..32b48f7eca 100644 --- a/src/MANYBODY/pair_gw.cpp +++ b/src/MANYBODY/pair_gw.cpp @@ -48,6 +48,7 @@ PairGW::PairGW(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); nelements = 0; elements = NULL; @@ -375,9 +376,14 @@ void PairGW::read_file(char *file) // open file on proc 0 if (comm->me == 0) { - PotentialFileReader reader(lmp, file, "GW"); + PotentialFileReader reader(lmp, file, "GW", 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); @@ -427,6 +433,11 @@ void PairGW::read_file(char *file) params[nparams].lam1 = values.next_double(); params[nparams].biga = 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_gw_zbl.cpp b/src/MANYBODY/pair_gw_zbl.cpp index 96fc742b42..ddc70174f6 100644 --- a/src/MANYBODY/pair_gw_zbl.cpp +++ b/src/MANYBODY/pair_gw_zbl.cpp @@ -70,9 +70,14 @@ void PairGWZBL::read_file(char *file) // open file on proc 0 if (comm->me == 0) { - PotentialFileReader reader(lmp, file, "GW/ZBL"); + PotentialFileReader reader(lmp, file, "GW/ZBL", 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); @@ -126,6 +131,11 @@ void PairGWZBL::read_file(char *file) params[nparams].ZBLcut = values.next_double(); params[nparams].ZBLexpscale = 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_nb3b_harmonic.cpp b/src/MANYBODY/pair_nb3b_harmonic.cpp index 7491c375db..6f1bfb8905 100644 --- a/src/MANYBODY/pair_nb3b_harmonic.cpp +++ b/src/MANYBODY/pair_nb3b_harmonic.cpp @@ -47,6 +47,7 @@ PairNb3bHarmonic::PairNb3bHarmonic(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); nelements = 0; elements = NULL; @@ -291,9 +292,14 @@ void PairNb3bHarmonic::read_file(char *file) // open file on proc 0 if (comm->me == 0) { - PotentialFileReader reader(lmp, file, "nb3b/harmonic"); + PotentialFileReader reader(lmp, file, "nb3b/harmonic", 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); @@ -331,6 +337,8 @@ void PairNb3bHarmonic::read_file(char *file) params[nparams].k_theta = values.next_double(); params[nparams].theta0 = values.next_double(); params[nparams].cutoff = values.next_double(); + + if (unit_convert) params[nparams].k_theta *= 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 b9e5431f6d..6c4b38397d 100644 --- a/unittest/formats/test_pair_unit_convert.cpp +++ b/unittest/formats/test_pair_unit_convert.cpp @@ -353,6 +353,178 @@ TEST_F(PairUnitConvertTest, eam_cd) EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error)); } +TEST_F(PairUnitConvertTest, eim) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "eim")) 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 eim"); + lmp->input->one("pair_coeff * * Na Cl ffield.eim Na Cl"); + 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 eim"); + lmp->input->one("pair_coeff * * Na Cl ffield.eim Na Cl"); + 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, gw) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "gw")) 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 gw"); + lmp->input->one("pair_coeff * * SiC.gw 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 gw"); + lmp->input->one("pair_coeff * * SiC.gw 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, gw_zbl) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "gw/zbl")) 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 gw/zbl"); + lmp->input->one("pair_coeff * * SiC.gw.zbl 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 gw/zbl"); + lmp->input->one("pair_coeff * * SiC.gw.zbl 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, nb3b_harmonic) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "nb3b/harmonic")) 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 nb3b/harmonic"); + lmp->input->one("pair_coeff * * MOH.nb3b.harmonic M O"); + 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 nb3b/harmonic"); + lmp->input->one("pair_coeff * * MOH.nb3b.harmonic M O"); + 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