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

Merge branch 'develop' into pair_style_tutorial

Browse Source
This commit is contained in:
Axel Kohlmeyer
2023-03-23 01:23:36 -04:00
parent 91eb544979 12e0131abf
commit 0007b7c9fe
35 changed files with 582 additions and 791 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) target_compile_definitions(lammps PRIVATE -DLMP_NO_SSL_CHECK)
endif() endif()
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) mark_as_advanced(KIM_EXTRA_UNITTESTS)
find_package(PkgConfig QUIET) find_package(PkgConfig QUIET)
set(DOWNLOAD_KIM_DEFAULT ON) set(DOWNLOAD_KIM_DEFAULT ON)

2
doc/src/Modify_gran_sub_mod.rst
View File

@ -25,7 +25,7 @@ specific implementation. For instance, from the GranSubModNormal class the
GranSubModNormalHooke, GranSubModNormalHertz, and GranSubModNormalJKR classes GranSubModNormalHooke, GranSubModNormalHertz, and GranSubModNormalJKR classes
are derived which calculate Hookean, Hertzian, or JKR normal forces, are derived which calculate Hookean, Hertzian, or JKR normal forces,
respectively. This modular structure simplifies the addition of new granular respectively. This modular structure simplifies the addition of new granular
contact models as as one only needs to create a new GranSubMod class without contact models as one only needs to create a new GranSubMod class without
having to modify the more complex PairGranular, FixGranWall, and GranularModel having to modify the more complex PairGranular, FixGranWall, and GranularModel
classes. Most GranSubMod methods are also already defined by the parent classes classes. Most GranSubMod methods are also already defined by the parent classes
so new contact models typically only require edits to a few relevant methods so new contact models typically only require edits to a few relevant methods

11
doc/src/bond_bpm_rotational.rst
View File

@ -10,7 +10,7 @@ Syntax
bond_style bpm/rotational keyword value attribute1 attribute2 ... bond_style bpm/rotational keyword value attribute1 attribute2 ...
* optional keyword = *overlay/pair* or *store/local* or *smooth* * optional keyword = *overlay/pair* or *store/local* or *smooth* or *break/no*
.. parsed-literal:: .. parsed-literal::
@ -30,6 +30,9 @@ Syntax
*smooth* value = *yes* or *no* *smooth* value = *yes* or *no*
smooths bond forces near the breaking point smooths bond forces near the breaking point
*break/no*
indicates that bonds should not break during a run
Examples Examples
"""""""" """"""""
@ -140,6 +143,12 @@ the *overlay/pair* keyword. These settings require specific
restrictions. Further details can be found in the `:doc: how to restrictions. Further details can be found in the `:doc: how to
<Howto_BPM>` page on BPMs. <Howto_BPM>` page on BPMs.
.. versionadded:: TBD
If the *break/no* keyword is used, then LAMMPS assumes bonds should not break
during a simulation run. This will prevent some unnecessary calculation.
However, if a bond does break, it will trigger an error.
If the *store/local* keyword is used, an internal fix will track bonds that If the *store/local* keyword is used, an internal fix will track bonds that
break during the simulation. Whenever a bond breaks, data is processed break during the simulation. Whenever a bond breaks, data is processed
and transferred to an internal fix labeled *fix_ID*. This allows the and transferred to an internal fix labeled *fix_ID*. This allows the

20
doc/src/bond_bpm_spring.rst
View File

@ -10,7 +10,7 @@ Syntax
bond_style bpm/spring keyword value attribute1 attribute2 ... bond_style bpm/spring keyword value attribute1 attribute2 ...
* optional keyword = *overlay/pair* or *store/local* or *smooth* * optional keyword = *overlay/pair* or *store/local* or *smooth* or *break/no*
.. parsed-literal:: .. parsed-literal::
@ -30,6 +30,9 @@ Syntax
*smooth* value = *yes* or *no* *smooth* value = *yes* or *no*
smooths bond forces near the breaking point smooths bond forces near the breaking point
*break/no*
indicates that bonds should not break during a run
Examples Examples
"""""""" """"""""
@ -47,7 +50,7 @@ Description
.. versionadded:: 4May2022 .. versionadded:: 4May2022
The *bpm/spring* bond style computes forces and torques based on The *bpm/spring* bond style computes forces based on
deviations from the initial reference state of the two atoms. The deviations from the initial reference state of the two atoms. The
reference state is stored by each bond when it is first computed in reference state is stored by each bond when it is first computed in
the setup of a run. Data is then preserved across run commands and is the setup of a run. Data is then preserved across run commands and is
@ -56,7 +59,8 @@ the system will not reset the reference state of a bond.
This bond style only applies central-body forces which conserve the This bond style only applies central-body forces which conserve the
translational and rotational degrees of freedom of a bonded set of translational and rotational degrees of freedom of a bonded set of
particles. The force has a magnitude of particles based on a model described by Clemmer and Robbins
:ref:`(Clemmer) <fragment-Clemmer>`. The force has a magnitude of
.. math:: .. math::
@ -105,6 +109,12 @@ the *overlay/pair* keyword. These settings require specific
restrictions. Further details can be found in the `:doc: how to restrictions. Further details can be found in the `:doc: how to
<Howto_BPM>` page on BPMs. <Howto_BPM>` page on BPMs.
.. versionadded:: TBD
If the *break/no* keyword is used, then LAMMPS assumes bonds should not break
during a simulation run. This will prevent some unnecessary calculation.
However, if a bond does break, it will trigger an error.
If the *store/local* keyword is used, an internal fix will track bonds that If the *store/local* keyword is used, an internal fix will track bonds that
break during the simulation. Whenever a bond breaks, data is processed break during the simulation. Whenever a bond breaks, data is processed
and transferred to an internal fix labeled *fix_ID*. This allows the and transferred to an internal fix labeled *fix_ID*. This allows the
@ -200,6 +210,10 @@ The option defaults are *smooth* = *yes*
---------- ----------
.. _fragment-Clemmer:
**(Clemmer)** Clemmer and Robbins, Phys. Rev. Lett. (2022).
.. _Groot4: .. _Groot4:
**(Groot)** Groot and Warren, J Chem Phys, 107, 4423-35 (1997). **(Groot)** Groot and Warren, J Chem Phys, 107, 4423-35 (1997).

1
doc/src/comm_modify.rst
View File

@ -69,6 +69,7 @@ For many systems this is an efficient algorithm, but for systems with
widely varying cutoffs for different type pairs, the *multi* or *multi/old* mode can widely varying cutoffs for different type pairs, the *multi* or *multi/old* mode can
be faster. In *multi*, each atom is assigned to a collection which should be faster. In *multi*, each atom is assigned to a collection which should
correspond to a set of atoms with similar interaction cutoffs. correspond to a set of atoms with similar interaction cutoffs.
See the :doc:`neighbor <neighbor>` command for a detailed description of collections.
In this case, each atom collection is assigned its own distance In this case, each atom collection is assigned its own distance
cutoff for communication purposes, and fewer atoms will be cutoff for communication purposes, and fewer atoms will be
communicated. in *multi/old*, a similar technique is used but atoms communicated. in *multi/old*, a similar technique is used but atoms

33
doc/src/neighbor.rst
View File

@ -59,9 +59,21 @@ long cutoff, but other type pairs have a much shorter cutoff. The
sized particles, where "size" may mean the physical size of the particle sized particles, where "size" may mean the physical size of the particle
or its cutoff distance for interacting with other particles. Different or its cutoff distance for interacting with other particles. Different
sets of bins are then used to construct the neighbor lists as as further sets of bins are then used to construct the neighbor lists as as further
described by Shire, Hanley, and Stratford :ref:`(Shire) <bytype-Shire>`. described by Shire, Hanley, and Stratford :ref:`(Shire) <multi-Shire>`
This imposes some extra setup overhead, but the searches themselves may and Monti et al. :ref:`(Monti) <multi-Monti>`. This imposes some extra
be much faster. By default, each atom type defines a separate collection setup overhead, but the searches themselves may be much faster.
For instance in a dense binary system in d-dimensions with a ratio of the size
of the largest to smallest collection bin :math:`\lambda`, the computational
costs of building a default neighbor list grows as :math:`\lambda^{2d}` while
the costs for *multi* grows as :math:`\lambda^d`, equivalent to the cost
of force evaluations, as argued in Monti et al. :ref:`(Monti) <multi-Monti>`.
In other words, the neighboring costs of *multi* are expected to scale the
same as force calculations, such that its relative cost is independent of
the particle size ratio. This is not the case for the default style which
becomes substantially more expensive with increasing size ratios.
By default in *multi*, each atom type defines a separate collection
of particles. For systems where two or more atom types have the same of particles. For systems where two or more atom types have the same
size (either physical size or cutoff distance), the definition of size (either physical size or cutoff distance), the definition of
collections can be customized, which can result in less overhead and collections can be customized, which can result in less overhead and
@ -75,8 +87,11 @@ An alternate style, *multi/old*, sets the bin size to 1/2 of the shortest
cutoff distance and multiple sets of bins are defined to search over for cutoff distance and multiple sets of bins are defined to search over for
different atom types. This algorithm used to be the default *multi* different atom types. This algorithm used to be the default *multi*
algorithm in LAMMPS but was found to be significantly slower than the new algorithm in LAMMPS but was found to be significantly slower than the new
approach. For now we are keeping the old option in case there are use cases approach. For the dense binary system, computational costs still grew as
where multi/old outperforms the new multi style. :math:`\lambda^{2d}` at large enough :math:`\lambda`. This is equivalent
to the default style, albeit with a smaller prefactor. For now we are
keeping the old option in case there are use cases where multi/old
outperforms the new multi style.
.. note:: .. note::
@ -118,6 +133,10 @@ Default
---------- ----------
.. _bytype-Shire: .. _multi-Shire:
**(Shire)** Shire, Hanley and Stratford, Comp Part Mech, (2020). **(Shire)** Shire, Hanley and Stratford, Comp. Part. Mech., (2020).
.. _multi-Monti:
**(Monti)** Monti, Clemmer, Srivastava, Silbert, Grest, and Lechman, Phys. Rev. E, (2022).

2
src/.gitignore vendored
View File

@ -56,6 +56,8 @@
/pair_lepton.cpp /pair_lepton.cpp
/pair_lepton.h /pair_lepton.h
/pair_lepton_coul.cpp
/pair_lepton_coul.h
/bond_lepton.cpp /bond_lepton.cpp
/bond_lepton.h /bond_lepton.h
/angle_lepton.cpp /angle_lepton.cpp

16
src/BPM/bond_bpm.cpp
View File

