diff --git a/cmake/Modules/Packages/KIM.cmake b/cmake/Modules/Packages/KIM.cmake index 2a2a1cde78..995d2d9490 100644 --- a/cmake/Modules/Packages/KIM.cmake +++ b/cmake/Modules/Packages/KIM.cmake @@ -19,7 +19,7 @@ if(CURL_FOUND) target_compile_definitions(lammps PRIVATE -DLMP_NO_SSL_CHECK) endif() endif() -set(KIM_EXTRA_UNITTESTS OFF CACHE STRING "Set extra unit tests verbose mode on/off. If on, extra tests are included.") +option(KIM_EXTRA_UNITTESTS "Enable extra unit tests for the KIM package." OFF) mark_as_advanced(KIM_EXTRA_UNITTESTS) find_package(PkgConfig QUIET) set(DOWNLOAD_KIM_DEFAULT ON) diff --git a/src/.gitignore b/src/.gitignore index ac4a776cfc..7c7155018a 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -56,6 +56,8 @@ /pair_lepton.cpp /pair_lepton.h +/pair_lepton_coul.cpp +/pair_lepton_coul.h /bond_lepton.cpp /bond_lepton.h /angle_lepton.cpp diff --git a/src/KIM/fix_store_kim.cpp b/src/KIM/fix_store_kim.cpp index 415ba4c410..44c2ccaa1d 100644 --- a/src/KIM/fix_store_kim.cpp +++ b/src/KIM/fix_store_kim.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -69,11 +68,11 @@ using namespace FixConst; /* ---------------------------------------------------------------------- */ -FixStoreKIM::FixStoreKIM(LAMMPS *lmp, int narg, char **arg) - : Fix(lmp, narg, arg), simulator_model(nullptr), model_name(nullptr), - model_units(nullptr), user_units(nullptr) +FixStoreKIM::FixStoreKIM(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), simulator_model(nullptr), model_name(nullptr), model_units(nullptr), + user_units(nullptr) { - if (narg != 3) error->all(FLERR,"Illegal fix STORE/KIM command"); + if (narg != 3) error->all(FLERR, "Illegal fix STORE/KIM command"); } /* ---------------------------------------------------------------------- */ @@ -83,25 +82,25 @@ FixStoreKIM::~FixStoreKIM() // free associated storage if (simulator_model) { - auto sm = (KIM_SimulatorModel *)simulator_model; + auto sm = (KIM_SimulatorModel *) simulator_model; KIM_SimulatorModel_Destroy(&sm); simulator_model = nullptr; } if (model_name) { - auto mn = (char *)model_name; + auto mn = (char *) model_name; delete[] mn; model_name = nullptr; } if (model_units) { - auto mu = (char *)model_units; + auto mu = (char *) model_units; delete[] mu; model_units = nullptr; } if (user_units) { - auto uu = (char *)user_units; + auto uu = (char *) user_units; delete[] uu; user_units = nullptr; } @@ -121,38 +120,44 @@ void FixStoreKIM::setptr(const std::string &name, void *ptr) { if (name == "simulator_model") { if (simulator_model) { - auto sm = (KIM_SimulatorModel *)simulator_model; + auto sm = (KIM_SimulatorModel *) simulator_model; KIM_SimulatorModel_Destroy(&sm); } simulator_model = ptr; } else if (name == "model_name") { if (model_name) { - auto mn = (char *)model_name; + auto mn = (char *) model_name; delete[] mn; } model_name = ptr; } else if (name == "model_units") { if (model_units) { - auto mu = (char *)model_units; + auto mu = (char *) model_units; delete[] mu; } model_units = ptr; } else if (name == "user_units") { if (user_units) { - auto uu = (char *)user_units; + auto uu = (char *) user_units; delete[] uu; } user_units = ptr; - } else error->all(FLERR,"Unknown property in fix STORE/KIM"); + } else + error->all(FLERR, "Unknown property in fix STORE/KIM"); } /* ---------------------------------------------------------------------- */ void *FixStoreKIM::getptr(const std::string &name) { - if (name == "simulator_model") return simulator_model; - else if (name == "model_name") return model_name; - else if (name == "model_units") return model_units; - else if (name == "user_units") return user_units; - else return nullptr; + if (name == "simulator_model") + return simulator_model; + else if (name == "model_name") + return model_name; + else if (name == "model_units") + return model_units; + else if (name == "user_units") + return user_units; + else + return nullptr; } diff --git a/src/KIM/kim_command.cpp b/src/KIM/kim_command.cpp index 820e4dcba7..57871e964b 100644 --- a/src/KIM/kim_command.cpp +++ b/src/KIM/kim_command.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -70,40 +69,40 @@ using namespace LAMMPS_NS; static constexpr const char *const cite_openkim = - "OpenKIM Project: doi:10.1007/s11837-011-0102-6\n\n" - "@Article{tadmor:elliott:2011,\n" - " author = {E. B. Tadmor and R. S. Elliott and J. P. Sethna and R. E. Miller " - "and C. A. Becker},\n" - " title = {The potential of atomistic simulations and the {K}nowledgebase of " - "{I}nteratomic {M}odels},\n" - " journal = {{JOM}},\n" - " year = 2011,\n" - " volume = 63,\n" - " number = 17,\n" - " pages = {17},\n" - " doi = {10.1007/s11837-011-0102-6}\n" - "}\n\n"; + "OpenKIM Project: doi:10.1007/s11837-011-0102-6\n\n" + "@Article{tadmor:elliott:2011,\n" + " author = {E. B. Tadmor and R. S. Elliott and J. P. Sethna and R. E. Miller " + "and C. A. Becker},\n" + " title = {The potential of atomistic simulations and the {K}nowledgebase of " + "{I}nteratomic {M}odels},\n" + " journal = {{JOM}},\n" + " year = 2011,\n" + " volume = 63,\n" + " number = 17,\n" + " pages = {17},\n" + " doi = {10.1007/s11837-011-0102-6}\n" + "}\n\n"; static constexpr const char *const cite_openkim_query = - "OpenKIM query: doi:10.1063/5.0014267\n\n" - "@Article{karls:bierbaum:2020,\n" - " author = {D. S. Karls and M. Bierbaum and A. A. Alemi and R. S. Elliott " - "and J. P. Sethna and E. B. Tadmor},\n" - " title = {The {O}pen{KIM} processing pipeline: {A} cloud-based automatic " - "material property computation engine},\n" - " journal = {{T}he {J}ournal of {C}hemical {P}hysics},\n" - " year = 2020,\n" - " volume = 153,\n" - " number = 6,\n" - " pages = {064104},\n" - " doi = {10.1063/5.0014267}\n" - "}\n\n"; + "OpenKIM query: doi:10.1063/5.0014267\n\n" + "@Article{karls:bierbaum:2020,\n" + " author = {D. S. Karls and M. Bierbaum and A. A. Alemi and R. S. Elliott " + "and J. P. Sethna and E. B. Tadmor},\n" + " title = {The {O}pen{KIM} processing pipeline: {A} cloud-based automatic " + "material property computation engine},\n" + " journal = {{T}he {J}ournal of {C}hemical {P}hysics},\n" + " year = 2020,\n" + " volume = 153,\n" + " number = 6,\n" + " pages = {064104},\n" + " doi = {10.1063/5.0014267}\n" + "}\n\n"; /* ---------------------------------------------------------------------- */ void KimCommand::command(int narg, char **arg) { - if (narg < 1) error->all(FLERR,"Illegal kim command"); + if (narg < 1) error->all(FLERR, "Illegal kim command"); const std::string subcmd(arg[0]); narg--; @@ -132,5 +131,6 @@ void KimCommand::command(int narg, char **arg) auto cmd = new KimQuery(lmp); cmd->command(narg, arg); delete cmd; - } else error->all(FLERR,"Unknown kim subcommand {}", subcmd); + } else + error->all(FLERR, "Unknown kim subcommand {}", subcmd); } diff --git a/src/KIM/kim_init.cpp b/src/KIM/kim_init.cpp index 037fb2857e..f7dd01be49 100644 --- a/src/KIM/kim_init.cpp +++ b/src/KIM/kim_init.cpp @@ -76,6 +76,8 @@ extern "C" { } using namespace LAMMPS_NS; +using kim_units::get_kim_unit_names; +using kim_units::lammps_unit_conversion; /* ---------------------------------------------------------------------- */ @@ -96,7 +98,7 @@ void KimInit::command(int narg, char **arg) error->all(FLERR, "Illegal 'kim init' command.\n" "The argument followed by unit_style {} is an optional argument and when " - "is used must be unit_conversion_mode", + "used must be unit_conversion_mode", user_units); } } else @@ -109,58 +111,10 @@ void KimInit::command(int narg, char **arg) if (universe->nprocs > 1) MPI_Barrier(universe->uworld); determine_model_type_and_units(model_name, user_units, &model_units, pkim); - write_log_cite(lmp, model_type, model_name); - do_init(model_name, user_units, model_units, pkim); } -/* ---------------------------------------------------------------------- */ - -namespace { -void get_kim_unit_names(char const *const system, KIM_LengthUnit &lengthUnit, - KIM_EnergyUnit &energyUnit, KIM_ChargeUnit &chargeUnit, - KIM_TemperatureUnit &temperatureUnit, KIM_TimeUnit &timeUnit, Error *error) -{ - const std::string system_str(system); - if (system_str == "real") { - lengthUnit = KIM_LENGTH_UNIT_A; - energyUnit = KIM_ENERGY_UNIT_kcal_mol; - chargeUnit = KIM_CHARGE_UNIT_e; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_fs; - } else if (system_str == "metal") { - lengthUnit = KIM_LENGTH_UNIT_A; - energyUnit = KIM_ENERGY_UNIT_eV; - chargeUnit = KIM_CHARGE_UNIT_e; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_ps; - } else if (system_str == "si") { - lengthUnit = KIM_LENGTH_UNIT_m; - energyUnit = KIM_ENERGY_UNIT_J; - chargeUnit = KIM_CHARGE_UNIT_C; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_s; - } else if (system_str == "cgs") { - lengthUnit = KIM_LENGTH_UNIT_cm; - energyUnit = KIM_ENERGY_UNIT_erg; - chargeUnit = KIM_CHARGE_UNIT_statC; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_s; - } else if (system_str == "electron") { - lengthUnit = KIM_LENGTH_UNIT_Bohr; - energyUnit = KIM_ENERGY_UNIT_Hartree; - chargeUnit = KIM_CHARGE_UNIT_e; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_fs; - } else if ((system_str == "lj") || (system_str == "micro") || (system_str == "nano")) { - error->all(FLERR, "LAMMPS unit_style {} not supported by KIM models", system_str); - } else { - error->all(FLERR, "Unknown unit_style"); - } -} -} // namespace - void KimInit::print_dirs(struct KIM_Collections *const collections) const { int kim_error = 0; @@ -312,12 +266,8 @@ void KimInit::do_init(char *model_name, char *user_units, char *model_units, KIM { // create storage proxy fix. delete existing fix, if needed. - int ifix = modify->find_fix("KIM_MODEL_STORE"); - if (ifix >= 0) modify->delete_fix(ifix); - modify->add_fix("KIM_MODEL_STORE all STORE/KIM"); - ifix = modify->find_fix("KIM_MODEL_STORE"); - - auto fix_store = dynamic_cast(modify->fix[ifix]); + if (modify->get_fix_by_id("KIM_MODEL_STORE")) modify->delete_fix("KIM_MODEL_STORE"); + auto fix_store = dynamic_cast(modify->add_fix("KIM_MODEL_STORE all STORE/KIM")); fix_store->setptr("model_name", (void *) model_name); fix_store->setptr("user_units", (void *) user_units); fix_store->setptr("model_units", (void *) model_units); @@ -471,9 +421,7 @@ void KimInit::do_variables(const std::string &from, const std::string &to) } ier = lammps_unit_conversion(units[i], from, to, conversion_factor); if (ier != 0) - error->all(FLERR, - "Unable to obtain conversion factor: " - "unit = {}; from = {}; to = {}", + error->all(FLERR, "Unable to obtain conversion factor: unit = {}; from = {}; to = {}", units[i], from, to); variable->internal_set(v_unit, conversion_factor); diff --git a/src/KIM/kim_interactions.cpp b/src/KIM/kim_interactions.cpp index 5a91109e53..1f4f84e648 100644 --- a/src/KIM/kim_interactions.cpp +++ b/src/KIM/kim_interactions.cpp @@ -85,11 +85,10 @@ using namespace LAMMPS_NS; void KimInteractions::command(int narg, char **arg) { - if (narg < 1) error->all(FLERR, "Illegal 'kim interactions' command"); + if (narg < 1) utils::missing_cmd_args(FLERR, "kim interactions", error); if (!domain->box_exist) - error->all(FLERR, "Must use 'kim interactions' command after " - "simulation box is defined"); + error->all(FLERR, "Use of 'kim interactions' before simulation box is defined"); do_setup(narg, arg); } @@ -103,11 +102,9 @@ void KimInteractions::do_setup(int narg, char **arg) if ((narg == 1) && (arg_str == "fixed_types")) { fixed_types = true; } else if (narg != atom->ntypes) { - error->all(FLERR, "Illegal 'kim interactions' command.\nThe " - "LAMMPS simulation has {} atom type(s), but " - "{} chemical species passed to the " - "'kim interactions' command", - atom->ntypes, narg); + error->all(FLERR, "Illegal 'kim interactions' command.\nThe LAMMPS simulation has {} atom " + "type(s), but {} chemical species passed to the 'kim interactions' command", + atom->ntypes, narg); } else { fixed_types = false; } @@ -119,16 +116,14 @@ void KimInteractions::do_setup(int narg, char **arg) // retrieve model name and pointer to simulator model class instance. // validate model name if not given as null pointer. - int ifix = modify->find_fix("KIM_MODEL_STORE"); - if (ifix >= 0) { - auto fix_store = dynamic_cast(modify->fix[ifix]); + auto fix_store = dynamic_cast(modify->get_fix_by_id("KIM_MODEL_STORE")); + if (fix_store) { model_name = (char *)fix_store->getptr("model_name"); simulatorModel = (KIM_SimulatorModel *)fix_store->getptr("simulator_model"); } else error->all(FLERR, "Must use 'kim init' before 'kim interactions'"); // Begin output to log file - input->write_echo("#=== BEGIN kim interactions ===========================" - "=======\n"); + input->write_echo("#=== BEGIN kim interactions ==================================\n"); if (simulatorModel) { auto first_visit = input->variable->find("kim_update"); @@ -165,8 +160,8 @@ void KimInteractions::do_setup(int narg, char **arg) if (atom_type_sym == sim_species) species_is_supported = true; } if (!species_is_supported) { - error->all(FLERR, "Species '{}' is not supported by this " - "KIM Simulator Model", atom_type_sym); + error->all(FLERR, "Species '{}' is not supported by this KIM Simulator Model", + atom_type_sym); } } } else { @@ -179,18 +174,15 @@ void KimInteractions::do_setup(int narg, char **arg) const char *sim_field, *sim_value; KIM_SimulatorModel_GetNumberOfSimulatorFields(simulatorModel, &sim_fields); for (int i = 0; i < sim_fields; ++i) { - KIM_SimulatorModel_GetSimulatorFieldMetadata( - simulatorModel, i, &sim_lines, &sim_field); + KIM_SimulatorModel_GetSimulatorFieldMetadata(simulatorModel, i, &sim_lines, &sim_field); const std::string sim_field_str(sim_field); if (sim_field_str == "units") { - KIM_SimulatorModel_GetSimulatorFieldLine( - simulatorModel, i, 0, &sim_value); + KIM_SimulatorModel_GetSimulatorFieldLine(simulatorModel, i, 0, &sim_value); - const std::string sim_value_str(sim_value); - const std::string unit_style_str(update->unit_style); - if (sim_value_str != unit_style_str) - error->all(FLERR, "Incompatible units for KIM Simulator Model"); + if (strcmp(sim_value, update->unit_style) != 0) + error->all(FLERR, "Incompatible units for KIM Simulator Model: {} vs {}", + sim_value, update->unit_style); } } @@ -217,8 +209,7 @@ void KimInteractions::do_setup(int narg, char **arg) no_model_definition = false; for (int j = 0; j < sim_lines; ++j) { - KIM_SimulatorModel_GetSimulatorFieldLine( - simulatorModel, i, j, &sim_value); + KIM_SimulatorModel_GetSimulatorFieldLine(simulatorModel, i, j, &sim_value); if (utils::strmatch(sim_value, "^KIM_SET_TYPE_PARAMETERS")) { // Notes regarding the KIM_SET_TYPE_PARAMETERS command // * This is an INTERNAL command. @@ -264,8 +255,7 @@ void KimInteractions::do_setup(int narg, char **arg) } // End output to log file - input->write_echo("#=== END kim interactions =============================" - "=======\n\n"); + input->write_echo("#=== END kim interactions ====================================\n\n"); } /* ---------------------------------------------------------------------- */ @@ -276,8 +266,7 @@ void KimInteractions::KIM_SET_TYPE_PARAMETERS(const std::string &input_line) con const std::string key = words[1]; if (key != "pair" && key != "charge") - error->one(FLERR, "Unrecognized KEY {} for " - "KIM_SET_TYPE_PARAMETERS command", key); + error->one(FLERR, "Unrecognized KEY {} for KIM_SET_TYPE_PARAMETERS command", key); std::string filename = words[2]; std::vector species(words.begin() + 3, words.end()); @@ -287,7 +276,7 @@ void KimInteractions::KIM_SET_TYPE_PARAMETERS(const std::string &input_line) con FILE *fp = nullptr; if (comm->me == 0) { fp = fopen(filename.c_str(), "r"); - if (fp == nullptr) error->one(FLERR, "Parameter file not found"); + if (fp == nullptr) error->one(FLERR, "Parameter file {} not found", filename); } char line[MAXLINE], *ptr; diff --git a/src/KIM/kim_param.cpp b/src/KIM/kim_param.cpp index 6473cb2ad0..f72df81989 100644 --- a/src/KIM/kim_param.cpp +++ b/src/KIM/kim_param.cpp @@ -63,6 +63,7 @@ #include "fix_store_kim.h" #include "force.h" #include "input.h" +#include "kim_units.h" #include "modify.h" #include "pair_kim.h" #include "variable.h" @@ -77,61 +78,7 @@ extern "C" } using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -namespace -{ -void get_kim_unit_names( - char const *const system, - KIM_LengthUnit &lengthUnit, - KIM_EnergyUnit &energyUnit, - KIM_ChargeUnit &chargeUnit, - KIM_TemperatureUnit &temperatureUnit, - KIM_TimeUnit &timeUnit, - Error *error) -{ - const std::string system_str(system); - if (system_str == "real") { - lengthUnit = KIM_LENGTH_UNIT_A; - energyUnit = KIM_ENERGY_UNIT_kcal_mol; - chargeUnit = KIM_CHARGE_UNIT_e; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_fs; - } else if (system_str == "metal") { - lengthUnit = KIM_LENGTH_UNIT_A; - energyUnit = KIM_ENERGY_UNIT_eV; - chargeUnit = KIM_CHARGE_UNIT_e; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_ps; - } else if (system_str == "si") { - lengthUnit = KIM_LENGTH_UNIT_m; - energyUnit = KIM_ENERGY_UNIT_J; - chargeUnit = KIM_CHARGE_UNIT_C; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_s; - } else if (system_str == "cgs") { - lengthUnit = KIM_LENGTH_UNIT_cm; - energyUnit = KIM_ENERGY_UNIT_erg; - chargeUnit = KIM_CHARGE_UNIT_statC; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_s; - } else if (system_str == "electron") { - lengthUnit = KIM_LENGTH_UNIT_Bohr; - energyUnit = KIM_ENERGY_UNIT_Hartree; - chargeUnit = KIM_CHARGE_UNIT_e; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_fs; - } else if ((system_str == "lj") || - (system_str == "micro") || - (system_str == "nano")) { - error->all(FLERR, "LAMMPS unit_style {} not supported " - "by KIM models", system_str); - } else { - error->all(FLERR, "Unknown unit_style"); - } -} -} // namespace +using kim_units::get_kim_unit_names; /* ---------------------------------------------------------------------- */ @@ -159,20 +106,17 @@ void KimParam::command(int narg, char **arg) error->all(FLERR, "Incorrect arguments in 'kim param' command.\n" "'kim param get/set' is mandatory"); - int const ifix = modify->find_fix("KIM_MODEL_STORE"); - if (ifix >= 0) { - auto fix_store = reinterpret_cast(modify->fix[ifix]); - - KIM_SimulatorModel *simulatorModel = - reinterpret_cast( - fix_store->getptr("simulator_model")); + auto fix_store = dynamic_cast(modify->get_fix_by_id("KIM_MODEL_STORE")); + if (fix_store) { + auto *simulatorModel = reinterpret_cast( + fix_store->getptr("simulator_model")); if (simulatorModel) error->all(FLERR, "'kim param' can only be used with a KIM Portable Model"); } - input->write_echo(fmt::format("#=== BEGIN kim param {} ===================" - "==================\n", kim_param_get_set)); + input->write_echo(fmt::format("#=== BEGIN kim param {} =====================================\n", + kim_param_get_set)); KIM_Model *pkim = nullptr; @@ -431,6 +375,6 @@ void KimParam::command(int narg, char **arg) } else error->all(FLERR, "This model has No mutable parameters"); - input->write_echo(fmt::format("#=== END kim param {} =====================" - "==================\n", kim_param_get_set)); + input->write_echo(fmt::format("#=== END kim param {} =======================================\n", + kim_param_get_set)); } diff --git a/src/KIM/kim_param.h b/src/KIM/kim_param.h index d160d2f1b1..daefb13017 100644 --- a/src/KIM/kim_param.h +++ b/src/KIM/kim_param.h @@ -67,7 +67,5 @@ class KimParam : protected Pointers { KimParam(class LAMMPS *lmp); void command(int, char **); }; - } // namespace LAMMPS_NS - #endif diff --git a/src/KIM/kim_property.cpp b/src/KIM/kim_property.cpp index 7d060ad8b1..6ddebeefc2 100644 --- a/src/KIM/kim_property.cpp +++ b/src/KIM/kim_property.cpp @@ -85,14 +85,13 @@ void KimProperty::command(int narg, char **arg) { #if LMP_PYTHON #if PY_MAJOR_VERSION >= 3 - if (narg < 2) error->all(FLERR, "Invalid 'kim property' command"); + if (narg < 2) utils::missing_cmd_args(FLERR, "kim property", error); const std::string subcmd(arg[0]); if ((subcmd != "create") && (subcmd != "destroy") && (subcmd != "modify") && (subcmd != "remove") && (subcmd != "dump")) { - std::string msg("Incorrect arguments in 'kim property' command.\n"); - msg += "'kim property create/destroy/modify/remove/dump' is mandatory"; - error->all(FLERR, msg); + error->all(FLERR, "Incorrect first argument {} to 'kim property' command.\n" + "One of create, destroy, modify, remove, or dump is mandatory", subcmd); } input->write_echo("#=== kim property =====================================" @@ -117,12 +116,11 @@ void KimProperty::command(int narg, char **arg) kim_property = PyImport_Import(obj); if (!kim_property) { PyGILState_Release(gstate); - std::string msg("Unable to import Python kim_property module!"); - msg += "\nkim-property Python package can be installed with pip:\n"; - msg += "'pip install kim-property'\nSee the installation instructions "; - msg += "at\nhttps://github.com/openkim/kim-property#installing-kim-"; - msg += "property\nfor detailed information"; - error->all(FLERR, msg); + error->all(FLERR, "Unable to import Python kim_property module!\n" + "kim-property Python package can be installed with pip:\n" + "'pip install kim-property'\nSee the installation instructions at\n" + "https://github.com/openkim/kim-property#installing-kim-property\n" + "for detailed information"); } // Decrementing of the reference count @@ -147,9 +145,8 @@ void KimProperty::command(int narg, char **arg) PyObject_GetAttrString(kim_property, "kim_property_create"); if (!pFunc) { PyGILState_Release(gstate); - std::string msg("Unable to get an attribute named "); - msg += "'kim_property_create' from a kim_property object"; - error->all(FLERR, msg); + error->all(FLERR, "Unable to get an attribute named 'kim_property_create' from " + "a kim_property object"); } // Decrementing of the reference count @@ -182,16 +179,13 @@ void KimProperty::command(int narg, char **arg) if (!pValue) { PyErr_Print(); PyGILState_Release(gstate); - std::string msg("Python 'kim_property_create' function "); - msg += "evaluation failed"; - error->one(FLERR, msg); + error->all(FLERR, "Python 'kim_property_create' function evaluation failed"); } // Python function returned a string value const char *pystr = PyUnicode_AsUTF8(pValue); if (kim_str) input->variable->set_string("kim_property_str", pystr); - else - input->variable->set(fmt::format("kim_property_str string '{}'", pystr)); + else input->variable->set(fmt::format("kim_property_str string '{}'", pystr)); Py_XDECREF(pArgs); Py_XDECREF(pFunc); @@ -216,9 +210,8 @@ void KimProperty::command(int narg, char **arg) PyObject_GetAttrString(kim_property, "kim_property_destroy"); if (!pFunc) { PyGILState_Release(gstate); - std::string msg("Unable to get an attribute named "); - msg += "'kim_property_destroy' from a kim_property object"; - error->all(FLERR, msg); + error->all(FLERR, "Unable to get an attribute named 'kim_property_destroy' " + "from a kim_property object"); } // Decrementing of the reference count @@ -244,9 +237,7 @@ void KimProperty::command(int narg, char **arg) if (!pValue) { PyErr_Print(); PyGILState_Release(gstate); - std::string msg("Python 'kim_property_destroy' function "); - msg += "evaluation failed"; - error->one(FLERR, msg); + error->all(FLERR, "Python 'kim_property_destroy' function evaluation failed"); } // Python function returned a string value @@ -276,9 +267,8 @@ void KimProperty::command(int narg, char **arg) PyObject_GetAttrString(kim_property, "kim_property_modify"); if (!pFunc) { PyGILState_Release(gstate); - std::string msg("Unable to get an attribute named "); - msg += "'kim_property_modify' from a kim_property object"; - error->all(FLERR, msg); + error->all(FLERR, "Unable to get an attribute named 'kim_property_modify' " + "from a kim_property object"); } // Decrementing of the reference count @@ -309,9 +299,7 @@ void KimProperty::command(int narg, char **arg) if (!pValue) { PyErr_Print(); PyGILState_Release(gstate); - std::string msg("Python 'kim_property_modify' function "); - msg += "evaluation failed"; - error->one(FLERR, msg); + error->all(FLERR, "Python 'kim_property_modify' function evaluation failed"); } // Python function returned a string value @@ -341,9 +329,8 @@ void KimProperty::command(int narg, char **arg) PyObject_GetAttrString(kim_property, "kim_property_remove"); if (!pFunc) { PyGILState_Release(gstate); - std::string msg("Unable to get an attribute named "); - msg += "'kim_property_remove' from a kim_property object"; - error->all(FLERR, msg); + error->all(FLERR, "Unable to get an attribute named " + "'kim_property_remove' from a kim_property object"); } // Decrementing of the reference count @@ -374,9 +361,7 @@ void KimProperty::command(int narg, char **arg) if (!pValue) { PyErr_Print(); PyGILState_Release(gstate); - std::string msg("Python 'kim_property_remove' function "); - msg += "evaluation failed"; - error->one(FLERR, msg); + error->all(FLERR, "Python 'kim_property_remove' function evaluation failed"); } // Python function returned a string value @@ -404,9 +389,8 @@ void KimProperty::command(int narg, char **arg) PyObject_GetAttrString(kim_property, "kim_property_dump"); if (!pFunc) { PyGILState_Release(gstate); - std::string msg("Unable to get an attribute named "); - msg += "'kim_property_dump' from a kim_property object"; - error->all(FLERR, msg); + error->all(FLERR, "Unable to get an attribute named " + "'kim_property_dump' from a kim_property object"); } // Decrementing of the reference count @@ -428,14 +412,12 @@ void KimProperty::command(int narg, char **arg) if (comm->me == 0) { // call the Python kim_property_dump function - // error check with one() since only some procs may fail + // error check with one() since only root process calls. pValue = PyObject_CallObject(pFunc, pArgs); if (!pValue) { PyErr_Print(); PyGILState_Release(gstate); - std::string msg("Python 'kim_property_dump' function "); - msg += "evaluation failed"; - error->one(FLERR, msg); + error->one(FLERR, "Python 'kim_property_dump' function evaluation failed"); } } else pValue = nullptr; diff --git a/src/KIM/kim_query.cpp b/src/KIM/kim_query.cpp index ebb1790e6e..610ec18fcc 100644 --- a/src/KIM/kim_query.cpp +++ b/src/KIM/kim_query.cpp @@ -136,11 +136,9 @@ void KimQuery::command(int narg, char **arg) // check the query_args format (a series of keyword=value pairs) for (int i = 2; i < narg; ++i) { if (!utils::strmatch(arg[i], "[=][\\[].*[\\]]")) - error->all(FLERR, "Illegal query format.\nInput argument " - "of `{}` to 'kim query' is wrong. The " - "query format is the keyword=[value], " - "where value is always an array of one or " - "more comma-separated items", arg[i]); + error->all(FLERR, "Illegal query format.\nInput argument of `{}` to 'kim query' is wrong. " + "The query format is the keyword=[value], where value is always an array of one " + "or more comma-separated items", arg[i]); } if (query_function != "get_available_models") { @@ -156,15 +154,13 @@ void KimQuery::command(int narg, char **arg) // if the model name is not provided by the user if (model_name.empty()) { // check if we had a kim init command by finding fix STORE/KIM - const int ifix = modify->find_fix("KIM_MODEL_STORE"); - if (ifix >= 0) { - auto fix_store = dynamic_cast(modify->fix[ifix]); + auto fix_store = dynamic_cast(modify->get_fix_by_id("KIM_MODEL_STORE")); + if (fix_store) { char *model_name_c = (char *) fix_store->getptr("model_name"); model_name = model_name_c; } else { - error->all(FLERR, "Illegal query format.\nMust use 'kim init' " - "before 'kim query' or must provide the model name " - "after query function with the format of " + error->all(FLERR, "Illegal query format.\nMust use 'kim init' before 'kim query' " + "or must provide the model name after query function with the format of " "'model=[model_name]'"); } } @@ -179,16 +175,14 @@ void KimQuery::command(int narg, char **arg) // and then the error message that was returned by the web server if (strlen(value) == 0) { - auto msg = fmt::format("OpenKIM query failed: {}", value + 1); - delete [] value; - error->all(FLERR, msg); + error->all(FLERR, "OpenKIM query failed: {}", value + 1); + delete[] value; } else if (strcmp(value, "EMPTY") == 0) { - delete [] value; + delete[] value; error->all(FLERR, "OpenKIM query returned no results"); } - input->write_echo("#=== BEGIN kim-query ==================================" - "=======\n"); + input->write_echo("#=== BEGIN kim-query =========================================\n"); // trim list of models to those that are installed on the system if (query_function == "get_available_models") { Tokenizer vals(value, ", \""); @@ -199,7 +193,7 @@ void KimQuery::command(int narg, char **arg) KIM_CollectionItemType typ; if (KIM_Collections_Create(&collections)) { - delete [] value; + delete[] value; error->all(FLERR, "Unable to access KIM Collections to find Model"); } @@ -219,7 +213,7 @@ void KimQuery::command(int narg, char **arg) fmt::format("# Missing OpenKIM models: {}\n\n", missing)); if (available.empty()) { - delete [] value; + delete[] value; error->all(FLERR,"There are no matching OpenKIM models installed on the system"); } @@ -262,13 +256,12 @@ void KimQuery::command(int narg, char **arg) input->variable->set(setcmd); input->write_echo(fmt::format("variable {}\n", setcmd)); } - input->write_echo("#=== END kim-query ====================================" - "=======\n\n"); + input->write_echo("#=== END kim-query ===========================================\n\n"); - delete [] value; + delete[] value; #else - error->all(FLERR, "Cannot use 'kim query' command when KIM package " - "is compiled without support for libcurl"); + error->all(FLERR, "Cannot use 'kim query' command when KIM package is compiled without " + "support for libcurl"); #endif } diff --git a/src/KIM/kim_units.cpp b/src/KIM/kim_units.cpp index 7c57e3b71a..7687481cb9 100644 --- a/src/KIM/kim_units.cpp +++ b/src/KIM/kim_units.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -54,179 +53,182 @@ Designed for use with the kim-api-2.0.2 (and newer) package ------------------------------------------------------------------------- */ +#include "kim_units.h" + +#include "error.h" + #include #include #include #include using namespace std; - -namespace -{ +using namespace LAMMPS_NS; +using namespace kim_units; // Constants of nature and basic conversion factors // Source: https://physics.nist.gov/cuu/Constants/Table/allascii.txt // Working with NIST values even when there are newer values for // compatibility with LAMMPS +// clang-format off + /*---------------------- Fundamental constants ------------------------ */ -double const boltz_si = 1.38064852e-23; // [J K^-1] Boltzmann's factor - // (NIST value) -double const Nav = 6.022140857e23; // [unitless] Avogadro's number - // (NIST value) -// double const Nav = 6.02214076e23; // [unitless] Avogadro's number - // (official value May 2019) -double const me_si = 9.10938356e-31; // [kg] electron rest mass - // (NIST value) -// double me_si = 9.10938291e-31; // [kg] electron rest mass -double const e_si = 1.6021766208e-19; // [C] elementary charge - // (charge of an electron/proton) - // (NIST value) +static double constexpr boltz_si = 1.38064852e-23; // [J K^-1] Boltzmann's factor + // (NIST value) +static double constexpr Nav = 6.022140857e23; // [unitless] Avogadro's number (NIST value) +// static double constexp Nav = 6.02214076e23; // [unitless] Avogadro's number + // (official value May 2019) +static double constexpr me_si = 9.10938356e-31; // [kg] electron rest mass (NIST value) +// static double constexpr me_si = 9.10938291e-31; // [kg] electron rest mass +static double constexpr e_si = 1.6021766208e-19; // [C] elementary charge + // (charge of an electron/proton) + // (NIST value) /*---------------------- Distance units ------------------------ */ -double const bohr_si = 5.2917721067e-11; // [m] Bohr unit (distance between - // nucleus and electron in H) - // (NIST value) -double const angstrom_si = 1e-10; // [m] Angstrom -double const centimeter_si = 1e-2; // [m] centimeter -double const micrometer_si = 1e-6; // [m] micrometer (micron) -double const nanometer_si = 1e-9; // [m] nanometer +static double constexpr bohr_si = 5.2917721067e-11; // [m] Bohr unit (distance between + // nucleus and electron in H) (NIST value) +static double constexpr angstrom_si = 1e-10; // [m] Angstrom +static double constexpr centimeter_si = 1e-2; // [m] centimeter +static double constexpr micrometer_si = 1e-6; // [m] micrometer (micron) +static double constexpr nanometer_si = 1e-9; // [m] nanometer /*---------------------- Mass units ------------------------ */ -double const gram_per_mole_si = 1e-3/Nav; // [kg] gram per mole -double const amu_si = 1e-3/Nav; // [kg] atomic mass unit (molecular - // weight) For example, the mean - // molecular weight of water - // is 18.015 atomic mass units - // (amu), so one mole of water - // weight 18.015 grams. -double const gram_si = 1e-3; // [kg] gram -double const picogram_si = 1e-15; // [kg] picogram -double const attogram_si = 1e-21; // [kg[ attogram +static double constexpr gram_per_mole_si = 1e-3/Nav; // [kg] gram per mole +static double constexpr amu_si = 1e-3/Nav; // [kg] atomic mass unit (molecular + // weight) For example, the mean + // molecular weight of water + // is 18.015 atomic mass units + // (amu), so one mole of water + // weight 18.015 grams. +static double constexpr gram_si = 1e-3; // [kg] gram +static double constexpr picogram_si = 1e-15; // [kg] picogram +static double constexpr attogram_si = 1e-21; // [kg[ attogram /*---------------------- Time units ------------------------ */ -double const atu_si = 2.418884326509e-17; // [s] atomic time unit - // ( = hbar/E_h where E_h is the - // Hartree energy) (NIST value) -double const atu_electron_si = atu_si*sqrt(amu_si/me_si); - // [s] atomic time unit - // used in electron system (see https://sourceforge.net/p/lammps/mailman/lammps-users/thread/BCA2BDB2-BA03-4280-896F-1E6120EF47B2%40caltech.edu/) -double const microsecond_si = 1e-6; // [s] microsecond -double const nanosecond_si = 1e-9; // [s] nanosecond -double const picosecond_si = 1e-12; // [s] picosecond -double const femtosecond_si = 1e-15; // [s] femtosecond +static double constexpr atu_si = 2.418884326509e-17; // [s] atomic time unit + // ( = hbar/E_h where E_h is the + // Hartree energy) (NIST value) +static double constexpr atu_electron_si = 5153034.567408186; // must not use sqrt() in constexpr +// static double constexpr atu_electron_si = atu_si*sqrt(amu_si/me_si); + // [s] atomic time unit + // used in electron system (see https://sourceforge.net/p/lammps/mailman/lammps-users/thread/BCA2BDB2-BA03-4280-896F-1E6120EF47B2%40caltech.edu/) +static double constexpr microsecond_si = 1e-6; // [s] microsecond +static double constexpr nanosecond_si = 1e-9; // [s] nanosecond +static double constexpr picosecond_si = 1e-12; // [s] picosecond +static double constexpr femtosecond_si = 1e-15; // [s] femtosecond /*---------------------- Density units ------------------------ */ -double const gram_per_centimetercu_si = - gram_si/pow(centimeter_si,3); // [kg/m^3] gram/centimeter^3 -double const amu_per_bohrcu_si = amu_si/pow(bohr_si,3); // [kg/m^3] amu/bohr^3 -double const picogram_per_micrometercu_si = - picogram_si/pow(micrometer_si,3); // [kg/m^3] picogram/micrometer^3 -double const attogram_per_nanometercu_si = - attogram_si/pow(nanometer_si,3); // [kg/m^3] attogram/nanometer^3 +static double constexpr gram_per_centimetercu_si = + gram_si/centimeter_si/centimeter_si/centimeter_si; // [kg/m^3] gram/centimeter^3 +static double constexpr amu_per_bohrcu_si = amu_si/bohr_si/bohr_si/bohr_si; // [kg/m^3] amu/bohr^3 +static double constexpr picogram_per_micrometercu_si = + picogram_si/micrometer_si/micrometer_si/micrometer_si; // [kg/m^3] picogram/micrometer^3 +static double constexpr attogram_per_nanometercu_si = + attogram_si/nanometer_si/nanometer_si/nanometer_si; // [kg/m^3] attogram/nanometer^3 /*---------------------- Energy/torque units ------------------------ */ -double const kcal_si = 4184.0; // [J] kilocalorie (heat energy - // involved in warming up one - // kilogram of water by one - // degree Kelvin) -double const ev_si = 1.6021766208e-19; // [J] electron volt (amount of - // energy gained or lost by the - // charge of a single electron - // moving across an electric - // potential difference of one - // volt.) (NIST value) -double const hartree_si = 4.359744650e-18; // [J] Hartree (approximately the - // electric potential energy of - // the hydrogen atom in its - // ground state) (NIST value) -double const kcal_per_mole_si = kcal_si/Nav;// [J] kcal/mole -double const erg_si = 1e-7; // [J] erg -double const dyne_centimeter_si = 1e-7; // [J[ dyne*centimeter -double const picogram_micrometersq_per_microsecondsq_si = - picogram_si*pow(micrometer_si,2)/pow(microsecond_si,2); - // [J] picogram*micrometer^2/ - // microsecond^2 -double const attogram_nanometersq_per_nanosecondsq_si = - attogram_si*pow(nanometer_si,2)/pow(nanosecond_si,2); - // [J] attogram*nanometer^2/ - // nanosecond^2 +static double constexpr kcal_si = 4184.0; // [J] kilocalorie (heat energy + // involved in warming up one + // kilogram of water by one + // degree Kelvin) +static double constexpr ev_si = 1.6021766208e-19; // [J] electron volt (amount of + // energy gained or lost by the + // charge of a single electron + // moving across an electric + // potential difference of one + // volt.) (NIST value) +static double constexpr hartree_si = 4.359744650e-18; // [J] Hartree (approximately the + // electric potential energy of + // the hydrogen atom in its + // ground state) (NIST value) +static double constexpr kcal_per_mole_si = kcal_si/Nav;// [J] kcal/mole +static double constexpr erg_si = 1e-7; // [J] erg +static double constexpr dyne_centimeter_si = 1e-7; // [J[ dyne*centimeter +static double constexpr picogram_micrometersq_per_microsecondsq_si = + picogram_si*micrometer_si/microsecond_si*micrometer_si/microsecond_si; + // [J] picogram*micrometer^2/ + // microsecond^2 +static double constexpr attogram_nanometersq_per_nanosecondsq_si = + attogram_si*nanometer_si/nanosecond_si*nanometer_si/nanosecond_si; + // [J] attogram*nanometer^2/ + // nanosecond^2 /*---------------------- Velocity units ------------------------ */ -double const angstrom_per_femtosecond_si = +static double constexpr angstrom_per_femtosecond_si = angstrom_si/femtosecond_si; // [m/s] Angstrom/femtosecond -double const angstrom_per_picosecond_si = +static double constexpr angstrom_per_picosecond_si = angstrom_si/picosecond_si; // [m/s] Angstrom/picosecond -double const micrometer_per_microsecond_si = +static double constexpr micrometer_per_microsecond_si = micrometer_si/microsecond_si; // [m/s] micrometer/microsecond -double const nanometer_per_nanosecond_si = +static double constexpr nanometer_per_nanosecond_si = nanometer_si/nanosecond_si; // [m/s] nanometer/nanosecond -double const centimeter_per_second_si = +static double constexpr centimeter_per_second_si = centimeter_si; // [m/s] centimeter/second -double const bohr_per_atu_electron_si = +static double constexpr bohr_per_atu_electron_si = bohr_si/atu_electron_si; // [m/s] bohr/atu /*---------------------- Force units ------------------------ */ -double const kcal_per_mole_angstrom_si = +static double constexpr kcal_per_mole_angstrom_si = kcal_per_mole_si/angstrom_si; // [N] kcal/(mole*Angstrom) -double const ev_per_angstrom_si = +static double constexpr ev_per_angstrom_si = ev_si/angstrom_si; // [N] eV/Angstrom -double const dyne_si = +static double constexpr dyne_si = dyne_centimeter_si/centimeter_si; // [N] dyne -double const hartree_per_bohr_si = +static double constexpr hartree_per_bohr_si = hartree_si/bohr_si; // [N] hartree/bohr -double const picogram_micrometer_per_microsecondsq_si = - picogram_si*micrometer_si/pow(microsecond_si,2); +static double constexpr picogram_micrometer_per_microsecondsq_si = + picogram_si*micrometer_si/microsecond_si/microsecond_si; // [N] picogram*micrometer/ // microsecond^2 -double const attogram_nanometer_per_nanosecondsq_si = - attogram_si*nanometer_si/pow(nanosecond_si,2); +static double constexpr attogram_nanometer_per_nanosecondsq_si = + attogram_si*nanometer_si/nanosecond_si/nanosecond_si; // [N] attogram*nanometer/ // nanosecond^2 /*---------------------- Pressure units ------------------------ */ -double const atmosphere_si = 101325.0; // [Pa] standard atmosphere (NIST value) -double const bar_si = 1e5; // [Pa] bar -double const dyne_per_centimetersq_si = - dyne_centimeter_si/pow(centimeter_si,3); +static double constexpr atmosphere_si = 101325.0; // [Pa] standard atmosphere (NIST value) +static double constexpr bar_si = 1e5; // [Pa] bar +static double constexpr dyne_per_centimetersq_si = + dyne_centimeter_si/centimeter_si/centimeter_si/centimeter_si; // [Pa] dyne/centimeter^2 -double const picogram_per_micrometer_microsecondsq_si = - picogram_si/(micrometer_si*pow(microsecond_si,2)); +static double constexpr picogram_per_micrometer_microsecondsq_si = + picogram_si/(micrometer_si*microsecond_si*microsecond_si); // [Pa] picogram/(micrometer* // microsecond^2) -double const attogram_per_nanometer_nanosecondsq_si = - attogram_si/(nanometer_si*pow(nanosecond_si,2)); +static double constexpr attogram_per_nanometer_nanosecondsq_si = + attogram_si/(nanometer_si*nanosecond_si*nanosecond_si); // [Pa] attogram/(nanometer*nanosecond^2) /*---------------------- Viscosity units ------------------------ */ -double const poise_si = 0.1; // [Pa*s] Poise -double const amu_per_bohr_femtosecond_si = +static double constexpr poise_si = 0.1; // [Pa*s] Poise +static double constexpr amu_per_bohr_femtosecond_si = amu_si/(bohr_si*femtosecond_si); // [Pa*s] amu/(bohr*femtosecond) -double const picogram_per_micrometer_microsecond_si = +static double constexpr picogram_per_micrometer_microsecond_si = picogram_si/(micrometer_si*microsecond_si); // [Pa*s] picogram/(micrometer* // microsecond) -double const attogram_per_nanometer_nanosecond_si = +static double constexpr attogram_per_nanometer_nanosecond_si = attogram_si/(nanometer_si*nanosecond_si); // [Pa*s] attogram/(nanometer* // nanosecond) @@ -234,39 +236,41 @@ double const attogram_per_nanometer_nanosecond_si = /*---------------------- Charge units ------------------------ */ -double const echarge_si = e_si; // [C] electron charge unit -double const statcoulomb_si = e_si/4.8032044e-10; // [C] Statcoulomb or esu +static double constexpr echarge_si = e_si; // [C] electron charge unit +static double constexpr statcoulomb_si = e_si/4.8032044e-10; // [C] Statcoulomb or esu // (value from LAMMPS units // documentation) -double const picocoulomb_si = 1e-12; // [C] picocoulomb +static double constexpr picocoulomb_si = 1e-12; // [C] picocoulomb /*---------------------- Dipole units ------------------------ */ -double const electron_angstrom_si = echarge_si*angstrom_si; +static double constexpr electron_angstrom_si = echarge_si*angstrom_si; // [C*m] electron*angstrom -double const statcoulomb_centimeter_si = statcoulomb_si*centimeter_si; +static double constexpr statcoulomb_centimeter_si = statcoulomb_si*centimeter_si; // [C*m] statcoulomb*centimeter -double const debye_si = 1e-18*statcoulomb_centimeter_si; +static double constexpr debye_si = 1e-18*statcoulomb_centimeter_si; // [C*m] Debye -double const picocoulomb_micrometer_si = picocoulomb_si*micrometer_si; +static double constexpr picocoulomb_micrometer_si = picocoulomb_si*micrometer_si; // [C*m] picocoulomb*micrometer -double const electron_nanometer_si = echarge_si*nanometer_si; +static double constexpr electron_nanometer_si = echarge_si*nanometer_si; // [C*m] electron*nanometer /*---------------------- Electric field units ------------------------ */ -double const volt_per_angstrom_si = 1.0/angstrom_si;// [V/m] volt/angstrom -double const statvolt_per_centimeter_si = +static double constexpr volt_per_angstrom_si = 1.0/angstrom_si;// [V/m] volt/angstrom +static double constexpr statvolt_per_centimeter_si = erg_si/(statcoulomb_si*centimeter_si); // [V/m] statvolt/centimeter -double const volt_per_centimeter_si = +static double constexpr volt_per_centimeter_si = 1.0/centimeter_si; // [V/m] volt/centimeter -double const volt_per_micrometer_si = +static double constexpr volt_per_micrometer_si = 1.0/micrometer_si; // [V/m] volt/micrometer -double const volt_per_nanometer_si = +static double constexpr volt_per_nanometer_si = 1.0/nanometer_si; // [V/m] volt/nanometer +namespace LAMMPS_NS { +namespace kim_units { // Define enumerations enum sys_type { @@ -388,20 +392,42 @@ enum units attogram_per_nanometercu = 1405 }; +void initialize_dictionaries(); +units get_lammps_system_unit(sys_type system_enum, unit_type unit_type_enum); +double get_mass_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_distance_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_time_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_energy_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_velocity_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_force_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_torque_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_temperature_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_pressure_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_viscosity_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_charge_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_dipole_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_efield_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_density_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_unit_conversion_factor(unit_type &unit_type_enum, sys_type from_system_enum, + sys_type to_system_enum); + +} +} + // Define dictionaries -map system_dic; -map unit_dic; -map units_real_dic; -map units_metal_dic; -map units_si_dic; -map units_cgs_dic; -map units_electron_dic; -map units_micro_dic; -map units_nano_dic; +static map system_dic; +static map unit_dic; +static map units_real_dic; +static map units_metal_dic; +static map units_si_dic; +static map units_cgs_dic; +static map units_electron_dic; +static map units_micro_dic; +static map units_nano_dic; /* ---------------------------------------------------------------------- */ -void initialize_dictionaries() +void kim_units::initialize_dictionaries() { system_dic["real"] = real; system_dic["metal"] = metal; @@ -537,7 +563,7 @@ void initialize_dictionaries() // Get the enumeration for the unit of type `unit_type_enum` // for LAMMPS system `system_enum`. -units get_lammps_system_unit(sys_type system_enum, unit_type unit_type_enum) +units kim_units::get_lammps_system_unit(sys_type system_enum, unit_type unit_type_enum) { switch(system_enum) { case real : @@ -561,7 +587,7 @@ units get_lammps_system_unit(sys_type system_enum, unit_type unit_type_enum) /* ---------------------------------------------------------------------- */ // Mass conversion -double get_mass_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_mass_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -619,7 +645,7 @@ double get_mass_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Distance conversion -double get_distance_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_distance_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -677,7 +703,7 @@ double get_distance_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Time conversion -double get_time_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_time_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -722,7 +748,7 @@ double get_time_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Energy conversion -double get_energy_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_energy_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; @@ -796,7 +822,7 @@ double get_energy_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Velocity conversion -double get_velocity_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_velocity_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -869,7 +895,7 @@ double get_velocity_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Force conversion -double get_force_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_force_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -942,7 +968,7 @@ double get_force_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Torque conversion -double get_torque_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_torque_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -1015,7 +1041,7 @@ double get_torque_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Temperature conversion -double get_temperature_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_temperature_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; @@ -1026,7 +1052,7 @@ double get_temperature_conversion_factor(units from_unit_enum, units to_unit_enu /* ---------------------------------------------------------------------- */ // Pressure conversion -double get_pressure_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_pressure_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -1084,7 +1110,7 @@ double get_pressure_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Viscosity conversion -double get_viscosity_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_viscosity_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -1129,7 +1155,7 @@ double get_viscosity_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Charge conversion -double get_charge_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_charge_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -1163,7 +1189,7 @@ double get_charge_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Dipole conversion -double get_dipole_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_dipole_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -1221,7 +1247,7 @@ double get_dipole_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Electric field conversion -double get_efield_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_efield_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -1278,8 +1304,8 @@ double get_efield_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ -// Demsity conversion -double get_density_conversion_factor(units from_unit_enum, units to_unit_enum) +// Density conversion +double kim_units::get_density_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -1325,9 +1351,9 @@ double get_density_conversion_factor(units from_unit_enum, units to_unit_enum) // This routine returns the unit conversion factor between the // `from_system_enum` to the `to_system_enum` for the `unit_type_enum`. -double get_unit_conversion_factor(unit_type &unit_type_enum, - sys_type from_system_enum, - sys_type to_system_enum) +double kim_units::get_unit_conversion_factor(unit_type &unit_type_enum, + sys_type from_system_enum, + sys_type to_system_enum) { units from_unit = get_lammps_system_unit(from_system_enum, unit_type_enum); units to_unit = get_lammps_system_unit(to_system_enum, unit_type_enum); @@ -1364,48 +1390,87 @@ double get_unit_conversion_factor(unit_type &unit_type_enum, } } -} // end of anonymous name space - /* ---------------------------------------------------------------------- */ - +// clang-format on // Wrapper to the routine that gets the unit conversion. Translates strings // to enumerations and then call get_unit_conversion_factor() -int lammps_unit_conversion(const string &unit_type_str, - const string &from_system_str, - const string &to_system_str, - double &conversion_factor) +int kim_units::lammps_unit_conversion(const string &unit_type_str, const string &from_system_str, + const string &to_system_str, double &conversion_factor) { - // initialize - conversion_factor = 0.0; - initialize_dictionaries(); + // initialize + conversion_factor = 0.0; + initialize_dictionaries(); - // convert input to enumeration - unit_type unit_type_enum; - { - map::const_iterator itr = unit_dic.find(unit_type_str); - if (itr != unit_dic.end()) unit_type_enum = itr->second; - else return 1; // error - } - sys_type from_system_enum; - { - map::const_iterator - itr = system_dic.find(from_system_str); - if (itr != system_dic.end()) from_system_enum = itr->second; - else return 1; // error - } - sys_type to_system_enum; - { - map::const_iterator - itr = system_dic.find(to_system_str); - if (itr != system_dic.end()) to_system_enum = itr->second; - else return 1; - } + // convert input to enumeration + unit_type unit_type_enum; + { + map::const_iterator itr = unit_dic.find(unit_type_str); + if (itr != unit_dic.end()) + unit_type_enum = itr->second; + else + return 1; // error + } + sys_type from_system_enum; + { + map::const_iterator itr = system_dic.find(from_system_str); + if (itr != system_dic.end()) + from_system_enum = itr->second; + else + return 1; // error + } + sys_type to_system_enum; + { + map::const_iterator itr = system_dic.find(to_system_str); + if (itr != system_dic.end()) + to_system_enum = itr->second; + else + return 1; + } - // process unit conversions - conversion_factor = get_unit_conversion_factor(unit_type_enum, - from_system_enum, - to_system_enum); - return 0; + // process unit conversions + conversion_factor = get_unit_conversion_factor(unit_type_enum, from_system_enum, to_system_enum); + return 0; } - +void kim_units::get_kim_unit_names(char const *const system, KIM_LengthUnit &lengthUnit, + KIM_EnergyUnit &energyUnit, KIM_ChargeUnit &chargeUnit, + KIM_TemperatureUnit &temperatureUnit, KIM_TimeUnit &timeUnit, + Error *error) +{ + const std::string system_str(system); + if (system_str == "real") { + lengthUnit = KIM_LENGTH_UNIT_A; + energyUnit = KIM_ENERGY_UNIT_kcal_mol; + chargeUnit = KIM_CHARGE_UNIT_e; + temperatureUnit = KIM_TEMPERATURE_UNIT_K; + timeUnit = KIM_TIME_UNIT_fs; + } else if (system_str == "metal") { + lengthUnit = KIM_LENGTH_UNIT_A; + energyUnit = KIM_ENERGY_UNIT_eV; + chargeUnit = KIM_CHARGE_UNIT_e; + temperatureUnit = KIM_TEMPERATURE_UNIT_K; + timeUnit = KIM_TIME_UNIT_ps; + } else if (system_str == "si") { + lengthUnit = KIM_LENGTH_UNIT_m; + energyUnit = KIM_ENERGY_UNIT_J; + chargeUnit = KIM_CHARGE_UNIT_C; + temperatureUnit = KIM_TEMPERATURE_UNIT_K; + timeUnit = KIM_TIME_UNIT_s; + } else if (system_str == "cgs") { + lengthUnit = KIM_LENGTH_UNIT_cm; + energyUnit = KIM_ENERGY_UNIT_erg; + chargeUnit = KIM_CHARGE_UNIT_statC; + temperatureUnit = KIM_TEMPERATURE_UNIT_K; + timeUnit = KIM_TIME_UNIT_s; + } else if (system_str == "electron") { + lengthUnit = KIM_LENGTH_UNIT_Bohr; + energyUnit = KIM_ENERGY_UNIT_Hartree; + chargeUnit = KIM_CHARGE_UNIT_e; + temperatureUnit = KIM_TEMPERATURE_UNIT_K; + timeUnit = KIM_TIME_UNIT_fs; + } else if ((system_str == "lj") || (system_str == "micro") || (system_str == "nano")) { + error->all(FLERR, "LAMMPS unit_style {} not supported by KIM models", system_str); + } else { + error->all(FLERR, "Unknown unit_style {}", system_str); + } +} diff --git a/src/KIM/kim_units.h b/src/KIM/kim_units.h index e85760ab82..dfb30c0a4c 100644 --- a/src/KIM/kim_units.h +++ b/src/KIM/kim_units.h @@ -56,7 +56,20 @@ #define LMP_KIM_UNITS_H #include +extern "C" { +#include "KIM_SimulatorHeaders.h" +} -int lammps_unit_conversion(const std::string &unit_type_str, const std::string &from_system_str, - const std::string &to_system_str, double &conversion_factor); +namespace LAMMPS_NS { +class Error; + +namespace kim_units { + int lammps_unit_conversion(const std::string &unit_type_str, const std::string &from_system_str, + const std::string &to_system_str, double &conversion_factor); + void get_kim_unit_names(char const *const system, KIM_LengthUnit &lengthUnit, + KIM_EnergyUnit &energyUnit, KIM_ChargeUnit &chargeUnit, + KIM_TemperatureUnit &temperatureUnit, KIM_TimeUnit &timeUnit, + Error *error); +} // namespace kim_units +} // namespace LAMMPS_NS #endif diff --git a/src/KIM/pair_kim.cpp b/src/KIM/pair_kim.cpp index 6a5a792797..14564162e2 100644 --- a/src/KIM/pair_kim.cpp +++ b/src/KIM/pair_kim.cpp @@ -218,8 +218,7 @@ void PairKIM::compute(int eflag, int vflag) kim_particleSpecies); memory->create(kim_particleContributing,lmps_maxalloc, "pair:kim_particleContributing"); - kimerror = kimerror || KIM_ComputeArguments_SetArgumentPointerInteger( - pargs, + kimerror = kimerror || KIM_ComputeArguments_SetArgumentPointerInteger(pargs, KIM_COMPUTE_ARGUMENT_NAME_particleContributing, kim_particleContributing); if (kimerror) @@ -248,7 +247,7 @@ void PairKIM::compute(int eflag, int vflag) // compute via KIM model int kimerror = KIM_Model_Compute(pkim, pargs); - if (kimerror) error->all(FLERR,"KIM Compute returned error"); + if (kimerror) error->all(FLERR,"KIM Compute returned error {}", kimerror); // compute virial before reverse comm! if (vflag_global) @@ -454,7 +453,7 @@ void PairKIM::coeff(int narg, char **arg) kimerror = KIM_Model_GetParameterMetadata(pkim, param_index, &kim_DataType, &extent, &str_name, &str_desc); if (kimerror) - error->all(FLERR,"KIM GetParameterMetadata returned error"); + error->all(FLERR,"KIM GetParameterMetadata returned error {}", kimerror); const std::string str_name_str(str_name); if (paramname == str_name_str) break; @@ -484,8 +483,8 @@ void PairKIM::coeff(int narg, char **arg) nubound = atoi(words[1].c_str()); if ((nubound < 1) || (nubound > extent) || (nlbound < 1) || (nlbound > nubound)) - error->all(FLERR,"Illegal index_range '{}-{}' for '{}' parameter with the extent of '{}'", - nlbound, nubound, paramname, extent); + error->all(FLERR,"Illegal index_range '{}-{}' for '{}' parameter with the extent " + "of '{}'", nlbound, nubound, paramname, extent); } else { nlbound = atoi(argtostr.c_str()); @@ -975,7 +974,7 @@ void PairKIM::set_lmps_flags() // determine if running with pair hybrid if (force->pair_match("hybrid",0)) - error->all(FLERR,"pair_kim does not support hybrid"); + error->all(FLERR,"Pair style must not be used as a hybrid sub-style"); const std::string unit_style_str(update->unit_style); @@ -1018,10 +1017,9 @@ void PairKIM::set_lmps_flags() } else if ((unit_style_str == "lj") || (unit_style_str == "micro") || (unit_style_str == "nano")) { - error->all(FLERR,"LAMMPS unit_style {} not supported " - "by KIM models", unit_style_str); + error->all(FLERR,"LAMMPS unit_style {} not supported by KIM models", unit_style_str); } else { - error->all(FLERR,"Unknown unit_style"); + error->all(FLERR,"Unknown unit_style {}", unit_style_str); } } @@ -1100,7 +1098,8 @@ void PairKIM::set_kim_model_has_flags() if (comm->me == 0) { if (KIM_SupportStatus_Equal(kim_model_support_for_energy, KIM_SUPPORT_STATUS_notSupported)) - error->warning(FLERR,"KIM Model does not provide 'partialEnergy'; Potential energy will be zero"); + error->warning(FLERR,"KIM Model does not provide 'partialEnergy'; " + "Potential energy will be zero"); if (KIM_SupportStatus_Equal(kim_model_support_for_forces, KIM_SUPPORT_STATUS_notSupported)) diff --git a/src/KOKKOS/fix_property_atom_kokkos.cpp b/src/KOKKOS/fix_property_atom_kokkos.cpp index 2c515510ee..1de07b39dc 100644 --- a/src/KOKKOS/fix_property_atom_kokkos.cpp +++ b/src/KOKKOS/fix_property_atom_kokkos.cpp @@ -24,8 +24,6 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum{MOLECULE,CHARGE,RMASS,INTEGER,DOUBLE}; - /* ---------------------------------------------------------------------- */ FixPropertyAtomKokkos::FixPropertyAtomKokkos(LAMMPS *lmp, int narg, char **arg) : @@ -44,33 +42,45 @@ FixPropertyAtomKokkos::FixPropertyAtomKokkos(LAMMPS *lmp, int narg, char **arg) void FixPropertyAtomKokkos::grow_arrays(int nmax) { - for (int m = 0; m < nvalue; m++) { - if (styles[m] == MOLECULE) { + for (int nv = 0; nv < nvalue; nv++) { + if (styles[nv] == MOLECULE) { memory->grow(atom->molecule,nmax,"atom:molecule"); size_t nbytes = (nmax-nmax_old) * sizeof(tagint); memset(&atom->molecule[nmax_old],0,nbytes); - } else if (styles[m] == CHARGE) { + } else if (styles[nv] == CHARGE) { memory->grow(atom->q,nmax,"atom:q"); size_t nbytes = (nmax-nmax_old) * sizeof(double); memset(&atom->q[nmax_old],0,nbytes); - } else if (styles[m] == RMASS) { + } else if (styles[nv] == RMASS) { memory->grow(atom->rmass,nmax,"atom:rmass"); size_t nbytes = (nmax-nmax_old) * sizeof(double); memset(&atom->rmass[nmax_old],0,nbytes); - } else if (styles[m] == INTEGER) { - memory->grow(atom->ivector[index[m]],nmax,"atom:ivector"); + } else if (styles[nv] == TEMPERATURE) { + memory->grow(atom->temperature, nmax, "atom:temperature"); + size_t nbytes = (nmax - nmax_old) * sizeof(double); + memset(&atom->temperature[nmax_old], 0, nbytes); + } else if (styles[nv] == HEATFLOW) { + memory->grow(atom->heatflow, nmax, "atom:heatflow"); + size_t nbytes = (nmax - nmax_old) * sizeof(double); + memset(&atom->heatflow[nmax_old], 0, nbytes); + } else if (styles[nv] == IVEC) { + memory->grow(atom->ivector[index[nv]],nmax,"atom:ivector"); size_t nbytes = (nmax-nmax_old) * sizeof(int); - memset(&atom->ivector[index[m]][nmax_old],0,nbytes); - } else if (styles[m] == DOUBLE) { + memset(&atom->ivector[index[nv]][nmax_old],0,nbytes); + } else if (styles[nv] == DVEC) { atomKK->sync(Device,DVECTOR_MASK); memoryKK->grow_kokkos(atomKK->k_dvector,atomKK->dvector,atomKK->k_dvector.extent(0),nmax, "atom:dvector"); atomKK->modified(Device,DVECTOR_MASK); - //memory->grow(atom->dvector[index[m]],nmax,"atom:dvector"); - //size_t nbytes = (nmax-nmax_old) * sizeof(double); - //memset(&atom->dvector[index[m]][nmax_old],0,nbytes); + } else if (styles[nv] == IARRAY) { + memory->grow(atom->iarray[index[nv]], nmax, cols[nv], "atom:iarray"); + size_t nbytes = (size_t) (nmax - nmax_old) * cols[nv] * sizeof(int); + if (nbytes) memset(&atom->iarray[index[nv]][nmax_old][0], 0, nbytes); + } else if (styles[nv] == DARRAY) { + memory->grow(atom->darray[index[nv]], nmax, cols[nv], "atom:darray"); + size_t nbytes = (size_t) (nmax - nmax_old) * cols[nv] * sizeof(double); + if (nbytes) memset(&atom->darray[index[nv]][nmax_old][0], 0, nbytes); } } - nmax_old = nmax; } diff --git a/src/KOKKOS/fix_rx_kokkos.cpp b/src/KOKKOS/fix_rx_kokkos.cpp index 963fb32a1a..ac1821e142 100644 --- a/src/KOKKOS/fix_rx_kokkos.cpp +++ b/src/KOKKOS/fix_rx_kokkos.cpp @@ -74,23 +74,17 @@ FixRxKokkos::FixRxKokkos(LAMMPS *lmp, int narg, char **arg) : datamask_modify = EMPTY_MASK; k_error_flag = DAT::tdual_int_scalar("FixRxKokkos::k_error_flag"); - - //printf("Inside FixRxKokkos::FixRxKokkos\n"); } template FixRxKokkos::~FixRxKokkos() { - //printf("Inside FixRxKokkos::~FixRxKokkos copymode= %d\n", copymode); if (copymode) return; if (localTempFlag) memoryKK->destroy_kokkos(k_dpdThetaLocal, dpdThetaLocal); memoryKK->destroy_kokkos(k_sumWeights, sumWeights); - //memoryKK->destroy_kokkos(k_sumWeights); - - //delete [] scratchSpace; memoryKK->destroy_kokkos(d_scratchSpace); memoryKK->destroy_kokkos(k_cutsq); @@ -147,7 +141,6 @@ void FixRxKokkos::init() template void FixRxKokkos::init_list(int, class NeighList* ptr) { - //printf("Inside FixRxKokkos::init_list\n"); this->list = ptr; } @@ -551,7 +544,6 @@ void FixRxKokkos::k_rkf45(const int neq, const double t_stop, Vector nfe += 6; if (maxIters && nit > maxIters) { - //fprintf(stderr,"atom[%d] took too many iterations in rkf45 %d %e %e\n", id, nit, t, t_stop); counter.nFails ++; break; // We should set an error here so that the solution is not used! @@ -562,8 +554,6 @@ void FixRxKokkos::k_rkf45(const int neq, const double t_stop, Vector counter.nSteps += nst; counter.nIters += nit; counter.nFuncs += nfe; - - //printf("id= %d nst= %d nit= %d\n", id, nst, nit); } /* ---------------------------------------------------------------------- */ @@ -733,9 +723,6 @@ int FixRxKokkos::rkf45_h0(const int neq, const double t, const doubl yddnrm = sqrt( yddnrm / double(neq) ); - //std::cout << "iter " << _iter << " hg " << hg << " y'' " << yddnrm << std::endl; - //std::cout << "ydot " << ydot[neq-1] << std::endl; - // should we accept this? if (hnew_is_ok || iter == max_iters) { hnew = hg; @@ -760,8 +747,6 @@ int FixRxKokkos::rkf45_h0(const int neq, const double t, const doubl hnew_is_ok = true; } - //printf("iter=%d, yddnrw=%e, hnew=%e, hmin=%e, hmax=%e\n", iter, yddnrm, hnew, hmin, hmax); - hg = hnew; iter ++; } @@ -805,13 +790,9 @@ void FixRxKokkos::rkf45(const int neq, const double t_stop, double * double t = 0.0; if (h < h_min) { - //fprintf(stderr,"hin not implemented yet\n"); - //exit(-1); nfe = rkf45_h0 (neq, t, t_stop, h_min, h_max, h, y, rwork, v_param); } - //printf("t= %e t_stop= %e h= %e\n", t, t_stop, h); - // Integrate until we reach the end time. while (fabs(t - t_stop) > tround) { double *yout = rwork; @@ -865,7 +846,6 @@ void FixRxKokkos::rkf45(const int neq, const double t_stop, double * nfe += 6; if (maxIters && nit > maxIters) { - //fprintf(stderr,"atom[%d] took too many iterations in rkf45 %d %e %e\n", id, nit, t, t_stop); counter.nFails ++; break; // We should set an error here so that the solution is not used! @@ -876,8 +856,6 @@ void FixRxKokkos::rkf45(const int neq, const double t_stop, double * counter.nSteps += nst; counter.nIters += nit; counter.nFuncs += nfe; - - //printf("id= %d nst= %d nit= %d\n", id, nst, nit); } /* ---------------------------------------------------------------------- */ @@ -902,9 +880,6 @@ int FixRxKokkos::rhs_dense(double /*t*/, const double *y, double *dy double *rxnRateLaw = userData->rxnRateLaw; double *kFor = userData->kFor; - //const double VDPD = domain->xprd * domain->yprd * domain->zprd / atom->natoms; - //const int nspecies = atom->nspecies_dpd; - for (int ispecies=0; ispecies::rhs_dense(double /*t*/, const double *y, double *dy for (int ispecies=0; ispecies::rhs_dense(double /*t*/, const double *y, double *dy for (int jrxn=0; jrxn::rhs_sparse(double /*t*/, const double *y, double *d { UserRHSData *userData = (UserRHSData *) v_params; - //const double VDPD = domain->xprd * domain->yprd * domain->zprd / atom->natoms; - #define kFor (userData->kFor) #define kRev (nullptr) #define rxnRateLaw (userData->rxnRateLaw) @@ -964,7 +935,6 @@ int FixRxKokkos::rhs_sparse(double /*t*/, const double *y, double *d for (int kk = 1; kk < maxReactants; ++kk) { const int k = nuk(i,kk); if (k == SparseKinetics_invalidIndex) break; - //if (k != SparseKinetics_invalidIndex) rxnRateLawForward *= powint( conc[k], inu(i,kk) ); } } else { @@ -972,7 +942,6 @@ int FixRxKokkos::rhs_sparse(double /*t*/, const double *y, double *d for (int kk = 1; kk < maxReactants; ++kk) { const int k = nuk(i,kk); if (k == SparseKinetics_invalidIndex) break; - //if (k != SparseKinetics_invalidIndex) rxnRateLawForward *= pow( conc[k], nu(i,kk) ); } } @@ -991,7 +960,6 @@ int FixRxKokkos::rhs_sparse(double /*t*/, const double *y, double *d for (int kk = 1; kk < maxReactants; ++kk) { const int k = nuk(i,kk); if (k == SparseKinetics_invalidIndex) break; - //if (k != SparseKinetics_invalidIndex) dydt[k] -= nu(i,kk) * rxnRateLaw[i]; } @@ -1000,7 +968,6 @@ int FixRxKokkos::rhs_sparse(double /*t*/, const double *y, double *d for (int kk = maxReactants+1; kk < maxSpecies; ++kk) { const int k = nuk(i,kk); if (k == SparseKinetics_invalidIndex) break; - //if (k != SparseKinetics_invalidIndex) dydt[k] += nu(i,kk) * rxnRateLaw[i]; } } @@ -1048,9 +1015,6 @@ int FixRxKokkos::k_rhs_dense(double /*t*/, const VectorType& y, Vect #define rxnRateLaw (userData.rxnRateLaw) #define kFor (userData.kFor ) - //const double VDPD = domain->xprd * domain->yprd * domain->zprd / atom->natoms; - //const int nspecies = atom->nspecies_dpd; - for (int ispecies=0; ispecies::k_rhs_dense(double /*t*/, const VectorType& y, Vect for (int ispecies=0; ispecies::k_rhs_dense(double /*t*/, const VectorType& y, Vect for (int jrxn=0; jrxn::k_rhs_sparse(double /*t*/, const VectorType& y, Vec for (int kk = 1; kk < maxReactants; ++kk) { const int k = nuk(i,kk); if (k == SparseKinetics_invalidIndex) break; - //if (k != SparseKinetics_invalidIndex) rxnRateLawForward *= powint( conc[k], inu(i,kk) ); } } else { @@ -1119,7 +1080,6 @@ int FixRxKokkos::k_rhs_sparse(double /*t*/, const VectorType& y, Vec for (int kk = 1; kk < maxReactants; ++kk) { const int k = nuk(i,kk); if (k == SparseKinetics_invalidIndex) break; - //if (k != SparseKinetics_invalidIndex) rxnRateLawForward *= pow( conc[k], nu(i,kk) ); } } @@ -1138,7 +1098,6 @@ int FixRxKokkos::k_rhs_sparse(double /*t*/, const VectorType& y, Vec for (int kk = 1; kk < maxReactants; ++kk) { const int k = nuk(i,kk); if (k == SparseKinetics_invalidIndex) break; - //if (k != SparseKinetics_invalidIndex) dydt[k] -= nu(i,kk) * rxnRateLaw[i]; } @@ -1147,7 +1106,6 @@ int FixRxKokkos::k_rhs_sparse(double /*t*/, const VectorType& y, Vec for (int kk = maxReactants+1; kk < maxSpecies; ++kk) { const int k = nuk(i,kk); if (k == SparseKinetics_invalidIndex) break; - //if (k != SparseKinetics_invalidIndex) dydt[k] += nu(i,kk) * rxnRateLaw[i]; } } @@ -1215,8 +1173,6 @@ void FixRxKokkos::operator()(SolverType, const int &i) const template void FixRxKokkos::create_kinetics_data() { - //printf("Inside FixRxKokkos::create_kinetics_data\n"); - memoryKK->create_kokkos( d_kineticsData.Arr, h_kineticsData.Arr, nreactions, "KineticsType::Arr"); memoryKK->create_kokkos( d_kineticsData.nArr, h_kineticsData.nArr, nreactions, "KineticsType::nArr"); memoryKK->create_kokkos( d_kineticsData.Ea, h_kineticsData.Ea, nreactions, "KineticsType::Ea"); @@ -1296,8 +1252,6 @@ void FixRxKokkos::create_kinetics_data() template void FixRxKokkos::setup_pre_force(int vflag) { - //printf("Inside FixRxKokkos::setup_pre_force restartFlag= %d\n", my_restartFlag); - if (my_restartFlag) my_restartFlag = 0; else @@ -1309,8 +1263,6 @@ void FixRxKokkos::setup_pre_force(int vflag) template void FixRxKokkos::pre_force(int vflag) { - //printf("Inside FixRxKokkos::pre_force localTempFlag= %d\n", localTempFlag); - this->solve_reactions( vflag, true ); } @@ -1407,8 +1359,6 @@ void FixRxKokkos::operator()(Tag_FixRxKokkos_solveSystems void FixRxKokkos::solve_reactions(const int /*vflag*/, const bool isPreForce) { - //printf("Inside FixRxKokkos::solve_reactions localTempFlag= %d isPreForce= %s\n", localTempFlag, isPreForce ? "True" : "false"); - copymode = 1; if (update_kinetics_data) @@ -1416,7 +1366,6 @@ void FixRxKokkos::solve_reactions(const int /*vflag*/, const bool is //TimerType timer_start = getTimeStamp(); - //const int nlocal = atom->nlocal; this->nlocal = atom->nlocal; const int nghost = atom->nghost; const int newton_pair = force->newton_pair; @@ -1477,9 +1426,6 @@ void FixRxKokkos::solve_reactions(const int /*vflag*/, const bool is // ... // Local references to the atomKK objects. - //typename ArrayTypes::t_efloat_1d d_dpdTheta = atomKK->k_dpdTheta.view(); - //typename ArrayTypes::t_float_2d d_dvector = atomKK->k_dvector.view(); - //typename ArrayTypes::t_int_1d d_mask = atomKK->k_mask.view(); this->d_dpdTheta = atomKK->k_dpdTheta.view(); this->d_dvector = atomKK->k_dvector.view(); this->d_mask = atomKK->k_mask.view(); @@ -1488,8 +1434,6 @@ void FixRxKokkos::solve_reactions(const int /*vflag*/, const bool is atomKK->sync( execution_space, MASK_MASK | DVECTOR_MASK | DPDTHETA_MASK ); // Set some constants outside of the parallel_for - //const double boltz = force->boltz; - //const double t_stop = update->dt; // DPD time-step and integration length. this->boltz = force->boltz; this->t_stop = update->dt; // DPD time-step and integration length. @@ -1505,13 +1449,6 @@ void FixRxKokkos::solve_reactions(const int /*vflag*/, const bool is d_diagnosticCounterPerODEnFuncs = k_diagnosticCounterPerODEnFuncs.template view(); Kokkos::parallel_for ( Kokkos::RangePolicy(0,nlocal), *this); - //Kokkos::parallel_for ( nlocal, - // LAMMPS_LAMBDA(const int i) - // { - // d_diagnosticCounterPerODEnSteps(i) = 0; - // d_diagnosticCounterPerODEnFuncs(i) = 0; - // } - // ); } // Error flag for any failures. @@ -1523,116 +1460,17 @@ void FixRxKokkos::solve_reactions(const int /*vflag*/, const bool is k_error_flag.template sync(); // Create scratch array space. - //const size_t scratchSpaceSize = (8*nspecies + 2*nreactions); this->scratchSpaceSize = (8*nspecies + 2*nreactions); - //double *scratchSpace = new double[ scratchSpaceSize * nlocal ]; - //typename ArrayTypes::t_double_1d d_scratchSpace("d_scratchSpace", scratchSpaceSize * nlocal); if (nlocal*scratchSpaceSize > d_scratchSpace.extent(0)) { memoryKK->destroy_kokkos (d_scratchSpace); memoryKK->create_kokkos (d_scratchSpace, nlocal*scratchSpaceSize, "FixRxKokkos::d_scratchSpace"); } -#if 0 - Kokkos::parallel_reduce( nlocal, LAMMPS_LAMBDA(int i, CounterType &counter) - { - if (d_mask(i) & groupbit) - { - //double *y = new double[8*nspecies]; - //double *rwork = y + nspecies; - - //StridedArrayType _y( y ); - //StridedArrayType _rwork( rwork ); - - StridedArrayType y( d_scratchSpace.data() + scratchSpaceSize * i ); - StridedArrayType rwork( &y[nspecies] ); - - //UserRHSData userData; - //userData.kFor = new double[nreactions]; - //userData.rxnRateLaw = new double[nreactions]; - - //UserRHSDataKokkos<1> userDataKokkos; - //userDataKokkos.kFor.m_data = userData.kFor; - //userDataKokkos.rxnRateLaw.m_data = userData.rxnRateLaw; - - UserRHSDataKokkos<1> userData; - userData.kFor.m_data = &( rwork[7*nspecies] ); - userData.rxnRateLaw.m_data = &( userData.kFor[ nreactions ] ); - - CounterType counter_i; - - const double theta = (localTempFlag) ? d_dpdThetaLocal(i) : d_dpdTheta(i); - - //Compute the reaction rate constants - for (int irxn = 0; irxn < nreactions; irxn++) - { - if (setRatesToZero) - userData.kFor[irxn] = 0.0; - else - { - userData.kFor[irxn] = d_kineticsData.Arr(irxn) * - pow(theta, d_kineticsData.nArr(irxn)) * - exp(-d_kineticsData.Ea(irxn) / boltz / theta); - } - } - - // Update ConcOld and initialize the ODE solution vector y[]. - for (int ispecies = 0; ispecies < nspecies; ispecies++) - { - const double tmp = d_dvector(ispecies, i); - d_dvector(ispecies+nspecies, i) = tmp; - y[ispecies] = tmp; - } - - // Solver the ODE system. - if (odeIntegrationFlag == ODE_LAMMPS_RK4) - { - k_rk4(t_stop, y, rwork, userData); - } - else if (odeIntegrationFlag == ODE_LAMMPS_RKF45) - { - k_rkf45(nspecies, t_stop, y, rwork, userData, counter_i); - - if (diagnosticFrequency == 1) - { - d_diagnosticCounterPerODEnSteps(i) = counter_i.nSteps; - d_diagnosticCounterPerODEnFuncs(i) = counter_i.nFuncs; - } - } - - // Store the solution back in dvector. - for (int ispecies = 0; ispecies < nspecies; ispecies++) - { - if (y[ispecies] < -1.0e-10) - { - //error->one(FLERR,"Computed concentration in RK solver is < -1.0e-10"); - k_error_flag.template view()() = 2; - // This should be an atomic update. - } - else if (y[ispecies] < MY_EPSILON) - y[ispecies] = 0.0; - - d_dvector(ispecies,i) = y[ispecies]; - } - - //delete [] y; - //delete [] userData.kFor; - //delete [] userData.rxnRateLaw; - - // Update the iteration statistics counter. Is this unique for each iteration? - counter += counter_i; - - } // if - } // parallel_for lambda-body - - , TotalCounters // reduction value for all iterations. - ); -#else if (setRatesToZero) Kokkos::parallel_reduce( Kokkos::RangePolicy >(0,nlocal), *this, TotalCounters); else Kokkos::parallel_reduce( Kokkos::RangePolicy >(0,nlocal), *this, TotalCounters); -#endif TimerType timer_ODE = getTimeStamp(); @@ -1652,16 +1490,8 @@ void FixRxKokkos::solve_reactions(const int /*vflag*/, const bool is atomKK->modified ( Host, DVECTOR_MASK ); - //TimerType timer_stop = getTimeStamp(); - double time_ODE = getElapsedTime(timer_localTemperature, timer_ODE); - //printf("me= %d kokkos total= %g temp= %g ode= %g comm= %g nlocal= %d nfc= %d %d\n", comm->me, - // getElapsedTime(timer_start, timer_stop), - // getElapsedTime(timer_start, timer_localTemperature), - // getElapsedTime(timer_localTemperature, timer_ODE), - // getElapsedTime(timer_ODE, timer_stop), nlocal, TotalCounters.nFuncs, TotalCounters.nSteps); - // Warn the user if a failure was detected in the ODE solver. if (TotalCounters.nFails > 0) { char sbuf[128]; @@ -1731,7 +1561,6 @@ void FixRxKokkos::odeDiagnostics() // Compute counters per dpd time-step. for (int i = 0; i < numCounters; ++i) { my_vals[i] = this->diagnosticCounter[i] / nTimes; - //printf("my sum[%d] = %f %d\n", i, my_vals[i], comm->me); } MPI_Allreduce (my_vals, sums, numCounters, MPI_DOUBLE, MPI_SUM, world); @@ -1770,7 +1599,6 @@ void FixRxKokkos::odeDiagnostics() double my_max[numCounters], my_min[numCounters]; - //const int nlocal = atom->nlocal; nlocal = atom->nlocal; HAT::t_int_1d h_mask = atomKK->k_mask.h_view; @@ -1980,9 +1808,6 @@ template template void FixRxKokkos::computeLocalTemperature() { - //typename ArrayTypes::t_x_array_randomread d_x = atomKK->k_x.view(); - //typename ArrayTypes::t_int_1d_randomread d_type = atomKK->k_type.view(); - //typename ArrayTypes::t_efloat_1d d_dpdTheta = atomKK->k_dpdTheta.view(); d_x = atomKK->k_x.view(); d_type = atomKK->k_type.view(); d_dpdTheta = atomKK->k_dpdTheta.view(); @@ -1993,21 +1818,9 @@ void FixRxKokkos::computeLocalTemperature() nlocal = atom->nlocal; const int nghost = atom->nghost; - //printf("Inside FixRxKokkos::computeLocalTemperature: %d %d %d %d %d %d %d\n", WT_FLAG, LOCAL_TEMP_FLAG, NEWTON_PAIR, (int)lmp->kokkos->neighflag, NEIGHFLAG, nlocal, nghost); - - // Pull from pairDPDE. The pairDPDEKK objects are protected so recreate here for now. - //pairDPDEKK->k_cutsq.template sync(); - //typename ArrayTypes::t_ffloat_2d d_cutsq = pairDPDEKK->k_cutsq.template view::tdual_ffloat_2d k_cutsq; - //typename ArrayTypes::t_ffloat_2d d_cutsq; - //double **h_cutsq; - { const int ntypes = atom->ntypes; - //memoryKK->create_kokkos (k_cutsq, h_cutsq, ntypes+1, ntypes+1, "pair:cutsq"); if (ntypes+1 > (int) k_cutsq.extent(0)) { memoryKK->destroy_kokkos (k_cutsq); memoryKK->create_kokkos (k_cutsq, ntypes+1, ntypes+1, "FixRxKokkos::k_cutsq"); @@ -2028,7 +1841,6 @@ void FixRxKokkos::computeLocalTemperature() // Initialize the local temperature weight array int sumWeightsCt = nlocal + (NEWTON_PAIR ? nghost : 0); - //memoryKK->create_kokkos (k_sumWeights, sumWeights, sumWeightsCt, "FixRxKokkos::sumWeights"); if (sumWeightsCt > (int)k_sumWeights.template view().extent(0)) { memoryKK->destroy_kokkos(k_sumWeights, sumWeights); memoryKK->create_kokkos (k_sumWeights, sumWeightsCt, "FixRxKokkos::sumWeights"); @@ -2036,14 +1848,6 @@ void FixRxKokkos::computeLocalTemperature() h_sumWeights = k_sumWeights.h_view; } - // Initialize the accumulator to zero ... - //Kokkos::parallel_for (sumWeightsCt, - // LAMMPS_LAMBDA(const int i) - // { - // d_sumWeights(i) = 0.0; - // } - // ); - Kokkos::parallel_for (Kokkos::RangePolicy(0, sumWeightsCt), *this); // Local list views. (This isn't working!) @@ -2051,9 +1855,6 @@ void FixRxKokkos::computeLocalTemperature() if (!list->kokkos) error->one(FLERR,"list is not a Kokkos list\n"); - //typename ArrayTypes::t_neighbors_2d d_neighbors = k_list->d_neighbors; - //typename ArrayTypes::t_int_1d d_ilist = k_list->d_ilist; - //typename ArrayTypes::t_int_1d d_numneigh = k_list->d_numneigh; d_neighbors = k_list->d_neighbors; d_ilist = k_list->d_ilist; d_numneigh = k_list->d_numneigh; @@ -2067,7 +1868,6 @@ void FixRxKokkos::computeLocalTemperature() { // Create an atomic view of sumWeights and dpdThetaLocal. Only needed // for Half/thread scenarios. - //typedef Kokkos::View< E_FLOAT*, typename DAT::t_efloat_1d::array_layout, typename KKDevice::value, Kokkos::MemoryTraits< AtomicF< NEIGHFLAG >::value> > AtomicViewType; typedef Kokkos::View< E_FLOAT*, typename DAT::t_efloat_1d::array_layout, typename KKDevice::value, Kokkos::MemoryTraits< AtomicF< NEIGHFLAG >::value> > AtomicViewType; AtomicViewType a_dpdThetaLocal = d_dpdThetaLocal; @@ -2174,8 +1974,6 @@ void FixRxKokkos::computeLocalTemperature() template int FixRxKokkos::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) { - //printf("inside FixRxKokkos::pack_forward_comm %d\n", comm->me); - HAT::t_float_2d h_dvector = atomKK->k_dvector.h_view; int m = 0; @@ -2186,9 +1984,6 @@ int FixRxKokkos::pack_forward_comm(int n, int *list, double *buf, in buf[m++] = h_dvector(ispecies+nspecies,jj); } } - - //printf("done with FixRxKokkos::pack_forward_comm %d\n", comm->me); - return m; } @@ -2197,8 +1992,6 @@ int FixRxKokkos::pack_forward_comm(int n, int *list, double *buf, in template void FixRxKokkos::unpack_forward_comm(int n, int first, double *buf) { - //printf("inside FixRxKokkos::unpack_forward_comm %d\n", comm->me); - HAT::t_float_2d h_dvector = atomKK->k_dvector.h_view; const int last = first + n ; @@ -2209,8 +2002,6 @@ void FixRxKokkos::unpack_forward_comm(int n, int first, double *buf) h_dvector(ispecies+nspecies,ii) = buf[m++]; } } - - //printf("done with FixRxKokkos::unpack_forward_comm %d\n", comm->me); } /* ---------------------------------------------------------------------- */ @@ -2218,7 +2009,6 @@ void FixRxKokkos::unpack_forward_comm(int n, int first, double *buf) template int FixRxKokkos::pack_reverse_comm(int n, int first, double *buf) { - //printf("inside FixRxKokkos::pack_reverse_comm %d %d %d\n", comm->me, first, n); // Sync the host view. k_dpdThetaLocal.template sync(); k_sumWeights. template sync(); @@ -2230,8 +2020,6 @@ int FixRxKokkos::pack_reverse_comm(int n, int first, double *buf) buf[m++] = h_dpdThetaLocal(i); buf[m++] = h_sumWeights(i); } - //printf("done with FixRxKokkos::pack_reverse_comm %d\n", comm->me); - return m; } @@ -2240,7 +2028,6 @@ int FixRxKokkos::pack_reverse_comm(int n, int first, double *buf) template void FixRxKokkos::unpack_reverse_comm(int n, int *list, double *buf) { - // printf("inside FixRxKokkos::unpack_reverse_comm %d\n", comm->me); int m = 0; for (int i = 0; i < n; i++) { const int j = list[i]; @@ -2252,8 +2039,6 @@ void FixRxKokkos::unpack_reverse_comm(int n, int *list, double *buf) // Signal that the host view has been modified. k_dpdThetaLocal.template modify(); k_sumWeights. template modify(); - - // printf("done with FixRxKokkos::unpack_reverse_comm %d\n", comm->me); } namespace LAMMPS_NS { diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp index ee824d291b..a3d4257cc2 100644 --- a/src/MANYBODY/pair_eam.cpp +++ b/src/MANYBODY/pair_eam.cpp @@ -852,8 +852,7 @@ double PairEAM::single(int i, int j, int itype, int jtype, /* ---------------------------------------------------------------------- */ -int PairEAM::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) +int PairEAM::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) { int i,j,m; diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp index 8dd97250a3..0b4ec9a638 100644 --- a/src/fix_adapt.cpp +++ b/src/fix_adapt.cpp @@ -48,9 +48,9 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), nadapt(0), anypair(0), anybond(0), anyangle(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr) { - if (narg < 5) error->all(FLERR,"Illegal fix adapt command"); + if (narg < 5) utils::missing_cmd_args(FLERR,"fix adapt", error); nevery = utils::inumeric(FLERR,arg[3],false,lmp); - if (nevery < 0) error->all(FLERR,"Illegal fix adapt command"); + if (nevery < 0) error->all(FLERR,"Illegal fix adapt every value {}", nevery); dynamic_group_allow = 1; create_attribute = 1; @@ -62,29 +62,29 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : int iarg = 4; while (iarg < narg) { if (strcmp(arg[iarg],"pair") == 0) { - if (iarg+6 > narg) error->all(FLERR,"Illegal fix adapt command"); + if (iarg+6 > narg) utils::missing_cmd_args(FLERR,"fix adapt pair", error); nadapt++; iarg += 6; } else if (strcmp(arg[iarg],"kspace") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR,"fix adapt kspace", error); nadapt++; iarg += 2; } else if (strcmp(arg[iarg],"atom") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command"); + if (iarg+3 > narg) utils::missing_cmd_args(FLERR,"fix adapt atom", error); nadapt++; iarg += 3; } else if (strcmp(arg[iarg],"bond") == 0) { - if (iarg+5 > narg) error->all(FLERR,"Illegal fix adapt command"); + if (iarg+5 > narg) utils::missing_cmd_args(FLERR,"fix adapt bond", error); nadapt++; iarg += 5; } else if (strcmp(arg[iarg],"angle") == 0) { - if (iarg+5 > narg) error->all(FLERR,"Illegal fix adapt command"); + if (iarg+5 > narg) utils::missing_cmd_args(FLERR,"fix adapt angle", error); nadapt++; iarg += 5; } else break; } - if (nadapt == 0) error->all(FLERR,"Illegal fix adapt command"); + if (nadapt == 0) error->all(FLERR,"Nothing to adapt in fix adapt command"); adapt = new Adapt[nadapt]; // parse keywords @@ -96,7 +96,6 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : iarg = 4; while (iarg < narg) { if (strcmp(arg[iarg],"pair") == 0) { - if (iarg+6 > narg) error->all(FLERR,"Illegal fix adapt command"); adapt[nadapt].which = PAIR; adapt[nadapt].pair = nullptr; adapt[nadapt].pstyle = utils::strdup(arg[iarg+1]); @@ -107,12 +106,11 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : adapt[nadapt].jlo,adapt[nadapt].jhi,error); if (utils::strmatch(arg[iarg+5],"^v_")) { adapt[nadapt].var = utils::strdup(arg[iarg+5]+2); - } else error->all(FLERR,"Illegal fix adapt command"); + } else error->all(FLERR,"Argument #{} must be variable not {}", iarg+6, arg[iarg+5]); nadapt++; iarg += 6; } else if (strcmp(arg[iarg],"bond") == 0) { - if (iarg+5 > narg) error->all(FLERR, "Illegal fix adapt command"); adapt[nadapt].which = BOND; adapt[nadapt].bond = nullptr; adapt[nadapt].bstyle = utils::strdup(arg[iarg+1]); @@ -121,12 +119,11 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : adapt[nadapt].ilo,adapt[nadapt].ihi,error); if (utils::strmatch(arg[iarg+4],"^v_")) { adapt[nadapt].var = utils::strdup(arg[iarg+4]+2); - } else error->all(FLERR,"Illegal fix adapt command"); + } else error->all(FLERR,"Argument #{} must be variable not {}", iarg+5, arg[iarg+4]); nadapt++; iarg += 5; } else if (strcmp(arg[iarg],"angle") == 0) { - if (iarg+5 > narg) error->all(FLERR, "Illegal fix adapt command"); adapt[nadapt].which = ANGLE; adapt[nadapt].angle = nullptr; adapt[nadapt].astyle = utils::strdup(arg[iarg+1]); @@ -135,21 +132,19 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : adapt[nadapt].ilo,adapt[nadapt].ihi,error); if (utils::strmatch(arg[iarg+4],"^v_")) { adapt[nadapt].var = utils::strdup(arg[iarg+4]+2); - } else error->all(FLERR,"Illegal fix adapt command"); + } else error->all(FLERR,"Argument #{} must be variable not {}", iarg+5, arg[iarg+4]); nadapt++; iarg += 5; } else if (strcmp(arg[iarg],"kspace") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); adapt[nadapt].which = KSPACE; if (utils::strmatch(arg[iarg+1],"^v_")) { adapt[nadapt].var = utils::strdup(arg[iarg+1]+2); - } else error->all(FLERR,"Illegal fix adapt command"); + } else error->all(FLERR,"Argument #{} must be variable not {}", iarg+2, arg[iarg+1]); nadapt++; iarg += 2; } else if (strcmp(arg[iarg],"atom") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command"); adapt[nadapt].which = ATOM; if (strcmp(arg[iarg+1],"diameter") == 0 || strcmp(arg[iarg+1],"diameter/disc") == 0) { @@ -160,10 +155,10 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg+1],"charge") == 0) { adapt[nadapt].atomparam = CHARGE; chgflag = 1; - } else error->all(FLERR,"Illegal fix adapt command"); + } else error->all(FLERR,"Unsupported per-atom property {} for fix adapt", arg[iarg+1]); if (utils::strmatch(arg[iarg+2],"^v_")) { adapt[nadapt].var = utils::strdup(arg[iarg+2]+2); - } else error->all(FLERR,"Illegal fix adapt command"); + } else error->all(FLERR,"Argument #{} must be variable not {}", iarg+3, arg[iarg+2]); nadapt++; iarg += 3; } else break; @@ -177,18 +172,18 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : while (iarg < narg) { if (strcmp(arg[iarg],"reset") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR,"fix adapt reset", error); resetflag = utils::logical(FLERR,arg[iarg+1],false,lmp); iarg += 2; } else if (strcmp(arg[iarg],"scale") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR,"fix adapt scale", error); scaleflag = utils::logical(FLERR,arg[iarg+1],false,lmp); iarg += 2; } else if (strcmp(arg[iarg],"mass") == 0) { - if (iarg+2 > narg)error->all(FLERR,"Illegal fix adapt command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR,"fix adapt mass", error); massflag = utils::logical(FLERR,arg[iarg+1],false,lmp); iarg += 2; - } else error->all(FLERR,"Illegal fix adapt command"); + } else error->all(FLERR,"Unknown fix adapt keyword {}", arg[iarg]); } // if scaleflag set with diameter or charge adaptation, @@ -202,22 +197,19 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : int n = atom->ntypes; for (int m = 0; m < nadapt; m++) - if (adapt[m].which == PAIR) - memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig"); + if (adapt[m].which == PAIR) memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig"); // allocate bond style arrays: n = atom->nbondtypes; for (int m = 0; m < nadapt; ++m) - if (adapt[m].which == BOND) - memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig"); + if (adapt[m].which == BOND) memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig"); // allocate angle style arrays: n = atom->nbondtypes; for (int m = 0; m < nadapt; ++m) - if (adapt[m].which == ANGLE) - memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig"); + if (adapt[m].which == ANGLE) memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig"); } /* ---------------------------------------------------------------------- */ @@ -324,7 +316,7 @@ void FixAdapt::init() if (group->dynamic[igroup]) for (i = 0; i < nadapt; i++) if (adapt[i].which == ATOM) - error->all(FLERR,"Cannot use dynamic group with fix adapt atom"); + error->all(FLERR,"Cannot use dynamic group {} with fix adapt atom", group->names[igroup]); // setup and error checks @@ -337,9 +329,9 @@ void FixAdapt::init() ad->ivar = input->variable->find(ad->var); if (ad->ivar < 0) - error->all(FLERR,"Variable name for fix adapt does not exist"); + error->all(FLERR,"Variable name {} for fix adapt does not exist", ad->var); if (!input->variable->equalstyle(ad->ivar)) - error->all(FLERR,"Variable for fix adapt is invalid style"); + error->all(FLERR,"Variable {} for fix adapt is invalid style", ad->var); if (ad->which == PAIR) { anypair = 1; @@ -368,13 +360,13 @@ void FixAdapt::init() void *ptr = ad->pair->extract(ad->pparam,ad->pdim); if (ptr == nullptr) - error->all(FLERR,"Fix adapt pair style param not supported"); + error->all(FLERR,"Fix adapt pair style {} param {} not supported", ad->pstyle, ad->pparam); // for pair styles only parameters that are 2-d arrays in atom types or // scalars are supported if (ad->pdim != 2 && ad->pdim != 0) - error->all(FLERR,"Fix adapt pair style param is not compatible"); + error->all(FLERR,"Pair style parameter {} is not compatible with fix adapt", ad->pparam); if (ad->pdim == 2) ad->array = (double **) ptr; if (ad->pdim == 0) ad->scalar = (double *) ptr; @@ -402,12 +394,12 @@ void FixAdapt::init() if (ad->bond == nullptr) ad->bond = force->bond_match(bstyle); if (ad->bond == nullptr ) - error->all(FLERR,"Fix adapt bond style does not exist"); + error->all(FLERR,"Fix adapt bond style {} does not exist", bstyle); void *ptr = ad->bond->extract(ad->bparam,ad->bdim); if (ptr == nullptr) - error->all(FLERR,"Fix adapt bond style param not supported"); + error->all(FLERR,"Fix adapt bond style parameter {} not supported", ad->bparam); // for bond styles, use a vector @@ -428,12 +420,12 @@ void FixAdapt::init() if (ad->angle == nullptr) ad->angle = force->angle_match(astyle); if (ad->angle == nullptr ) - error->all(FLERR,"Fix adapt angle style does not exist"); + error->all(FLERR,"Fix adapt angle style {} does not exist", astyle); void *ptr = ad->angle->extract(ad->aparam,ad->adim); if (ptr == nullptr) - error->all(FLERR,"Fix adapt angle style param not supported"); + error->all(FLERR,"Fix adapt angle style parameter {} not supported", ad->aparam); // for angle styles, use a vector @@ -446,7 +438,7 @@ void FixAdapt::init() } else if (ad->which == KSPACE) { if (force->kspace == nullptr) - error->all(FLERR,"Fix adapt kspace style does not exist"); + error->all(FLERR,"Fix adapt expected a kspace style but none was defined"); kspace_scale = (double *) force->kspace->extract("scale"); } else if (ad->which == ATOM) { diff --git a/src/fix_property_atom.cpp b/src/fix_property_atom.cpp index 4b43feca74..994b4f0f19 100644 --- a/src/fix_property_atom.cpp +++ b/src/fix_property_atom.cpp @@ -25,8 +25,6 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum { MOLECULE, CHARGE, RMASS, TEMPERATURE, HEATFLOW, IVEC, DVEC, IARRAY, DARRAY }; - /* ---------------------------------------------------------------------- */ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) : diff --git a/src/fix_property_atom.h b/src/fix_property_atom.h index 7d9f50b75e..92497d6188 100644 --- a/src/fix_property_atom.h +++ b/src/fix_property_atom.h @@ -27,11 +27,12 @@ namespace LAMMPS_NS { class FixPropertyAtom : public Fix { public: FixPropertyAtom(class LAMMPS *, int, char **); - ~FixPropertyAtom() override; int setmask() override; void init() override; + enum { MOLECULE, CHARGE, RMASS, TEMPERATURE, HEATFLOW, IVEC, DVEC, IARRAY, DARRAY }; + void read_data_section(char *, int, char *, tagint) override; bigint read_data_skip_lines(char *) override; void write_data_section_size(int, int &, int &) override; @@ -55,10 +56,10 @@ class FixPropertyAtom : public Fix { int nvalue, border; int molecule_flag, q_flag, rmass_flag; // flags for specific fields int temperature_flag, heatflow_flag; - int *styles; // style of each value, see enum - int *index; // indices into atom custom data structs - int *cols; // columns per value, for arrays - char *astyle; // atom style at instantiation + int *styles; // style of each value, see enum + int *index; // indices into atom custom data structs + int *cols; // columns per value, for arrays + char *astyle; // atom style at instantiation int values_peratom; // # of values per atom, including multiple for arrays int nmax_old; // length of peratom arrays the last time they grew diff --git a/tools/fep/fdti.py b/tools/fep/fdti.py index 1f3f86f001..42f0a5fd2a 100755 --- a/tools/fep/fdti.py +++ b/tools/fep/fdti.py @@ -10,7 +10,7 @@ if len(sys.argv) < 3: print("usage: fdti.py temperature hderiv < out.fep") sys.exit() -rt = 0.008314 / 4.184 * float(sys.argv[1]) +rt = 0.008314 / 4.184 * float(sys.argv[1]) # in kcal/mol hderiv = float(sys.argv[2]) line = sys.stdin.readline() @@ -23,7 +23,7 @@ if len(tok) == 4: v = float(tok[3]) lo = -rt * math.log(float(tok[2]) / v) -i = 0 +i = 1 sum = 0.0 for line in sys.stdin: tok = line.split() diff --git a/tools/fep/fep.py b/tools/fep/fep.py index 22ebaf38c6..c7656c5ef5 100755 --- a/tools/fep/fep.py +++ b/tools/fep/fep.py @@ -9,7 +9,7 @@ if len(sys.argv) < 2: print("usage: fep.py temperature < out.fep") sys.exit() -rt = 0.008314 / 4.184 * float(sys.argv[1]) +rt = 0.008314 / 4.184 * float(sys.argv[1]) # in kcal/mol v = 1.0 sum = 0.0 diff --git a/tools/fep/nti.py b/tools/fep/nti.py index ffe47fb908..cd426163bc 100755 --- a/tools/fep/nti.py +++ b/tools/fep/nti.py @@ -19,7 +19,7 @@ while line.startswith("#"): tok = line.split() lo = float(tok[1]) -i = 0 +i = 1 sum = 0.0 for line in sys.stdin: tok = line.split() diff --git a/tools/python/pizza/dump.py b/tools/python/pizza/dump.py index 3dabc9b1b4..d57bdf63a3 100644 --- a/tools/python/pizza/dump.py +++ b/tools/python/pizza/dump.py @@ -657,11 +657,14 @@ class dump: atoms = snap.atoms nvalues = len(atoms[0]) + keys = dict() + for pair in self.names.items(): + keys[pair[1]] = pair[0] for i in range(snap.natoms): if not snap.aselect[i]: continue line = "" for j in range(nvalues): - if (j < 2): + if keys[j] == 'id' or keys[j] == 'type' or keys[j] == 'mol': line += str(int(atoms[i][j])) + " " else: line += str(atoms[i][j]) + " " diff --git a/unittest/commands/test_kim_commands.cpp b/unittest/commands/test_kim_commands.cpp index 988b0db69a..2e6e758c76 100644 --- a/unittest/commands/test_kim_commands.cpp +++ b/unittest/commands/test_kim_commands.cpp @@ -92,14 +92,14 @@ TEST_F(KimCommandsTest, kim_interactions) { if (!LAMMPS::is_installed_pkg("KIM")) GTEST_SKIP(); - TEST_FAILURE(".*ERROR: Illegal 'kim interactions' command.*", command("kim interactions");); + TEST_FAILURE(".*ERROR: Illegal kim interactions command: missing argument.*", + command("kim interactions");); BEGIN_HIDE_OUTPUT(); command("kim init LennardJones_Ar real"); END_HIDE_OUTPUT(); - TEST_FAILURE(".*ERROR: Must use 'kim interactions' command " - "after simulation box is defined.*", + TEST_FAILURE(".*ERROR: Use of 'kim interactions' before simulation box is defined.*", command("kim interactions Ar");); BEGIN_HIDE_OUTPUT(); @@ -410,11 +410,11 @@ TEST_F(KimCommandsTest, kim_property) "3 >= 3.6 support.*", command("kim property");); } else { - TEST_FAILURE(".*ERROR: Invalid 'kim property' command.*", command("kim property");); - TEST_FAILURE(".*ERROR: Invalid 'kim property' command.*", command("kim property create");); - TEST_FAILURE(".*ERROR: Incorrect arguments in 'kim property' command." - "\n'kim property create/destroy/modify/remove/dump' " - "is mandatory.*", + TEST_FAILURE(".*ERROR: Illegal kim property command: missing argument.*", + command("kim property");); + TEST_FAILURE(".*ERROR: Illegal kim property command: missing argument.*", + command("kim property create");); + TEST_FAILURE(".*ERROR: Incorrect first argument unknown to 'kim property' command.*", command("kim property unknown 1 atomic-mass");); } #if defined(KIM_EXTRA_UNITTESTS)