ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #3693 from akohlmey/collected-small-changes

Browse Source 
Collected small changes and fixes
This commit is contained in:
Axel Kohlmeyer
2023-03-22 21:56:22 -04:00
committed by GitHub
parent 8149386aec 872a4e2999
commit 12e0131abf
24 changed files with 481 additions and 755 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
cmake/Modules/Packages/KIM.cmake
View File

@ -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)

2
src/.gitignore vendored
View File

@ -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

43
src/KIM/fix_store_kim.cpp
View File

@ -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;
}

58
src/KIM/kim_command.cpp
View File

@ -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);
}

64
src/KIM/kim_init.cpp
View File

@ -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<FixStoreKIM *>(modify->fix[ifix]);
if (modify->get_fix_by_id("KIM_MODEL_STORE")) modify->delete_fix("KIM_MODEL_STORE");
auto fix_store = dynamic_cast<FixStoreKIM *>(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);

49
src/KIM/kim_interactions.cpp
View File

@ -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<FixStoreKIM *>(modify->fix[ifix]);
auto fix_store = dynamic_cast<FixStoreKIM *>(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<std::string> 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;

76
src/KIM/kim_param.cpp
View File

@ -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<FixStoreKIM *>(modify->fix[ifix]);
KIM_SimulatorModel *simulatorModel =
reinterpret_cast<KIM_SimulatorModel *>(
fix_store->getptr("simulator_model"));
auto fix_store = dynamic_cast<FixStoreKIM *>(modify->get_fix_by_id("KIM_MODEL_STORE"));
if (fix_store) {
auto *simulatorModel = reinterpret_cast<KIM_SimulatorModel *>(
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));
}

2
src/KIM/kim_param.h
View File

@ -67,7 +67,5 @@ class KimParam : protected Pointers {
KimParam(class LAMMPS *lmp);
void command(int, char **);
};
} // namespace LAMMPS_NS
#endif

68
src/KIM/kim_property.cpp
View File

@ -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;

41
src/KIM/kim_query.cpp
View File

@ -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<FixStoreKIM *>(modify->fix[ifix]);
auto fix_store = dynamic_cast<FixStoreKIM *>(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
}

425
src/KIM/kim_units.cpp
View File

@ -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 <cmath>
#include <map>
#include <string>
#include <utility>
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<string, sys_type> system_dic;
map<string, unit_type> unit_dic;
map<unit_type, units> units_real_dic;
map<unit_type, units> units_metal_dic;
map<unit_type, units> units_si_dic;
map<unit_type, units> units_cgs_dic;
map<unit_type, units> units_electron_dic;
map<unit_type, units> units_micro_dic;
map<unit_type, units> units_nano_dic;
static map<string, sys_type> system_dic;
static map<string, unit_type> unit_dic;
static map<unit_type, units> units_real_dic;
static map<unit_type, units> units_metal_dic;
static map<unit_type, units> units_si_dic;
static map<unit_type, units> units_cgs_dic;
static map<unit_type, units> units_electron_dic;
static map<unit_type, units> units_micro_dic;
static map<unit_type, units> 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<units, map<units, double> > 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<units, map<units, double> > 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<units, map<units, double> > 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<units, map<units, double> > 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<units, map<units, double> > 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<units, map<units, double> > 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<units, map<units, double> > 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<units, map<units, double> > 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<units, map<units, double> > 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<units, map<units, double> > 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<units, map<units, double> > 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<units, map<units, double> > 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<units, map<units, double> > 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<units, map<units, double> > 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<string, unit_type>::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<string, sys_type>::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<string, sys_type>::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<string, unit_type>::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<string, sys_type>::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<string, sys_type>::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);
}
}

17
src/KIM/kim_units.h
View File

@ -56,7 +56,20 @@
#define LMP_KIM_UNITS_H
#include <string>
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

21
src/KIM/pair_kim.cpp
View File

@ -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))

38
src/KOKKOS/fix_property_atom_kokkos.cpp
View File

@ -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;
}

215
src/KOKKOS/fix_rx_kokkos.cpp
View File