@ -40,6 +40,7 @@ BondBPM::BondBPM(LAMMPS *_lmp) :
{ {
overlay_flag = 0; overlay_flag = 0;
prop_atom_flag = 0; prop_atom_flag = 0;
break_flag = 1;
nvalues = 0; nvalues = 0;
r0_max_estimate = 0.0; r0_max_estimate = 0.0;
@ -104,21 +105,21 @@ void BondBPM::init_style()
} }
} else { } else {
// Require atoms know about all of their bonds and if they break // Require atoms know about all of their bonds and if they break
if (force->newton_bond) if (force->newton_bond && break_flag)
error->all(FLERR, "Without overlay/pair, BPM bond styles require Newton bond off"); error->all(FLERR, "Without overlay/pair or break/no, BPM bond styles require Newton bond off");
// special lj must be 0 1 1 to censor pair forces between bonded particles // special lj must be 0 1 1 to censor pair forces between bonded particles
// special coulomb must be 1 1 1 to ensure all pairs are included in the // special coulomb must be 1 1 1 to ensure all pairs are included in the
// neighbor list and 1-3 and 1-4 special bond lists are skipped // neighbor list and 1-3 and 1-4 special bond lists are skipped
if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0) if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0)
error->all(FLERR, error->all(FLERR,
"Without overlay/pair, BPM bond sytles requires special LJ weights = 0,1,1"); "Without overlay/pair, BPM bond styles requires special LJ weights = 0,1,1");
if (force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 || if (force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 ||
force->special_coul[3] != 1.0) force->special_coul[3] != 1.0)
error->all(FLERR, error->all(FLERR,
"Without overlay/pair, BPM bond sytles requires special Coulomb weights = 1,1,1"); "Without overlay/pair, BPM bond styles requires special Coulomb weights = 1,1,1");
if (id_fix_dummy) { if (id_fix_dummy && break_flag) {
id_fix_update = utils::strdup("BPM_UPDATE_SPECIAL_BONDS"); id_fix_update = utils::strdup("BPM_UPDATE_SPECIAL_BONDS");
fix_update_special_bonds = dynamic_cast<FixUpdateSpecialBonds *>(modify->replace_fix( fix_update_special_bonds = dynamic_cast<FixUpdateSpecialBonds *>(modify->replace_fix(
id_fix_dummy, fmt::format("{} all UPDATE_SPECIAL_BONDS", id_fix_update), 1)); id_fix_dummy, fmt::format("{} all UPDATE_SPECIAL_BONDS", id_fix_update), 1));
@ -188,6 +189,9 @@ void BondBPM::settings(int narg, char **arg)
} else if (strcmp(arg[iarg], "overlay/pair") == 0) { } else if (strcmp(arg[iarg], "overlay/pair") == 0) {
overlay_flag = 1; overlay_flag = 1;
iarg++; iarg++;
} else if (strcmp(arg[iarg], "break/no") == 0) {
break_flag = 0;
iarg++;
} else { } else {
leftover_iarg.push_back(iarg); leftover_iarg.push_back(iarg);
iarg++; iarg++;
@ -328,6 +332,8 @@ void BondBPM::read_restart(FILE *fp)
void BondBPM::process_broken(int i, int j) void BondBPM::process_broken(int i, int j)
{ {
if (break_flag)
error->one(FLERR, "BPM bond broke with break/no option");
if (fix_store_local) { if (fix_store_local) {
for (int n = 0; n < nvalues; n++) (this->*pack_choice[n])(n, i, j); for (int n = 0; n < nvalues; n++) (this->*pack_choice[n])(n, i, j);

2
src/BPM/bond_bpm.h
View File

@ -50,7 +50,7 @@ class BondBPM : public Bond {
FnPtrPack *pack_choice; // ptrs to pack functions FnPtrPack *pack_choice; // ptrs to pack functions
double *output_data; double *output_data;
int prop_atom_flag, nvalues, overlay_flag; int prop_atom_flag, nvalues, overlay_flag, break_flag;
int index_x_ref, index_y_ref, index_z_ref; int index_x_ref, index_y_ref, index_z_ref;
void pack_id1(int, int, int); void pack_id1(int, int, int);

43
src/KIM/fix_store_kim.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -69,11 +68,11 @@ using namespace FixConst;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixStoreKIM::FixStoreKIM(LAMMPS *lmp, int narg, char **arg) FixStoreKIM::FixStoreKIM(LAMMPS *lmp, int narg, char **arg) :
: Fix(lmp, narg, arg), simulator_model(nullptr), model_name(nullptr), Fix(lmp, narg, arg), simulator_model(nullptr), model_name(nullptr), model_units(nullptr),
model_units(nullptr), user_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 // free associated storage
if (simulator_model) { if (simulator_model) {
auto sm = (KIM_SimulatorModel *)simulator_model; auto sm = (KIM_SimulatorModel *) simulator_model;
KIM_SimulatorModel_Destroy(&sm); KIM_SimulatorModel_Destroy(&sm);
simulator_model = nullptr; simulator_model = nullptr;
} }
if (model_name) { if (model_name) {
auto mn = (char *)model_name; auto mn = (char *) model_name;
delete[] mn; delete[] mn;
model_name = nullptr; model_name = nullptr;
} }
if (model_units) { if (model_units) {
auto mu = (char *)model_units; auto mu = (char *) model_units;
delete[] mu; delete[] mu;
model_units = nullptr; model_units = nullptr;
} }
if (user_units) { if (user_units) {
auto uu = (char *)user_units; auto uu = (char *) user_units;
delete[] uu; delete[] uu;
user_units = nullptr; user_units = nullptr;
} }
@ -121,38 +120,44 @@ void FixStoreKIM::setptr(const std::string &name, void *ptr)
{ {
if (name == "simulator_model") { if (name == "simulator_model") {
if (simulator_model) { if (simulator_model) {
auto sm = (KIM_SimulatorModel *)simulator_model; auto sm = (KIM_SimulatorModel *) simulator_model;
KIM_SimulatorModel_Destroy(&sm); KIM_SimulatorModel_Destroy(&sm);
} }
simulator_model = ptr; simulator_model = ptr;
} else if (name == "model_name") { } else if (name == "model_name") {
if (model_name) { if (model_name) {
auto mn = (char *)model_name; auto mn = (char *) model_name;
delete[] mn; delete[] mn;
} }
model_name = ptr; model_name = ptr;
} else if (name == "model_units") { } else if (name == "model_units") {
if (model_units) { if (model_units) {
auto mu = (char *)model_units; auto mu = (char *) model_units;
delete[] mu; delete[] mu;
} }
model_units = ptr; model_units = ptr;
} else if (name == "user_units") { } else if (name == "user_units") {
if (user_units) { if (user_units) {
auto uu = (char *)user_units; auto uu = (char *) user_units;
delete[] uu; delete[] uu;
} }
user_units = ptr; 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) void *FixStoreKIM::getptr(const std::string &name)
{ {
if (name == "simulator_model") return simulator_model; if (name == "simulator_model")
else if (name == "model_name") return model_name; return simulator_model;
else if (name == "model_units") return model_units; else if (name == "model_name")
else if (name == "user_units") return user_units; return model_name;
else return nullptr; 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 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -70,40 +69,40 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
static constexpr const char *const cite_openkim = static constexpr const char *const cite_openkim =
"OpenKIM Project: doi:10.1007/s11837-011-0102-6\n\n" "OpenKIM Project: doi:10.1007/s11837-011-0102-6\n\n"
"@Article{tadmor:elliott:2011,\n" "@Article{tadmor:elliott:2011,\n"
" author = {E. B. Tadmor and R. S. Elliott and J. P. Sethna and R. E. Miller " " author = {E. B. Tadmor and R. S. Elliott and J. P. Sethna and R. E. Miller "
"and C. A. Becker},\n" "and C. A. Becker},\n"
" title = {The potential of atomistic simulations and the {K}nowledgebase of " " title = {The potential of atomistic simulations and the {K}nowledgebase of "
"{I}nteratomic {M}odels},\n" "{I}nteratomic {M}odels},\n"
" journal = {{JOM}},\n" " journal = {{JOM}},\n"
" year = 2011,\n" " year = 2011,\n"
" volume = 63,\n" " volume = 63,\n"
" number = 17,\n" " number = 17,\n"
" pages = {17},\n" " pages = {17},\n"
" doi = {10.1007/s11837-011-0102-6}\n" " doi = {10.1007/s11837-011-0102-6}\n"
"}\n\n"; "}\n\n";
static constexpr const char *const cite_openkim_query = static constexpr const char *const cite_openkim_query =
"OpenKIM query: doi:10.1063/5.0014267\n\n" "OpenKIM query: doi:10.1063/5.0014267\n\n"
"@Article{karls:bierbaum:2020,\n" "@Article{karls:bierbaum:2020,\n"
" author = {D. S. Karls and M. Bierbaum and A. A. Alemi and R. S. Elliott " " author = {D. S. Karls and M. Bierbaum and A. A. Alemi and R. S. Elliott "
"and J. P. Sethna and E. B. Tadmor},\n" "and J. P. Sethna and E. B. Tadmor},\n"
" title = {The {O}pen{KIM} processing pipeline: {A} cloud-based automatic " " title = {The {O}pen{KIM} processing pipeline: {A} cloud-based automatic "
"material property computation engine},\n" "material property computation engine},\n"
" journal = {{T}he {J}ournal of {C}hemical {P}hysics},\n" " journal = {{T}he {J}ournal of {C}hemical {P}hysics},\n"
" year = 2020,\n" " year = 2020,\n"
" volume = 153,\n" " volume = 153,\n"
" number = 6,\n" " number = 6,\n"
" pages = {064104},\n" " pages = {064104},\n"
" doi = {10.1063/5.0014267}\n" " doi = {10.1063/5.0014267}\n"
"}\n\n"; "}\n\n";
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void KimCommand::command(int narg, char **arg) 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]); const std::string subcmd(arg[0]);
narg--; narg--;
@ -132,5 +131,6 @@ void KimCommand::command(int narg, char **arg)
auto cmd = new KimQuery(lmp); auto cmd = new KimQuery(lmp);
cmd->command(narg, arg); cmd->command(narg, arg);
delete cmd; 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 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, error->all(FLERR,
"Illegal 'kim init' command.\n" "Illegal 'kim init' command.\n"
"The argument followed by unit_style {} is an optional argument and when " "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); user_units);
} }
} else } else
@ -109,58 +111,10 @@ void KimInit::command(int narg, char **arg)
if (universe->nprocs > 1) MPI_Barrier(universe->uworld); if (universe->nprocs > 1) MPI_Barrier(universe->uworld);
determine_model_type_and_units(model_name, user_units, &model_units, pkim); determine_model_type_and_units(model_name, user_units, &model_units, pkim);
write_log_cite(lmp, model_type, model_name); write_log_cite(lmp, model_type, model_name);
do_init(model_name, user_units, model_units, pkim); 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 void KimInit::print_dirs(struct KIM_Collections *const collections) const
{ {
int kim_error = 0; 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. // create storage proxy fix. delete existing fix, if needed.
int ifix = modify->find_fix("KIM_MODEL_STORE"); if (modify->get_fix_by_id("KIM_MODEL_STORE")) modify->delete_fix("KIM_MODEL_STORE");
if (ifix >= 0) modify->delete_fix(ifix); auto fix_store = dynamic_cast<FixStoreKIM *>(modify->add_fix("KIM_MODEL_STORE all STORE/KIM"));
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]);
fix_store->setptr("model_name", (void *) model_name); fix_store->setptr("model_name", (void *) model_name);
fix_store->setptr("user_units", (void *) user_units); fix_store->setptr("user_units", (void *) user_units);
fix_store->setptr("model_units", (void *) model_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); ier = lammps_unit_conversion(units[i], from, to, conversion_factor);
if (ier != 0) if (ier != 0)
error->all(FLERR, error->all(FLERR, "Unable to obtain conversion factor: unit = {}; from = {}; to = {}",
"Unable to obtain conversion factor: "
"unit = {}; from = {}; to = {}",
units[i], from, to); units[i], from, to);
variable->internal_set(v_unit, conversion_factor); 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) 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) if (!domain->box_exist)
error->all(FLERR, "Must use 'kim interactions' command after " error->all(FLERR, "Use of 'kim interactions' before simulation box is defined");
"simulation box is defined");
do_setup(narg, arg); do_setup(narg, arg);
} }
@ -103,11 +102,9 @@ void KimInteractions::do_setup(int narg, char **arg)
if ((narg == 1) && (arg_str == "fixed_types")) { if ((narg == 1) && (arg_str == "fixed_types")) {
fixed_types = true; fixed_types = true;
} else if (narg != atom->ntypes) { } else if (narg != atom->ntypes) {
error->all(FLERR, "Illegal 'kim interactions' command.\nThe " error->all(FLERR, "Illegal 'kim interactions' command.\nThe LAMMPS simulation has {} atom "
"LAMMPS simulation has {} atom type(s), but " "type(s), but {} chemical species passed to the 'kim interactions' command",
"{} chemical species passed to the " atom->ntypes, narg);
"'kim interactions' command",
atom->ntypes, narg);
} else { } else {
fixed_types = false; 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. // retrieve model name and pointer to simulator model class instance.
// validate model name if not given as null pointer. // validate model name if not given as null pointer.
int ifix = modify->find_fix("KIM_MODEL_STORE"); auto fix_store = dynamic_cast<FixStoreKIM *>(modify->get_fix_by_id("KIM_MODEL_STORE"));
if (ifix >= 0) { if (fix_store) {
auto fix_store = dynamic_cast<FixStoreKIM *>(modify->fix[ifix]);
model_name = (char *)fix_store->getptr("model_name"); model_name = (char *)fix_store->getptr("model_name");
simulatorModel = (KIM_SimulatorModel *)fix_store->getptr("simulator_model"); simulatorModel = (KIM_SimulatorModel *)fix_store->getptr("simulator_model");
} else error->all(FLERR, "Must use 'kim init' before 'kim interactions'"); } else error->all(FLERR, "Must use 'kim init' before 'kim interactions'");
// Begin output to log file // Begin output to log file
input->write_echo("#=== BEGIN kim interactions ===========================" input->write_echo("#=== BEGIN kim interactions ==================================\n");
"=======\n");
if (simulatorModel) { if (simulatorModel) {
auto first_visit = input->variable->find("kim_update"); 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 (atom_type_sym == sim_species) species_is_supported = true;
} }
if (!species_is_supported) { if (!species_is_supported) {
error->all(FLERR, "Species '{}' is not supported by this " error->all(FLERR, "Species '{}' is not supported by this KIM Simulator Model",
"KIM Simulator Model", atom_type_sym); atom_type_sym);
} }
} }
} else { } else {
@ -179,18 +174,15 @@ void KimInteractions::do_setup(int narg, char **arg)
const char *sim_field, *sim_value; const char *sim_field, *sim_value;
KIM_SimulatorModel_GetNumberOfSimulatorFields(simulatorModel, &sim_fields); KIM_SimulatorModel_GetNumberOfSimulatorFields(simulatorModel, &sim_fields);
for (int i = 0; i < sim_fields; ++i) { for (int i = 0; i < sim_fields; ++i) {
KIM_SimulatorModel_GetSimulatorFieldMetadata( KIM_SimulatorModel_GetSimulatorFieldMetadata(simulatorModel, i, &sim_lines, &sim_field);
simulatorModel, i, &sim_lines, &sim_field);
const std::string sim_field_str(sim_field); const std::string sim_field_str(sim_field);
if (sim_field_str == "units") { if (sim_field_str == "units") {
KIM_SimulatorModel_GetSimulatorFieldLine( KIM_SimulatorModel_GetSimulatorFieldLine(simulatorModel, i, 0, &sim_value);
simulatorModel, i, 0, &sim_value);
const std::string sim_value_str(sim_value); if (strcmp(sim_value, update->unit_style) != 0)
const std::string unit_style_str(update->unit_style); error->all(FLERR, "Incompatible units for KIM Simulator Model: {} vs {}",
if (sim_value_str != unit_style_str) sim_value, update->unit_style);
error->all(FLERR, "Incompatible units for KIM Simulator Model");
} }
} }
@ -217,8 +209,7 @@ void KimInteractions::do_setup(int narg, char **arg)
no_model_definition = false; no_model_definition = false;
for (int j = 0; j < sim_lines; ++j) { for (int j = 0; j < sim_lines; ++j) {
KIM_SimulatorModel_GetSimulatorFieldLine( KIM_SimulatorModel_GetSimulatorFieldLine(simulatorModel, i, j, &sim_value);
simulatorModel, i, j, &sim_value);
if (utils::strmatch(sim_value, "^KIM_SET_TYPE_PARAMETERS")) { if (utils::strmatch(sim_value, "^KIM_SET_TYPE_PARAMETERS")) {
// Notes regarding the KIM_SET_TYPE_PARAMETERS command // Notes regarding the KIM_SET_TYPE_PARAMETERS command
// * This is an INTERNAL command. // * This is an INTERNAL command.
@ -264,8 +255,7 @@ void KimInteractions::do_setup(int narg, char **arg)
} }
// End output to log file // End output to log file
input->write_echo("#=== END kim interactions =============================" input->write_echo("#=== END kim interactions ====================================\n\n");
"=======\n\n");
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -276,8 +266,7 @@ void KimInteractions::KIM_SET_TYPE_PARAMETERS(const std::string &input_line) con
const std::string key = words[1]; const std::string key = words[1];
if (key != "pair" && key != "charge") if (key != "pair" && key != "charge")
error->one(FLERR, "Unrecognized KEY {} for " error->one(FLERR, "Unrecognized KEY {} for KIM_SET_TYPE_PARAMETERS command", key);
"KIM_SET_TYPE_PARAMETERS command", key);
std::string filename = words[2]; std::string filename = words[2];
std::vector<std::string> species(words.begin() + 3, words.end()); 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; FILE *fp = nullptr;
if (comm->me == 0) { if (comm->me == 0) {
fp = fopen(filename.c_str(), "r"); 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; char line[MAXLINE], *ptr;

76
src/KIM/kim_param.cpp
View File

@ -63,6 +63,7 @@
#include "fix_store_kim.h" #include "fix_store_kim.h"
#include "force.h" #include "force.h"
#include "input.h" #include "input.h"
#include "kim_units.h"
#include "modify.h" #include "modify.h"
#include "pair_kim.h" #include "pair_kim.h"
#include "variable.h" #include "variable.h"
@ -77,61 +78,7 @@ extern "C"
} }
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using kim_units::get_kim_unit_names;
/* ---------------------------------------------------------------------- */
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
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -159,20 +106,17 @@ void KimParam::command(int narg, char **arg)
error->all(FLERR, "Incorrect arguments in 'kim param' command.\n" error->all(FLERR, "Incorrect arguments in 'kim param' command.\n"
"'kim param get/set' is mandatory"); "'kim param get/set' is mandatory");
int const ifix = modify->find_fix("KIM_MODEL_STORE"); auto fix_store = dynamic_cast<FixStoreKIM *>(modify->get_fix_by_id("KIM_MODEL_STORE"));
if (ifix >= 0) { if (fix_store) {
auto fix_store = reinterpret_cast<FixStoreKIM *>(modify->fix[ifix]); auto *simulatorModel = reinterpret_cast<KIM_SimulatorModel *>(
fix_store->getptr("simulator_model"));
KIM_SimulatorModel *simulatorModel =
reinterpret_cast<KIM_SimulatorModel *>(
fix_store->getptr("simulator_model"));
if (simulatorModel) if (simulatorModel)
error->all(FLERR, "'kim param' can only be used with a KIM Portable Model"); error->all(FLERR, "'kim param' can only be used with a KIM Portable Model");
} }
input->write_echo(fmt::format("#=== BEGIN kim param {} ===================" input->write_echo(fmt::format("#=== BEGIN kim param {} =====================================\n",
"==================\n", kim_param_get_set)); kim_param_get_set));
KIM_Model *pkim = nullptr; KIM_Model *pkim = nullptr;
@ -431,6 +375,6 @@ void KimParam::command(int narg, char **arg)
} else } else
error->all(FLERR, "This model has No mutable parameters"); error->all(FLERR, "This model has No mutable parameters");
input->write_echo(fmt::format("#=== END kim param {} =====================" input->write_echo(fmt::format("#=== END kim param {} =======================================\n",
"==================\n", kim_param_get_set)); kim_param_get_set));
} }

2
src/KIM/kim_param.h
View File

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

68
src/KIM/kim_property.cpp
View File

@ -85,14 +85,13 @@ void KimProperty::command(int narg, char **arg)
{ {
#if LMP_PYTHON #if LMP_PYTHON
#if PY_MAJOR_VERSION >= 3 #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]); const std::string subcmd(arg[0]);
if ((subcmd != "create") && (subcmd != "destroy") && (subcmd != "modify") && if ((subcmd != "create") && (subcmd != "destroy") && (subcmd != "modify") &&
(subcmd != "remove") && (subcmd != "dump")) { (subcmd != "remove") && (subcmd != "dump")) {
std::string msg("Incorrect arguments in 'kim property' command.\n"); error->all(FLERR, "Incorrect first argument {} to 'kim property' command.\n"
msg += "'kim property create/destroy/modify/remove/dump' is mandatory"; "One of create, destroy, modify, remove, or dump is mandatory", subcmd);
error->all(FLERR, msg);
} }
input->write_echo("#=== kim property =====================================" input->write_echo("#=== kim property ====================================="
@ -117,12 +116,11 @@ void KimProperty::command(int narg, char **arg)
kim_property = PyImport_Import(obj); kim_property = PyImport_Import(obj);
if (!kim_property) { if (!kim_property) {
PyGILState_Release(gstate); PyGILState_Release(gstate);
std::string msg("Unable to import Python kim_property module!"); error->all(FLERR, "Unable to import Python kim_property module!\n"
msg += "\nkim-property Python package can be installed with pip:\n"; "kim-property Python package can be installed with pip:\n"
msg += "'pip install kim-property'\nSee the installation instructions "; "'pip install kim-property'\nSee the installation instructions at\n"
msg += "at\nhttps://github.com/openkim/kim-property#installing-kim-"; "https://github.com/openkim/kim-property#installing-kim-property\n"
msg += "property\nfor detailed information"; "for detailed information");
error->all(FLERR, msg);
} }
// Decrementing of the reference count // Decrementing of the reference count
@ -147,9 +145,8 @@ void KimProperty::command(int narg, char **arg)
PyObject_GetAttrString(kim_property, "kim_property_create"); PyObject_GetAttrString(kim_property, "kim_property_create");
if (!pFunc) { if (!pFunc) {
PyGILState_Release(gstate); PyGILState_Release(gstate);
std::string msg("Unable to get an attribute named "); error->all(FLERR, "Unable to get an attribute named 'kim_property_create' from "
msg += "'kim_property_create' from a kim_property object"; "a kim_property object");
error->all(FLERR, msg);
} }
// Decrementing of the reference count // Decrementing of the reference count
@ -182,16 +179,13 @@ void KimProperty::command(int narg, char **arg)
if (!pValue) { if (!pValue) {
PyErr_Print(); PyErr_Print();
PyGILState_Release(gstate); PyGILState_Release(gstate);
std::string msg("Python 'kim_property_create' function "); error->all(FLERR, "Python 'kim_property_create' function evaluation failed");
msg += "evaluation failed";
error->one(FLERR, msg);
} }
// Python function returned a string value // Python function returned a string value
const char *pystr = PyUnicode_AsUTF8(pValue); const char *pystr = PyUnicode_AsUTF8(pValue);
if (kim_str) input->variable->set_string("kim_property_str", pystr); if (kim_str) input->variable->set_string("kim_property_str", pystr);
else else input->variable->set(fmt::format("kim_property_str string '{}'", pystr));
input->variable->set(fmt::format("kim_property_str string '{}'", pystr));
Py_XDECREF(pArgs); Py_XDECREF(pArgs);
Py_XDECREF(pFunc); Py_XDECREF(pFunc);
@ -216,9 +210,8 @@ void KimProperty::command(int narg, char **arg)
PyObject_GetAttrString(kim_property, "kim_property_destroy"); PyObject_GetAttrString(kim_property, "kim_property_destroy");
if (!pFunc) { if (!pFunc) {
PyGILState_Release(gstate); PyGILState_Release(gstate);
std::string msg("Unable to get an attribute named "); error->all(FLERR, "Unable to get an attribute named 'kim_property_destroy' "
msg += "'kim_property_destroy' from a kim_property object"; "from a kim_property object");
error->all(FLERR, msg);
} }
// Decrementing of the reference count // Decrementing of the reference count
@ -244,9 +237,7 @@ void KimProperty::command(int narg, char **arg)
if (!pValue) { if (!pValue) {
PyErr_Print(); PyErr_Print();
PyGILState_Release(gstate); PyGILState_Release(gstate);
std::string msg("Python 'kim_property_destroy' function "); error->all(FLERR, "Python 'kim_property_destroy' function evaluation failed");
msg += "evaluation failed";
error->one(FLERR, msg);
} }
// Python function returned a string value // Python function returned a string value
@ -276,9 +267,8 @@ void KimProperty::command(int narg, char **arg)
PyObject_GetAttrString(kim_property, "kim_property_modify"); PyObject_GetAttrString(kim_property, "kim_property_modify");
if (!pFunc) { if (!pFunc) {
PyGILState_Release(gstate); PyGILState_Release(gstate);
std::string msg("Unable to get an attribute named "); error->all(FLERR, "Unable to get an attribute named 'kim_property_modify' "
msg += "'kim_property_modify' from a kim_property object"; "from a kim_property object");
error->all(FLERR, msg);
} }
// Decrementing of the reference count // Decrementing of the reference count
@ -309,9 +299,7 @@ void KimProperty::command(int narg, char **arg)
if (!pValue) { if (!pValue) {
PyErr_Print(); PyErr_Print();
PyGILState_Release(gstate); PyGILState_Release(gstate);
std::string msg("Python 'kim_property_modify' function "); error->all(FLERR, "Python 'kim_property_modify' function evaluation failed");
msg += "evaluation failed";
error->one(FLERR, msg);
} }
// Python function returned a string value // Python function returned a string value
@ -341,9 +329,8 @@ void KimProperty::command(int narg, char **arg)
PyObject_GetAttrString(kim_property, "kim_property_remove"); PyObject_GetAttrString(kim_property, "kim_property_remove");
if (!pFunc) { if (!pFunc) {
PyGILState_Release(gstate); PyGILState_Release(gstate);
std::string msg("Unable to get an attribute named "); error->all(FLERR, "Unable to get an attribute named "
msg += "'kim_property_remove' from a kim_property object"; "'kim_property_remove' from a kim_property object");
error->all(FLERR, msg);
} }
// Decrementing of the reference count // Decrementing of the reference count
@ -374,9 +361,7 @@ void KimProperty::command(int narg, char **arg)
if (!pValue) { if (!pValue) {
PyErr_Print(); PyErr_Print();
PyGILState_Release(gstate); PyGILState_Release(gstate);
std::string msg("Python 'kim_property_remove' function "); error->all(FLERR, "Python 'kim_property_remove' function evaluation failed");
msg += "evaluation failed";
error->one(FLERR, msg);
} }
// Python function returned a string value // Python function returned a string value
@ -404,9 +389,8 @@ void KimProperty::command(int narg, char **arg)
PyObject_GetAttrString(kim_property, "kim_property_dump"); PyObject_GetAttrString(kim_property, "kim_property_dump");
if (!pFunc) { if (!pFunc) {
PyGILState_Release(gstate); PyGILState_Release(gstate);
std::string msg("Unable to get an attribute named "); error->all(FLERR, "Unable to get an attribute named "
msg += "'kim_property_dump' from a kim_property object"; "'kim_property_dump' from a kim_property object");
error->all(FLERR, msg);
} }
// Decrementing of the reference count // Decrementing of the reference count
@ -428,14 +412,12 @@ void KimProperty::command(int narg, char **arg)
if (comm->me == 0) { if (comm->me == 0) {
// call the Python kim_property_dump function // 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); pValue = PyObject_CallObject(pFunc, pArgs);
if (!pValue) { if (!pValue) {
PyErr_Print(); PyErr_Print();
PyGILState_Release(gstate); PyGILState_Release(gstate);
std::string msg("Python 'kim_property_dump' function "); error->one(FLERR, "Python 'kim_property_dump' function evaluation failed");
msg += "evaluation failed";
error->one(FLERR, msg);
} }
} else } else
pValue = nullptr; 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) // check the query_args format (a series of keyword=value pairs)
for (int i = 2; i < narg; ++i) { for (int i = 2; i < narg; ++i) {
if (!utils::strmatch(arg[i], "[=][\\[].*[\\]]")) if (!utils::strmatch(arg[i], "[=][\\[].*[\\]]"))
error->all(FLERR, "Illegal query format.\nInput argument " error->all(FLERR, "Illegal query format.\nInput argument of `{}` to 'kim query' is wrong. "
"of `{}` to 'kim query' is wrong. The " "The query format is the keyword=[value], where value is always an array of one "
"query format is the keyword=[value], " "or more comma-separated items", arg[i]);
"where value is always an array of one or "
"more comma-separated items", arg[i]);
} }
if (query_function != "get_available_models") { 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 the model name is not provided by the user
if (model_name.empty()) { if (model_name.empty()) {
// check if we had a kim init command by finding fix STORE/KIM // check if we had a kim init command by finding fix STORE/KIM
const int ifix = modify->find_fix("KIM_MODEL_STORE"); auto fix_store = dynamic_cast<FixStoreKIM *>(modify->get_fix_by_id("KIM_MODEL_STORE"));
if (ifix >= 0) { if (fix_store) {
auto fix_store = dynamic_cast<FixStoreKIM *>(modify->fix[ifix]);
char *model_name_c = (char *) fix_store->getptr("model_name"); char *model_name_c = (char *) fix_store->getptr("model_name");
model_name = model_name_c; model_name = model_name_c;
} else { } else {
error->all(FLERR, "Illegal query format.\nMust use 'kim init' " error->all(FLERR, "Illegal query format.\nMust use 'kim init' before 'kim query' "
"before 'kim query' or must provide the model name " "or must provide the model name after query function with the format of "
"after query function with the format of "
"'model=[model_name]'"); "'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 // and then the error message that was returned by the web server
if (strlen(value) == 0) { if (strlen(value) == 0) {
auto msg = fmt::format("OpenKIM query failed: {}", value + 1); error->all(FLERR, "OpenKIM query failed: {}", value + 1);
delete [] value; delete[] value;
error->all(FLERR, msg);
} else if (strcmp(value, "EMPTY") == 0) { } else if (strcmp(value, "EMPTY") == 0) {
delete [] value; delete[] value;
error->all(FLERR, "OpenKIM query returned no results"); error->all(FLERR, "OpenKIM query returned no results");
} }
input->write_echo("#=== BEGIN kim-query ==================================" input->write_echo("#=== BEGIN kim-query =========================================\n");
"=======\n");
// trim list of models to those that are installed on the system // trim list of models to those that are installed on the system
if (query_function == "get_available_models") { if (query_function == "get_available_models") {
Tokenizer vals(value, ", \""); Tokenizer vals(value, ", \"");
@ -199,7 +193,7 @@ void KimQuery::command(int narg, char **arg)
KIM_CollectionItemType typ; KIM_CollectionItemType typ;
if (KIM_Collections_Create(&collections)) { if (KIM_Collections_Create(&collections)) {
delete [] value; delete[] value;
error->all(FLERR, "Unable to access KIM Collections to find Model"); 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)); fmt::format("# Missing OpenKIM models: {}\n\n", missing));
if (available.empty()) { if (available.empty()) {
delete [] value; delete[] value;
error->all(FLERR,"There are no matching OpenKIM models installed on the system"); 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->variable->set(setcmd);
input->write_echo(fmt::format("variable {}\n", setcmd)); input->write_echo(fmt::format("variable {}\n", setcmd));
} }
input->write_echo("#=== END kim-query ====================================" input->write_echo("#=== END kim-query ===========================================\n\n");
"=======\n\n");
delete [] value; delete[] value;
#else #else
error->all(FLERR, "Cannot use 'kim query' command when KIM package " error->all(FLERR, "Cannot use 'kim query' command when KIM package is compiled without "
"is compiled without support for libcurl"); "support for libcurl");
#endif #endif
} }

425
src/KIM/kim_units.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -54,179 +53,182 @@
Designed for use with the kim-api-2.0.2 (and newer) package Designed for use with the kim-api-2.0.2 (and newer) package
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#include "kim_units.h"
#include "error.h"
#include <cmath> #include <cmath>
#include <map> #include <map>
#include <string> #include <string>
#include <utility> #include <utility>
using namespace std; using namespace std;
using namespace LAMMPS_NS;
namespace using namespace kim_units;
{
// Constants of nature and basic conversion factors // Constants of nature and basic conversion factors
// Source: https://physics.nist.gov/cuu/Constants/Table/allascii.txt // Source: https://physics.nist.gov/cuu/Constants/Table/allascii.txt
// Working with NIST values even when there are newer values for // Working with NIST values even when there are newer values for
// compatibility with LAMMPS // compatibility with LAMMPS
// clang-format off
/*---------------------- /*----------------------
Fundamental constants Fundamental constants
------------------------ */ ------------------------ */
double const boltz_si = 1.38064852e-23; // [J K^-1] Boltzmann's factor static double constexpr boltz_si = 1.38064852e-23; // [J K^-1] Boltzmann's factor
// (NIST value) // (NIST value)
double const Nav = 6.022140857e23; // [unitless] Avogadro's number static double constexpr Nav = 6.022140857e23; // [unitless] Avogadro's number (NIST value)
// (NIST value) // static double constexp Nav = 6.02214076e23; // [unitless] Avogadro's number
// double const Nav = 6.02214076e23; // [unitless] Avogadro's number // (official value May 2019)
// (official value May 2019) static double constexpr me_si = 9.10938356e-31; // [kg] electron rest mass (NIST value)
double const me_si = 9.10938356e-31; // [kg] electron rest mass // static double constexpr me_si = 9.10938291e-31; // [kg] electron rest mass
// (NIST value) static double constexpr e_si = 1.6021766208e-19; // [C] elementary charge
// double me_si = 9.10938291e-31; // [kg] electron rest mass // (charge of an electron/proton)
double const e_si = 1.6021766208e-19; // [C] elementary charge // (NIST value)
// (charge of an electron/proton)
// (NIST value)
/*---------------------- /*----------------------
Distance units Distance units
------------------------ */ ------------------------ */
double const bohr_si = 5.2917721067e-11; // [m] Bohr unit (distance between static double constexpr bohr_si = 5.2917721067e-11; // [m] Bohr unit (distance between
// nucleus and electron in H) // nucleus and electron in H) (NIST value)
// (NIST value) static double constexpr angstrom_si = 1e-10; // [m] Angstrom
double const angstrom_si = 1e-10; // [m] Angstrom static double constexpr centimeter_si = 1e-2; // [m] centimeter
double const centimeter_si = 1e-2; // [m] centimeter static double constexpr micrometer_si = 1e-6; // [m] micrometer (micron)
double const micrometer_si = 1e-6; // [m] micrometer (micron) static double constexpr nanometer_si = 1e-9; // [m] nanometer
double const nanometer_si = 1e-9; // [m] nanometer
/*---------------------- /*----------------------
Mass units Mass units
------------------------ */ ------------------------ */
double const gram_per_mole_si = 1e-3/Nav; // [kg] gram per mole static double constexpr gram_per_mole_si = 1e-3/Nav; // [kg] gram per mole
double const amu_si = 1e-3/Nav; // [kg] atomic mass unit (molecular static double constexpr amu_si = 1e-3/Nav; // [kg] atomic mass unit (molecular
// weight) For example, the mean // weight) For example, the mean
// molecular weight of water // molecular weight of water
// is 18.015 atomic mass units // is 18.015 atomic mass units
// (amu), so one mole of water // (amu), so one mole of water
// weight 18.015 grams. // weight 18.015 grams.
double const gram_si = 1e-3; // [kg] gram static double constexpr gram_si = 1e-3; // [kg] gram
double const picogram_si = 1e-15; // [kg] picogram static double constexpr picogram_si = 1e-15; // [kg] picogram
double const attogram_si = 1e-21; // [kg[ attogram static double constexpr attogram_si = 1e-21; // [kg[ attogram
/*---------------------- /*----------------------
Time units Time units
------------------------ */ ------------------------ */
double const atu_si = 2.418884326509e-17; // [s] atomic time unit static double constexpr atu_si = 2.418884326509e-17; // [s] atomic time unit
// ( = hbar/E_h where E_h is the // ( = hbar/E_h where E_h is the
// Hartree energy) (NIST value) // Hartree energy) (NIST value)
double const atu_electron_si = atu_si*sqrt(amu_si/me_si); static double constexpr atu_electron_si = 5153034.567408186; // must not use sqrt() in constexpr
// [s] atomic time unit // static double constexpr atu_electron_si = atu_si*sqrt(amu_si/me_si);
// used in electron system (see https://sourceforge.net/p/lammps/mailman/lammps-users/thread/BCA2BDB2-BA03-4280-896F-1E6120EF47B2%40caltech.edu/) // [s] atomic time unit
double const microsecond_si = 1e-6; // [s] microsecond // used in electron system (see https://sourceforge.net/p/lammps/mailman/lammps-users/thread/BCA2BDB2-BA03-4280-896F-1E6120EF47B2%40caltech.edu/)
double const nanosecond_si = 1e-9; // [s] nanosecond static double constexpr microsecond_si = 1e-6; // [s] microsecond
double const picosecond_si = 1e-12; // [s] picosecond static double constexpr nanosecond_si = 1e-9; // [s] nanosecond
double const femtosecond_si = 1e-15; // [s] femtosecond static double constexpr picosecond_si = 1e-12; // [s] picosecond
static double constexpr femtosecond_si = 1e-15; // [s] femtosecond
/*---------------------- /*----------------------
Density units Density units
------------------------ */ ------------------------ */
double const gram_per_centimetercu_si = static double constexpr gram_per_centimetercu_si =
gram_si/pow(centimeter_si,3); // [kg/m^3] gram/centimeter^3 gram_si/centimeter_si/centimeter_si/centimeter_si; // [kg/m^3] gram/centimeter^3
double const amu_per_bohrcu_si = amu_si/pow(bohr_si,3); // [kg/m^3] amu/bohr^3 static double constexpr amu_per_bohrcu_si = amu_si/bohr_si/bohr_si/bohr_si; // [kg/m^3] amu/bohr^3
double const picogram_per_micrometercu_si = static double constexpr picogram_per_micrometercu_si =
picogram_si/pow(micrometer_si,3); // [kg/m^3] picogram/micrometer^3 picogram_si/micrometer_si/micrometer_si/micrometer_si; // [kg/m^3] picogram/micrometer^3
double const attogram_per_nanometercu_si = static double constexpr attogram_per_nanometercu_si =
attogram_si/pow(nanometer_si,3); // [kg/m^3] attogram/nanometer^3 attogram_si/nanometer_si/nanometer_si/nanometer_si; // [kg/m^3] attogram/nanometer^3
/*---------------------- /*----------------------
Energy/torque units Energy/torque units
------------------------ */ ------------------------ */
double const kcal_si = 4184.0; // [J] kilocalorie (heat energy static double constexpr kcal_si = 4184.0; // [J] kilocalorie (heat energy
// involved in warming up one // involved in warming up one
// kilogram of water by one // kilogram of water by one
// degree Kelvin) // degree Kelvin)
double const ev_si = 1.6021766208e-19; // [J] electron volt (amount of static double constexpr ev_si = 1.6021766208e-19; // [J] electron volt (amount of
// energy gained or lost by the // energy gained or lost by the
// charge of a single electron // charge of a single electron
// moving across an electric // moving across an electric
// potential difference of one // potential difference of one
// volt.) (NIST value) // volt.) (NIST value)
double const hartree_si = 4.359744650e-18; // [J] Hartree (approximately the static double constexpr hartree_si = 4.359744650e-18; // [J] Hartree (approximately the
// electric potential energy of // electric potential energy of
// the hydrogen atom in its // the hydrogen atom in its
// ground state) (NIST value) // ground state) (NIST value)
double const kcal_per_mole_si = kcal_si/Nav;// [J] kcal/mole static double constexpr kcal_per_mole_si = kcal_si/Nav;// [J] kcal/mole
double const erg_si = 1e-7; // [J] erg static double constexpr erg_si = 1e-7; // [J] erg
double const dyne_centimeter_si = 1e-7; // [J[ dyne*centimeter static double constexpr dyne_centimeter_si = 1e-7; // [J[ dyne*centimeter
double const picogram_micrometersq_per_microsecondsq_si = static double constexpr picogram_micrometersq_per_microsecondsq_si =
picogram_si*pow(micrometer_si,2)/pow(microsecond_si,2); picogram_si*micrometer_si/microsecond_si*micrometer_si/microsecond_si;
// [J] picogram*micrometer^2/ // [J] picogram*micrometer^2/
// microsecond^2 // microsecond^2
double const attogram_nanometersq_per_nanosecondsq_si = static double constexpr attogram_nanometersq_per_nanosecondsq_si =
attogram_si*pow(nanometer_si,2)/pow(nanosecond_si,2); attogram_si*nanometer_si/nanosecond_si*nanometer_si/nanosecond_si;
// [J] attogram*nanometer^2/ // [J] attogram*nanometer^2/
// nanosecond^2 // nanosecond^2
/*---------------------- /*----------------------
Velocity units Velocity units
------------------------ */ ------------------------ */
double const angstrom_per_femtosecond_si = static double constexpr angstrom_per_femtosecond_si =
angstrom_si/femtosecond_si; // [m/s] Angstrom/femtosecond 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 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 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 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 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 bohr_si/atu_electron_si; // [m/s] bohr/atu
/*---------------------- /*----------------------
Force units 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) 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 ev_si/angstrom_si; // [N] eV/Angstrom
double const dyne_si = static double constexpr dyne_si =
dyne_centimeter_si/centimeter_si; // [N] dyne 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 hartree_si/bohr_si; // [N] hartree/bohr
double const picogram_micrometer_per_microsecondsq_si = static double constexpr picogram_micrometer_per_microsecondsq_si =
picogram_si*micrometer_si/pow(microsecond_si,2); picogram_si*micrometer_si/microsecond_si/microsecond_si;
// [N] picogram*micrometer/ // [N] picogram*micrometer/
// microsecond^2 // microsecond^2
double const attogram_nanometer_per_nanosecondsq_si = static double constexpr attogram_nanometer_per_nanosecondsq_si =
attogram_si*nanometer_si/pow(nanosecond_si,2); attogram_si*nanometer_si/nanosecond_si/nanosecond_si;
// [N] attogram*nanometer/ // [N] attogram*nanometer/
// nanosecond^2 // nanosecond^2
/*---------------------- /*----------------------
Pressure units Pressure units
------------------------ */ ------------------------ */
double const atmosphere_si = 101325.0; // [Pa] standard atmosphere (NIST value) static double constexpr atmosphere_si = 101325.0; // [Pa] standard atmosphere (NIST value)
double const bar_si = 1e5; // [Pa] bar static double constexpr bar_si = 1e5; // [Pa] bar
double const dyne_per_centimetersq_si = static double constexpr dyne_per_centimetersq_si =
dyne_centimeter_si/pow(centimeter_si,3); dyne_centimeter_si/centimeter_si/centimeter_si/centimeter_si;
// [Pa] dyne/centimeter^2 // [Pa] dyne/centimeter^2
double const picogram_per_micrometer_microsecondsq_si = static double constexpr picogram_per_micrometer_microsecondsq_si =
picogram_si/(micrometer_si*pow(microsecond_si,2)); picogram_si/(micrometer_si*microsecond_si*microsecond_si);
// [Pa] picogram/(micrometer* // [Pa] picogram/(micrometer*
// microsecond^2) // microsecond^2)
double const attogram_per_nanometer_nanosecondsq_si = static double constexpr attogram_per_nanometer_nanosecondsq_si =
attogram_si/(nanometer_si*pow(nanosecond_si,2)); attogram_si/(nanometer_si*nanosecond_si*nanosecond_si);
// [Pa] attogram/(nanometer*nanosecond^2) // [Pa] attogram/(nanometer*nanosecond^2)
/*---------------------- /*----------------------
Viscosity units Viscosity units
------------------------ */ ------------------------ */
double const poise_si = 0.1; // [Pa*s] Poise static double constexpr poise_si = 0.1; // [Pa*s] Poise
double const amu_per_bohr_femtosecond_si = static double constexpr amu_per_bohr_femtosecond_si =
amu_si/(bohr_si*femtosecond_si); // [Pa*s] amu/(bohr*femtosecond) 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); picogram_si/(micrometer_si*microsecond_si);
// [Pa*s] picogram/(micrometer* // [Pa*s] picogram/(micrometer*
// microsecond) // microsecond)
double const attogram_per_nanometer_nanosecond_si = static double constexpr attogram_per_nanometer_nanosecond_si =
attogram_si/(nanometer_si*nanosecond_si); attogram_si/(nanometer_si*nanosecond_si);
// [Pa*s] attogram/(nanometer* // [Pa*s] attogram/(nanometer*
// nanosecond) // nanosecond)
@ -234,39 +236,41 @@ double const attogram_per_nanometer_nanosecond_si =
/*---------------------- /*----------------------
Charge units Charge units
------------------------ */ ------------------------ */
double const echarge_si = e_si; // [C] electron charge unit static double constexpr echarge_si = e_si; // [C] electron charge unit
double const statcoulomb_si = e_si/4.8032044e-10; // [C] Statcoulomb or esu static double constexpr statcoulomb_si = e_si/4.8032044e-10; // [C] Statcoulomb or esu
// (value from LAMMPS units // (value from LAMMPS units
// documentation) // documentation)
double const picocoulomb_si = 1e-12; // [C] picocoulomb static double constexpr picocoulomb_si = 1e-12; // [C] picocoulomb
/*---------------------- /*----------------------
Dipole units 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 // [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 // [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 // [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 // [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 // [C*m] electron*nanometer
/*---------------------- /*----------------------
Electric field units Electric field units
------------------------ */ ------------------------ */
double const volt_per_angstrom_si = 1.0/angstrom_si;// [V/m] volt/angstrom static double constexpr volt_per_angstrom_si = 1.0/angstrom_si;// [V/m] volt/angstrom
double const statvolt_per_centimeter_si = static double constexpr statvolt_per_centimeter_si =
erg_si/(statcoulomb_si*centimeter_si); // [V/m] statvolt/centimeter 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 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 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 1.0/nanometer_si; // [V/m] volt/nanometer
namespace LAMMPS_NS {
namespace kim_units {
// Define enumerations // Define enumerations
enum sys_type enum sys_type
{ {
@ -388,20 +392,42 @@ enum units
attogram_per_nanometercu = 1405 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 // Define dictionaries
map<string, sys_type> system_dic; static map<string, sys_type> system_dic;
map<string, unit_type> unit_dic; static map<string, unit_type> unit_dic;
map<unit_type, units> units_real_dic; static map<unit_type, units> units_real_dic;
map<unit_type, units> units_metal_dic; static map<unit_type, units> units_metal_dic;
map<unit_type, units> units_si_dic; static map<unit_type, units> units_si_dic;
map<unit_type, units> units_cgs_dic; static map<unit_type, units> units_cgs_dic;
map<unit_type, units> units_electron_dic; static map<unit_type, units> units_electron_dic;
map<unit_type, units> units_micro_dic; static map<unit_type, units> units_micro_dic;
map<unit_type, units> units_nano_dic; static map<unit_type, units> units_nano_dic;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void initialize_dictionaries() void kim_units::initialize_dictionaries()
{ {
system_dic["real"] = real; system_dic["real"] = real;
system_dic["metal"] = metal; system_dic["metal"] = metal;
@ -537,7 +563,7 @@ void initialize_dictionaries()
// Get the enumeration for the unit of type `unit_type_enum` // Get the enumeration for the unit of type `unit_type_enum`
// for LAMMPS system `system_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) { switch(system_enum) {
case real : case real :
@ -561,7 +587,7 @@ units get_lammps_system_unit(sys_type system_enum, unit_type unit_type_enum)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
// Mass conversion // 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; map<units, map<units, double> > conv;
double to_si; double to_si;
@ -619,7 +645,7 @@ double get_mass_conversion_factor(units from_unit_enum, units to_unit_enum)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
// Distance conversion // 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; map<units, map<units, double> > conv;
double to_si; double to_si;
@ -677,7 +703,7 @@ double get_distance_conversion_factor(units from_unit_enum, units to_unit_enum)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
// Time conversion // 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; map<units, map<units, double> > conv;
double to_si; double to_si;
@ -722,7 +748,7 @@ double get_time_conversion_factor(units from_unit_enum, units to_unit_enum)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
// Energy conversion // 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; 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 // 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; map<units, map<units, double> > conv;
double to_si; double to_si;
@ -869,7 +895,7 @@ double get_velocity_conversion_factor(units from_unit_enum, units to_unit_enum)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
// Force conversion // 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; map<units, map<units, double> > conv;
double to_si; double to_si;
@ -942,7 +968,7 @@ double get_force_conversion_factor(units from_unit_enum, units to_unit_enum)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
// Torque conversion // 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; map<units, map<units, double> > conv;
double to_si; double to_si;
@ -1015,7 +1041,7 @@ double get_torque_conversion_factor(units from_unit_enum, units to_unit_enum)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
// Temperature conversion // 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; 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 // 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; map<units, map<units, double> > conv;
double to_si; double to_si;
@ -1084,7 +1110,7 @@ double get_pressure_conversion_factor(units from_unit_enum, units to_unit_enum)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
// Viscosity conversion // 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; map<units, map<units, double> > conv;
double to_si; double to_si;
@ -1129,7 +1155,7 @@ double get_viscosity_conversion_factor(units from_unit_enum, units to_unit_enum)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
// Charge conversion // 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; map<units, map<units, double> > conv;
double to_si; double to_si;
@ -1163,7 +1189,7 @@ double get_charge_conversion_factor(units from_unit_enum, units to_unit_enum)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
// Dipole conversion // 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; map<units, map<units, double> > conv;
double to_si; double to_si;
@ -1221,7 +1247,7 @@ double get_dipole_conversion_factor(units from_unit_enum, units to_unit_enum)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
// Electric field conversion // 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; map<units, map<units, double> > conv;
double to_si; double to_si;
@ -1278,8 +1304,8 @@ double get_efield_conversion_factor(units from_unit_enum, units to_unit_enum)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
// Demsity conversion // Density conversion
double get_density_conversion_factor(units from_unit_enum, units to_unit_enum) double kim_units::get_density_conversion_factor(units from_unit_enum, units to_unit_enum)
{ {
map<units, map<units, double> > conv; map<units, map<units, double> > conv;
double to_si; 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 // This routine returns the unit conversion factor between the
// `from_system_enum` to the `to_system_enum` for the `unit_type_enum`. // `from_system_enum` to the `to_system_enum` for the `unit_type_enum`.
double get_unit_conversion_factor(unit_type &unit_type_enum, double kim_units::get_unit_conversion_factor(unit_type &unit_type_enum,
sys_type from_system_enum, sys_type from_system_enum,
sys_type to_system_enum) sys_type to_system_enum)
{ {
units from_unit = get_lammps_system_unit(from_system_enum, unit_type_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); 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 // Wrapper to the routine that gets the unit conversion. Translates strings
// to enumerations and then call get_unit_conversion_factor() // to enumerations and then call get_unit_conversion_factor()
int lammps_unit_conversion(const string &unit_type_str, int kim_units::lammps_unit_conversion(const string &unit_type_str, const string &from_system_str,
const string &from_system_str, const string &to_system_str, double &conversion_factor)
const string &to_system_str,
double &conversion_factor)
{ {
// initialize // initialize
conversion_factor = 0.0; conversion_factor = 0.0;
initialize_dictionaries(); initialize_dictionaries();
// convert input to enumeration // convert input to enumeration
unit_type unit_type_enum; unit_type unit_type_enum;
{ {
map<string, unit_type>::const_iterator itr = unit_dic.find(unit_type_str); map<string, unit_type>::const_iterator itr = unit_dic.find(unit_type_str);
if (itr != unit_dic.end()) unit_type_enum = itr->second; if (itr != unit_dic.end())
else return 1; // error unit_type_enum = itr->second;
} else
sys_type from_system_enum; return 1; // error
{ }
map<string, sys_type>::const_iterator sys_type from_system_enum;
itr = system_dic.find(from_system_str); {
if (itr != system_dic.end()) from_system_enum = itr->second; map<string, sys_type>::const_iterator itr = system_dic.find(from_system_str);
else return 1; // error if (itr != system_dic.end())
} from_system_enum = itr->second;
sys_type to_system_enum; else
{ return 1; // error
map<string, sys_type>::const_iterator }
itr = system_dic.find(to_system_str); sys_type to_system_enum;
if (itr != system_dic.end()) to_system_enum = itr->second; {
else return 1; 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 // process unit conversions
conversion_factor = get_unit_conversion_factor(unit_type_enum, conversion_factor = get_unit_conversion_factor(unit_type_enum, from_system_enum, to_system_enum);
from_system_enum, return 0;
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 #define LMP_KIM_UNITS_H
#include <string> #include <string>
extern "C" {
#include "KIM_SimulatorHeaders.h"
}
int lammps_unit_conversion(const std::string &unit_type_str, const std::string &from_system_str, namespace LAMMPS_NS {
const std::string &to_system_str, double &conversion_factor); 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 #endif

21
src/KIM/pair_kim.cpp
View File

@ -218,8 +218,7 @@ void PairKIM::compute(int eflag, int vflag)
kim_particleSpecies); kim_particleSpecies);
memory->create(kim_particleContributing,lmps_maxalloc, memory->create(kim_particleContributing,lmps_maxalloc,
"pair:kim_particleContributing"); "pair:kim_particleContributing");
kimerror = kimerror || KIM_ComputeArguments_SetArgumentPointerInteger( kimerror = kimerror || KIM_ComputeArguments_SetArgumentPointerInteger(pargs,
pargs,
KIM_COMPUTE_ARGUMENT_NAME_particleContributing, KIM_COMPUTE_ARGUMENT_NAME_particleContributing,
kim_particleContributing); kim_particleContributing);
if (kimerror) if (kimerror)
@ -248,7 +247,7 @@ void PairKIM::compute(int eflag, int vflag)
// compute via KIM model // compute via KIM model
int kimerror = KIM_Model_Compute(pkim, pargs); 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! // compute virial before reverse comm!
if (vflag_global) if (vflag_global)
@ -454,7 +453,7 @@ void PairKIM::coeff(int narg, char **arg)
kimerror = KIM_Model_GetParameterMetadata(pkim, param_index, &kim_DataType, kimerror = KIM_Model_GetParameterMetadata(pkim, param_index, &kim_DataType,
&extent, &str_name, &str_desc); &extent, &str_name, &str_desc);
if (kimerror) 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); const std::string str_name_str(str_name);
if (paramname == str_name_str) break; if (paramname == str_name_str) break;
@ -484,8 +483,8 @@ void PairKIM::coeff(int narg, char **arg)
nubound = atoi(words[1].c_str()); nubound = atoi(words[1].c_str());
if ((nubound < 1) || (nubound > extent) || (nlbound < 1) || (nlbound > nubound)) if ((nubound < 1) || (nubound > extent) || (nlbound < 1) || (nlbound > nubound))
error->all(FLERR,"Illegal index_range '{}-{}' for '{}' parameter with the extent of '{}'", error->all(FLERR,"Illegal index_range '{}-{}' for '{}' parameter with the extent "
nlbound, nubound, paramname, extent); "of '{}'", nlbound, nubound, paramname, extent);
} else { } else {
nlbound = atoi(argtostr.c_str()); nlbound = atoi(argtostr.c_str());
@ -975,7 +974,7 @@ void PairKIM::set_lmps_flags()
// determine if running with pair hybrid // determine if running with pair hybrid
if (force->pair_match("hybrid",0)) 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); const std::string unit_style_str(update->unit_style);
@ -1018,10 +1017,9 @@ void PairKIM::set_lmps_flags()
} else if ((unit_style_str == "lj") || } else if ((unit_style_str == "lj") ||
(unit_style_str == "micro") || (unit_style_str == "micro") ||
(unit_style_str == "nano")) { (unit_style_str == "nano")) {
error->all(FLERR,"LAMMPS unit_style {} not supported " error->all(FLERR,"LAMMPS unit_style {} not supported by KIM models", unit_style_str);
"by KIM models", unit_style_str);
} else { } 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 (comm->me == 0) {
if (KIM_SupportStatus_Equal(kim_model_support_for_energy, KIM_SUPPORT_STATUS_notSupported)) 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, if (KIM_SupportStatus_Equal(kim_model_support_for_forces,
KIM_SUPPORT_STATUS_notSupported)) KIM_SUPPORT_STATUS_notSupported))

38
src/KOKKOS/fix_property_atom_kokkos.cpp
View File

@ -24,8 +24,6 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
enum{MOLECULE,CHARGE,RMASS,INTEGER,DOUBLE};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixPropertyAtomKokkos::FixPropertyAtomKokkos(LAMMPS *lmp, int narg, char **arg) : 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) void FixPropertyAtomKokkos::grow_arrays(int nmax)
{ {
for (int m = 0; m < nvalue; m++) { for (int nv = 0; nv < nvalue; nv++) {
if (styles[m] == MOLECULE) { if (styles[nv] == MOLECULE) {
memory->grow(atom->molecule,nmax,"atom:molecule"); memory->grow(atom->molecule,nmax,"atom:molecule");
size_t nbytes = (nmax-nmax_old) * sizeof(tagint); size_t nbytes = (nmax-nmax_old) * sizeof(tagint);
memset(&atom->molecule[nmax_old],0,nbytes); memset(&atom->molecule[nmax_old],0,nbytes);
} else if (styles[m] == CHARGE) { } else if (styles[nv] == CHARGE) {
memory->grow(atom->q,nmax,"atom:q"); memory->grow(atom->q,nmax,"atom:q");
size_t nbytes = (nmax-nmax_old) * sizeof(double); size_t nbytes = (nmax-nmax_old) * sizeof(double);
memset(&atom->q[nmax_old],0,nbytes); memset(&atom->q[nmax_old],0,nbytes);
} else if (styles[m] == RMASS) { } else if (styles[nv] == RMASS) {
memory->grow(atom->rmass,nmax,"atom:rmass"); memory->grow(atom->rmass,nmax,"atom:rmass");
size_t nbytes = (nmax-nmax_old) * sizeof(double); size_t nbytes = (nmax-nmax_old) * sizeof(double);
memset(&atom->rmass[nmax_old],0,nbytes); memset(&atom->rmass[nmax_old],0,nbytes);
} else if (styles[m] == INTEGER) { } else if (styles[nv] == TEMPERATURE) {
memory->grow(atom->ivector[index[m]],nmax,"atom:ivector"); 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); size_t nbytes = (nmax-nmax_old) * sizeof(int);
memset(&atom->ivector[index[m]][nmax_old],0,nbytes); memset(&atom->ivector[index[nv]][nmax_old],0,nbytes);
} else if (styles[m] == DOUBLE) { } else if (styles[nv] == DVEC) {
atomKK->sync(Device,DVECTOR_MASK); atomKK->sync(Device,DVECTOR_MASK);
memoryKK->grow_kokkos(atomKK->k_dvector,atomKK->dvector,atomKK->k_dvector.extent(0),nmax, memoryKK->grow_kokkos(atomKK->k_dvector,atomKK->dvector,atomKK->k_dvector.extent(0),nmax,
"atom:dvector"); "atom:dvector");
atomKK->modified(Device,DVECTOR_MASK); atomKK->modified(Device,DVECTOR_MASK);
//memory->grow(atom->dvector[index[m]],nmax,"atom:dvector"); } else if (styles[nv] == IARRAY) {
//size_t nbytes = (nmax-nmax_old) * sizeof(double); memory->grow(atom->iarray[index[nv]], nmax, cols[nv], "atom:iarray");
//memset(&atom->dvector[index[m]][nmax_old],0,nbytes); 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; 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; datamask_modify = EMPTY_MASK;
k_error_flag = DAT::tdual_int_scalar("FixRxKokkos::k_error_flag"); k_error_flag = DAT::tdual_int_scalar("FixRxKokkos::k_error_flag");
//printf("Inside FixRxKokkos::FixRxKokkos\n");
} }
template <typename DeviceType> template <typename DeviceType>
FixRxKokkos<DeviceType>::~FixRxKokkos() FixRxKokkos<DeviceType>::~FixRxKokkos()
{ {
//printf("Inside FixRxKokkos::~FixRxKokkos copymode= %d\n", copymode);
if (copymode) return; if (copymode) return;
if (localTempFlag) if (localTempFlag)
memoryKK->destroy_kokkos(k_dpdThetaLocal, dpdThetaLocal); memoryKK->destroy_kokkos(k_dpdThetaLocal, dpdThetaLocal);
memoryKK->destroy_kokkos(k_sumWeights, sumWeights); memoryKK->destroy_kokkos(k_sumWeights, sumWeights);
//memoryKK->destroy_kokkos(k_sumWeights);
//delete [] scratchSpace;
memoryKK->destroy_kokkos(d_scratchSpace); memoryKK->destroy_kokkos(d_scratchSpace);
memoryKK->destroy_kokkos(k_cutsq); memoryKK->destroy_kokkos(k_cutsq);
@ -147,7 +141,6 @@ void FixRxKokkos<DeviceType>::init()
template <typename DeviceType> template <typename DeviceType>
void FixRxKokkos<DeviceType>::init_list(int, class NeighList* ptr) void FixRxKokkos<DeviceType>::init_list(int, class NeighList* ptr)
{ {
//printf("Inside FixRxKokkos::init_list\n");
this->list = ptr; this->list = ptr;
} }
@ -551,7 +544,6 @@ void FixRxKokkos<DeviceType>::k_rkf45(const int neq, const double t_stop, Vector
nfe += 6; nfe += 6;
if (maxIters && nit > maxIters) { 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 ++; counter.nFails ++;
break; break;
// We should set an error here so that the solution is not used! // 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.nSteps += nst;
counter.nIters += nit; counter.nIters += nit;
counter.nFuncs += nfe; 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) ); 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? // should we accept this?
if (hnew_is_ok || iter == max_iters) { if (hnew_is_ok || iter == max_iters) {
hnew = hg; hnew = hg;
@ -760,8 +747,6 @@ int FixRxKokkos<DeviceType>::rkf45_h0(const int neq, const double t, const doubl
hnew_is_ok = true; hnew_is_ok = true;
} }
//printf("iter=%d, yddnrw=%e, hnew=%e, hmin=%e, hmax=%e\n", iter, yddnrm, hnew, hmin, hmax);
hg = hnew; hg = hnew;
iter ++; iter ++;
} }
@ -805,13 +790,9 @@ void FixRxKokkos<DeviceType>::rkf45(const int neq, const double t_stop, double *
double t = 0.0; double t = 0.0;
if (h < h_min) { 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); 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. // Integrate until we reach the end time.
while (fabs(t - t_stop) > tround) { while (fabs(t - t_stop) > tround) {
double *yout = rwork; double *yout = rwork;
@ -865,7 +846,6 @@ void FixRxKokkos<DeviceType>::rkf45(const int neq, const double t_stop, double *
nfe += 6; nfe += 6;
if (maxIters && nit > maxIters) { 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 ++; counter.nFails ++;
break; break;
// We should set an error here so that the solution is not used! // 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.nSteps += nst;
counter.nIters += nit; counter.nIters += nit;
counter.nFuncs += nfe; 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 *rxnRateLaw = userData->rxnRateLaw;
double *kFor = userData->kFor; 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++) for (int ispecies=0; ispecies<nspecies; ispecies++)
dydt[ispecies] = 0.0; 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++) { for (int ispecies=0; ispecies<nspecies; ispecies++) {
const double concentration = y[ispecies]/VDPD; const double concentration = y[ispecies]/VDPD;
rxnRateLawForward *= pow( concentration, d_kineticsData.stoichReactants(jrxn,ispecies) ); rxnRateLawForward *= pow( concentration, d_kineticsData.stoichReactants(jrxn,ispecies) );
//rxnRateLawForward *= pow(concentration,stoichReactants[jrxn][ispecies]);
} }
rxnRateLaw[jrxn] = rxnRateLawForward; 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++) for (int jrxn=0; jrxn<nreactions; jrxn++)
{ {
dydt[ispecies] += d_kineticsData.stoich(jrxn,ispecies) *VDPD*rxnRateLaw[jrxn]; dydt[ispecies] += d_kineticsData.stoich(jrxn,ispecies) *VDPD*rxnRateLaw[jrxn];
//dydt[ispecies] += stoich[jrxn][ispecies]*VDPD*rxnRateLaw[jrxn];
} }
return 0; return 0;
@ -938,8 +911,6 @@ int FixRxKokkos<DeviceType>::rhs_sparse(double /*t*/, const double *y, double *d
{ {
UserRHSData *userData = (UserRHSData *) v_params; UserRHSData *userData = (UserRHSData *) v_params;
//const double VDPD = domain->xprd * domain->yprd * domain->zprd / atom->natoms;
#define kFor (userData->kFor) #define kFor (userData->kFor)
#define kRev (nullptr) #define kRev (nullptr)
#define rxnRateLaw (userData->rxnRateLaw) #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) { for (int kk = 1; kk < maxReactants; ++kk) {
const int k = nuk(i,kk); const int k = nuk(i,kk);
if (k == SparseKinetics_invalidIndex) break; if (k == SparseKinetics_invalidIndex) break;
//if (k != SparseKinetics_invalidIndex)
rxnRateLawForward *= powint( conc[k], inu(i,kk) ); rxnRateLawForward *= powint( conc[k], inu(i,kk) );
} }
} else { } else {
@ -972,7 +942,6 @@ int FixRxKokkos<DeviceType>::rhs_sparse(double /*t*/, const double *y, double *d
for (int kk = 1; kk < maxReactants; ++kk) { for (int kk = 1; kk < maxReactants; ++kk) {
const int k = nuk(i,kk); const int k = nuk(i,kk);
if (k == SparseKinetics_invalidIndex) break; if (k == SparseKinetics_invalidIndex) break;
//if (k != SparseKinetics_invalidIndex)
rxnRateLawForward *= pow( conc[k], nu(i,kk) ); 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) { for (int kk = 1; kk < maxReactants; ++kk) {
const int k = nuk(i,kk); const int k = nuk(i,kk);
if (k == SparseKinetics_invalidIndex) break; if (k == SparseKinetics_invalidIndex) break;
//if (k != SparseKinetics_invalidIndex)
dydt[k] -= nu(i,kk) * rxnRateLaw[i]; 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) { for (int kk = maxReactants+1; kk < maxSpecies; ++kk) {
const int k = nuk(i,kk); const int k = nuk(i,kk);
if (k == SparseKinetics_invalidIndex) break; if (k == SparseKinetics_invalidIndex) break;
//if (k != SparseKinetics_invalidIndex)
dydt[k] += nu(i,kk) * rxnRateLaw[i]; 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 rxnRateLaw (userData.rxnRateLaw)
#define kFor (userData.kFor ) #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++) for (int ispecies=0; ispecies<nspecies; ispecies++)
dydt[ispecies] = 0.0; 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++) { for (int ispecies=0; ispecies<nspecies; ispecies++) {
const double concentration = y[ispecies]/VDPD; const double concentration = y[ispecies]/VDPD;
rxnRateLawForward *= pow( concentration, d_kineticsData.stoichReactants(jrxn,ispecies) ); rxnRateLawForward *= pow( concentration, d_kineticsData.stoichReactants(jrxn,ispecies) );
//rxnRateLawForward *= pow(concentration,stoichReactants[jrxn][ispecies]);
} }
rxnRateLaw[jrxn] = rxnRateLawForward; 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++) for (int jrxn=0; jrxn<nreactions; jrxn++)
{ {
dydt[ispecies] += d_kineticsData.stoich(jrxn,ispecies) *VDPD*rxnRateLaw[jrxn]; dydt[ispecies] += d_kineticsData.stoich(jrxn,ispecies) *VDPD*rxnRateLaw[jrxn];
//dydt[ispecies] += stoich[jrxn][ispecies]*VDPD*rxnRateLaw[jrxn];
} }
#undef rxnRateLaw #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) { for (int kk = 1; kk < maxReactants; ++kk) {
const int k = nuk(i,kk); const int k = nuk(i,kk);
if (k == SparseKinetics_invalidIndex) break; if (k == SparseKinetics_invalidIndex) break;
//if (k != SparseKinetics_invalidIndex)
rxnRateLawForward *= powint( conc[k], inu(i,kk) ); rxnRateLawForward *= powint( conc[k], inu(i,kk) );
} }
} else { } else {
@ -1119,7 +1080,6 @@ int FixRxKokkos<DeviceType>::k_rhs_sparse(double /*t*/, const VectorType& y, Vec
for (int kk = 1; kk < maxReactants; ++kk) { for (int kk = 1; kk < maxReactants; ++kk) {
const int k = nuk(i,kk); const int k = nuk(i,kk);
if (k == SparseKinetics_invalidIndex) break; if (k == SparseKinetics_invalidIndex) break;
//if (k != SparseKinetics_invalidIndex)
rxnRateLawForward *= pow( conc[k], nu(i,kk) ); 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) { for (int kk = 1; kk < maxReactants; ++kk) {
const int k = nuk(i,kk); const int k = nuk(i,kk);
if (k == SparseKinetics_invalidIndex) break; if (k == SparseKinetics_invalidIndex) break;
//if (k != SparseKinetics_invalidIndex)
dydt[k] -= nu(i,kk) * rxnRateLaw[i]; 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) { for (int kk = maxReactants+1; kk < maxSpecies; ++kk) {
const int k = nuk(i,kk); const int k = nuk(i,kk);
if (k == SparseKinetics_invalidIndex) break; if (k == SparseKinetics_invalidIndex) break;
//if (k != SparseKinetics_invalidIndex)
dydt[k] += nu(i,kk) * rxnRateLaw[i]; dydt[k] += nu(i,kk) * rxnRateLaw[i];
} }
} }
@ -1215,8 +1173,6 @@ void FixRxKokkos<DeviceType>::operator()(SolverType, const int &i) const
template <typename DeviceType> template <typename DeviceType>
void FixRxKokkos<DeviceType>::create_kinetics_data() 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.Arr, h_kineticsData.Arr, nreactions, "KineticsType::Arr");
memoryKK->create_kokkos( d_kineticsData.nArr, h_kineticsData.nArr, nreactions, "KineticsType::nArr"); memoryKK->create_kokkos( d_kineticsData.nArr, h_kineticsData.nArr, nreactions, "KineticsType::nArr");
memoryKK->create_kokkos( d_kineticsData.Ea, h_kineticsData.Ea, nreactions, "KineticsType::Ea"); 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> template <typename DeviceType>
void FixRxKokkos<DeviceType>::setup_pre_force(int vflag) void FixRxKokkos<DeviceType>::setup_pre_force(int vflag)
{ {
//printf("Inside FixRxKokkos<DeviceType>::setup_pre_force restartFlag= %d\n", my_restartFlag);
if (my_restartFlag) if (my_restartFlag)
my_restartFlag = 0; my_restartFlag = 0;
else else
@ -1309,8 +1263,6 @@ void FixRxKokkos<DeviceType>::setup_pre_force(int vflag)
template <typename DeviceType> template <typename DeviceType>
void FixRxKokkos<DeviceType>::pre_force(int vflag) void FixRxKokkos<DeviceType>::pre_force(int vflag)
{ {
//printf("Inside FixRxKokkos<DeviceType>::pre_force localTempFlag= %d\n", localTempFlag);
this->solve_reactions( vflag, true ); this->solve_reactions( vflag, true );
} }
@ -1407,8 +1359,6 @@ void FixRxKokkos<DeviceType>::operator()(Tag_FixRxKokkos_solveSystems<ZERO_RATES
template <typename DeviceType> template <typename DeviceType>
void FixRxKokkos<DeviceType>::solve_reactions(const int /*vflag*/, const bool isPreForce) 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; copymode = 1;
if (update_kinetics_data) if (update_kinetics_data)
@ -1416,7 +1366,6 @@ void FixRxKokkos<DeviceType>::solve_reactions(const int /*vflag*/, const bool is
//TimerType timer_start = getTimeStamp(); //TimerType timer_start = getTimeStamp();
//const int nlocal = atom->nlocal;
this->nlocal = atom->nlocal; this->nlocal = atom->nlocal;
const int nghost = atom->nghost; const int nghost = atom->nghost;
const int newton_pair = force->newton_pair; 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. // 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_dpdTheta = atomKK->k_dpdTheta.view<DeviceType>();
this->d_dvector = atomKK->k_dvector.view<DeviceType>(); this->d_dvector = atomKK->k_dvector.view<DeviceType>();
this->d_mask = atomKK->k_mask.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 ); atomKK->sync( execution_space, MASK_MASK | DVECTOR_MASK | DPDTHETA_MASK );
// Set some constants outside of the parallel_for // 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->boltz = force->boltz;
this->t_stop = update->dt; // DPD time-step and integration length. 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>(); d_diagnosticCounterPerODEnFuncs = k_diagnosticCounterPerODEnFuncs.template view<DeviceType>();
Kokkos::parallel_for ( Kokkos::RangePolicy<DeviceType, Tag_FixRxKokkos_zeroCounterViews>(0,nlocal), *this); 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. // 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>(); k_error_flag.template sync<DeviceType>();
// Create scratch array space. // Create scratch array space.
//const size_t scratchSpaceSize = (8*nspecies + 2*nreactions);
this->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)) { if (nlocal*scratchSpaceSize > d_scratchSpace.extent(0)) {
memoryKK->destroy_kokkos (d_scratchSpace); memoryKK->destroy_kokkos (d_scratchSpace);
memoryKK->create_kokkos (d_scratchSpace, nlocal*scratchSpaceSize, "FixRxKokkos::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) if (setRatesToZero)
Kokkos::parallel_reduce( Kokkos::RangePolicy<DeviceType, Tag_FixRxKokkos_solveSystems<true > >(0,nlocal), *this, TotalCounters); Kokkos::parallel_reduce( Kokkos::RangePolicy<DeviceType, Tag_FixRxKokkos_solveSystems<true > >(0,nlocal), *this, TotalCounters);
else else
Kokkos::parallel_reduce( Kokkos::RangePolicy<DeviceType, Tag_FixRxKokkos_solveSystems<false> >(0,nlocal), *this, TotalCounters); Kokkos::parallel_reduce( Kokkos::RangePolicy<DeviceType, Tag_FixRxKokkos_solveSystems<false> >(0,nlocal), *this, TotalCounters);
#endif
TimerType timer_ODE = getTimeStamp(); TimerType timer_ODE = getTimeStamp();
@ -1652,16 +1490,8 @@ void FixRxKokkos<DeviceType>::solve_reactions(const int /*vflag*/, const bool is
atomKK->modified ( Host, DVECTOR_MASK ); atomKK->modified ( Host, DVECTOR_MASK );
//TimerType timer_stop = getTimeStamp();
double time_ODE = getElapsedTime(timer_localTemperature, timer_ODE); 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. // Warn the user if a failure was detected in the ODE solver.
if (TotalCounters.nFails > 0) { if (TotalCounters.nFails > 0) {
char sbuf[128]; char sbuf[128];
@ -1731,7 +1561,6 @@ void FixRxKokkos<DeviceType>::odeDiagnostics()
// Compute counters per dpd time-step. // Compute counters per dpd time-step.
for (int i = 0; i < numCounters; ++i) { for (int i = 0; i < numCounters; ++i) {
my_vals[i] = this->diagnosticCounter[i] / nTimes; 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); 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]; double my_max[numCounters], my_min[numCounters];
//const int nlocal = atom->nlocal;
nlocal = atom->nlocal; nlocal = atom->nlocal;
HAT::t_int_1d h_mask = atomKK->k_mask.h_view; 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> template <int WT_FLAG, int LOCAL_TEMP_FLAG, bool NEWTON_PAIR, int NEIGHFLAG>
void FixRxKokkos<DeviceType>::computeLocalTemperature() 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_x = atomKK->k_x.view<DeviceType>();
d_type = atomKK->k_type.view<DeviceType>(); d_type = atomKK->k_type.view<DeviceType>();
d_dpdTheta = atomKK->k_dpdTheta.view<DeviceType>(); d_dpdTheta = atomKK->k_dpdTheta.view<DeviceType>();
@ -1993,21 +1818,9 @@ void FixRxKokkos<DeviceType>::computeLocalTemperature()
nlocal = atom->nlocal; nlocal = atom->nlocal;
const int nghost = atom->nghost; 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; 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)) { if (ntypes+1 > (int) k_cutsq.extent(0)) {
memoryKK->destroy_kokkos (k_cutsq); memoryKK->destroy_kokkos (k_cutsq);
memoryKK->create_kokkos (k_cutsq, ntypes+1, ntypes+1, "FixRxKokkos::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 // Initialize the local temperature weight array
int sumWeightsCt = nlocal + (NEWTON_PAIR ? nghost : 0); 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)) { if (sumWeightsCt > (int)k_sumWeights.template view<DeviceType>().extent(0)) {
memoryKK->destroy_kokkos(k_sumWeights, sumWeights); memoryKK->destroy_kokkos(k_sumWeights, sumWeights);
memoryKK->create_kokkos (k_sumWeights, sumWeightsCt, "FixRxKokkos::sumWeights"); memoryKK->create_kokkos (k_sumWeights, sumWeightsCt, "FixRxKokkos::sumWeights");
@ -2036,14 +1848,6 @@ void FixRxKokkos<DeviceType>::computeLocalTemperature()
h_sumWeights = k_sumWeights.h_view; 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); Kokkos::parallel_for (Kokkos::RangePolicy<DeviceType, Tag_FixRxKokkos_zeroTemperatureViews>(0, sumWeightsCt), *this);
// Local list views. (This isn't working!) // Local list views. (This isn't working!)
@ -2051,9 +1855,6 @@ void FixRxKokkos<DeviceType>::computeLocalTemperature()
if (!list->kokkos) if (!list->kokkos)
error->one(FLERR,"list is not a Kokkos list\n"); 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_neighbors = k_list->d_neighbors;
d_ilist = k_list->d_ilist; d_ilist = k_list->d_ilist;
d_numneigh = k_list->d_numneigh; d_numneigh = k_list->d_numneigh;
@ -2067,7 +1868,6 @@ void FixRxKokkos<DeviceType>::computeLocalTemperature()
{ {
// Create an atomic view of sumWeights and dpdThetaLocal. Only needed // Create an atomic view of sumWeights and dpdThetaLocal. Only needed
// for Half/thread scenarios. // 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; 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; AtomicViewType a_dpdThetaLocal = d_dpdThetaLocal;
@ -2174,8 +1974,6 @@ void FixRxKokkos<DeviceType>::computeLocalTemperature()
template <typename DeviceType> template <typename DeviceType>
int FixRxKokkos<DeviceType>::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) 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; HAT::t_float_2d h_dvector = atomKK->k_dvector.h_view;
int m = 0; 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); buf[m++] = h_dvector(ispecies+nspecies,jj);
} }
} }
//printf("done with FixRxKokkos::pack_forward_comm %d\n", comm->me);
return m; return m;
} }
@ -2197,8 +1992,6 @@ int FixRxKokkos<DeviceType>::pack_forward_comm(int n, int *list, double *buf, in
template <typename DeviceType> template <typename DeviceType>
void FixRxKokkos<DeviceType>::unpack_forward_comm(int n, int first, double *buf) 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; HAT::t_float_2d h_dvector = atomKK->k_dvector.h_view;
const int last = first + n ; 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++]; 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> template <typename DeviceType>
int FixRxKokkos<DeviceType>::pack_reverse_comm(int n, int first, double *buf) 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. // Sync the host view.
k_dpdThetaLocal.template sync<LMPHostType>(); k_dpdThetaLocal.template sync<LMPHostType>();
k_sumWeights. 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_dpdThetaLocal(i);
buf[m++] = h_sumWeights(i); buf[m++] = h_sumWeights(i);
} }
//printf("done with FixRxKokkos::pack_reverse_comm %d\n", comm->me);
return m; return m;
} }
@ -2240,7 +2028,6 @@ int FixRxKokkos<DeviceType>::pack_reverse_comm(int n, int first, double *buf)
template <typename DeviceType> template <typename DeviceType>
void FixRxKokkos<DeviceType>::unpack_reverse_comm(int n, int *list, double *buf) 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; int m = 0;
for (int i = 0; i < n; i++) { for (int i = 0; i < n; i++) {
const int j = list[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. // Signal that the host view has been modified.
k_dpdThetaLocal.template modify<LMPHostType>(); k_dpdThetaLocal.template modify<LMPHostType>();
k_sumWeights. template modify<LMPHostType>(); k_sumWeights. template modify<LMPHostType>();
// printf("done with FixRxKokkos::unpack_reverse_comm %d\n", comm->me);
} }
namespace LAMMPS_NS { 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 PairEAM::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
int /*pbc_flag*/, int * /*pbc*/)
{ {
int i,j,m; int i,j,m;

34
src/RIGID/fix_rigid.cpp
View File

@ -2523,9 +2523,17 @@ int FixRigid::pack_exchange(int i, double *buf)
buf[2] = displace[i][0]; buf[2] = displace[i][0];
buf[3] = displace[i][1]; buf[3] = displace[i][1];
buf[4] = displace[i][2]; buf[4] = displace[i][2];
if (!extended) return 5;
// must also pack vatom if per-atom virial calculated on this timestep
// since vatom is calculated before and after atom migration
int m = 5; int m = 5;
if (vflag_atom)
for (int k = 0; k < 6; k++)
buf[m++] = vatom[i][k];
if (!extended) return m;
buf[m++] = eflags[i]; buf[m++] = eflags[i];
for (int j = 0; j < orientflag; j++) for (int j = 0; j < orientflag; j++)
buf[m++] = orient[i][j]; buf[m++] = orient[i][j];
@ -2535,13 +2543,6 @@ int FixRigid::pack_exchange(int i, double *buf)
buf[m++] = dorient[i][2]; buf[m++] = dorient[i][2];
} }
// must also pack vatom if per-atom virial calculated on this timestep
// since vatom is calculated before and after atom migration
if (vflag_atom)
for (int k = 0; k < 6; k++)
buf[m++] = vatom[i][k];
return m; return m;
} }
@ -2556,9 +2557,17 @@ int FixRigid::unpack_exchange(int nlocal, double *buf)
displace[nlocal][0] = buf[2]; displace[nlocal][0] = buf[2];
displace[nlocal][1] = buf[3]; displace[nlocal][1] = buf[3];
displace[nlocal][2] = buf[4]; displace[nlocal][2] = buf[4];
if (!extended) return 5;
// must also unpack vatom if per-atom virial calculated on this timestep
// since vatom is calculated before and after atom migration
int m = 5; int m = 5;
if (vflag_atom)
for (int k = 0; k < 6; k++)
vatom[nlocal][k] = buf[m++];
if (!extended) return m;
eflags[nlocal] = static_cast<int> (buf[m++]); eflags[nlocal] = static_cast<int> (buf[m++]);
for (int j = 0; j < orientflag; j++) for (int j = 0; j < orientflag; j++)
orient[nlocal][j] = buf[m++]; orient[nlocal][j] = buf[m++];
@ -2568,13 +2577,6 @@ int FixRigid::unpack_exchange(int nlocal, double *buf)
dorient[nlocal][2] = buf[m++]; dorient[nlocal][2] = buf[m++];
} }
// must also unpack vatom if per-atom virial calculated on this timestep
// since vatom is calculated before and after atom migration
if (vflag_atom)
for (int k = 0; k < 6; k++)
vatom[nlocal][k] = buf[m++];
return m; return m;
} }

5
src/compute_chunk_atom.cpp
View File

@ -691,9 +691,10 @@ void ComputeChunkAtom::compute_ichunk()
if (invoked_ichunk == update->ntimestep) return; if (invoked_ichunk == update->ntimestep) return;
// if old IDs persist via storage in fixstore, then just retrieve them // if old IDs persist via storage in fixstore, then just retrieve them
// yes if idsflag = ONCE, and already done once // restore = yes if idsflag = ONCE, and already done once
// or if idsflag = NFREQ and lock is in place and are on later timestep // or if idsflag = NFREQ and lock is in place and are on later timestep
// else proceed to recalculate per-atom chunk assignments // else proceed to recalculate per-atom chunk assignments
// if restoring, update invoked_ichunk only for NFREQ case
const int nlocal = atom->nlocal; const int nlocal = atom->nlocal;
int restore = 0; int restore = 0;
@ -701,7 +702,7 @@ void ComputeChunkAtom::compute_ichunk()
if (idsflag == NFREQ && lockfix && update->ntimestep > lockstart) restore = 1; if (idsflag == NFREQ && lockfix && update->ntimestep > lockstart) restore = 1;
if (restore) { if (restore) {
invoked_ichunk = update->ntimestep; if (idsflag == NFREQ) invoked_ichunk = update->ntimestep;
double *vstore = fixstore->vstore; double *vstore = fixstore->vstore;
for (i = 0; i < nlocal; i++) ichunk[i] = static_cast<int>(vstore[i]); for (i = 0; i < nlocal; i++) ichunk[i] = static_cast<int>(vstore[i]);
return; return;

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), Fix(lmp, narg, arg), nadapt(0), anypair(0), anybond(0), anyangle(0),
id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr) 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); 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; dynamic_group_allow = 1;
create_attribute = 1; create_attribute = 1;
@ -62,29 +62,29 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) :
int iarg = 4; int iarg = 4;
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg],"pair") == 0) { 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++; nadapt++;
iarg += 6; iarg += 6;
} else if (strcmp(arg[iarg],"kspace") == 0) { } 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++; nadapt++;
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"atom") == 0) { } 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++; nadapt++;
iarg += 3; iarg += 3;
} else if (strcmp(arg[iarg],"bond") == 0) { } 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++; nadapt++;
iarg += 5; iarg += 5;
} else if (strcmp(arg[iarg],"angle") == 0) { } 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++; nadapt++;
iarg += 5; iarg += 5;
} else break; } 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]; adapt = new Adapt[nadapt];
// parse keywords // parse keywords
@ -96,7 +96,6 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) :
iarg = 4; iarg = 4;
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg],"pair") == 0) { if (strcmp(arg[iarg],"pair") == 0) {
if (iarg+6 > narg) error->all(FLERR,"Illegal fix adapt command");
adapt[nadapt].which = PAIR; adapt[nadapt].which = PAIR;
adapt[nadapt].pair = nullptr; adapt[nadapt].pair = nullptr;
adapt[nadapt].pstyle = utils::strdup(arg[iarg+1]); 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); adapt[nadapt].jlo,adapt[nadapt].jhi,error);
if (utils::strmatch(arg[iarg+5],"^v_")) { if (utils::strmatch(arg[iarg+5],"^v_")) {
adapt[nadapt].var = utils::strdup(arg[iarg+5]+2); 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++; nadapt++;
iarg += 6; iarg += 6;
} else if (strcmp(arg[iarg],"bond") == 0) { } else if (strcmp(arg[iarg],"bond") == 0) {
if (iarg+5 > narg) error->all(FLERR, "Illegal fix adapt command");
adapt[nadapt].which = BOND; adapt[nadapt].which = BOND;
adapt[nadapt].bond = nullptr; adapt[nadapt].bond = nullptr;
adapt[nadapt].bstyle = utils::strdup(arg[iarg+1]); 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); adapt[nadapt].ilo,adapt[nadapt].ihi,error);
if (utils::strmatch(arg[iarg+4],"^v_")) { if (utils::strmatch(arg[iarg+4],"^v_")) {
adapt[nadapt].var = utils::strdup(arg[iarg+4]+2); 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++; nadapt++;
iarg += 5; iarg += 5;
} else if (strcmp(arg[iarg],"angle") == 0) { } else if (strcmp(arg[iarg],"angle") == 0) {
if (iarg+5 > narg) error->all(FLERR, "Illegal fix adapt command");
adapt[nadapt].which = ANGLE; adapt[nadapt].which = ANGLE;
adapt[nadapt].angle = nullptr; adapt[nadapt].angle = nullptr;
adapt[nadapt].astyle = utils::strdup(arg[iarg+1]); 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); adapt[nadapt].ilo,adapt[nadapt].ihi,error);
if (utils::strmatch(arg[iarg+4],"^v_")) { if (utils::strmatch(arg[iarg+4],"^v_")) {
adapt[nadapt].var = utils::strdup(arg[iarg+4]+2); 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++; nadapt++;
iarg += 5; iarg += 5;
} else if (strcmp(arg[iarg],"kspace") == 0) { } else if (strcmp(arg[iarg],"kspace") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command");
adapt[nadapt].which = KSPACE; adapt[nadapt].which = KSPACE;
if (utils::strmatch(arg[iarg+1],"^v_")) { if (utils::strmatch(arg[iarg+1],"^v_")) {
adapt[nadapt].var = utils::strdup(arg[iarg+1]+2); 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++; nadapt++;
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"atom") == 0) { } else if (strcmp(arg[iarg],"atom") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command");
adapt[nadapt].which = ATOM; adapt[nadapt].which = ATOM;
if (strcmp(arg[iarg+1],"diameter") == 0 || if (strcmp(arg[iarg+1],"diameter") == 0 ||
strcmp(arg[iarg+1],"diameter/disc") == 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) { } else if (strcmp(arg[iarg+1],"charge") == 0) {
adapt[nadapt].atomparam = CHARGE; adapt[nadapt].atomparam = CHARGE;
chgflag = 1; 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_")) { if (utils::strmatch(arg[iarg+2],"^v_")) {
adapt[nadapt].var = utils::strdup(arg[iarg+2]+2); 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++; nadapt++;
iarg += 3; iarg += 3;
} else break; } else break;
@ -177,18 +172,18 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) :
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg],"reset") == 0) { 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); resetflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"scale") == 0) { } 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); scaleflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"mass") == 0) { } 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); massflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2; 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, // if scaleflag set with diameter or charge adaptation,
@ -202,22 +197,19 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) :
int n = atom->ntypes; int n = atom->ntypes;
for (int m = 0; m < nadapt; m++) for (int m = 0; m < nadapt; m++)
if (adapt[m].which == PAIR) if (adapt[m].which == PAIR) memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig");
memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig");
// allocate bond style arrays: // allocate bond style arrays:
n = atom->nbondtypes; n = atom->nbondtypes;
for (int m = 0; m < nadapt; ++m) for (int m = 0; m < nadapt; ++m)
if (adapt[m].which == BOND) if (adapt[m].which == BOND) memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig");
memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig");
// allocate angle style arrays: // allocate angle style arrays:
n = atom->nbondtypes; n = atom->nbondtypes;
for (int m = 0; m < nadapt; ++m) for (int m = 0; m < nadapt; ++m)
if (adapt[m].which == ANGLE) if (adapt[m].which == ANGLE) memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig");
memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig");
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -324,7 +316,7 @@ void FixAdapt::init()
if (group->dynamic[igroup]) if (group->dynamic[igroup])
for (i = 0; i < nadapt; i++) for (i = 0; i < nadapt; i++)
if (adapt[i].which == ATOM) 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 // setup and error checks
@ -337,9 +329,9 @@ void FixAdapt::init()
ad->ivar = input->variable->find(ad->var); ad->ivar = input->variable->find(ad->var);
if (ad->ivar < 0) 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)) 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) { if (ad->which == PAIR) {
anypair = 1; anypair = 1;
@ -368,13 +360,13 @@ void FixAdapt::init()
void *ptr = ad->pair->extract(ad->pparam,ad->pdim); void *ptr = ad->pair->extract(ad->pparam,ad->pdim);
if (ptr == nullptr) 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 // for pair styles only parameters that are 2-d arrays in atom types or
// scalars are supported // scalars are supported
if (ad->pdim != 2 && ad->pdim != 0) 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 == 2) ad->array = (double **) ptr;
if (ad->pdim == 0) ad->scalar = (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) ad->bond = force->bond_match(bstyle);
if (ad->bond == nullptr ) 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); void *ptr = ad->bond->extract(ad->bparam,ad->bdim);
if (ptr == nullptr) 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 // 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) ad->angle = force->angle_match(astyle);
if (ad->angle == nullptr ) 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); void *ptr = ad->angle->extract(ad->aparam,ad->adim);
if (ptr == nullptr) 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 // for angle styles, use a vector
@ -446,7 +438,7 @@ void FixAdapt::init()
} else if (ad->which == KSPACE) { } else if (ad->which == KSPACE) {
if (force->kspace == nullptr) 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"); kspace_scale = (double *) force->kspace->extract("scale");
} else if (ad->which == ATOM) { } else if (ad->which == ATOM) {

2
src/fix_property_atom.cpp
View File

@ -25,8 +25,6 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
enum { MOLECULE, CHARGE, RMASS, TEMPERATURE, HEATFLOW, IVEC, DVEC, IARRAY, DARRAY };
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) : 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 { class FixPropertyAtom : public Fix {
public: public:
FixPropertyAtom(class LAMMPS *, int, char **); FixPropertyAtom(class LAMMPS *, int, char **);
~FixPropertyAtom() override; ~FixPropertyAtom() override;
int setmask() override; int setmask() override;
void init() override; void init() override;
enum { MOLECULE, CHARGE, RMASS, TEMPERATURE, HEATFLOW, IVEC, DVEC, IARRAY, DARRAY };
void read_data_section(char *, int, char *, tagint) override; void read_data_section(char *, int, char *, tagint) override;
bigint read_data_skip_lines(char *) override; bigint read_data_skip_lines(char *) override;
void write_data_section_size(int, int &, int &) override; void write_data_section_size(int, int &, int &) override;
@ -55,10 +56,10 @@ class FixPropertyAtom : public Fix {
int nvalue, border; int nvalue, border;
int molecule_flag, q_flag, rmass_flag; // flags for specific fields int molecule_flag, q_flag, rmass_flag; // flags for specific fields
int temperature_flag, heatflow_flag; int temperature_flag, heatflow_flag;
int *styles; // style of each value, see enum int *styles; // style of each value, see enum
int *index; // indices into atom custom data structs int *index; // indices into atom custom data structs
int *cols; // columns per value, for arrays int *cols; // columns per value, for arrays
char *astyle; // atom style at instantiation char *astyle; // atom style at instantiation
int values_peratom; // # of values per atom, including multiple for arrays int values_peratom; // # of values per atom, including multiple for arrays
int nmax_old; // length of peratom arrays the last time they grew int nmax_old; // length of peratom arrays the last time they grew

3
src/fix_update_special_bonds.cpp
View File

@ -72,6 +72,9 @@ void FixUpdateSpecialBonds::setup(int /*vflag*/)
force->special_coul[3] != 1.0) force->special_coul[3] != 1.0)
error->all(FLERR, "Fix update/special/bonds requires special Coulomb weights = 1,1,1"); error->all(FLERR, "Fix update/special/bonds requires special Coulomb weights = 1,1,1");
// Implies neighbor->special_flag = [X, 2, 1, 1] // Implies neighbor->special_flag = [X, 2, 1, 1]
if (utils::strmatch(force->pair_style, "^hybrid"))
error->all(FLERR, "Cannot use fix update/special/bonds with hybrid pair styles");
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

10
src/neighbor.cpp
View File

@ -95,6 +95,16 @@ static const char cite_neigh_multi[] =
" Detection Applied to Investigate the Quasi-Static Limit},\n" " Detection Applied to Investigate the Quasi-Static Limit},\n"
" journal = {Computational Particle Mechanics},\n" " journal = {Computational Particle Mechanics},\n"
" year = {2020}\n" " year = {2020}\n"
"@article{Monti2022,\n"
" author = {Monti, Joseph M. and Clemmer, Joel T. and Srivastava, \n"
" Ishan and Silbert, Leonardo E. and Grest, Gary S. \n"
" and Lechman, Jeremy B.},\n"
" title = {Large-scale frictionless jamming with power-law particle \n"
" size distributions},\n"
" journal = {Phys. Rev. E},\n"
" volume = {106}\n"
" issue = {3}\n"
" year = {2022}\n"
"}\n\n"; "}\n\n";
// template for factory functions: // template for factory functions:

4
tools/fep/fdti.py
View File

@ -10,7 +10,7 @@ if len(sys.argv) < 3:
print("usage: fdti.py temperature hderiv < out.fep") print("usage: fdti.py temperature hderiv < out.fep")
sys.exit() 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]) hderiv = float(sys.argv[2])
line = sys.stdin.readline() line = sys.stdin.readline()
@ -23,7 +23,7 @@ if len(tok) == 4:
v = float(tok[3]) v = float(tok[3])
lo = -rt * math.log(float(tok[2]) / v) lo = -rt * math.log(float(tok[2]) / v)
i = 0 i = 1
sum = 0.0 sum = 0.0
for line in sys.stdin: for line in sys.stdin:
tok = line.split() 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") print("usage: fep.py temperature < out.fep")
sys.exit() 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 v = 1.0
sum = 0.0 sum = 0.0

