diff --git a/src/MANYBODY/pair_eam_cd.cpp b/src/MANYBODY/pair_eam_cd.cpp index 4eeeb173ec..c8ebea9b61 100644 --- a/src/MANYBODY/pair_eam_cd.cpp +++ b/src/MANYBODY/pair_eam_cd.cpp @@ -29,6 +29,8 @@ #include "memory.h" #include "error.h" #include "tokenizer.h" +#include "utils.h" +#include "fmt/format.h" using namespace LAMMPS_NS; @@ -40,6 +42,7 @@ PairEAMCD::PairEAMCD(LAMMPS *lmp, int _cdeamVersion) { single_enable = 0; restartinfo = 0; + unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); rhoB = NULL; D_values = NULL; @@ -500,12 +503,11 @@ void PairEAMCD::read_h_coeff(char *filename) FILE *fptr; char line[MAXLINE]; char nextline[MAXLINE]; - fptr = force->open_potential(filename); - if (fptr == NULL) { - char str[128]; - snprintf(str,128,"Cannot open EAM potential file %s", filename); - error->one(FLERR,str); - } + int convert_flag = unit_convert_flag; + fptr = force->open_potential(filename, &convert_flag); + if (fptr == NULL) + error->one(FLERR,fmt::format("Cannot open EAMCD potential file {}", + filename)); // h coefficients are stored at the end of the file. // Skip to last line of file. diff --git a/src/force.cpp b/src/force.cpp index 5824e35648..0a1a4ec4fd 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -1011,7 +1011,7 @@ tagint Force::tnumeric(const char *file, int line, char *str) if fails, search in dir specified by env variable LAMMPS_POTENTIALS ------------------------------------------------------------------------- */ -FILE *Force::open_potential(const char *name) +FILE *Force::open_potential(const char *name, int *auto_convert) { std::string filepath = utils::get_potential_file_path(name); @@ -1024,9 +1024,35 @@ FILE *Force::open_potential(const char *name) utils::logmesg(lmp, fmt::format("Reading potential file {} " "with DATE: {}\n", name, date)); } - if (!units.empty() && (units != unit_style)) { - error->one(FLERR, fmt::format("Potential file {} requires {} units " - "but {} units are in use", name, units, unit_style)); + + if (auto_convert == nullptr) { + if (units != unit_style) { + error->one(FLERR, fmt::format("Potential file {} requires {} units " + "but {} units are in use", name, units, + unit_style)); + return nullptr; + } + } else { + if (units == unit_style) { + *auto_convert = utils::NOCONVERT; + } else { + if ((units == "metal") && (unit_style == "real") + && (*auto_convert & utils::METAL2REAL)) { + *auto_convert = utils::METAL2REAL; + } else if ((units == "real") && (unit_style == "metal") + && (*auto_convert & utils::REAL2METAL)) { + *auto_convert = utils::REAL2METAL; + } else { + error->one(FLERR, fmt::format("Potential file {} requires {} units " + "but {} units are in use", name, + units, unit_style)); + return nullptr; + } + } + if (*auto_convert != utils::NOCONVERT) + lmp->error->warning(FLERR, fmt::format("Converting potential file in " + "{} units to {} units", + units, unit_style)); } return fopen(filepath.c_str(), "r"); } diff --git a/src/force.h b/src/force.h index fa0c0f6877..b566f02fb7 100644 --- a/src/force.h +++ b/src/force.h @@ -135,7 +135,7 @@ class Force : protected Pointers { bigint bnumeric(const char *, int, char *); tagint tnumeric(const char *, int, char *); - FILE *open_potential(const char *); + FILE *open_potential(const char *, int *auto_convert = nullptr); bigint memory_usage(); diff --git a/unittest/formats/test_pair_unit_convert.cpp b/unittest/formats/test_pair_unit_convert.cpp index f397846ffd..b9e5431f6d 100644 --- a/unittest/formats/test_pair_unit_convert.cpp +++ b/unittest/formats/test_pair_unit_convert.cpp @@ -224,6 +224,135 @@ TEST_F(PairUnitConvertTest, eam) EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error)); } +TEST_F(PairUnitConvertTest, eam_alloy) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "eam/alloy")) 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 eam/alloy"); + lmp->input->one("pair_coeff * * AlCu.eam.alloy Al Cu"); + 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 eam/alloy"); + lmp->input->one("pair_coeff * * AlCu.eam.alloy Al Cu"); + 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, eam_fs) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "eam/fs")) 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 eam/fs"); + lmp->input->one("pair_coeff * * FeP_mm.eam.fs Fe P"); + 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 eam/fs"); + lmp->input->one("pair_coeff * * FeP_mm.eam.fs Fe P"); + 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, eam_cd) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "eam/cd")) 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 eam/cd"); + lmp->input->one("pair_coeff * * FeCr.cdeam Cr Fe"); + 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 eam/cd"); + lmp->input->one("pair_coeff * * FeCr.cdeam Cr Fe"); + 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