@ -74,23 +74,17 @@ FixRxKokkos<DeviceType>::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 <typename DeviceType>
FixRxKokkos<DeviceType>::~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<DeviceType>::init()
template <typename DeviceType>
void FixRxKokkos<DeviceType>::init_list(int, class NeighList* ptr)
{
//printf("Inside FixRxKokkos::init_list\n");
this->list = ptr;
}
@ -551,7 +544,6 @@ void FixRxKokkos<DeviceType>::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<DeviceType>::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<DeviceType>::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<DeviceType>::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<DeviceType>::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<DeviceType>::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<DeviceType>::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<DeviceType>::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<nspecies; ispecies++)
dydt[ispecies] = 0.0;
@ -915,7 +890,6 @@ int FixRxKokkos<DeviceType>::rhs_dense(double /*t*/, const double *y, double *dy
for (int ispecies=0; ispecies<nspecies; ispecies++) {
const double concentration = y[ispecies]/VDPD;
rxnRateLawForward *= pow( concentration, d_kineticsData.stoichReactants(jrxn,ispecies) );
//rxnRateLawForward *= pow(concentration,stoichReactants[jrxn][ispecies]);
}
rxnRateLaw[jrxn] = rxnRateLawForward;
}
@ -925,7 +899,6 @@ int FixRxKokkos<DeviceType>::rhs_dense(double /*t*/, const double *y, double *dy
for (int jrxn=0; jrxn<nreactions; jrxn++)
{
dydt[ispecies] += d_kineticsData.stoich(jrxn,ispecies) *VDPD*rxnRateLaw[jrxn];
//dydt[ispecies] += stoich[jrxn][ispecies]*VDPD*rxnRateLaw[jrxn];
}
return 0;
@ -938,8 +911,6 @@ int FixRxKokkos<DeviceType>::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<DeviceType>::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<DeviceType>::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<DeviceType>::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<DeviceType>::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<DeviceType>::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<nspecies; ispecies++)
dydt[ispecies] = 0.0;
@ -1061,7 +1025,6 @@ int FixRxKokkos<DeviceType>::k_rhs_dense(double /*t*/, const VectorType& y, Vect
for (int ispecies=0; ispecies<nspecies; ispecies++) {
const double concentration = y[ispecies]/VDPD;
rxnRateLawForward *= pow( concentration, d_kineticsData.stoichReactants(jrxn,ispecies) );
//rxnRateLawForward *= pow(concentration,stoichReactants[jrxn][ispecies]);
}
rxnRateLaw[jrxn] = rxnRateLawForward;
}
@ -1071,7 +1034,6 @@ int FixRxKokkos<DeviceType>::k_rhs_dense(double /*t*/, const VectorType& y, Vect
for (int jrxn=0; jrxn<nreactions; jrxn++)
{
dydt[ispecies] += d_kineticsData.stoich(jrxn,ispecies) *VDPD*rxnRateLaw[jrxn];
//dydt[ispecies] += stoich[jrxn][ispecies]*VDPD*rxnRateLaw[jrxn];
}
#undef rxnRateLaw
@ -1111,7 +1073,6 @@ int FixRxKokkos<DeviceType>::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<DeviceType>::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<DeviceType>::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<DeviceType>::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<DeviceType>::operator()(SolverType, const int &i) const
template <typename DeviceType>
void FixRxKokkos<DeviceType>::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<DeviceType>::create_kinetics_data()
template <typename DeviceType>
void FixRxKokkos<DeviceType>::setup_pre_force(int vflag)
{
//printf("Inside FixRxKokkos<DeviceType>::setup_pre_force restartFlag= %d\n", my_restartFlag);
if (my_restartFlag)
my_restartFlag = 0;
else
@ -1309,8 +1263,6 @@ void FixRxKokkos<DeviceType>::setup_pre_force(int vflag)
template <typename DeviceType>
void FixRxKokkos<DeviceType>::pre_force(int vflag)
{
//printf("Inside FixRxKokkos<DeviceType>::pre_force localTempFlag= %d\n", localTempFlag);
this->solve_reactions( vflag, true );
}
@ -1407,8 +1359,6 @@ void FixRxKokkos<DeviceType>::operator()(Tag_FixRxKokkos_solveSystems<ZERO_RATES
template <typename DeviceType>
void FixRxKokkos<DeviceType>::solve_reactions(const int /*vflag*/, const bool isPreForce)
{
//printf("Inside FixRxKokkos<DeviceType>::solve_reactions localTempFlag= %d isPreForce= %s\n", localTempFlag, isPreForce ? "True" : "false");
copymode = 1;
if (update_kinetics_data)
@ -1416,7 +1366,6 @@ void FixRxKokkos<DeviceType>::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<DeviceType>::solve_reactions(const int /*vflag*/, const bool is
// ...
// Local references to the atomKK objects.
//typename ArrayTypes<DeviceType>::t_efloat_1d d_dpdTheta = atomKK->k_dpdTheta.view<DeviceType>();
//typename ArrayTypes<DeviceType>::t_float_2d d_dvector = atomKK->k_dvector.view<DeviceType>();
//typename ArrayTypes<DeviceType>::t_int_1d d_mask = atomKK->k_mask.view<DeviceType>();
this->d_dpdTheta = atomKK->k_dpdTheta.view<DeviceType>();
this->d_dvector = atomKK->k_dvector.view<DeviceType>();
this->d_mask = atomKK->k_mask.view<DeviceType>();
@ -1488,8 +1434,6 @@ void FixRxKokkos<DeviceType>::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<DeviceType>::solve_reactions(const int /*vflag*/, const bool is
d_diagnosticCounterPerODEnFuncs = k_diagnosticCounterPerODEnFuncs.template view<DeviceType>();
Kokkos::parallel_for ( Kokkos::RangePolicy<DeviceType, Tag_FixRxKokkos_zeroCounterViews>(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<DeviceType>::solve_reactions(const int /*vflag*/, const bool is
k_error_flag.template sync<DeviceType>();
// 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<DeviceType>::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<double,1> _y( y );
//StridedArrayType<double,1> _rwork( rwork );
StridedArrayType<double,1> y( d_scratchSpace.data() + scratchSpaceSize * i );
StridedArrayType<double,1> 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<DeviceType>()() = 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<DeviceType, Tag_FixRxKokkos_solveSystems<true > >(0,nlocal), *this, TotalCounters);
else
Kokkos::parallel_reduce( Kokkos::RangePolicy<DeviceType, Tag_FixRxKokkos_solveSystems<false> >(0,nlocal), *this, TotalCounters);
#endif
TimerType timer_ODE = getTimeStamp();
@ -1652,16 +1490,8 @@ void FixRxKokkos<DeviceType>::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<DeviceType>::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<DeviceType>::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 <typename DeviceType>
template <int WT_FLAG, int LOCAL_TEMP_FLAG, bool NEWTON_PAIR, int NEIGHFLAG>
void FixRxKokkos<DeviceType>::computeLocalTemperature()
{
//typename ArrayTypes<DeviceType>::t_x_array_randomread d_x = atomKK->k_x.view<DeviceType>();
//typename ArrayTypes<DeviceType>::t_int_1d_randomread d_type = atomKK->k_type.view<DeviceType>();
//typename ArrayTypes<DeviceType>::t_efloat_1d d_dpdTheta = atomKK->k_dpdTheta.view<DeviceType>();
d_x = atomKK->k_x.view<DeviceType>();
d_type = atomKK->k_type.view<DeviceType>();
d_dpdTheta = atomKK->k_dpdTheta.view<DeviceType>();
@ -1993,21 +1818,9 @@ void FixRxKokkos<DeviceType>::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<DeviceType>();
//typename ArrayTypes<DeviceType>::t_ffloat_2d d_cutsq = pairDPDEKK->k_cutsq.template view<DeviceType();
//!< Copies pulled from pairDPDE for local use since pairDPDEKK's objects are protected.
//typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cutsq;
//typename ArrayTypes<DeviceType>::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<DeviceType>::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<DeviceType>().extent(0)) {
memoryKK->destroy_kokkos(k_sumWeights, sumWeights);
memoryKK->create_kokkos (k_sumWeights, sumWeightsCt, "FixRxKokkos::sumWeights");
@ -2036,14 +1848,6 @@ void FixRxKokkos<DeviceType>::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<DeviceType, Tag_FixRxKokkos_zeroTemperatureViews>(0, sumWeightsCt), *this);
// Local list views. (This isn't working!)
@ -2051,9 +1855,6 @@ void FixRxKokkos<DeviceType>::computeLocalTemperature()
if (!list->kokkos)
error->one(FLERR,"list is not a Kokkos list\n");
//typename ArrayTypes<DeviceType>::t_neighbors_2d d_neighbors = k_list->d_neighbors;
//typename ArrayTypes<DeviceType>::t_int_1d d_ilist = k_list->d_ilist;
//typename ArrayTypes<DeviceType>::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<DeviceType>::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<DeviceType>::value, Kokkos::MemoryTraits< AtomicF< NEIGHFLAG >::value> > AtomicViewType;
typedef Kokkos::View< E_FLOAT*, typename DAT::t_efloat_1d::array_layout, typename KKDevice<DeviceType>::value, Kokkos::MemoryTraits< AtomicF< NEIGHFLAG >::value> > AtomicViewType;
AtomicViewType a_dpdThetaLocal = d_dpdThetaLocal;
@ -2174,8 +1974,6 @@ void FixRxKokkos<DeviceType>::computeLocalTemperature()
template <typename DeviceType>
int FixRxKokkos<DeviceType>::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<DeviceType>::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<DeviceType>::pack_forward_comm(int n, int *list, double *buf, in
template <typename DeviceType>
void FixRxKokkos<DeviceType>::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<DeviceType>::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<DeviceType>::unpack_forward_comm(int n, int first, double *buf)
template <typename DeviceType>
int FixRxKokkos<DeviceType>::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<LMPHostType>();
k_sumWeights. template sync<LMPHostType>();
@ -2230,8 +2020,6 @@ int FixRxKokkos<DeviceType>::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<DeviceType>::pack_reverse_comm(int n, int first, double *buf)
template <typename DeviceType>
void FixRxKokkos<DeviceType>::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<DeviceType>::unpack_reverse_comm(int n, int *list, double *buf)
// Signal that the host view has been modified.
k_dpdThetaLocal.template modify<LMPHostType>();
k_sumWeights. template modify<LMPHostType>();
// printf("done with FixRxKokkos::unpack_reverse_comm %d\n", comm->me);
}
namespace LAMMPS_NS {

3
src/MANYBODY/pair_eam.cpp
View File

@ -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;

70
src/fix_adapt.cpp
View File

@ -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) {

2
src/fix_property_atom.cpp
View File

@ -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) :

11
src/fix_property_atom.h
View File

@ -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

4
tools/fep/fdti.py
View File

@ -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()

2
tools/fep/fep.py
View File

@ -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

2
tools/fep/nti.py
View File

@ -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()

5
tools/python/pizza/dump.py
View File

@ -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]) + " "

16
unittest/commands/test_kim_commands.cpp
View File

@ -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)