2
tools/fep/nti.py
View File

@ -19,7 +19,7 @@ while line.startswith("#"):
tok = line.split() tok = line.split()
lo = float(tok[1]) lo = float(tok[1])
i = 0 i = 1
sum = 0.0 sum = 0.0
for line in sys.stdin: for line in sys.stdin:
tok = line.split() tok = line.split()

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

@ -657,11 +657,14 @@ class dump:
atoms = snap.atoms atoms = snap.atoms
nvalues = len(atoms[0]) nvalues = len(atoms[0])
keys = dict()
for pair in self.names.items():
keys[pair[1]] = pair[0]
for i in range(snap.natoms): for i in range(snap.natoms):
if not snap.aselect[i]: continue if not snap.aselect[i]: continue
line = "" line = ""
for j in range(nvalues): 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])) + " " line += str(int(atoms[i][j])) + " "
else: else:
line += str(atoms[i][j]) + " " 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(); 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(); BEGIN_HIDE_OUTPUT();
command("kim init LennardJones_Ar real"); command("kim init LennardJones_Ar real");
END_HIDE_OUTPUT(); END_HIDE_OUTPUT();
TEST_FAILURE(".*ERROR: Must use 'kim interactions' command " TEST_FAILURE(".*ERROR: Use of 'kim interactions' before simulation box is defined.*",
"after simulation box is defined.*",
command("kim interactions Ar");); command("kim interactions Ar"););
BEGIN_HIDE_OUTPUT(); BEGIN_HIDE_OUTPUT();
@ -410,11 +410,11 @@ TEST_F(KimCommandsTest, kim_property)
"3 >= 3.6 support.*", "3 >= 3.6 support.*",
command("kim property");); command("kim property"););
} else { } else {
TEST_FAILURE(".*ERROR: Invalid 'kim property' command.*", command("kim property");); TEST_FAILURE(".*ERROR: Illegal kim property command: missing argument.*",
TEST_FAILURE(".*ERROR: Invalid 'kim property' command.*", command("kim property create");); command("kim property"););
TEST_FAILURE(".*ERROR: Incorrect arguments in 'kim property' command." TEST_FAILURE(".*ERROR: Illegal kim property command: missing argument.*",
"\n'kim property create/destroy/modify/remove/dump' " command("kim property create"););
"is mandatory.*", TEST_FAILURE(".*ERROR: Incorrect first argument unknown to 'kim property' command.*",
command("kim property unknown 1 atomic-mass");); command("kim property unknown 1 atomic-mass"););
} }
#if defined(KIM_EXTRA_UNITTESTS) #if defined(KIM_EXTRA_UNITTESTS)