diff --git a/src/MANYBODY/pair_tersoff.cpp b/src/MANYBODY/pair_tersoff.cpp index 4ba7471e60..1f2593e188 100644 --- a/src/MANYBODY/pair_tersoff.cpp +++ b/src/MANYBODY/pair_tersoff.cpp @@ -49,6 +49,7 @@ PairTersoff::PairTersoff(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; @@ -400,8 +401,16 @@ void PairTersoff::read_file(char *file) // open file on proc 0 if (comm->me == 0) { - PotentialFileReader reader(lmp, file, "Tersoff"); - char * line; + PotentialFileReader reader(lmp, file, "Tersoff", 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); + printf("unit_conver=%d\n",unit_convert); + printf("conversion_factor=%g\n",conversion_factor); while((line = reader.next_line(NPARAMS_PER_LINE))) { try { @@ -453,6 +462,11 @@ void PairTersoff::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/pair.cpp b/src/pair.cpp index c18b84a8ed..432680b324 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -66,6 +66,7 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp) no_virial_fdotr_compute = 0; writedata = 0; ghostneigh = 0; + unit_convert_flag = utils::NOCONVERT; nextra = 0; pvector = NULL; diff --git a/src/pair.h b/src/pair.h index d0aee8cd7f..0664e97462 100644 --- a/src/pair.h +++ b/src/pair.h @@ -53,6 +53,7 @@ class Pair : protected Pointers { int respa_enable; // 1 if inner/middle/outer rRESPA routines int one_coeff; // 1 if allows only one coeff * * call int manybody_flag; // 1 if a manybody potential + int unit_convert_flag; // value != 0 indicates support for unit conversion. int no_virial_fdotr_compute; // 1 if does not invoke virial_fdotr_compute() int writedata; // 1 if writes coeffs to data file int ghostneigh; // 1 if pair style needs neighbors of ghosts diff --git a/src/potential_file_reader.cpp b/src/potential_file_reader.cpp index 1cc77351a9..924af94eb3 100644 --- a/src/potential_file_reader.cpp +++ b/src/potential_file_reader.cpp @@ -31,11 +31,13 @@ using namespace LAMMPS_NS; PotentialFileReader::PotentialFileReader(LAMMPS *lmp, const std::string &filename, - const std::string &potential_name) : + const std::string &potential_name, + const int auto_convert) : Pointers(lmp), reader(nullptr), filename(filename), - filetype(potential_name + " potential") + filetype(potential_name + " potential"), + unit_convert(auto_convert) { if (comm->me != 0) { error->one(FLERR, "FileReader should only be called by proc 0!"); @@ -155,9 +157,21 @@ TextFileReader * PotentialFileReader::open_potential(const std::string& path) { utils::logmesg(lmp, fmt::format("Reading potential file {} with DATE: {}\n", filename, date)); } - if (!units.empty() && (units != unit_style)) { - lmp->error->one(FLERR, fmt::format("Potential file {} requires {} units " - "but {} units are in use",filename, units, unit_style)); + if (units.empty()) { + unit_convert == utils::NOCONVERT; + } else { + if (units == unit_style) { + unit_convert = utils::NOCONVERT; + } else { + if ((units == "metal") && (unit_style == "real") && (unit_convert | utils::METAL2REAL)) { + unit_convert = utils::METAL2REAL; + } else if ((units == "real") && (unit_style == "metal") && (unit_convert | utils::REAL2METAL)) { + unit_convert = utils::REAL2METAL; + } else { + lmp->error->one(FLERR, fmt::format("Potential file {} requires {} units " + "but {} units are in use",filename, units, unit_style)); + } + } } return new TextFileReader(filepath, filetype); diff --git a/src/potential_file_reader.h b/src/potential_file_reader.h index a73f5fdbaa..9094d3957a 100644 --- a/src/potential_file_reader.h +++ b/src/potential_file_reader.h @@ -28,22 +28,25 @@ namespace LAMMPS_NS { class PotentialFileReader : protected Pointers { protected: - TextFileReader * reader; + TextFileReader *reader; std::string filename; std::string filetype; + int unit_convert; - TextFileReader * open_potential(const std::string& path); + TextFileReader *open_potential(const std::string& path); public: - PotentialFileReader(class LAMMPS *lmp, const std::string &filename, const std::string &potential_name); + PotentialFileReader(class LAMMPS *lmp, const std::string &filename, + const std::string &potential_name, + const int auto_convert = 0); virtual ~PotentialFileReader(); void ignore_comments(bool value); void skip_line(); - char * next_line(int nparams = 0); - void next_dvector(double * list, int n); - ValueTokenizer next_values(int nparams, const std::string & separators = TOKENIZER_DEFAULT_SEPARATORS); + char *next_line(int nparams = 0); + void next_dvector(double *list, int n); + ValueTokenizer next_values(int nparams, const std::string &separators = TOKENIZER_DEFAULT_SEPARATORS); // convenience functions double next_double(); @@ -51,6 +54,9 @@ namespace LAMMPS_NS tagint next_tagint(); bigint next_bigint(); std::string next_string(); + + // unit conversion info + int get_unit_convert() const { return unit_convert; } }; } // namespace LAMMPS_NS diff --git a/src/utils.cpp b/src/utils.cpp index 1e1ded3c0b..42d69d26a2 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -583,6 +583,36 @@ std::string utils::get_potential_units(const std::string & path, const std::stri return ""; } +/* ---------------------------------------------------------------------- + return bitmask of supported conversions for a given property +------------------------------------------------------------------------- */ +int utils::get_supported_conversions(const int property) +{ + if (property == ENERGY) { + return METAL2REAL | REAL2METAL; + } + return NOCONVERT; +} + +/* ---------------------------------------------------------------------- + return conversion factor for a given property and conversion setting + return 0.0 if unknown. +------------------------------------------------------------------------- */ + +double utils::get_conversion_factor(const int property, const int conversion) +{ + if (property == ENERGY) { + if (conversion == NOCONVERT) { + return 1.0; + } else if (conversion == METAL2REAL) { + return 23.060549; + } else if (conversion == REAL2METAL) { + return 0.04336410204284381954; + } + } + return 0.0; +} + /* ------------------------------------------------------------------ */ extern "C" { diff --git a/src/utils.h b/src/utils.h index 91e5024121..c83764391f 100644 --- a/src/utils.h +++ b/src/utils.h @@ -240,6 +240,24 @@ namespace LAMMPS_NS { * \return UNITS field if present */ std::string get_potential_units(const std::string & path, const std::string & potential_name); + + enum { NOCONVERT = 0, METAL2REAL = 1, REAL2METAL = 1<<1 }; + enum { UNKNOWN = 0, ENERGY }; + + /** + * \brief Return bitmask of available conversion factors for a given propert + * \param property property to be converted + * \return bitmask indicating available conversions + */ + int get_supported_conversions(const int property); + + /** + * \brief Return unit conversion factor for given property and selected from/to units + * \param property property to be converted + * \param conversion constant indicating the conversion + * \return conversion factor + */ + double get_conversion_factor(const int property, const int conversion); } } diff --git a/unittest/formats/CMakeLists.txt b/unittest/formats/CMakeLists.txt index aaf6731855..bf4300ff2e 100644 --- a/unittest/formats/CMakeLists.txt +++ b/unittest/formats/CMakeLists.txt @@ -3,6 +3,11 @@ add_executable(test_atom_styles test_atom_styles.cpp) target_link_libraries(test_atom_styles PRIVATE lammps GTest::GMock GTest::GTest) add_test(NAME AtomStyles COMMAND test_atom_styles WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +add_executable(test_pair_unit_convert test_pair_unit_convert.cpp) +target_link_libraries(test_pair_unit_convert PRIVATE lammps GTest::GMock GTest::GTest) +add_test(NAME PairUnitConvert COMMAND test_pair_unit_convert WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +set_tests_properties(PairUnitConvert PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}") + add_executable(test_potential_file_reader test_potential_file_reader.cpp) target_link_libraries(test_potential_file_reader PRIVATE lammps GTest::GMock GTest::GTest) add_test(NAME PotentialFileReader COMMAND test_potential_file_reader WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/unittest/formats/test_pair_unit_convert.cpp b/unittest/formats/test_pair_unit_convert.cpp new file mode 100644 index 0000000000..0dc4f1ec6b --- /dev/null +++ b/unittest/formats/test_pair_unit_convert.cpp @@ -0,0 +1,206 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "atom.h" +#include "force.h" +#include "info.h" +#include "input.h" +#include "lammps.h" +#include "pair.h" +#include "utils.h" +#include "gmock/gmock.h" +#include "gtest/gtest.h" + +#include +#include +#include + +// whether to print verbose output (i.e. not capturing LAMMPS screen output). +bool verbose = false; + +namespace LAMMPS_NS { +using ::testing::Eq; + +// eV to kcal/mol conversion constant (CODATA 2018) +const double ev_convert = utils::get_conversion_factor(utils::ENERGY, utils::METAL2REAL); +const double rel_error = 1.0e-13; + +class PairUnitConvertTest : public ::testing::Test { +protected: + LAMMPS *lmp; + Info *info; + double fold[4][3]; + + void SetUp() override + { + const char *args[] = {"PairUnitConvertTest", "-log", "none", "-echo", "screen", "-nocite"}; + char **argv = (char **)args; + int argc = sizeof(args) / sizeof(char *); + if (!verbose) ::testing::internal::CaptureStdout(); + lmp = new LAMMPS(argc, argv, MPI_COMM_WORLD); + if (!verbose) ::testing::internal::GetCapturedStdout(); + ASSERT_NE(lmp, nullptr); + if (!verbose) ::testing::internal::CaptureStdout(); + info = new Info(lmp); + lmp->input->one("units metal"); + lmp->input->one("dimension 3"); + lmp->input->one("region box block -4 4 -4 4 -4 4"); + lmp->input->one("create_box 2 box"); + lmp->input->one("create_atoms 1 single -1.1 1.2 0.0 units box"); + lmp->input->one("create_atoms 1 single -1.2 -1.1 0.0 units box"); + lmp->input->one("create_atoms 2 single 0.9 1.0 0.0 units box"); + lmp->input->one("create_atoms 2 single 1.0 -0.9 0.0 units box"); + lmp->input->one("pair_style zero 4.0"); + lmp->input->one("pair_coeff * *"); + lmp->input->one("mass * 1.0"); + lmp->input->one("write_data test_pair_unit_convert.data nocoeff"); + lmp->input->one("clear"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + } + + void TearDown() override + { + if (!verbose) ::testing::internal::CaptureStdout(); + delete info; + delete lmp; + if (!verbose) ::testing::internal::GetCapturedStdout(); + remove("test_pair_unit_convert.data"); + } +}; + +TEST_F(PairUnitConvertTest, zero) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "zero")) 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 zero 6.0"); + lmp->input->one("pair_coeff * *"); + lmp->input->one("run 0 post no"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + // copy energy and force from first step + 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 zero 6.0"); + lmp->input->one("pair_coeff * *"); + lmp->input->one("run 0 post no"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + 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, lj_cut) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "lj/cut")) 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 lj/cut 6.0"); + lmp->input->one("pair_coeff * * 0.01014286346782117 3.504"); + lmp->input->one("run 0 post no"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + // copy energy and force from first step + 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 lj/cut 6.0"); + lmp->input->one("pair_coeff * * 0.2339 3.504"); + lmp->input->one("run 0 post no"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + 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) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "tersoff")) 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"); + lmp->input->one("pair_coeff * * SiC.tersoff Si C"); + lmp->input->one("run 0 post no"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + // copy energy and force from first step + 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"); + lmp->input->one("pair_coeff * * SiC.tersoff Si C"); + lmp->input->one("run 0 post no"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + 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)); +} + +} // namespace LAMMPS_NS + +int main(int argc, char **argv) +{ + MPI_Init(&argc, &argv); + ::testing::InitGoogleMock(&argc, argv); + if ((argc > 1) && (strcmp(argv[1], "-v") == 0)) verbose = true; + + int rv = RUN_ALL_TESTS(); + MPI_Finalize(); + return rv; +}