From 120f5cf7f120f4aa30c677ddaef1cee77129dc9d Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 16 Mar 2023 15:20:18 -0600 Subject: [PATCH 01/22] Minor patches to BPM, multi, and rigid --- doc/src/bond_bpm_rotational.rst | 9 ++++++++- doc/src/bond_bpm_spring.rst | 28 ++++++++++++++++++-------- doc/src/comm_modify.rst | 1 + doc/src/neighbor.rst | 21 +++++++++++++++----- src/BPM/bond_bpm.cpp | 16 ++++++++++----- src/BPM/bond_bpm.h | 2 +- src/RIGID/fix_rigid.cpp | 34 +++++++++++++++++--------------- src/fix_update_special_bonds.cpp | 3 +++ src/neighbor.cpp | 10 ++++++++++ 9 files changed, 88 insertions(+), 36 deletions(-) diff --git a/doc/src/bond_bpm_rotational.rst b/doc/src/bond_bpm_rotational.rst index bc2f383828..0975fb35f6 100644 --- a/doc/src/bond_bpm_rotational.rst +++ b/doc/src/bond_bpm_rotational.rst @@ -10,7 +10,7 @@ Syntax 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:: @@ -30,6 +30,9 @@ Syntax *smooth* value = *yes* or *no* smooths bond forces near the breaking point + *break/no* + indicates that bonds should not break during a run + Examples """""""" @@ -138,6 +141,10 @@ the *overlay/pair* keyword. These settings require specific restrictions. Further details can be found in the `:doc: how to ` page on BPMs. +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 break during the simulation. Whenever a bond breaks, data is processed and transferred to an internal fix labeled *fix_ID*. This allows the diff --git a/doc/src/bond_bpm_spring.rst b/doc/src/bond_bpm_spring.rst index de31357afe..fe5c63216d 100644 --- a/doc/src/bond_bpm_spring.rst +++ b/doc/src/bond_bpm_spring.rst @@ -10,7 +10,7 @@ Syntax 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:: @@ -30,6 +30,9 @@ Syntax *smooth* value = *yes* or *no* smooths bond forces near the breaking point + *break/no* + indicates that bonds should not break during a run + Examples """""""" @@ -45,16 +48,17 @@ Examples Description """"""""""" -The *bpm/spring* bond style computes forces and torques based on -deviations from the initial reference state of the two atoms. The -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 -written to :doc:`binary restart files ` such that restarting -the system will not reset the reference state of a bond. +The *bpm/spring* bond style computes forces based on deviations from +the initial reference state of the two atoms. The 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 written to +:doc:`binary restart files ` such that restarting the system +will not reset the reference state of a bond. This bond style only applies central-body forces which conserve the 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) `. The force has a magnitude of .. math:: @@ -103,6 +107,10 @@ the *overlay/pair* keyword. These settings require specific restrictions. Further details can be found in the `:doc: how to ` page on BPMs. +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 break during the simulation. Whenever a bond breaks, data is processed and transferred to an internal fix labeled *fix_ID*. This allows the @@ -198,6 +206,10 @@ The option defaults are *smooth* = *yes* ---------- +.. _fragment-Clemmer: + +**(Clemmer)** Clemmer and Robbins, Phys. Rev. Lett. (2022). + .. _Groot4: **(Groot)** Groot and Warren, J Chem Phys, 107, 4423-35 (1997). diff --git a/doc/src/comm_modify.rst b/doc/src/comm_modify.rst index d0526b792b..14f55378c2 100644 --- a/doc/src/comm_modify.rst +++ b/doc/src/comm_modify.rst @@ -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 be faster. In *multi*, each atom is assigned to a collection which should correspond to a set of atoms with similar interaction cutoffs. +See the :doc:`neighbor ` command for a detailed description of collections. In this case, each atom collection is assigned its own distance cutoff for communication purposes, and fewer atoms will be communicated. in *multi/old*, a similar technique is used but atoms diff --git a/doc/src/neighbor.rst b/doc/src/neighbor.rst index 663170ef47..bf9a13923c 100644 --- a/doc/src/neighbor.rst +++ b/doc/src/neighbor.rst @@ -59,9 +59,16 @@ long cutoff, but other type pairs have a much shorter cutoff. The sized particles, where "size" may mean the physical size of the particle or its cutoff distance for interacting with other particles. Different sets of bins are then used to construct the neighbor lists as as further -described by Shire, Hanley, and Stratford :ref:`(Shire) `. -This imposes some extra setup overhead, but the searches themselves may -be much faster. By default, each atom type defines a separate collection +described by Shire, Hanley, and Stratford :ref:`(Shire) ` +and Monti et al. :ref:`(Monti) `. This imposes some extra +setup overhead, but the searches themselves may be much faster. For +instance in a dense binary system with a ratio of the size of the largest +to smallest collection bin :math:`\lamda`, the computational costs of +building a default neighbor list grows as :math:`\lamda^6` while the costs +for *multi* grows as :math:`\lamda^3`, equivalent to the cost of force +evaluations, as identified in Monti et al. :ref:`(Monti) `. + +By default in *multi*, each atom type defines a separate collection of particles. For systems where two or more atom types have the same size (either physical size or cutoff distance), the definition of collections can be customized, which can result in less overhead and @@ -118,6 +125,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). diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp index 5a20059860..724ec0a7bc 100644 --- a/src/BPM/bond_bpm.cpp +++ b/src/BPM/bond_bpm.cpp @@ -39,6 +39,7 @@ BondBPM::BondBPM(LAMMPS *_lmp) : { overlay_flag = 0; prop_atom_flag = 0; + break_flag = 1; nvalues = 0; r0_max_estimate = 0.0; @@ -103,21 +104,21 @@ void BondBPM::init_style() } } else { // Require atoms know about all of their bonds and if they break - if (force->newton_bond) - error->all(FLERR, "Without overlay/pair, BPM bond styles require Newton bond off"); + if (force->newton_bond && break_flag) + 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 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 if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0) 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 || force->special_coul[3] != 1.0) 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"); fix_update_special_bonds = dynamic_cast(modify->replace_fix( id_fix_dummy, fmt::format("{} all UPDATE_SPECIAL_BONDS", id_fix_update), 1)); @@ -187,6 +188,9 @@ void BondBPM::settings(int narg, char **arg) } else if (strcmp(arg[iarg], "overlay/pair") == 0) { overlay_flag = 1; iarg++; + } else if (strcmp(arg[iarg], "break/no") == 0) { + break_flag = 0; + iarg++; } else { leftover_iarg.push_back(iarg); iarg++; @@ -328,6 +332,8 @@ void BondBPM::read_restart(FILE *fp) 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) { for (int n = 0; n < nvalues; n++) (this->*pack_choice[n])(n, i, j); diff --git a/src/BPM/bond_bpm.h b/src/BPM/bond_bpm.h index 7a32ec2d0d..470307add1 100644 --- a/src/BPM/bond_bpm.h +++ b/src/BPM/bond_bpm.h @@ -52,7 +52,7 @@ class BondBPM : public Bond { FnPtrPack *pack_choice; // ptrs to pack functions 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; void pack_id1(int, int, int); diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index f0bcf630b7..caae812468 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/src/RIGID/fix_rigid.cpp @@ -2523,9 +2523,17 @@ int FixRigid::pack_exchange(int i, double *buf) buf[2] = displace[i][0]; buf[3] = displace[i][1]; 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; + if (vflag_atom) + for (int k = 0; k < 6; k++) + buf[m++] = vatom[i][k]; + + if (!extended) return m; + buf[m++] = eflags[i]; for (int j = 0; j < orientflag; j++) buf[m++] = orient[i][j]; @@ -2535,13 +2543,6 @@ int FixRigid::pack_exchange(int i, double *buf) 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; } @@ -2556,9 +2557,17 @@ int FixRigid::unpack_exchange(int nlocal, double *buf) displace[nlocal][0] = buf[2]; displace[nlocal][1] = buf[3]; 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; + if (vflag_atom) + for (int k = 0; k < 6; k++) + vatom[nlocal][k] = buf[m++]; + + if (!extended) return m; + eflags[nlocal] = static_cast (buf[m++]); for (int j = 0; j < orientflag; j++) orient[nlocal][j] = buf[m++]; @@ -2568,13 +2577,6 @@ int FixRigid::unpack_exchange(int nlocal, double *buf) 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; } diff --git a/src/fix_update_special_bonds.cpp b/src/fix_update_special_bonds.cpp index 581918dbd3..9ffd84c7de 100644 --- a/src/fix_update_special_bonds.cpp +++ b/src/fix_update_special_bonds.cpp @@ -73,6 +73,9 @@ void FixUpdateSpecialBonds::setup(int /*vflag*/) force->special_coul[3] != 1.0) error->all(FLERR, "Fix update/special/bonds requires special Coulomb weights = 1,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"); } /* ---------------------------------------------------------------------- diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 298ec29994..6696971e7d 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -95,6 +95,16 @@ static const char cite_neigh_multi[] = " Detection Applied to Investigate the Quasi-Static Limit},\n" " journal = {Computational Particle Mechanics},\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"; // template for factory functions: From c9af040be96b5b5e607b8b6494ff79c4a0dc5bf1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 18 Mar 2023 20:00:58 -0400 Subject: [PATCH 02/22] improve error message --- src/fix_adapt.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp index 8dd97250a3..c90aa7c5b6 100644 --- a/src/fix_adapt.cpp +++ b/src/fix_adapt.cpp @@ -368,7 +368,7 @@ void FixAdapt::init() void *ptr = ad->pair->extract(ad->pparam,ad->pdim); if (ptr == nullptr) - error->all(FLERR,"Fix adapt pair style param not supported"); + error->all(FLERR,"Fix adapt pair style {} param {} not supported", ad->pstyle, ad->pparam); // for pair styles only parameters that are 2-d arrays in atom types or // scalars are supported From ece7697f6a3da0e8a7593321ba77fe846cbb98bf Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 18 Mar 2023 23:16:59 -0400 Subject: [PATCH 03/22] cosmetic --- src/MANYBODY/pair_eam.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp index ee824d291b..a3d4257cc2 100644 --- a/src/MANYBODY/pair_eam.cpp +++ b/src/MANYBODY/pair_eam.cpp @@ -852,8 +852,7 @@ double PairEAM::single(int i, int j, int itype, int jtype, /* ---------------------------------------------------------------------- */ -int PairEAM::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) +int PairEAM::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) { int i,j,m; From 03b63de58808697f29effb5be0d823972a85c4e7 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 19 Mar 2023 16:05:53 -0400 Subject: [PATCH 04/22] update .gitignore for recently added styles --- src/.gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/.gitignore b/src/.gitignore index ac4a776cfc..7c7155018a 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -56,6 +56,8 @@ /pair_lepton.cpp /pair_lepton.h +/pair_lepton_coul.cpp +/pair_lepton_coul.h /bond_lepton.cpp /bond_lepton.h /angle_lepton.cpp From 57e86346a6936a9c1d23e5f24d1d79b8f12e190a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 19 Mar 2023 18:20:31 -0400 Subject: [PATCH 05/22] fix bug where floating point data was formatted as integer since the second column was assumed to be the type --- tools/python/pizza/dump.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tools/python/pizza/dump.py b/tools/python/pizza/dump.py index 3dabc9b1b4..d57bdf63a3 100644 --- a/tools/python/pizza/dump.py +++ b/tools/python/pizza/dump.py @@ -657,11 +657,14 @@ class dump: atoms = snap.atoms nvalues = len(atoms[0]) + keys = dict() + for pair in self.names.items(): + keys[pair[1]] = pair[0] for i in range(snap.natoms): if not snap.aselect[i]: continue line = "" for j in range(nvalues): - if (j < 2): + if keys[j] == 'id' or keys[j] == 'type' or keys[j] == 'mol': line += str(int(atoms[i][j])) + " " else: line += str(atoms[i][j]) + " " From b0e7d9702bfc83ea35056fdb8278033eea9005c2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 20 Mar 2023 01:02:22 -0400 Subject: [PATCH 06/22] modernize and simplify --- src/KIM/kim_init.cpp | 8 ++----- src/KIM/kim_interactions.cpp | 39 +++++++++++++--------------------- src/KIM/kim_param.cpp | 19 +++++++---------- src/KIM/kim_query.cpp | 41 +++++++++++++++--------------------- 4 files changed, 42 insertions(+), 65 deletions(-) diff --git a/src/KIM/kim_init.cpp b/src/KIM/kim_init.cpp index 037fb2857e..b5b744d71b 100644 --- a/src/KIM/kim_init.cpp +++ b/src/KIM/kim_init.cpp @@ -312,12 +312,8 @@ void KimInit::do_init(char *model_name, char *user_units, char *model_units, KIM { // create storage proxy fix. delete existing fix, if needed. - int ifix = modify->find_fix("KIM_MODEL_STORE"); - if (ifix >= 0) modify->delete_fix(ifix); - modify->add_fix("KIM_MODEL_STORE all STORE/KIM"); - ifix = modify->find_fix("KIM_MODEL_STORE"); - - auto fix_store = dynamic_cast(modify->fix[ifix]); + if (modify->get_fix_by_id("KIM_MODEL_STORE")) modify->delete_fix("KIM_MODEL_STORE"); + auto fix_store = dynamic_cast(modify->add_fix("KIM_MODEL_STORE all STORE/KIM")); fix_store->setptr("model_name", (void *) model_name); fix_store->setptr("user_units", (void *) user_units); fix_store->setptr("model_units", (void *) model_units); diff --git a/src/KIM/kim_interactions.cpp b/src/KIM/kim_interactions.cpp index 5a91109e53..024586d6b3 100644 --- a/src/KIM/kim_interactions.cpp +++ b/src/KIM/kim_interactions.cpp @@ -103,11 +103,9 @@ void KimInteractions::do_setup(int narg, char **arg) if ((narg == 1) && (arg_str == "fixed_types")) { fixed_types = true; } else if (narg != atom->ntypes) { - error->all(FLERR, "Illegal 'kim interactions' command.\nThe " - "LAMMPS simulation has {} atom type(s), but " - "{} chemical species passed to the " - "'kim interactions' command", - atom->ntypes, narg); + error->all(FLERR, "Illegal 'kim interactions' command.\nThe LAMMPS simulation has {} atom " + "type(s), but {} chemical species passed to the 'kim interactions' command", + atom->ntypes, narg); } else { fixed_types = false; } @@ -119,16 +117,14 @@ void KimInteractions::do_setup(int narg, char **arg) // retrieve model name and pointer to simulator model class instance. // validate model name if not given as null pointer. - int ifix = modify->find_fix("KIM_MODEL_STORE"); - if (ifix >= 0) { - auto fix_store = dynamic_cast(modify->fix[ifix]); + auto fix_store = dynamic_cast(modify->get_fix_by_id("KIM_MODEL_STORE")); + if (fix_store) { model_name = (char *)fix_store->getptr("model_name"); simulatorModel = (KIM_SimulatorModel *)fix_store->getptr("simulator_model"); } else error->all(FLERR, "Must use 'kim init' before 'kim interactions'"); // Begin output to log file - input->write_echo("#=== BEGIN kim interactions ===========================" - "=======\n"); + input->write_echo("#=== BEGIN kim interactions ==================================\n"); if (simulatorModel) { auto first_visit = input->variable->find("kim_update"); @@ -165,8 +161,8 @@ void KimInteractions::do_setup(int narg, char **arg) if (atom_type_sym == sim_species) species_is_supported = true; } if (!species_is_supported) { - error->all(FLERR, "Species '{}' is not supported by this " - "KIM Simulator Model", atom_type_sym); + error->all(FLERR, "Species '{}' is not supported by this KIM Simulator Model", + atom_type_sym); } } } else { @@ -179,18 +175,15 @@ void KimInteractions::do_setup(int narg, char **arg) const char *sim_field, *sim_value; KIM_SimulatorModel_GetNumberOfSimulatorFields(simulatorModel, &sim_fields); for (int i = 0; i < sim_fields; ++i) { - KIM_SimulatorModel_GetSimulatorFieldMetadata( - simulatorModel, i, &sim_lines, &sim_field); + KIM_SimulatorModel_GetSimulatorFieldMetadata(simulatorModel, i, &sim_lines, &sim_field); const std::string sim_field_str(sim_field); if (sim_field_str == "units") { - KIM_SimulatorModel_GetSimulatorFieldLine( - simulatorModel, i, 0, &sim_value); + KIM_SimulatorModel_GetSimulatorFieldLine(simulatorModel, i, 0, &sim_value); - const std::string sim_value_str(sim_value); - const std::string unit_style_str(update->unit_style); - if (sim_value_str != unit_style_str) - error->all(FLERR, "Incompatible units for KIM Simulator Model"); + if (strcmp(sim_value, update->unit_style) != 0) + error->all(FLERR, "Incompatible units for KIM Simulator Model: {} vs {}", + sim_value, update->unit_style); } } @@ -217,8 +210,7 @@ void KimInteractions::do_setup(int narg, char **arg) no_model_definition = false; for (int j = 0; j < sim_lines; ++j) { - KIM_SimulatorModel_GetSimulatorFieldLine( - simulatorModel, i, j, &sim_value); + KIM_SimulatorModel_GetSimulatorFieldLine(simulatorModel, i, j, &sim_value); if (utils::strmatch(sim_value, "^KIM_SET_TYPE_PARAMETERS")) { // Notes regarding the KIM_SET_TYPE_PARAMETERS command // * This is an INTERNAL command. @@ -264,8 +256,7 @@ void KimInteractions::do_setup(int narg, char **arg) } // End output to log file - input->write_echo("#=== END kim interactions =============================" - "=======\n\n"); + input->write_echo("#=== END kim interactions ====================================\n\n"); } /* ---------------------------------------------------------------------- */ diff --git a/src/KIM/kim_param.cpp b/src/KIM/kim_param.cpp index 6473cb2ad0..daf761873d 100644 --- a/src/KIM/kim_param.cpp +++ b/src/KIM/kim_param.cpp @@ -159,20 +159,17 @@ void KimParam::command(int narg, char **arg) error->all(FLERR, "Incorrect arguments in 'kim param' command.\n" "'kim param get/set' is mandatory"); - int const ifix = modify->find_fix("KIM_MODEL_STORE"); - if (ifix >= 0) { - auto fix_store = reinterpret_cast(modify->fix[ifix]); - - KIM_SimulatorModel *simulatorModel = - reinterpret_cast( - fix_store->getptr("simulator_model")); + auto fix_store = dynamic_cast(modify->get_fix_by_id("KIM_MODEL_STORE")); + if (fix_store) { + auto *simulatorModel = reinterpret_cast( + fix_store->getptr("simulator_model")); if (simulatorModel) error->all(FLERR, "'kim param' can only be used with a KIM Portable Model"); } - input->write_echo(fmt::format("#=== BEGIN kim param {} ===================" - "==================\n", kim_param_get_set)); + input->write_echo(fmt::format("#=== BEGIN kim param {} =====================================\n", + kim_param_get_set)); KIM_Model *pkim = nullptr; @@ -431,6 +428,6 @@ void KimParam::command(int narg, char **arg) } else error->all(FLERR, "This model has No mutable parameters"); - input->write_echo(fmt::format("#=== END kim param {} =====================" - "==================\n", kim_param_get_set)); + input->write_echo(fmt::format("#=== END kim param {} =======================================\n", + kim_param_get_set)); } diff --git a/src/KIM/kim_query.cpp b/src/KIM/kim_query.cpp index ebb1790e6e..610ec18fcc 100644 --- a/src/KIM/kim_query.cpp +++ b/src/KIM/kim_query.cpp @@ -136,11 +136,9 @@ void KimQuery::command(int narg, char **arg) // check the query_args format (a series of keyword=value pairs) for (int i = 2; i < narg; ++i) { if (!utils::strmatch(arg[i], "[=][\\[].*[\\]]")) - error->all(FLERR, "Illegal query format.\nInput argument " - "of `{}` to 'kim query' is wrong. The " - "query format is the keyword=[value], " - "where value is always an array of one or " - "more comma-separated items", arg[i]); + error->all(FLERR, "Illegal query format.\nInput argument of `{}` to 'kim query' is wrong. " + "The query format is the keyword=[value], where value is always an array of one " + "or more comma-separated items", arg[i]); } if (query_function != "get_available_models") { @@ -156,15 +154,13 @@ void KimQuery::command(int narg, char **arg) // if the model name is not provided by the user if (model_name.empty()) { // check if we had a kim init command by finding fix STORE/KIM - const int ifix = modify->find_fix("KIM_MODEL_STORE"); - if (ifix >= 0) { - auto fix_store = dynamic_cast(modify->fix[ifix]); + auto fix_store = dynamic_cast(modify->get_fix_by_id("KIM_MODEL_STORE")); + if (fix_store) { char *model_name_c = (char *) fix_store->getptr("model_name"); model_name = model_name_c; } else { - error->all(FLERR, "Illegal query format.\nMust use 'kim init' " - "before 'kim query' or must provide the model name " - "after query function with the format of " + error->all(FLERR, "Illegal query format.\nMust use 'kim init' before 'kim query' " + "or must provide the model name after query function with the format of " "'model=[model_name]'"); } } @@ -179,16 +175,14 @@ void KimQuery::command(int narg, char **arg) // and then the error message that was returned by the web server if (strlen(value) == 0) { - auto msg = fmt::format("OpenKIM query failed: {}", value + 1); - delete [] value; - error->all(FLERR, msg); + error->all(FLERR, "OpenKIM query failed: {}", value + 1); + delete[] value; } else if (strcmp(value, "EMPTY") == 0) { - delete [] value; + delete[] value; error->all(FLERR, "OpenKIM query returned no results"); } - input->write_echo("#=== BEGIN kim-query ==================================" - "=======\n"); + input->write_echo("#=== BEGIN kim-query =========================================\n"); // trim list of models to those that are installed on the system if (query_function == "get_available_models") { Tokenizer vals(value, ", \""); @@ -199,7 +193,7 @@ void KimQuery::command(int narg, char **arg) KIM_CollectionItemType typ; if (KIM_Collections_Create(&collections)) { - delete [] value; + delete[] value; error->all(FLERR, "Unable to access KIM Collections to find Model"); } @@ -219,7 +213,7 @@ void KimQuery::command(int narg, char **arg) fmt::format("# Missing OpenKIM models: {}\n\n", missing)); if (available.empty()) { - delete [] value; + delete[] value; error->all(FLERR,"There are no matching OpenKIM models installed on the system"); } @@ -262,13 +256,12 @@ void KimQuery::command(int narg, char **arg) input->variable->set(setcmd); input->write_echo(fmt::format("variable {}\n", setcmd)); } - input->write_echo("#=== END kim-query ====================================" - "=======\n\n"); + input->write_echo("#=== END kim-query ===========================================\n\n"); - delete [] value; + delete[] value; #else - error->all(FLERR, "Cannot use 'kim query' command when KIM package " - "is compiled without support for libcurl"); + error->all(FLERR, "Cannot use 'kim query' command when KIM package is compiled without " + "support for libcurl"); #endif } From f4314076ae20cc2e0d9cbb19e2cea81714c4727a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 20 Mar 2023 06:17:37 -0400 Subject: [PATCH 07/22] address scoping issues, use constexpr for unit constants, modernize --- src/KIM/fix_store_kim.cpp | 43 ++-- src/KIM/kim_command.cpp | 58 ++--- src/KIM/kim_init.cpp | 56 +---- src/KIM/kim_interactions.cpp | 10 +- src/KIM/kim_param.cpp | 57 +---- src/KIM/kim_param.h | 2 - src/KIM/kim_property.cpp | 68 +++--- src/KIM/kim_units.cpp | 425 ++++++++++++++++++++--------------- src/KIM/kim_units.h | 17 +- src/KIM/pair_kim.cpp | 21 +- 10 files changed, 358 insertions(+), 399 deletions(-) diff --git a/src/KIM/fix_store_kim.cpp b/src/KIM/fix_store_kim.cpp index 415ba4c410..44c2ccaa1d 100644 --- a/src/KIM/fix_store_kim.cpp +++ b/src/KIM/fix_store_kim.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -69,11 +68,11 @@ using namespace FixConst; /* ---------------------------------------------------------------------- */ -FixStoreKIM::FixStoreKIM(LAMMPS *lmp, int narg, char **arg) - : Fix(lmp, narg, arg), simulator_model(nullptr), model_name(nullptr), - model_units(nullptr), user_units(nullptr) +FixStoreKIM::FixStoreKIM(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), simulator_model(nullptr), model_name(nullptr), model_units(nullptr), + user_units(nullptr) { - if (narg != 3) error->all(FLERR,"Illegal fix STORE/KIM command"); + if (narg != 3) error->all(FLERR, "Illegal fix STORE/KIM command"); } /* ---------------------------------------------------------------------- */ @@ -83,25 +82,25 @@ FixStoreKIM::~FixStoreKIM() // free associated storage if (simulator_model) { - auto sm = (KIM_SimulatorModel *)simulator_model; + auto sm = (KIM_SimulatorModel *) simulator_model; KIM_SimulatorModel_Destroy(&sm); simulator_model = nullptr; } if (model_name) { - auto mn = (char *)model_name; + auto mn = (char *) model_name; delete[] mn; model_name = nullptr; } if (model_units) { - auto mu = (char *)model_units; + auto mu = (char *) model_units; delete[] mu; model_units = nullptr; } if (user_units) { - auto uu = (char *)user_units; + auto uu = (char *) user_units; delete[] uu; user_units = nullptr; } @@ -121,38 +120,44 @@ void FixStoreKIM::setptr(const std::string &name, void *ptr) { if (name == "simulator_model") { if (simulator_model) { - auto sm = (KIM_SimulatorModel *)simulator_model; + auto sm = (KIM_SimulatorModel *) simulator_model; KIM_SimulatorModel_Destroy(&sm); } simulator_model = ptr; } else if (name == "model_name") { if (model_name) { - auto mn = (char *)model_name; + auto mn = (char *) model_name; delete[] mn; } model_name = ptr; } else if (name == "model_units") { if (model_units) { - auto mu = (char *)model_units; + auto mu = (char *) model_units; delete[] mu; } model_units = ptr; } else if (name == "user_units") { if (user_units) { - auto uu = (char *)user_units; + auto uu = (char *) user_units; delete[] uu; } user_units = ptr; - } else error->all(FLERR,"Unknown property in fix STORE/KIM"); + } else + error->all(FLERR, "Unknown property in fix STORE/KIM"); } /* ---------------------------------------------------------------------- */ void *FixStoreKIM::getptr(const std::string &name) { - if (name == "simulator_model") return simulator_model; - else if (name == "model_name") return model_name; - else if (name == "model_units") return model_units; - else if (name == "user_units") return user_units; - else return nullptr; + if (name == "simulator_model") + return simulator_model; + else if (name == "model_name") + return model_name; + else if (name == "model_units") + return model_units; + else if (name == "user_units") + return user_units; + else + return nullptr; } diff --git a/src/KIM/kim_command.cpp b/src/KIM/kim_command.cpp index 820e4dcba7..57871e964b 100644 --- a/src/KIM/kim_command.cpp +++ b/src/KIM/kim_command.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -70,40 +69,40 @@ using namespace LAMMPS_NS; static constexpr const char *const cite_openkim = - "OpenKIM Project: doi:10.1007/s11837-011-0102-6\n\n" - "@Article{tadmor:elliott:2011,\n" - " author = {E. B. Tadmor and R. S. Elliott and J. P. Sethna and R. E. Miller " - "and C. A. Becker},\n" - " title = {The potential of atomistic simulations and the {K}nowledgebase of " - "{I}nteratomic {M}odels},\n" - " journal = {{JOM}},\n" - " year = 2011,\n" - " volume = 63,\n" - " number = 17,\n" - " pages = {17},\n" - " doi = {10.1007/s11837-011-0102-6}\n" - "}\n\n"; + "OpenKIM Project: doi:10.1007/s11837-011-0102-6\n\n" + "@Article{tadmor:elliott:2011,\n" + " author = {E. B. Tadmor and R. S. Elliott and J. P. Sethna and R. E. Miller " + "and C. A. Becker},\n" + " title = {The potential of atomistic simulations and the {K}nowledgebase of " + "{I}nteratomic {M}odels},\n" + " journal = {{JOM}},\n" + " year = 2011,\n" + " volume = 63,\n" + " number = 17,\n" + " pages = {17},\n" + " doi = {10.1007/s11837-011-0102-6}\n" + "}\n\n"; static constexpr const char *const cite_openkim_query = - "OpenKIM query: doi:10.1063/5.0014267\n\n" - "@Article{karls:bierbaum:2020,\n" - " author = {D. S. Karls and M. Bierbaum and A. A. Alemi and R. S. Elliott " - "and J. P. Sethna and E. B. Tadmor},\n" - " title = {The {O}pen{KIM} processing pipeline: {A} cloud-based automatic " - "material property computation engine},\n" - " journal = {{T}he {J}ournal of {C}hemical {P}hysics},\n" - " year = 2020,\n" - " volume = 153,\n" - " number = 6,\n" - " pages = {064104},\n" - " doi = {10.1063/5.0014267}\n" - "}\n\n"; + "OpenKIM query: doi:10.1063/5.0014267\n\n" + "@Article{karls:bierbaum:2020,\n" + " author = {D. S. Karls and M. Bierbaum and A. A. Alemi and R. S. Elliott " + "and J. P. Sethna and E. B. Tadmor},\n" + " title = {The {O}pen{KIM} processing pipeline: {A} cloud-based automatic " + "material property computation engine},\n" + " journal = {{T}he {J}ournal of {C}hemical {P}hysics},\n" + " year = 2020,\n" + " volume = 153,\n" + " number = 6,\n" + " pages = {064104},\n" + " doi = {10.1063/5.0014267}\n" + "}\n\n"; /* ---------------------------------------------------------------------- */ void KimCommand::command(int narg, char **arg) { - if (narg < 1) error->all(FLERR,"Illegal kim command"); + if (narg < 1) error->all(FLERR, "Illegal kim command"); const std::string subcmd(arg[0]); narg--; @@ -132,5 +131,6 @@ void KimCommand::command(int narg, char **arg) auto cmd = new KimQuery(lmp); cmd->command(narg, arg); delete cmd; - } else error->all(FLERR,"Unknown kim subcommand {}", subcmd); + } else + error->all(FLERR, "Unknown kim subcommand {}", subcmd); } diff --git a/src/KIM/kim_init.cpp b/src/KIM/kim_init.cpp index b5b744d71b..f7dd01be49 100644 --- a/src/KIM/kim_init.cpp +++ b/src/KIM/kim_init.cpp @@ -76,6 +76,8 @@ extern "C" { } using namespace LAMMPS_NS; +using kim_units::get_kim_unit_names; +using kim_units::lammps_unit_conversion; /* ---------------------------------------------------------------------- */ @@ -96,7 +98,7 @@ void KimInit::command(int narg, char **arg) error->all(FLERR, "Illegal 'kim init' command.\n" "The argument followed by unit_style {} is an optional argument and when " - "is used must be unit_conversion_mode", + "used must be unit_conversion_mode", user_units); } } else @@ -109,58 +111,10 @@ void KimInit::command(int narg, char **arg) if (universe->nprocs > 1) MPI_Barrier(universe->uworld); determine_model_type_and_units(model_name, user_units, &model_units, pkim); - write_log_cite(lmp, model_type, model_name); - do_init(model_name, user_units, model_units, pkim); } -/* ---------------------------------------------------------------------- */ - -namespace { -void get_kim_unit_names(char const *const system, KIM_LengthUnit &lengthUnit, - KIM_EnergyUnit &energyUnit, KIM_ChargeUnit &chargeUnit, - KIM_TemperatureUnit &temperatureUnit, KIM_TimeUnit &timeUnit, Error *error) -{ - const std::string system_str(system); - if (system_str == "real") { - lengthUnit = KIM_LENGTH_UNIT_A; - energyUnit = KIM_ENERGY_UNIT_kcal_mol; - chargeUnit = KIM_CHARGE_UNIT_e; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_fs; - } else if (system_str == "metal") { - lengthUnit = KIM_LENGTH_UNIT_A; - energyUnit = KIM_ENERGY_UNIT_eV; - chargeUnit = KIM_CHARGE_UNIT_e; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_ps; - } else if (system_str == "si") { - lengthUnit = KIM_LENGTH_UNIT_m; - energyUnit = KIM_ENERGY_UNIT_J; - chargeUnit = KIM_CHARGE_UNIT_C; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_s; - } else if (system_str == "cgs") { - lengthUnit = KIM_LENGTH_UNIT_cm; - energyUnit = KIM_ENERGY_UNIT_erg; - chargeUnit = KIM_CHARGE_UNIT_statC; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_s; - } else if (system_str == "electron") { - lengthUnit = KIM_LENGTH_UNIT_Bohr; - energyUnit = KIM_ENERGY_UNIT_Hartree; - chargeUnit = KIM_CHARGE_UNIT_e; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_fs; - } else if ((system_str == "lj") || (system_str == "micro") || (system_str == "nano")) { - error->all(FLERR, "LAMMPS unit_style {} not supported by KIM models", system_str); - } else { - error->all(FLERR, "Unknown unit_style"); - } -} -} // namespace - void KimInit::print_dirs(struct KIM_Collections *const collections) const { int kim_error = 0; @@ -467,9 +421,7 @@ void KimInit::do_variables(const std::string &from, const std::string &to) } ier = lammps_unit_conversion(units[i], from, to, conversion_factor); if (ier != 0) - error->all(FLERR, - "Unable to obtain conversion factor: " - "unit = {}; from = {}; to = {}", + error->all(FLERR, "Unable to obtain conversion factor: unit = {}; from = {}; to = {}", units[i], from, to); variable->internal_set(v_unit, conversion_factor); diff --git a/src/KIM/kim_interactions.cpp b/src/KIM/kim_interactions.cpp index 024586d6b3..1f4f84e648 100644 --- a/src/KIM/kim_interactions.cpp +++ b/src/KIM/kim_interactions.cpp @@ -85,11 +85,10 @@ using namespace LAMMPS_NS; void KimInteractions::command(int narg, char **arg) { - if (narg < 1) error->all(FLERR, "Illegal 'kim interactions' command"); + if (narg < 1) utils::missing_cmd_args(FLERR, "kim interactions", error); if (!domain->box_exist) - error->all(FLERR, "Must use 'kim interactions' command after " - "simulation box is defined"); + error->all(FLERR, "Use of 'kim interactions' before simulation box is defined"); do_setup(narg, arg); } @@ -267,8 +266,7 @@ void KimInteractions::KIM_SET_TYPE_PARAMETERS(const std::string &input_line) con const std::string key = words[1]; if (key != "pair" && key != "charge") - error->one(FLERR, "Unrecognized KEY {} for " - "KIM_SET_TYPE_PARAMETERS command", key); + error->one(FLERR, "Unrecognized KEY {} for KIM_SET_TYPE_PARAMETERS command", key); std::string filename = words[2]; std::vector species(words.begin() + 3, words.end()); @@ -278,7 +276,7 @@ void KimInteractions::KIM_SET_TYPE_PARAMETERS(const std::string &input_line) con FILE *fp = nullptr; if (comm->me == 0) { fp = fopen(filename.c_str(), "r"); - if (fp == nullptr) error->one(FLERR, "Parameter file not found"); + if (fp == nullptr) error->one(FLERR, "Parameter file {} not found", filename); } char line[MAXLINE], *ptr; diff --git a/src/KIM/kim_param.cpp b/src/KIM/kim_param.cpp index daf761873d..f72df81989 100644 --- a/src/KIM/kim_param.cpp +++ b/src/KIM/kim_param.cpp @@ -63,6 +63,7 @@ #include "fix_store_kim.h" #include "force.h" #include "input.h" +#include "kim_units.h" #include "modify.h" #include "pair_kim.h" #include "variable.h" @@ -77,61 +78,7 @@ extern "C" } using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -namespace -{ -void get_kim_unit_names( - char const *const system, - KIM_LengthUnit &lengthUnit, - KIM_EnergyUnit &energyUnit, - KIM_ChargeUnit &chargeUnit, - KIM_TemperatureUnit &temperatureUnit, - KIM_TimeUnit &timeUnit, - Error *error) -{ - const std::string system_str(system); - if (system_str == "real") { - lengthUnit = KIM_LENGTH_UNIT_A; - energyUnit = KIM_ENERGY_UNIT_kcal_mol; - chargeUnit = KIM_CHARGE_UNIT_e; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_fs; - } else if (system_str == "metal") { - lengthUnit = KIM_LENGTH_UNIT_A; - energyUnit = KIM_ENERGY_UNIT_eV; - chargeUnit = KIM_CHARGE_UNIT_e; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_ps; - } else if (system_str == "si") { - lengthUnit = KIM_LENGTH_UNIT_m; - energyUnit = KIM_ENERGY_UNIT_J; - chargeUnit = KIM_CHARGE_UNIT_C; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_s; - } else if (system_str == "cgs") { - lengthUnit = KIM_LENGTH_UNIT_cm; - energyUnit = KIM_ENERGY_UNIT_erg; - chargeUnit = KIM_CHARGE_UNIT_statC; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_s; - } else if (system_str == "electron") { - lengthUnit = KIM_LENGTH_UNIT_Bohr; - energyUnit = KIM_ENERGY_UNIT_Hartree; - chargeUnit = KIM_CHARGE_UNIT_e; - temperatureUnit = KIM_TEMPERATURE_UNIT_K; - timeUnit = KIM_TIME_UNIT_fs; - } else if ((system_str == "lj") || - (system_str == "micro") || - (system_str == "nano")) { - error->all(FLERR, "LAMMPS unit_style {} not supported " - "by KIM models", system_str); - } else { - error->all(FLERR, "Unknown unit_style"); - } -} -} // namespace +using kim_units::get_kim_unit_names; /* ---------------------------------------------------------------------- */ diff --git a/src/KIM/kim_param.h b/src/KIM/kim_param.h index d160d2f1b1..daefb13017 100644 --- a/src/KIM/kim_param.h +++ b/src/KIM/kim_param.h @@ -67,7 +67,5 @@ class KimParam : protected Pointers { KimParam(class LAMMPS *lmp); void command(int, char **); }; - } // namespace LAMMPS_NS - #endif diff --git a/src/KIM/kim_property.cpp b/src/KIM/kim_property.cpp index 7d060ad8b1..6ddebeefc2 100644 --- a/src/KIM/kim_property.cpp +++ b/src/KIM/kim_property.cpp @@ -85,14 +85,13 @@ void KimProperty::command(int narg, char **arg) { #if LMP_PYTHON #if PY_MAJOR_VERSION >= 3 - if (narg < 2) error->all(FLERR, "Invalid 'kim property' command"); + if (narg < 2) utils::missing_cmd_args(FLERR, "kim property", error); const std::string subcmd(arg[0]); if ((subcmd != "create") && (subcmd != "destroy") && (subcmd != "modify") && (subcmd != "remove") && (subcmd != "dump")) { - std::string msg("Incorrect arguments in 'kim property' command.\n"); - msg += "'kim property create/destroy/modify/remove/dump' is mandatory"; - error->all(FLERR, msg); + error->all(FLERR, "Incorrect first argument {} to 'kim property' command.\n" + "One of create, destroy, modify, remove, or dump is mandatory", subcmd); } input->write_echo("#=== kim property =====================================" @@ -117,12 +116,11 @@ void KimProperty::command(int narg, char **arg) kim_property = PyImport_Import(obj); if (!kim_property) { PyGILState_Release(gstate); - std::string msg("Unable to import Python kim_property module!"); - msg += "\nkim-property Python package can be installed with pip:\n"; - msg += "'pip install kim-property'\nSee the installation instructions "; - msg += "at\nhttps://github.com/openkim/kim-property#installing-kim-"; - msg += "property\nfor detailed information"; - error->all(FLERR, msg); + error->all(FLERR, "Unable to import Python kim_property module!\n" + "kim-property Python package can be installed with pip:\n" + "'pip install kim-property'\nSee the installation instructions at\n" + "https://github.com/openkim/kim-property#installing-kim-property\n" + "for detailed information"); } // Decrementing of the reference count @@ -147,9 +145,8 @@ void KimProperty::command(int narg, char **arg) PyObject_GetAttrString(kim_property, "kim_property_create"); if (!pFunc) { PyGILState_Release(gstate); - std::string msg("Unable to get an attribute named "); - msg += "'kim_property_create' from a kim_property object"; - error->all(FLERR, msg); + error->all(FLERR, "Unable to get an attribute named 'kim_property_create' from " + "a kim_property object"); } // Decrementing of the reference count @@ -182,16 +179,13 @@ void KimProperty::command(int narg, char **arg) if (!pValue) { PyErr_Print(); PyGILState_Release(gstate); - std::string msg("Python 'kim_property_create' function "); - msg += "evaluation failed"; - error->one(FLERR, msg); + error->all(FLERR, "Python 'kim_property_create' function evaluation failed"); } // Python function returned a string value const char *pystr = PyUnicode_AsUTF8(pValue); if (kim_str) input->variable->set_string("kim_property_str", pystr); - else - input->variable->set(fmt::format("kim_property_str string '{}'", pystr)); + else input->variable->set(fmt::format("kim_property_str string '{}'", pystr)); Py_XDECREF(pArgs); Py_XDECREF(pFunc); @@ -216,9 +210,8 @@ void KimProperty::command(int narg, char **arg) PyObject_GetAttrString(kim_property, "kim_property_destroy"); if (!pFunc) { PyGILState_Release(gstate); - std::string msg("Unable to get an attribute named "); - msg += "'kim_property_destroy' from a kim_property object"; - error->all(FLERR, msg); + error->all(FLERR, "Unable to get an attribute named 'kim_property_destroy' " + "from a kim_property object"); } // Decrementing of the reference count @@ -244,9 +237,7 @@ void KimProperty::command(int narg, char **arg) if (!pValue) { PyErr_Print(); PyGILState_Release(gstate); - std::string msg("Python 'kim_property_destroy' function "); - msg += "evaluation failed"; - error->one(FLERR, msg); + error->all(FLERR, "Python 'kim_property_destroy' function evaluation failed"); } // Python function returned a string value @@ -276,9 +267,8 @@ void KimProperty::command(int narg, char **arg) PyObject_GetAttrString(kim_property, "kim_property_modify"); if (!pFunc) { PyGILState_Release(gstate); - std::string msg("Unable to get an attribute named "); - msg += "'kim_property_modify' from a kim_property object"; - error->all(FLERR, msg); + error->all(FLERR, "Unable to get an attribute named 'kim_property_modify' " + "from a kim_property object"); } // Decrementing of the reference count @@ -309,9 +299,7 @@ void KimProperty::command(int narg, char **arg) if (!pValue) { PyErr_Print(); PyGILState_Release(gstate); - std::string msg("Python 'kim_property_modify' function "); - msg += "evaluation failed"; - error->one(FLERR, msg); + error->all(FLERR, "Python 'kim_property_modify' function evaluation failed"); } // Python function returned a string value @@ -341,9 +329,8 @@ void KimProperty::command(int narg, char **arg) PyObject_GetAttrString(kim_property, "kim_property_remove"); if (!pFunc) { PyGILState_Release(gstate); - std::string msg("Unable to get an attribute named "); - msg += "'kim_property_remove' from a kim_property object"; - error->all(FLERR, msg); + error->all(FLERR, "Unable to get an attribute named " + "'kim_property_remove' from a kim_property object"); } // Decrementing of the reference count @@ -374,9 +361,7 @@ void KimProperty::command(int narg, char **arg) if (!pValue) { PyErr_Print(); PyGILState_Release(gstate); - std::string msg("Python 'kim_property_remove' function "); - msg += "evaluation failed"; - error->one(FLERR, msg); + error->all(FLERR, "Python 'kim_property_remove' function evaluation failed"); } // Python function returned a string value @@ -404,9 +389,8 @@ void KimProperty::command(int narg, char **arg) PyObject_GetAttrString(kim_property, "kim_property_dump"); if (!pFunc) { PyGILState_Release(gstate); - std::string msg("Unable to get an attribute named "); - msg += "'kim_property_dump' from a kim_property object"; - error->all(FLERR, msg); + error->all(FLERR, "Unable to get an attribute named " + "'kim_property_dump' from a kim_property object"); } // Decrementing of the reference count @@ -428,14 +412,12 @@ void KimProperty::command(int narg, char **arg) if (comm->me == 0) { // call the Python kim_property_dump function - // error check with one() since only some procs may fail + // error check with one() since only root process calls. pValue = PyObject_CallObject(pFunc, pArgs); if (!pValue) { PyErr_Print(); PyGILState_Release(gstate); - std::string msg("Python 'kim_property_dump' function "); - msg += "evaluation failed"; - error->one(FLERR, msg); + error->one(FLERR, "Python 'kim_property_dump' function evaluation failed"); } } else pValue = nullptr; diff --git a/src/KIM/kim_units.cpp b/src/KIM/kim_units.cpp index 7c57e3b71a..7687481cb9 100644 --- a/src/KIM/kim_units.cpp +++ b/src/KIM/kim_units.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -54,179 +53,182 @@ Designed for use with the kim-api-2.0.2 (and newer) package ------------------------------------------------------------------------- */ +#include "kim_units.h" + +#include "error.h" + #include #include #include #include using namespace std; - -namespace -{ +using namespace LAMMPS_NS; +using namespace kim_units; // Constants of nature and basic conversion factors // Source: https://physics.nist.gov/cuu/Constants/Table/allascii.txt // Working with NIST values even when there are newer values for // compatibility with LAMMPS +// clang-format off + /*---------------------- Fundamental constants ------------------------ */ -double const boltz_si = 1.38064852e-23; // [J K^-1] Boltzmann's factor - // (NIST value) -double const Nav = 6.022140857e23; // [unitless] Avogadro's number - // (NIST value) -// double const Nav = 6.02214076e23; // [unitless] Avogadro's number - // (official value May 2019) -double const me_si = 9.10938356e-31; // [kg] electron rest mass - // (NIST value) -// double me_si = 9.10938291e-31; // [kg] electron rest mass -double const e_si = 1.6021766208e-19; // [C] elementary charge - // (charge of an electron/proton) - // (NIST value) +static double constexpr boltz_si = 1.38064852e-23; // [J K^-1] Boltzmann's factor + // (NIST value) +static double constexpr Nav = 6.022140857e23; // [unitless] Avogadro's number (NIST value) +// static double constexp Nav = 6.02214076e23; // [unitless] Avogadro's number + // (official value May 2019) +static double constexpr me_si = 9.10938356e-31; // [kg] electron rest mass (NIST value) +// static double constexpr me_si = 9.10938291e-31; // [kg] electron rest mass +static double constexpr e_si = 1.6021766208e-19; // [C] elementary charge + // (charge of an electron/proton) + // (NIST value) /*---------------------- Distance units ------------------------ */ -double const bohr_si = 5.2917721067e-11; // [m] Bohr unit (distance between - // nucleus and electron in H) - // (NIST value) -double const angstrom_si = 1e-10; // [m] Angstrom -double const centimeter_si = 1e-2; // [m] centimeter -double const micrometer_si = 1e-6; // [m] micrometer (micron) -double const nanometer_si = 1e-9; // [m] nanometer +static double constexpr bohr_si = 5.2917721067e-11; // [m] Bohr unit (distance between + // nucleus and electron in H) (NIST value) +static double constexpr angstrom_si = 1e-10; // [m] Angstrom +static double constexpr centimeter_si = 1e-2; // [m] centimeter +static double constexpr micrometer_si = 1e-6; // [m] micrometer (micron) +static double constexpr nanometer_si = 1e-9; // [m] nanometer /*---------------------- Mass units ------------------------ */ -double const gram_per_mole_si = 1e-3/Nav; // [kg] gram per mole -double const amu_si = 1e-3/Nav; // [kg] atomic mass unit (molecular - // weight) For example, the mean - // molecular weight of water - // is 18.015 atomic mass units - // (amu), so one mole of water - // weight 18.015 grams. -double const gram_si = 1e-3; // [kg] gram -double const picogram_si = 1e-15; // [kg] picogram -double const attogram_si = 1e-21; // [kg[ attogram +static double constexpr gram_per_mole_si = 1e-3/Nav; // [kg] gram per mole +static double constexpr amu_si = 1e-3/Nav; // [kg] atomic mass unit (molecular + // weight) For example, the mean + // molecular weight of water + // is 18.015 atomic mass units + // (amu), so one mole of water + // weight 18.015 grams. +static double constexpr gram_si = 1e-3; // [kg] gram +static double constexpr picogram_si = 1e-15; // [kg] picogram +static double constexpr attogram_si = 1e-21; // [kg[ attogram /*---------------------- Time units ------------------------ */ -double const atu_si = 2.418884326509e-17; // [s] atomic time unit - // ( = hbar/E_h where E_h is the - // Hartree energy) (NIST value) -double const atu_electron_si = atu_si*sqrt(amu_si/me_si); - // [s] atomic time unit - // used in electron system (see https://sourceforge.net/p/lammps/mailman/lammps-users/thread/BCA2BDB2-BA03-4280-896F-1E6120EF47B2%40caltech.edu/) -double const microsecond_si = 1e-6; // [s] microsecond -double const nanosecond_si = 1e-9; // [s] nanosecond -double const picosecond_si = 1e-12; // [s] picosecond -double const femtosecond_si = 1e-15; // [s] femtosecond +static double constexpr atu_si = 2.418884326509e-17; // [s] atomic time unit + // ( = hbar/E_h where E_h is the + // Hartree energy) (NIST value) +static double constexpr atu_electron_si = 5153034.567408186; // must not use sqrt() in constexpr +// static double constexpr atu_electron_si = atu_si*sqrt(amu_si/me_si); + // [s] atomic time unit + // used in electron system (see https://sourceforge.net/p/lammps/mailman/lammps-users/thread/BCA2BDB2-BA03-4280-896F-1E6120EF47B2%40caltech.edu/) +static double constexpr microsecond_si = 1e-6; // [s] microsecond +static double constexpr nanosecond_si = 1e-9; // [s] nanosecond +static double constexpr picosecond_si = 1e-12; // [s] picosecond +static double constexpr femtosecond_si = 1e-15; // [s] femtosecond /*---------------------- Density units ------------------------ */ -double const gram_per_centimetercu_si = - gram_si/pow(centimeter_si,3); // [kg/m^3] gram/centimeter^3 -double const amu_per_bohrcu_si = amu_si/pow(bohr_si,3); // [kg/m^3] amu/bohr^3 -double const picogram_per_micrometercu_si = - picogram_si/pow(micrometer_si,3); // [kg/m^3] picogram/micrometer^3 -double const attogram_per_nanometercu_si = - attogram_si/pow(nanometer_si,3); // [kg/m^3] attogram/nanometer^3 +static double constexpr gram_per_centimetercu_si = + gram_si/centimeter_si/centimeter_si/centimeter_si; // [kg/m^3] gram/centimeter^3 +static double constexpr amu_per_bohrcu_si = amu_si/bohr_si/bohr_si/bohr_si; // [kg/m^3] amu/bohr^3 +static double constexpr picogram_per_micrometercu_si = + picogram_si/micrometer_si/micrometer_si/micrometer_si; // [kg/m^3] picogram/micrometer^3 +static double constexpr attogram_per_nanometercu_si = + attogram_si/nanometer_si/nanometer_si/nanometer_si; // [kg/m^3] attogram/nanometer^3 /*---------------------- Energy/torque units ------------------------ */ -double const kcal_si = 4184.0; // [J] kilocalorie (heat energy - // involved in warming up one - // kilogram of water by one - // degree Kelvin) -double const ev_si = 1.6021766208e-19; // [J] electron volt (amount of - // energy gained or lost by the - // charge of a single electron - // moving across an electric - // potential difference of one - // volt.) (NIST value) -double const hartree_si = 4.359744650e-18; // [J] Hartree (approximately the - // electric potential energy of - // the hydrogen atom in its - // ground state) (NIST value) -double const kcal_per_mole_si = kcal_si/Nav;// [J] kcal/mole -double const erg_si = 1e-7; // [J] erg -double const dyne_centimeter_si = 1e-7; // [J[ dyne*centimeter -double const picogram_micrometersq_per_microsecondsq_si = - picogram_si*pow(micrometer_si,2)/pow(microsecond_si,2); - // [J] picogram*micrometer^2/ - // microsecond^2 -double const attogram_nanometersq_per_nanosecondsq_si = - attogram_si*pow(nanometer_si,2)/pow(nanosecond_si,2); - // [J] attogram*nanometer^2/ - // nanosecond^2 +static double constexpr kcal_si = 4184.0; // [J] kilocalorie (heat energy + // involved in warming up one + // kilogram of water by one + // degree Kelvin) +static double constexpr ev_si = 1.6021766208e-19; // [J] electron volt (amount of + // energy gained or lost by the + // charge of a single electron + // moving across an electric + // potential difference of one + // volt.) (NIST value) +static double constexpr hartree_si = 4.359744650e-18; // [J] Hartree (approximately the + // electric potential energy of + // the hydrogen atom in its + // ground state) (NIST value) +static double constexpr kcal_per_mole_si = kcal_si/Nav;// [J] kcal/mole +static double constexpr erg_si = 1e-7; // [J] erg +static double constexpr dyne_centimeter_si = 1e-7; // [J[ dyne*centimeter +static double constexpr picogram_micrometersq_per_microsecondsq_si = + picogram_si*micrometer_si/microsecond_si*micrometer_si/microsecond_si; + // [J] picogram*micrometer^2/ + // microsecond^2 +static double constexpr attogram_nanometersq_per_nanosecondsq_si = + attogram_si*nanometer_si/nanosecond_si*nanometer_si/nanosecond_si; + // [J] attogram*nanometer^2/ + // nanosecond^2 /*---------------------- Velocity units ------------------------ */ -double const angstrom_per_femtosecond_si = +static double constexpr angstrom_per_femtosecond_si = angstrom_si/femtosecond_si; // [m/s] Angstrom/femtosecond -double const angstrom_per_picosecond_si = +static double constexpr angstrom_per_picosecond_si = angstrom_si/picosecond_si; // [m/s] Angstrom/picosecond -double const micrometer_per_microsecond_si = +static double constexpr micrometer_per_microsecond_si = micrometer_si/microsecond_si; // [m/s] micrometer/microsecond -double const nanometer_per_nanosecond_si = +static double constexpr nanometer_per_nanosecond_si = nanometer_si/nanosecond_si; // [m/s] nanometer/nanosecond -double const centimeter_per_second_si = +static double constexpr centimeter_per_second_si = centimeter_si; // [m/s] centimeter/second -double const bohr_per_atu_electron_si = +static double constexpr bohr_per_atu_electron_si = bohr_si/atu_electron_si; // [m/s] bohr/atu /*---------------------- Force units ------------------------ */ -double const kcal_per_mole_angstrom_si = +static double constexpr kcal_per_mole_angstrom_si = kcal_per_mole_si/angstrom_si; // [N] kcal/(mole*Angstrom) -double const ev_per_angstrom_si = +static double constexpr ev_per_angstrom_si = ev_si/angstrom_si; // [N] eV/Angstrom -double const dyne_si = +static double constexpr dyne_si = dyne_centimeter_si/centimeter_si; // [N] dyne -double const hartree_per_bohr_si = +static double constexpr hartree_per_bohr_si = hartree_si/bohr_si; // [N] hartree/bohr -double const picogram_micrometer_per_microsecondsq_si = - picogram_si*micrometer_si/pow(microsecond_si,2); +static double constexpr picogram_micrometer_per_microsecondsq_si = + picogram_si*micrometer_si/microsecond_si/microsecond_si; // [N] picogram*micrometer/ // microsecond^2 -double const attogram_nanometer_per_nanosecondsq_si = - attogram_si*nanometer_si/pow(nanosecond_si,2); +static double constexpr attogram_nanometer_per_nanosecondsq_si = + attogram_si*nanometer_si/nanosecond_si/nanosecond_si; // [N] attogram*nanometer/ // nanosecond^2 /*---------------------- Pressure units ------------------------ */ -double const atmosphere_si = 101325.0; // [Pa] standard atmosphere (NIST value) -double const bar_si = 1e5; // [Pa] bar -double const dyne_per_centimetersq_si = - dyne_centimeter_si/pow(centimeter_si,3); +static double constexpr atmosphere_si = 101325.0; // [Pa] standard atmosphere (NIST value) +static double constexpr bar_si = 1e5; // [Pa] bar +static double constexpr dyne_per_centimetersq_si = + dyne_centimeter_si/centimeter_si/centimeter_si/centimeter_si; // [Pa] dyne/centimeter^2 -double const picogram_per_micrometer_microsecondsq_si = - picogram_si/(micrometer_si*pow(microsecond_si,2)); +static double constexpr picogram_per_micrometer_microsecondsq_si = + picogram_si/(micrometer_si*microsecond_si*microsecond_si); // [Pa] picogram/(micrometer* // microsecond^2) -double const attogram_per_nanometer_nanosecondsq_si = - attogram_si/(nanometer_si*pow(nanosecond_si,2)); +static double constexpr attogram_per_nanometer_nanosecondsq_si = + attogram_si/(nanometer_si*nanosecond_si*nanosecond_si); // [Pa] attogram/(nanometer*nanosecond^2) /*---------------------- Viscosity units ------------------------ */ -double const poise_si = 0.1; // [Pa*s] Poise -double const amu_per_bohr_femtosecond_si = +static double constexpr poise_si = 0.1; // [Pa*s] Poise +static double constexpr amu_per_bohr_femtosecond_si = amu_si/(bohr_si*femtosecond_si); // [Pa*s] amu/(bohr*femtosecond) -double const picogram_per_micrometer_microsecond_si = +static double constexpr picogram_per_micrometer_microsecond_si = picogram_si/(micrometer_si*microsecond_si); // [Pa*s] picogram/(micrometer* // microsecond) -double const attogram_per_nanometer_nanosecond_si = +static double constexpr attogram_per_nanometer_nanosecond_si = attogram_si/(nanometer_si*nanosecond_si); // [Pa*s] attogram/(nanometer* // nanosecond) @@ -234,39 +236,41 @@ double const attogram_per_nanometer_nanosecond_si = /*---------------------- Charge units ------------------------ */ -double const echarge_si = e_si; // [C] electron charge unit -double const statcoulomb_si = e_si/4.8032044e-10; // [C] Statcoulomb or esu +static double constexpr echarge_si = e_si; // [C] electron charge unit +static double constexpr statcoulomb_si = e_si/4.8032044e-10; // [C] Statcoulomb or esu // (value from LAMMPS units // documentation) -double const picocoulomb_si = 1e-12; // [C] picocoulomb +static double constexpr picocoulomb_si = 1e-12; // [C] picocoulomb /*---------------------- Dipole units ------------------------ */ -double const electron_angstrom_si = echarge_si*angstrom_si; +static double constexpr electron_angstrom_si = echarge_si*angstrom_si; // [C*m] electron*angstrom -double const statcoulomb_centimeter_si = statcoulomb_si*centimeter_si; +static double constexpr statcoulomb_centimeter_si = statcoulomb_si*centimeter_si; // [C*m] statcoulomb*centimeter -double const debye_si = 1e-18*statcoulomb_centimeter_si; +static double constexpr debye_si = 1e-18*statcoulomb_centimeter_si; // [C*m] Debye -double const picocoulomb_micrometer_si = picocoulomb_si*micrometer_si; +static double constexpr picocoulomb_micrometer_si = picocoulomb_si*micrometer_si; // [C*m] picocoulomb*micrometer -double const electron_nanometer_si = echarge_si*nanometer_si; +static double constexpr electron_nanometer_si = echarge_si*nanometer_si; // [C*m] electron*nanometer /*---------------------- Electric field units ------------------------ */ -double const volt_per_angstrom_si = 1.0/angstrom_si;// [V/m] volt/angstrom -double const statvolt_per_centimeter_si = +static double constexpr volt_per_angstrom_si = 1.0/angstrom_si;// [V/m] volt/angstrom +static double constexpr statvolt_per_centimeter_si = erg_si/(statcoulomb_si*centimeter_si); // [V/m] statvolt/centimeter -double const volt_per_centimeter_si = +static double constexpr volt_per_centimeter_si = 1.0/centimeter_si; // [V/m] volt/centimeter -double const volt_per_micrometer_si = +static double constexpr volt_per_micrometer_si = 1.0/micrometer_si; // [V/m] volt/micrometer -double const volt_per_nanometer_si = +static double constexpr volt_per_nanometer_si = 1.0/nanometer_si; // [V/m] volt/nanometer +namespace LAMMPS_NS { +namespace kim_units { // Define enumerations enum sys_type { @@ -388,20 +392,42 @@ enum units attogram_per_nanometercu = 1405 }; +void initialize_dictionaries(); +units get_lammps_system_unit(sys_type system_enum, unit_type unit_type_enum); +double get_mass_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_distance_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_time_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_energy_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_velocity_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_force_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_torque_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_temperature_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_pressure_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_viscosity_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_charge_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_dipole_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_efield_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_density_conversion_factor(units from_unit_enum, units to_unit_enum); +double get_unit_conversion_factor(unit_type &unit_type_enum, sys_type from_system_enum, + sys_type to_system_enum); + +} +} + // Define dictionaries -map system_dic; -map unit_dic; -map units_real_dic; -map units_metal_dic; -map units_si_dic; -map units_cgs_dic; -map units_electron_dic; -map units_micro_dic; -map units_nano_dic; +static map system_dic; +static map unit_dic; +static map units_real_dic; +static map units_metal_dic; +static map units_si_dic; +static map units_cgs_dic; +static map units_electron_dic; +static map units_micro_dic; +static map units_nano_dic; /* ---------------------------------------------------------------------- */ -void initialize_dictionaries() +void kim_units::initialize_dictionaries() { system_dic["real"] = real; system_dic["metal"] = metal; @@ -537,7 +563,7 @@ void initialize_dictionaries() // Get the enumeration for the unit of type `unit_type_enum` // for LAMMPS system `system_enum`. -units get_lammps_system_unit(sys_type system_enum, unit_type unit_type_enum) +units kim_units::get_lammps_system_unit(sys_type system_enum, unit_type unit_type_enum) { switch(system_enum) { case real : @@ -561,7 +587,7 @@ units get_lammps_system_unit(sys_type system_enum, unit_type unit_type_enum) /* ---------------------------------------------------------------------- */ // Mass conversion -double get_mass_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_mass_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -619,7 +645,7 @@ double get_mass_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Distance conversion -double get_distance_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_distance_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -677,7 +703,7 @@ double get_distance_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Time conversion -double get_time_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_time_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -722,7 +748,7 @@ double get_time_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Energy conversion -double get_energy_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_energy_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; @@ -796,7 +822,7 @@ double get_energy_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Velocity conversion -double get_velocity_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_velocity_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -869,7 +895,7 @@ double get_velocity_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Force conversion -double get_force_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_force_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -942,7 +968,7 @@ double get_force_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Torque conversion -double get_torque_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_torque_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -1015,7 +1041,7 @@ double get_torque_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Temperature conversion -double get_temperature_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_temperature_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; @@ -1026,7 +1052,7 @@ double get_temperature_conversion_factor(units from_unit_enum, units to_unit_enu /* ---------------------------------------------------------------------- */ // Pressure conversion -double get_pressure_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_pressure_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -1084,7 +1110,7 @@ double get_pressure_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Viscosity conversion -double get_viscosity_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_viscosity_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -1129,7 +1155,7 @@ double get_viscosity_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Charge conversion -double get_charge_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_charge_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -1163,7 +1189,7 @@ double get_charge_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Dipole conversion -double get_dipole_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_dipole_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -1221,7 +1247,7 @@ double get_dipole_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ // Electric field conversion -double get_efield_conversion_factor(units from_unit_enum, units to_unit_enum) +double kim_units::get_efield_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -1278,8 +1304,8 @@ double get_efield_conversion_factor(units from_unit_enum, units to_unit_enum) /* ---------------------------------------------------------------------- */ -// Demsity conversion -double get_density_conversion_factor(units from_unit_enum, units to_unit_enum) +// Density conversion +double kim_units::get_density_conversion_factor(units from_unit_enum, units to_unit_enum) { map > conv; double to_si; @@ -1325,9 +1351,9 @@ double get_density_conversion_factor(units from_unit_enum, units to_unit_enum) // This routine returns the unit conversion factor between the // `from_system_enum` to the `to_system_enum` for the `unit_type_enum`. -double get_unit_conversion_factor(unit_type &unit_type_enum, - sys_type from_system_enum, - sys_type to_system_enum) +double kim_units::get_unit_conversion_factor(unit_type &unit_type_enum, + sys_type from_system_enum, + sys_type to_system_enum) { units from_unit = get_lammps_system_unit(from_system_enum, unit_type_enum); units to_unit = get_lammps_system_unit(to_system_enum, unit_type_enum); @@ -1364,48 +1390,87 @@ double get_unit_conversion_factor(unit_type &unit_type_enum, } } -} // end of anonymous name space - /* ---------------------------------------------------------------------- */ - +// clang-format on // Wrapper to the routine that gets the unit conversion. Translates strings // to enumerations and then call get_unit_conversion_factor() -int lammps_unit_conversion(const string &unit_type_str, - const string &from_system_str, - const string &to_system_str, - double &conversion_factor) +int kim_units::lammps_unit_conversion(const string &unit_type_str, const string &from_system_str, + const string &to_system_str, double &conversion_factor) { - // initialize - conversion_factor = 0.0; - initialize_dictionaries(); + // initialize + conversion_factor = 0.0; + initialize_dictionaries(); - // convert input to enumeration - unit_type unit_type_enum; - { - map::const_iterator itr = unit_dic.find(unit_type_str); - if (itr != unit_dic.end()) unit_type_enum = itr->second; - else return 1; // error - } - sys_type from_system_enum; - { - map::const_iterator - itr = system_dic.find(from_system_str); - if (itr != system_dic.end()) from_system_enum = itr->second; - else return 1; // error - } - sys_type to_system_enum; - { - map::const_iterator - itr = system_dic.find(to_system_str); - if (itr != system_dic.end()) to_system_enum = itr->second; - else return 1; - } + // convert input to enumeration + unit_type unit_type_enum; + { + map::const_iterator itr = unit_dic.find(unit_type_str); + if (itr != unit_dic.end()) + unit_type_enum = itr->second; + else + return 1; // error + } + sys_type from_system_enum; + { + map::const_iterator itr = system_dic.find(from_system_str); + if (itr != system_dic.end()) + from_system_enum = itr->second; + else + return 1; // error + } + sys_type to_system_enum; + { + map::const_iterator itr = system_dic.find(to_system_str); + if (itr != system_dic.end()) + to_system_enum = itr->second; + else + return 1; + } - // process unit conversions - conversion_factor = get_unit_conversion_factor(unit_type_enum, - from_system_enum, - to_system_enum); - return 0; + // process unit conversions + conversion_factor = get_unit_conversion_factor(unit_type_enum, from_system_enum, to_system_enum); + return 0; } - +void kim_units::get_kim_unit_names(char const *const system, KIM_LengthUnit &lengthUnit, + KIM_EnergyUnit &energyUnit, KIM_ChargeUnit &chargeUnit, + KIM_TemperatureUnit &temperatureUnit, KIM_TimeUnit &timeUnit, + Error *error) +{ + const std::string system_str(system); + if (system_str == "real") { + lengthUnit = KIM_LENGTH_UNIT_A; + energyUnit = KIM_ENERGY_UNIT_kcal_mol; + chargeUnit = KIM_CHARGE_UNIT_e; + temperatureUnit = KIM_TEMPERATURE_UNIT_K; + timeUnit = KIM_TIME_UNIT_fs; + } else if (system_str == "metal") { + lengthUnit = KIM_LENGTH_UNIT_A; + energyUnit = KIM_ENERGY_UNIT_eV; + chargeUnit = KIM_CHARGE_UNIT_e; + temperatureUnit = KIM_TEMPERATURE_UNIT_K; + timeUnit = KIM_TIME_UNIT_ps; + } else if (system_str == "si") { + lengthUnit = KIM_LENGTH_UNIT_m; + energyUnit = KIM_ENERGY_UNIT_J; + chargeUnit = KIM_CHARGE_UNIT_C; + temperatureUnit = KIM_TEMPERATURE_UNIT_K; + timeUnit = KIM_TIME_UNIT_s; + } else if (system_str == "cgs") { + lengthUnit = KIM_LENGTH_UNIT_cm; + energyUnit = KIM_ENERGY_UNIT_erg; + chargeUnit = KIM_CHARGE_UNIT_statC; + temperatureUnit = KIM_TEMPERATURE_UNIT_K; + timeUnit = KIM_TIME_UNIT_s; + } else if (system_str == "electron") { + lengthUnit = KIM_LENGTH_UNIT_Bohr; + energyUnit = KIM_ENERGY_UNIT_Hartree; + chargeUnit = KIM_CHARGE_UNIT_e; + temperatureUnit = KIM_TEMPERATURE_UNIT_K; + timeUnit = KIM_TIME_UNIT_fs; + } else if ((system_str == "lj") || (system_str == "micro") || (system_str == "nano")) { + error->all(FLERR, "LAMMPS unit_style {} not supported by KIM models", system_str); + } else { + error->all(FLERR, "Unknown unit_style {}", system_str); + } +} diff --git a/src/KIM/kim_units.h b/src/KIM/kim_units.h index e85760ab82..dfb30c0a4c 100644 --- a/src/KIM/kim_units.h +++ b/src/KIM/kim_units.h @@ -56,7 +56,20 @@ #define LMP_KIM_UNITS_H #include +extern "C" { +#include "KIM_SimulatorHeaders.h" +} -int lammps_unit_conversion(const std::string &unit_type_str, const std::string &from_system_str, - const std::string &to_system_str, double &conversion_factor); +namespace LAMMPS_NS { +class Error; + +namespace kim_units { + int lammps_unit_conversion(const std::string &unit_type_str, const std::string &from_system_str, + const std::string &to_system_str, double &conversion_factor); + void get_kim_unit_names(char const *const system, KIM_LengthUnit &lengthUnit, + KIM_EnergyUnit &energyUnit, KIM_ChargeUnit &chargeUnit, + KIM_TemperatureUnit &temperatureUnit, KIM_TimeUnit &timeUnit, + Error *error); +} // namespace kim_units +} // namespace LAMMPS_NS #endif diff --git a/src/KIM/pair_kim.cpp b/src/KIM/pair_kim.cpp index 6a5a792797..14564162e2 100644 --- a/src/KIM/pair_kim.cpp +++ b/src/KIM/pair_kim.cpp @@ -218,8 +218,7 @@ void PairKIM::compute(int eflag, int vflag) kim_particleSpecies); memory->create(kim_particleContributing,lmps_maxalloc, "pair:kim_particleContributing"); - kimerror = kimerror || KIM_ComputeArguments_SetArgumentPointerInteger( - pargs, + kimerror = kimerror || KIM_ComputeArguments_SetArgumentPointerInteger(pargs, KIM_COMPUTE_ARGUMENT_NAME_particleContributing, kim_particleContributing); if (kimerror) @@ -248,7 +247,7 @@ void PairKIM::compute(int eflag, int vflag) // compute via KIM model int kimerror = KIM_Model_Compute(pkim, pargs); - if (kimerror) error->all(FLERR,"KIM Compute returned error"); + if (kimerror) error->all(FLERR,"KIM Compute returned error {}", kimerror); // compute virial before reverse comm! if (vflag_global) @@ -454,7 +453,7 @@ void PairKIM::coeff(int narg, char **arg) kimerror = KIM_Model_GetParameterMetadata(pkim, param_index, &kim_DataType, &extent, &str_name, &str_desc); if (kimerror) - error->all(FLERR,"KIM GetParameterMetadata returned error"); + error->all(FLERR,"KIM GetParameterMetadata returned error {}", kimerror); const std::string str_name_str(str_name); if (paramname == str_name_str) break; @@ -484,8 +483,8 @@ void PairKIM::coeff(int narg, char **arg) nubound = atoi(words[1].c_str()); if ((nubound < 1) || (nubound > extent) || (nlbound < 1) || (nlbound > nubound)) - error->all(FLERR,"Illegal index_range '{}-{}' for '{}' parameter with the extent of '{}'", - nlbound, nubound, paramname, extent); + error->all(FLERR,"Illegal index_range '{}-{}' for '{}' parameter with the extent " + "of '{}'", nlbound, nubound, paramname, extent); } else { nlbound = atoi(argtostr.c_str()); @@ -975,7 +974,7 @@ void PairKIM::set_lmps_flags() // determine if running with pair hybrid if (force->pair_match("hybrid",0)) - error->all(FLERR,"pair_kim does not support hybrid"); + error->all(FLERR,"Pair style must not be used as a hybrid sub-style"); const std::string unit_style_str(update->unit_style); @@ -1018,10 +1017,9 @@ void PairKIM::set_lmps_flags() } else if ((unit_style_str == "lj") || (unit_style_str == "micro") || (unit_style_str == "nano")) { - error->all(FLERR,"LAMMPS unit_style {} not supported " - "by KIM models", unit_style_str); + error->all(FLERR,"LAMMPS unit_style {} not supported by KIM models", unit_style_str); } else { - error->all(FLERR,"Unknown unit_style"); + error->all(FLERR,"Unknown unit_style {}", unit_style_str); } } @@ -1100,7 +1098,8 @@ void PairKIM::set_kim_model_has_flags() if (comm->me == 0) { if (KIM_SupportStatus_Equal(kim_model_support_for_energy, KIM_SUPPORT_STATUS_notSupported)) - error->warning(FLERR,"KIM Model does not provide 'partialEnergy'; Potential energy will be zero"); + error->warning(FLERR,"KIM Model does not provide 'partialEnergy'; " + "Potential energy will be zero"); if (KIM_SupportStatus_Equal(kim_model_support_for_forces, KIM_SUPPORT_STATUS_notSupported)) From e317b0eb7e60ff8e8a48f7f4649fc3db6a1d037e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 20 Mar 2023 07:11:47 -0400 Subject: [PATCH 08/22] update kim unit tests for recent changes, make variable for extended test an option --- cmake/Modules/Packages/KIM.cmake | 2 +- unittest/commands/test_kim_commands.cpp | 16 ++++++++-------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/cmake/Modules/Packages/KIM.cmake b/cmake/Modules/Packages/KIM.cmake index 2a2a1cde78..995d2d9490 100644 --- a/cmake/Modules/Packages/KIM.cmake +++ b/cmake/Modules/Packages/KIM.cmake @@ -19,7 +19,7 @@ if(CURL_FOUND) target_compile_definitions(lammps PRIVATE -DLMP_NO_SSL_CHECK) endif() endif() -set(KIM_EXTRA_UNITTESTS OFF CACHE STRING "Set extra unit tests verbose mode on/off. If on, extra tests are included.") +option(KIM_EXTRA_UNITTESTS "Enable extra unit tests for the KIM package." OFF) mark_as_advanced(KIM_EXTRA_UNITTESTS) find_package(PkgConfig QUIET) set(DOWNLOAD_KIM_DEFAULT ON) diff --git a/unittest/commands/test_kim_commands.cpp b/unittest/commands/test_kim_commands.cpp index 988b0db69a..2e6e758c76 100644 --- a/unittest/commands/test_kim_commands.cpp +++ b/unittest/commands/test_kim_commands.cpp @@ -92,14 +92,14 @@ TEST_F(KimCommandsTest, kim_interactions) { if (!LAMMPS::is_installed_pkg("KIM")) GTEST_SKIP(); - TEST_FAILURE(".*ERROR: Illegal 'kim interactions' command.*", command("kim interactions");); + TEST_FAILURE(".*ERROR: Illegal kim interactions command: missing argument.*", + command("kim interactions");); BEGIN_HIDE_OUTPUT(); command("kim init LennardJones_Ar real"); END_HIDE_OUTPUT(); - TEST_FAILURE(".*ERROR: Must use 'kim interactions' command " - "after simulation box is defined.*", + TEST_FAILURE(".*ERROR: Use of 'kim interactions' before simulation box is defined.*", command("kim interactions Ar");); BEGIN_HIDE_OUTPUT(); @@ -410,11 +410,11 @@ TEST_F(KimCommandsTest, kim_property) "3 >= 3.6 support.*", command("kim property");); } else { - TEST_FAILURE(".*ERROR: Invalid 'kim property' command.*", command("kim property");); - TEST_FAILURE(".*ERROR: Invalid 'kim property' command.*", command("kim property create");); - TEST_FAILURE(".*ERROR: Incorrect arguments in 'kim property' command." - "\n'kim property create/destroy/modify/remove/dump' " - "is mandatory.*", + TEST_FAILURE(".*ERROR: Illegal kim property command: missing argument.*", + command("kim property");); + TEST_FAILURE(".*ERROR: Illegal kim property command: missing argument.*", + command("kim property create");); + TEST_FAILURE(".*ERROR: Incorrect first argument unknown to 'kim property' command.*", command("kim property unknown 1 atomic-mass");); } #if defined(KIM_EXTRA_UNITTESTS) From a320f242477e10dc1df4fe4073c8550898d9591a Mon Sep 17 00:00:00 2001 From: agiliopadua Date: Mon, 20 Mar 2023 14:46:01 +0100 Subject: [PATCH 09/22] Fixed bug in fep tools --- tools/fep/fdti.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/fep/fdti.py b/tools/fep/fdti.py index 1f3f86f001..c439f72eb0 100755 --- a/tools/fep/fdti.py +++ b/tools/fep/fdti.py @@ -23,7 +23,7 @@ if len(tok) == 4: v = float(tok[3]) lo = -rt * math.log(float(tok[2]) / v) -i = 0 +i = 1 sum = 0.0 for line in sys.stdin: tok = line.split() From 3280d6b7864a5c15e65e66e268f4e00bf08e0fa5 Mon Sep 17 00:00:00 2001 From: agiliopadua Date: Mon, 20 Mar 2023 14:47:40 +0100 Subject: [PATCH 10/22] Fixed bug in fep tools --- tools/fep/nti.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/fep/nti.py b/tools/fep/nti.py index ffe47fb908..cd426163bc 100755 --- a/tools/fep/nti.py +++ b/tools/fep/nti.py @@ -19,7 +19,7 @@ while line.startswith("#"): tok = line.split() lo = float(tok[1]) -i = 0 +i = 1 sum = 0.0 for line in sys.stdin: tok = line.split() From ecc626964806c5fb3a99dc2a69ac5dd512b3e13e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 20 Mar 2023 10:52:20 -0400 Subject: [PATCH 11/22] add useful comments --- tools/fep/fdti.py | 2 +- tools/fep/fep.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tools/fep/fdti.py b/tools/fep/fdti.py index c439f72eb0..42f0a5fd2a 100755 --- a/tools/fep/fdti.py +++ b/tools/fep/fdti.py @@ -10,7 +10,7 @@ if len(sys.argv) < 3: print("usage: fdti.py temperature hderiv < out.fep") sys.exit() -rt = 0.008314 / 4.184 * float(sys.argv[1]) +rt = 0.008314 / 4.184 * float(sys.argv[1]) # in kcal/mol hderiv = float(sys.argv[2]) line = sys.stdin.readline() diff --git a/tools/fep/fep.py b/tools/fep/fep.py index 22ebaf38c6..c7656c5ef5 100755 --- a/tools/fep/fep.py +++ b/tools/fep/fep.py @@ -9,7 +9,7 @@ if len(sys.argv) < 2: print("usage: fep.py temperature < out.fep") sys.exit() -rt = 0.008314 / 4.184 * float(sys.argv[1]) +rt = 0.008314 / 4.184 * float(sys.argv[1]) # in kcal/mol v = 1.0 sum = 0.0 From ac9389f5cbdaa2d0e16bb805c432498fb391fcf7 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 20 Mar 2023 14:24:21 -0600 Subject: [PATCH 12/22] Slight rewrite --- doc/src/neighbor.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/neighbor.rst b/doc/src/neighbor.rst index bf9a13923c..b5ca4d15a8 100644 --- a/doc/src/neighbor.rst +++ b/doc/src/neighbor.rst @@ -66,7 +66,7 @@ instance in a dense binary system with a ratio of the size of the largest to smallest collection bin :math:`\lamda`, the computational costs of building a default neighbor list grows as :math:`\lamda^6` while the costs for *multi* grows as :math:`\lamda^3`, equivalent to the cost of force -evaluations, as identified in Monti et al. :ref:`(Monti) `. +evaluations, as argued in Monti et al. :ref:`(Monti) `. By default in *multi*, each atom type defines a separate collection of particles. For systems where two or more atom types have the same From d15e13d4756270b42fb7762282028cd69916fddf Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 20 Mar 2023 14:56:39 -0600 Subject: [PATCH 13/22] Reverting mistakenly deleted line, fixing duplicated text in granular doc --- doc/src/Modify_gran_sub_mod.rst | 2 +- doc/src/bond_bpm_spring.rst | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/doc/src/Modify_gran_sub_mod.rst b/doc/src/Modify_gran_sub_mod.rst index e1d559e6bf..b3b4a13cac 100644 --- a/doc/src/Modify_gran_sub_mod.rst +++ b/doc/src/Modify_gran_sub_mod.rst @@ -25,7 +25,7 @@ specific implementation. For instance, from the GranSubModNormal class the GranSubModNormalHooke, GranSubModNormalHertz, and GranSubModNormalJKR classes are derived which calculate Hookean, Hertzian, or JKR normal forces, 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 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 diff --git a/doc/src/bond_bpm_spring.rst b/doc/src/bond_bpm_spring.rst index 1de814b488..93908710d1 100644 --- a/doc/src/bond_bpm_spring.rst +++ b/doc/src/bond_bpm_spring.rst @@ -48,6 +48,8 @@ Examples Description """"""""""" +.. versionadded:: 4May2022 + The *bpm/spring* bond style computes forces based on deviations from the initial reference state of the two atoms. The reference state is stored by each bond when it is first computed in From badfd0bc40072d2da1a4de14eb0c098682b360bb Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 20 Mar 2023 15:34:52 -0600 Subject: [PATCH 14/22] Specifying dimensions, lamda->lambda --- doc/src/neighbor.rst | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/src/neighbor.rst b/doc/src/neighbor.rst index b5ca4d15a8..c4a2caac81 100644 --- a/doc/src/neighbor.rst +++ b/doc/src/neighbor.rst @@ -62,11 +62,11 @@ sets of bins are then used to construct the neighbor lists as as further described by Shire, Hanley, and Stratford :ref:`(Shire) ` and Monti et al. :ref:`(Monti) `. This imposes some extra setup overhead, but the searches themselves may be much faster. For -instance in a dense binary system with a ratio of the size of the largest -to smallest collection bin :math:`\lamda`, the computational costs of -building a default neighbor list grows as :math:`\lamda^6` while the costs -for *multi* grows as :math:`\lamda^3`, equivalent to the cost of force -evaluations, as argued in Monti et al. :ref:`(Monti) `. +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^d` 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) `. By default in *multi*, each atom type defines a separate collection of particles. For systems where two or more atom types have the same From 4351ada794a5d5709c6bdbbde3a19467edd290eb Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 20 Mar 2023 21:05:11 -0400 Subject: [PATCH 15/22] improve error message throughout fix adap --- src/fix_adapt.cpp | 68 +++++++++++++++++++++-------------------------- 1 file changed, 30 insertions(+), 38 deletions(-) diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp index c90aa7c5b6..0b4ec9a638 100644 --- a/src/fix_adapt.cpp +++ b/src/fix_adapt.cpp @@ -48,9 +48,9 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), nadapt(0), anypair(0), anybond(0), anyangle(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr) { - if (narg < 5) error->all(FLERR,"Illegal fix adapt command"); + if (narg < 5) utils::missing_cmd_args(FLERR,"fix adapt", error); nevery = utils::inumeric(FLERR,arg[3],false,lmp); - if (nevery < 0) error->all(FLERR,"Illegal fix adapt command"); + if (nevery < 0) error->all(FLERR,"Illegal fix adapt every value {}", nevery); dynamic_group_allow = 1; create_attribute = 1; @@ -62,29 +62,29 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : int iarg = 4; while (iarg < narg) { if (strcmp(arg[iarg],"pair") == 0) { - if (iarg+6 > narg) error->all(FLERR,"Illegal fix adapt command"); + if (iarg+6 > narg) utils::missing_cmd_args(FLERR,"fix adapt pair", error); nadapt++; iarg += 6; } else if (strcmp(arg[iarg],"kspace") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR,"fix adapt kspace", error); nadapt++; iarg += 2; } else if (strcmp(arg[iarg],"atom") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command"); + if (iarg+3 > narg) utils::missing_cmd_args(FLERR,"fix adapt atom", error); nadapt++; iarg += 3; } else if (strcmp(arg[iarg],"bond") == 0) { - if (iarg+5 > narg) error->all(FLERR,"Illegal fix adapt command"); + if (iarg+5 > narg) utils::missing_cmd_args(FLERR,"fix adapt bond", error); nadapt++; iarg += 5; } else if (strcmp(arg[iarg],"angle") == 0) { - if (iarg+5 > narg) error->all(FLERR,"Illegal fix adapt command"); + if (iarg+5 > narg) utils::missing_cmd_args(FLERR,"fix adapt angle", error); nadapt++; iarg += 5; } else break; } - if (nadapt == 0) error->all(FLERR,"Illegal fix adapt command"); + if (nadapt == 0) error->all(FLERR,"Nothing to adapt in fix adapt command"); adapt = new Adapt[nadapt]; // parse keywords @@ -96,7 +96,6 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : iarg = 4; while (iarg < narg) { if (strcmp(arg[iarg],"pair") == 0) { - if (iarg+6 > narg) error->all(FLERR,"Illegal fix adapt command"); adapt[nadapt].which = PAIR; adapt[nadapt].pair = nullptr; adapt[nadapt].pstyle = utils::strdup(arg[iarg+1]); @@ -107,12 +106,11 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : adapt[nadapt].jlo,adapt[nadapt].jhi,error); if (utils::strmatch(arg[iarg+5],"^v_")) { adapt[nadapt].var = utils::strdup(arg[iarg+5]+2); - } else error->all(FLERR,"Illegal fix adapt command"); + } else error->all(FLERR,"Argument #{} must be variable not {}", iarg+6, arg[iarg+5]); nadapt++; iarg += 6; } else if (strcmp(arg[iarg],"bond") == 0) { - if (iarg+5 > narg) error->all(FLERR, "Illegal fix adapt command"); adapt[nadapt].which = BOND; adapt[nadapt].bond = nullptr; adapt[nadapt].bstyle = utils::strdup(arg[iarg+1]); @@ -121,12 +119,11 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : adapt[nadapt].ilo,adapt[nadapt].ihi,error); if (utils::strmatch(arg[iarg+4],"^v_")) { adapt[nadapt].var = utils::strdup(arg[iarg+4]+2); - } else error->all(FLERR,"Illegal fix adapt command"); + } else error->all(FLERR,"Argument #{} must be variable not {}", iarg+5, arg[iarg+4]); nadapt++; iarg += 5; } else if (strcmp(arg[iarg],"angle") == 0) { - if (iarg+5 > narg) error->all(FLERR, "Illegal fix adapt command"); adapt[nadapt].which = ANGLE; adapt[nadapt].angle = nullptr; adapt[nadapt].astyle = utils::strdup(arg[iarg+1]); @@ -135,21 +132,19 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : adapt[nadapt].ilo,adapt[nadapt].ihi,error); if (utils::strmatch(arg[iarg+4],"^v_")) { adapt[nadapt].var = utils::strdup(arg[iarg+4]+2); - } else error->all(FLERR,"Illegal fix adapt command"); + } else error->all(FLERR,"Argument #{} must be variable not {}", iarg+5, arg[iarg+4]); nadapt++; iarg += 5; } else if (strcmp(arg[iarg],"kspace") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); adapt[nadapt].which = KSPACE; if (utils::strmatch(arg[iarg+1],"^v_")) { adapt[nadapt].var = utils::strdup(arg[iarg+1]+2); - } else error->all(FLERR,"Illegal fix adapt command"); + } else error->all(FLERR,"Argument #{} must be variable not {}", iarg+2, arg[iarg+1]); nadapt++; iarg += 2; } else if (strcmp(arg[iarg],"atom") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command"); adapt[nadapt].which = ATOM; if (strcmp(arg[iarg+1],"diameter") == 0 || strcmp(arg[iarg+1],"diameter/disc") == 0) { @@ -160,10 +155,10 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg+1],"charge") == 0) { adapt[nadapt].atomparam = CHARGE; chgflag = 1; - } else error->all(FLERR,"Illegal fix adapt command"); + } else error->all(FLERR,"Unsupported per-atom property {} for fix adapt", arg[iarg+1]); if (utils::strmatch(arg[iarg+2],"^v_")) { adapt[nadapt].var = utils::strdup(arg[iarg+2]+2); - } else error->all(FLERR,"Illegal fix adapt command"); + } else error->all(FLERR,"Argument #{} must be variable not {}", iarg+3, arg[iarg+2]); nadapt++; iarg += 3; } else break; @@ -177,18 +172,18 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : while (iarg < narg) { if (strcmp(arg[iarg],"reset") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR,"fix adapt reset", error); resetflag = utils::logical(FLERR,arg[iarg+1],false,lmp); iarg += 2; } else if (strcmp(arg[iarg],"scale") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR,"fix adapt scale", error); scaleflag = utils::logical(FLERR,arg[iarg+1],false,lmp); iarg += 2; } else if (strcmp(arg[iarg],"mass") == 0) { - if (iarg+2 > narg)error->all(FLERR,"Illegal fix adapt command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR,"fix adapt mass", error); massflag = utils::logical(FLERR,arg[iarg+1],false,lmp); iarg += 2; - } else error->all(FLERR,"Illegal fix adapt command"); + } else error->all(FLERR,"Unknown fix adapt keyword {}", arg[iarg]); } // if scaleflag set with diameter or charge adaptation, @@ -202,22 +197,19 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : int n = atom->ntypes; for (int m = 0; m < nadapt; m++) - if (adapt[m].which == PAIR) - memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig"); + if (adapt[m].which == PAIR) memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig"); // allocate bond style arrays: n = atom->nbondtypes; for (int m = 0; m < nadapt; ++m) - if (adapt[m].which == BOND) - memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig"); + if (adapt[m].which == BOND) memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig"); // allocate angle style arrays: n = atom->nbondtypes; for (int m = 0; m < nadapt; ++m) - if (adapt[m].which == ANGLE) - memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig"); + if (adapt[m].which == ANGLE) memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig"); } /* ---------------------------------------------------------------------- */ @@ -324,7 +316,7 @@ void FixAdapt::init() if (group->dynamic[igroup]) for (i = 0; i < nadapt; i++) if (adapt[i].which == ATOM) - error->all(FLERR,"Cannot use dynamic group with fix adapt atom"); + error->all(FLERR,"Cannot use dynamic group {} with fix adapt atom", group->names[igroup]); // setup and error checks @@ -337,9 +329,9 @@ void FixAdapt::init() ad->ivar = input->variable->find(ad->var); if (ad->ivar < 0) - error->all(FLERR,"Variable name for fix adapt does not exist"); + error->all(FLERR,"Variable name {} for fix adapt does not exist", ad->var); if (!input->variable->equalstyle(ad->ivar)) - error->all(FLERR,"Variable for fix adapt is invalid style"); + error->all(FLERR,"Variable {} for fix adapt is invalid style", ad->var); if (ad->which == PAIR) { anypair = 1; @@ -374,7 +366,7 @@ void FixAdapt::init() // scalars are supported if (ad->pdim != 2 && ad->pdim != 0) - error->all(FLERR,"Fix adapt pair style param is not compatible"); + error->all(FLERR,"Pair style parameter {} is not compatible with fix adapt", ad->pparam); if (ad->pdim == 2) ad->array = (double **) ptr; if (ad->pdim == 0) ad->scalar = (double *) ptr; @@ -402,12 +394,12 @@ void FixAdapt::init() if (ad->bond == nullptr) ad->bond = force->bond_match(bstyle); if (ad->bond == nullptr ) - error->all(FLERR,"Fix adapt bond style does not exist"); + error->all(FLERR,"Fix adapt bond style {} does not exist", bstyle); void *ptr = ad->bond->extract(ad->bparam,ad->bdim); if (ptr == nullptr) - error->all(FLERR,"Fix adapt bond style param not supported"); + error->all(FLERR,"Fix adapt bond style parameter {} not supported", ad->bparam); // for bond styles, use a vector @@ -428,12 +420,12 @@ void FixAdapt::init() if (ad->angle == nullptr) ad->angle = force->angle_match(astyle); if (ad->angle == nullptr ) - error->all(FLERR,"Fix adapt angle style does not exist"); + error->all(FLERR,"Fix adapt angle style {} does not exist", astyle); void *ptr = ad->angle->extract(ad->aparam,ad->adim); if (ptr == nullptr) - error->all(FLERR,"Fix adapt angle style param not supported"); + error->all(FLERR,"Fix adapt angle style parameter {} not supported", ad->aparam); // for angle styles, use a vector @@ -446,7 +438,7 @@ void FixAdapt::init() } else if (ad->which == KSPACE) { if (force->kspace == nullptr) - error->all(FLERR,"Fix adapt kspace style does not exist"); + error->all(FLERR,"Fix adapt expected a kspace style but none was defined"); kspace_scale = (double *) force->kspace->extract("scale"); } else if (ad->which == ATOM) { From 1101383b51b966fc59b354e6c6f28cef84ee1db5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 20 Mar 2023 22:07:02 -0400 Subject: [PATCH 16/22] add versionadded tags --- doc/src/bond_bpm_rotational.rst | 2 ++ doc/src/bond_bpm_spring.rst | 2 ++ 2 files changed, 4 insertions(+) diff --git a/doc/src/bond_bpm_rotational.rst b/doc/src/bond_bpm_rotational.rst index 19fd607169..0b3fa4442c 100644 --- a/doc/src/bond_bpm_rotational.rst +++ b/doc/src/bond_bpm_rotational.rst @@ -143,6 +143,8 @@ the *overlay/pair* keyword. These settings require specific restrictions. Further details can be found in the `:doc: how to ` 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. diff --git a/doc/src/bond_bpm_spring.rst b/doc/src/bond_bpm_spring.rst index 93908710d1..5391d81420 100644 --- a/doc/src/bond_bpm_spring.rst +++ b/doc/src/bond_bpm_spring.rst @@ -109,6 +109,8 @@ the *overlay/pair* keyword. These settings require specific restrictions. Further details can be found in the `:doc: how to ` 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. From 669397b0929024a5f82088502f0c6ce4cdddaae4 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Wed, 22 Mar 2023 09:36:29 -0600 Subject: [PATCH 17/22] fixing exponent in multi documentation --- doc/src/neighbor.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/neighbor.rst b/doc/src/neighbor.rst index c4a2caac81..4ee3d859f8 100644 --- a/doc/src/neighbor.rst +++ b/doc/src/neighbor.rst @@ -64,7 +64,7 @@ and Monti et al. :ref:`(Monti) `. This imposes some extra 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^d` while +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) `. From cc2106397af7f2cd3f3f07d12091350d7cb2f01e Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Wed, 22 Mar 2023 10:18:50 -0600 Subject: [PATCH 18/22] fix ids once bug in compute chunk/atom --- src/compute_chunk_atom.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/compute_chunk_atom.cpp b/src/compute_chunk_atom.cpp index f5c68c8f25..fc70a3246f 100644 --- a/src/compute_chunk_atom.cpp +++ b/src/compute_chunk_atom.cpp @@ -691,9 +691,10 @@ void ComputeChunkAtom::compute_ichunk() if (invoked_ichunk == update->ntimestep) return; // 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 // else proceed to recalculate per-atom chunk assignments + // if restoring, update invoked_ichunk only for NFREQ case const int nlocal = atom->nlocal; int restore = 0; @@ -701,7 +702,7 @@ void ComputeChunkAtom::compute_ichunk() if (idsflag == NFREQ && lockfix && update->ntimestep > lockstart) restore = 1; if (restore) { - invoked_ichunk = update->ntimestep; + if (idsflag == NFREQ) invoked_ichunk = update->ntimestep; double *vstore = fixstore->vstore; for (i = 0; i < nlocal; i++) ichunk[i] = static_cast(vstore[i]); return; From 1370f0571486880179310483b1c6977dd517403d Mon Sep 17 00:00:00 2001 From: jtclemm Date: Wed, 22 Mar 2023 11:33:17 -0600 Subject: [PATCH 19/22] Elaborating on the scaling of multi --- doc/src/neighbor.rst | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/doc/src/neighbor.rst b/doc/src/neighbor.rst index 4ee3d859f8..9dcdd1daf2 100644 --- a/doc/src/neighbor.rst +++ b/doc/src/neighbor.rst @@ -61,12 +61,17 @@ or its cutoff distance for interacting with other particles. Different sets of bins are then used to construct the neighbor lists as as further described by Shire, Hanley, and Stratford :ref:`(Shire) ` and Monti et al. :ref:`(Monti) `. This imposes some extra -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 +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) `. +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 @@ -82,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 different atom types. This algorithm used to be the default *multi* 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 -where multi/old outperforms the new multi style. +approach. For the dense binary system, computational costs still grew as +: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:: From b5e1bbfa6f9294d837df9f26c3296e6b19cf04fe Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 22 Mar 2023 15:02:14 -0400 Subject: [PATCH 20/22] move fix property/atom property type enumerator to class definition in header --- src/KOKKOS/fix_property_atom_kokkos.cpp | 2 -- src/fix_property_atom.cpp | 2 -- src/fix_property_atom.h | 3 ++- 3 files changed, 2 insertions(+), 5 deletions(-) diff --git a/src/KOKKOS/fix_property_atom_kokkos.cpp b/src/KOKKOS/fix_property_atom_kokkos.cpp index 2c515510ee..6c7c7c6d37 100644 --- a/src/KOKKOS/fix_property_atom_kokkos.cpp +++ b/src/KOKKOS/fix_property_atom_kokkos.cpp @@ -24,8 +24,6 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum{MOLECULE,CHARGE,RMASS,INTEGER,DOUBLE}; - /* ---------------------------------------------------------------------- */ FixPropertyAtomKokkos::FixPropertyAtomKokkos(LAMMPS *lmp, int narg, char **arg) : diff --git a/src/fix_property_atom.cpp b/src/fix_property_atom.cpp index 4b43feca74..994b4f0f19 100644 --- a/src/fix_property_atom.cpp +++ b/src/fix_property_atom.cpp @@ -25,8 +25,6 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum { MOLECULE, CHARGE, RMASS, TEMPERATURE, HEATFLOW, IVEC, DVEC, IARRAY, DARRAY }; - /* ---------------------------------------------------------------------- */ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) : diff --git a/src/fix_property_atom.h b/src/fix_property_atom.h index 7d9f50b75e..f2f26d511b 100644 --- a/src/fix_property_atom.h +++ b/src/fix_property_atom.h @@ -27,11 +27,12 @@ namespace LAMMPS_NS { class FixPropertyAtom : public Fix { public: FixPropertyAtom(class LAMMPS *, int, char **); - ~FixPropertyAtom() override; int setmask() override; void init() override; + enum { MOLECULE, CHARGE, RMASS, TEMPERATURE, HEATFLOW, IVEC, DVEC, IARRAY, DARRAY }; + void read_data_section(char *, int, char *, tagint) override; bigint read_data_skip_lines(char *) override; void write_data_section_size(int, int &, int &) override; From 56ac387e6f4e0c0b619da85c7366c67e18a41dbb Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 22 Mar 2023 15:02:49 -0400 Subject: [PATCH 21/22] synchronize Kokkos version of grow_atoms() with base version --- src/KOKKOS/fix_property_atom_kokkos.cpp | 36 ++++++++++++++++--------- 1 file changed, 24 insertions(+), 12 deletions(-) diff --git a/src/KOKKOS/fix_property_atom_kokkos.cpp b/src/KOKKOS/fix_property_atom_kokkos.cpp index 6c7c7c6d37..1de07b39dc 100644 --- a/src/KOKKOS/fix_property_atom_kokkos.cpp +++ b/src/KOKKOS/fix_property_atom_kokkos.cpp @@ -42,33 +42,45 @@ FixPropertyAtomKokkos::FixPropertyAtomKokkos(LAMMPS *lmp, int narg, char **arg) void FixPropertyAtomKokkos::grow_arrays(int nmax) { - for (int m = 0; m < nvalue; m++) { - if (styles[m] == MOLECULE) { + for (int nv = 0; nv < nvalue; nv++) { + if (styles[nv] == MOLECULE) { memory->grow(atom->molecule,nmax,"atom:molecule"); size_t nbytes = (nmax-nmax_old) * sizeof(tagint); memset(&atom->molecule[nmax_old],0,nbytes); - } else if (styles[m] == CHARGE) { + } else if (styles[nv] == CHARGE) { memory->grow(atom->q,nmax,"atom:q"); size_t nbytes = (nmax-nmax_old) * sizeof(double); memset(&atom->q[nmax_old],0,nbytes); - } else if (styles[m] == RMASS) { + } else if (styles[nv] == RMASS) { memory->grow(atom->rmass,nmax,"atom:rmass"); size_t nbytes = (nmax-nmax_old) * sizeof(double); memset(&atom->rmass[nmax_old],0,nbytes); - } else if (styles[m] == INTEGER) { - memory->grow(atom->ivector[index[m]],nmax,"atom:ivector"); + } else if (styles[nv] == TEMPERATURE) { + memory->grow(atom->temperature, nmax, "atom:temperature"); + size_t nbytes = (nmax - nmax_old) * sizeof(double); + memset(&atom->temperature[nmax_old], 0, nbytes); + } else if (styles[nv] == HEATFLOW) { + memory->grow(atom->heatflow, nmax, "atom:heatflow"); + size_t nbytes = (nmax - nmax_old) * sizeof(double); + memset(&atom->heatflow[nmax_old], 0, nbytes); + } else if (styles[nv] == IVEC) { + memory->grow(atom->ivector[index[nv]],nmax,"atom:ivector"); size_t nbytes = (nmax-nmax_old) * sizeof(int); - memset(&atom->ivector[index[m]][nmax_old],0,nbytes); - } else if (styles[m] == DOUBLE) { + memset(&atom->ivector[index[nv]][nmax_old],0,nbytes); + } else if (styles[nv] == DVEC) { atomKK->sync(Device,DVECTOR_MASK); memoryKK->grow_kokkos(atomKK->k_dvector,atomKK->dvector,atomKK->k_dvector.extent(0),nmax, "atom:dvector"); atomKK->modified(Device,DVECTOR_MASK); - //memory->grow(atom->dvector[index[m]],nmax,"atom:dvector"); - //size_t nbytes = (nmax-nmax_old) * sizeof(double); - //memset(&atom->dvector[index[m]][nmax_old],0,nbytes); + } else if (styles[nv] == IARRAY) { + memory->grow(atom->iarray[index[nv]], nmax, cols[nv], "atom:iarray"); + size_t nbytes = (size_t) (nmax - nmax_old) * cols[nv] * sizeof(int); + if (nbytes) memset(&atom->iarray[index[nv]][nmax_old][0], 0, nbytes); + } else if (styles[nv] == DARRAY) { + memory->grow(atom->darray[index[nv]], nmax, cols[nv], "atom:darray"); + size_t nbytes = (size_t) (nmax - nmax_old) * cols[nv] * sizeof(double); + if (nbytes) memset(&atom->darray[index[nv]][nmax_old][0], 0, nbytes); } } - nmax_old = nmax; } From 872a4e299924e91679158ab5bf6417d48f30c331 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 22 Mar 2023 15:03:52 -0400 Subject: [PATCH 22/22] cosmetic. remove commented out debug code. --- src/KOKKOS/fix_rx_kokkos.cpp | 215 ----------------------------------- src/fix_property_atom.h | 8 +- 2 files changed, 4 insertions(+), 219 deletions(-) diff --git a/src/KOKKOS/fix_rx_kokkos.cpp b/src/KOKKOS/fix_rx_kokkos.cpp index 963fb32a1a..ac1821e142 100644 --- a/src/KOKKOS/fix_rx_kokkos.cpp +++ b/src/KOKKOS/fix_rx_kokkos.cpp @@ -74,23 +74,17 @@ FixRxKokkos::FixRxKokkos(LAMMPS *lmp, int narg, char **arg) : datamask_modify = EMPTY_MASK; k_error_flag = DAT::tdual_int_scalar("FixRxKokkos::k_error_flag"); - - //printf("Inside FixRxKokkos::FixRxKokkos\n"); } template FixRxKokkos::~FixRxKokkos() { - //printf("Inside FixRxKokkos::~FixRxKokkos copymode= %d\n", copymode); if (copymode) return; if (localTempFlag) memoryKK->destroy_kokkos(k_dpdThetaLocal, dpdThetaLocal); memoryKK->destroy_kokkos(k_sumWeights, sumWeights); - //memoryKK->destroy_kokkos(k_sumWeights); - - //delete [] scratchSpace; memoryKK->destroy_kokkos(d_scratchSpace); memoryKK->destroy_kokkos(k_cutsq); @@ -147,7 +141,6 @@ void FixRxKokkos::init() template void FixRxKokkos::init_list(int, class NeighList* ptr) { - //printf("Inside FixRxKokkos::init_list\n"); this->list = ptr; } @@ -551,7 +544,6 @@ void FixRxKokkos::k_rkf45(const int neq, const double t_stop, Vector nfe += 6; if (maxIters && nit > maxIters) { - //fprintf(stderr,"atom[%d] took too many iterations in rkf45 %d %e %e\n", id, nit, t, t_stop); counter.nFails ++; break; // We should set an error here so that the solution is not used! @@ -562,8 +554,6 @@ void FixRxKokkos::k_rkf45(const int neq, const double t_stop, Vector counter.nSteps += nst; counter.nIters += nit; counter.nFuncs += nfe; - - //printf("id= %d nst= %d nit= %d\n", id, nst, nit); } /* ---------------------------------------------------------------------- */ @@ -733,9 +723,6 @@ int FixRxKokkos::rkf45_h0(const int neq, const double t, const doubl yddnrm = sqrt( yddnrm / double(neq) ); - //std::cout << "iter " << _iter << " hg " << hg << " y'' " << yddnrm << std::endl; - //std::cout << "ydot " << ydot[neq-1] << std::endl; - // should we accept this? if (hnew_is_ok || iter == max_iters) { hnew = hg; @@ -760,8 +747,6 @@ int FixRxKokkos::rkf45_h0(const int neq, const double t, const doubl hnew_is_ok = true; } - //printf("iter=%d, yddnrw=%e, hnew=%e, hmin=%e, hmax=%e\n", iter, yddnrm, hnew, hmin, hmax); - hg = hnew; iter ++; } @@ -805,13 +790,9 @@ void FixRxKokkos::rkf45(const int neq, const double t_stop, double * double t = 0.0; if (h < h_min) { - //fprintf(stderr,"hin not implemented yet\n"); - //exit(-1); nfe = rkf45_h0 (neq, t, t_stop, h_min, h_max, h, y, rwork, v_param); } - //printf("t= %e t_stop= %e h= %e\n", t, t_stop, h); - // Integrate until we reach the end time. while (fabs(t - t_stop) > tround) { double *yout = rwork; @@ -865,7 +846,6 @@ void FixRxKokkos::rkf45(const int neq, const double t_stop, double * nfe += 6; if (maxIters && nit > maxIters) { - //fprintf(stderr,"atom[%d] took too many iterations in rkf45 %d %e %e\n", id, nit, t, t_stop); counter.nFails ++; break; // We should set an error here so that the solution is not used! @@ -876,8 +856,6 @@ void FixRxKokkos::rkf45(const int neq, const double t_stop, double * counter.nSteps += nst; counter.nIters += nit; counter.nFuncs += nfe; - - //printf("id= %d nst= %d nit= %d\n", id, nst, nit); } /* ---------------------------------------------------------------------- */ @@ -902,9 +880,6 @@ int FixRxKokkos::rhs_dense(double /*t*/, const double *y, double *dy double *rxnRateLaw = userData->rxnRateLaw; double *kFor = userData->kFor; - //const double VDPD = domain->xprd * domain->yprd * domain->zprd / atom->natoms; - //const int nspecies = atom->nspecies_dpd; - for (int ispecies=0; ispecies::rhs_dense(double /*t*/, const double *y, double *dy for (int ispecies=0; ispecies::rhs_dense(double /*t*/, const double *y, double *dy for (int jrxn=0; jrxn::rhs_sparse(double /*t*/, const double *y, double *d { UserRHSData *userData = (UserRHSData *) v_params; - //const double VDPD = domain->xprd * domain->yprd * domain->zprd / atom->natoms; - #define kFor (userData->kFor) #define kRev (nullptr) #define rxnRateLaw (userData->rxnRateLaw) @@ -964,7 +935,6 @@ int FixRxKokkos::rhs_sparse(double /*t*/, const double *y, double *d for (int kk = 1; kk < maxReactants; ++kk) { const int k = nuk(i,kk); if (k == SparseKinetics_invalidIndex) break; - //if (k != SparseKinetics_invalidIndex) rxnRateLawForward *= powint( conc[k], inu(i,kk) ); } } else { @@ -972,7 +942,6 @@ int FixRxKokkos::rhs_sparse(double /*t*/, const double *y, double *d for (int kk = 1; kk < maxReactants; ++kk) { const int k = nuk(i,kk); if (k == SparseKinetics_invalidIndex) break; - //if (k != SparseKinetics_invalidIndex) rxnRateLawForward *= pow( conc[k], nu(i,kk) ); } } @@ -991,7 +960,6 @@ int FixRxKokkos::rhs_sparse(double /*t*/, const double *y, double *d for (int kk = 1; kk < maxReactants; ++kk) { const int k = nuk(i,kk); if (k == SparseKinetics_invalidIndex) break; - //if (k != SparseKinetics_invalidIndex) dydt[k] -= nu(i,kk) * rxnRateLaw[i]; } @@ -1000,7 +968,6 @@ int FixRxKokkos::rhs_sparse(double /*t*/, const double *y, double *d for (int kk = maxReactants+1; kk < maxSpecies; ++kk) { const int k = nuk(i,kk); if (k == SparseKinetics_invalidIndex) break; - //if (k != SparseKinetics_invalidIndex) dydt[k] += nu(i,kk) * rxnRateLaw[i]; } } @@ -1048,9 +1015,6 @@ int FixRxKokkos::k_rhs_dense(double /*t*/, const VectorType& y, Vect #define rxnRateLaw (userData.rxnRateLaw) #define kFor (userData.kFor ) - //const double VDPD = domain->xprd * domain->yprd * domain->zprd / atom->natoms; - //const int nspecies = atom->nspecies_dpd; - for (int ispecies=0; ispecies::k_rhs_dense(double /*t*/, const VectorType& y, Vect for (int ispecies=0; ispecies::k_rhs_dense(double /*t*/, const VectorType& y, Vect for (int jrxn=0; jrxn::k_rhs_sparse(double /*t*/, const VectorType& y, Vec for (int kk = 1; kk < maxReactants; ++kk) { const int k = nuk(i,kk); if (k == SparseKinetics_invalidIndex) break; - //if (k != SparseKinetics_invalidIndex) rxnRateLawForward *= powint( conc[k], inu(i,kk) ); } } else { @@ -1119,7 +1080,6 @@ int FixRxKokkos::k_rhs_sparse(double /*t*/, const VectorType& y, Vec for (int kk = 1; kk < maxReactants; ++kk) { const int k = nuk(i,kk); if (k == SparseKinetics_invalidIndex) break; - //if (k != SparseKinetics_invalidIndex) rxnRateLawForward *= pow( conc[k], nu(i,kk) ); } } @@ -1138,7 +1098,6 @@ int FixRxKokkos::k_rhs_sparse(double /*t*/, const VectorType& y, Vec for (int kk = 1; kk < maxReactants; ++kk) { const int k = nuk(i,kk); if (k == SparseKinetics_invalidIndex) break; - //if (k != SparseKinetics_invalidIndex) dydt[k] -= nu(i,kk) * rxnRateLaw[i]; } @@ -1147,7 +1106,6 @@ int FixRxKokkos::k_rhs_sparse(double /*t*/, const VectorType& y, Vec for (int kk = maxReactants+1; kk < maxSpecies; ++kk) { const int k = nuk(i,kk); if (k == SparseKinetics_invalidIndex) break; - //if (k != SparseKinetics_invalidIndex) dydt[k] += nu(i,kk) * rxnRateLaw[i]; } } @@ -1215,8 +1173,6 @@ void FixRxKokkos::operator()(SolverType, const int &i) const template void FixRxKokkos::create_kinetics_data() { - //printf("Inside FixRxKokkos::create_kinetics_data\n"); - memoryKK->create_kokkos( d_kineticsData.Arr, h_kineticsData.Arr, nreactions, "KineticsType::Arr"); memoryKK->create_kokkos( d_kineticsData.nArr, h_kineticsData.nArr, nreactions, "KineticsType::nArr"); memoryKK->create_kokkos( d_kineticsData.Ea, h_kineticsData.Ea, nreactions, "KineticsType::Ea"); @@ -1296,8 +1252,6 @@ void FixRxKokkos::create_kinetics_data() template void FixRxKokkos::setup_pre_force(int vflag) { - //printf("Inside FixRxKokkos::setup_pre_force restartFlag= %d\n", my_restartFlag); - if (my_restartFlag) my_restartFlag = 0; else @@ -1309,8 +1263,6 @@ void FixRxKokkos::setup_pre_force(int vflag) template void FixRxKokkos::pre_force(int vflag) { - //printf("Inside FixRxKokkos::pre_force localTempFlag= %d\n", localTempFlag); - this->solve_reactions( vflag, true ); } @@ -1407,8 +1359,6 @@ void FixRxKokkos::operator()(Tag_FixRxKokkos_solveSystems void FixRxKokkos::solve_reactions(const int /*vflag*/, const bool isPreForce) { - //printf("Inside FixRxKokkos::solve_reactions localTempFlag= %d isPreForce= %s\n", localTempFlag, isPreForce ? "True" : "false"); - copymode = 1; if (update_kinetics_data) @@ -1416,7 +1366,6 @@ void FixRxKokkos::solve_reactions(const int /*vflag*/, const bool is //TimerType timer_start = getTimeStamp(); - //const int nlocal = atom->nlocal; this->nlocal = atom->nlocal; const int nghost = atom->nghost; const int newton_pair = force->newton_pair; @@ -1477,9 +1426,6 @@ void FixRxKokkos::solve_reactions(const int /*vflag*/, const bool is // ... // Local references to the atomKK objects. - //typename ArrayTypes::t_efloat_1d d_dpdTheta = atomKK->k_dpdTheta.view(); - //typename ArrayTypes::t_float_2d d_dvector = atomKK->k_dvector.view(); - //typename ArrayTypes::t_int_1d d_mask = atomKK->k_mask.view(); this->d_dpdTheta = atomKK->k_dpdTheta.view(); this->d_dvector = atomKK->k_dvector.view(); this->d_mask = atomKK->k_mask.view(); @@ -1488,8 +1434,6 @@ void FixRxKokkos::solve_reactions(const int /*vflag*/, const bool is atomKK->sync( execution_space, MASK_MASK | DVECTOR_MASK | DPDTHETA_MASK ); // Set some constants outside of the parallel_for - //const double boltz = force->boltz; - //const double t_stop = update->dt; // DPD time-step and integration length. this->boltz = force->boltz; this->t_stop = update->dt; // DPD time-step and integration length. @@ -1505,13 +1449,6 @@ void FixRxKokkos::solve_reactions(const int /*vflag*/, const bool is d_diagnosticCounterPerODEnFuncs = k_diagnosticCounterPerODEnFuncs.template view(); Kokkos::parallel_for ( Kokkos::RangePolicy(0,nlocal), *this); - //Kokkos::parallel_for ( nlocal, - // LAMMPS_LAMBDA(const int i) - // { - // d_diagnosticCounterPerODEnSteps(i) = 0; - // d_diagnosticCounterPerODEnFuncs(i) = 0; - // } - // ); } // Error flag for any failures. @@ -1523,116 +1460,17 @@ void FixRxKokkos::solve_reactions(const int /*vflag*/, const bool is k_error_flag.template sync(); // Create scratch array space. - //const size_t scratchSpaceSize = (8*nspecies + 2*nreactions); this->scratchSpaceSize = (8*nspecies + 2*nreactions); - //double *scratchSpace = new double[ scratchSpaceSize * nlocal ]; - //typename ArrayTypes::t_double_1d d_scratchSpace("d_scratchSpace", scratchSpaceSize * nlocal); if (nlocal*scratchSpaceSize > d_scratchSpace.extent(0)) { memoryKK->destroy_kokkos (d_scratchSpace); memoryKK->create_kokkos (d_scratchSpace, nlocal*scratchSpaceSize, "FixRxKokkos::d_scratchSpace"); } -#if 0 - Kokkos::parallel_reduce( nlocal, LAMMPS_LAMBDA(int i, CounterType &counter) - { - if (d_mask(i) & groupbit) - { - //double *y = new double[8*nspecies]; - //double *rwork = y + nspecies; - - //StridedArrayType _y( y ); - //StridedArrayType _rwork( rwork ); - - StridedArrayType y( d_scratchSpace.data() + scratchSpaceSize * i ); - StridedArrayType rwork( &y[nspecies] ); - - //UserRHSData userData; - //userData.kFor = new double[nreactions]; - //userData.rxnRateLaw = new double[nreactions]; - - //UserRHSDataKokkos<1> userDataKokkos; - //userDataKokkos.kFor.m_data = userData.kFor; - //userDataKokkos.rxnRateLaw.m_data = userData.rxnRateLaw; - - UserRHSDataKokkos<1> userData; - userData.kFor.m_data = &( rwork[7*nspecies] ); - userData.rxnRateLaw.m_data = &( userData.kFor[ nreactions ] ); - - CounterType counter_i; - - const double theta = (localTempFlag) ? d_dpdThetaLocal(i) : d_dpdTheta(i); - - //Compute the reaction rate constants - for (int irxn = 0; irxn < nreactions; irxn++) - { - if (setRatesToZero) - userData.kFor[irxn] = 0.0; - else - { - userData.kFor[irxn] = d_kineticsData.Arr(irxn) * - pow(theta, d_kineticsData.nArr(irxn)) * - exp(-d_kineticsData.Ea(irxn) / boltz / theta); - } - } - - // Update ConcOld and initialize the ODE solution vector y[]. - for (int ispecies = 0; ispecies < nspecies; ispecies++) - { - const double tmp = d_dvector(ispecies, i); - d_dvector(ispecies+nspecies, i) = tmp; - y[ispecies] = tmp; - } - - // Solver the ODE system. - if (odeIntegrationFlag == ODE_LAMMPS_RK4) - { - k_rk4(t_stop, y, rwork, userData); - } - else if (odeIntegrationFlag == ODE_LAMMPS_RKF45) - { - k_rkf45(nspecies, t_stop, y, rwork, userData, counter_i); - - if (diagnosticFrequency == 1) - { - d_diagnosticCounterPerODEnSteps(i) = counter_i.nSteps; - d_diagnosticCounterPerODEnFuncs(i) = counter_i.nFuncs; - } - } - - // Store the solution back in dvector. - for (int ispecies = 0; ispecies < nspecies; ispecies++) - { - if (y[ispecies] < -1.0e-10) - { - //error->one(FLERR,"Computed concentration in RK solver is < -1.0e-10"); - k_error_flag.template view()() = 2; - // This should be an atomic update. - } - else if (y[ispecies] < MY_EPSILON) - y[ispecies] = 0.0; - - d_dvector(ispecies,i) = y[ispecies]; - } - - //delete [] y; - //delete [] userData.kFor; - //delete [] userData.rxnRateLaw; - - // Update the iteration statistics counter. Is this unique for each iteration? - counter += counter_i; - - } // if - } // parallel_for lambda-body - - , TotalCounters // reduction value for all iterations. - ); -#else if (setRatesToZero) Kokkos::parallel_reduce( Kokkos::RangePolicy >(0,nlocal), *this, TotalCounters); else Kokkos::parallel_reduce( Kokkos::RangePolicy >(0,nlocal), *this, TotalCounters); -#endif TimerType timer_ODE = getTimeStamp(); @@ -1652,16 +1490,8 @@ void FixRxKokkos::solve_reactions(const int /*vflag*/, const bool is atomKK->modified ( Host, DVECTOR_MASK ); - //TimerType timer_stop = getTimeStamp(); - double time_ODE = getElapsedTime(timer_localTemperature, timer_ODE); - //printf("me= %d kokkos total= %g temp= %g ode= %g comm= %g nlocal= %d nfc= %d %d\n", comm->me, - // getElapsedTime(timer_start, timer_stop), - // getElapsedTime(timer_start, timer_localTemperature), - // getElapsedTime(timer_localTemperature, timer_ODE), - // getElapsedTime(timer_ODE, timer_stop), nlocal, TotalCounters.nFuncs, TotalCounters.nSteps); - // Warn the user if a failure was detected in the ODE solver. if (TotalCounters.nFails > 0) { char sbuf[128]; @@ -1731,7 +1561,6 @@ void FixRxKokkos::odeDiagnostics() // Compute counters per dpd time-step. for (int i = 0; i < numCounters; ++i) { my_vals[i] = this->diagnosticCounter[i] / nTimes; - //printf("my sum[%d] = %f %d\n", i, my_vals[i], comm->me); } MPI_Allreduce (my_vals, sums, numCounters, MPI_DOUBLE, MPI_SUM, world); @@ -1770,7 +1599,6 @@ void FixRxKokkos::odeDiagnostics() double my_max[numCounters], my_min[numCounters]; - //const int nlocal = atom->nlocal; nlocal = atom->nlocal; HAT::t_int_1d h_mask = atomKK->k_mask.h_view; @@ -1980,9 +1808,6 @@ template template void FixRxKokkos::computeLocalTemperature() { - //typename ArrayTypes::t_x_array_randomread d_x = atomKK->k_x.view(); - //typename ArrayTypes::t_int_1d_randomread d_type = atomKK->k_type.view(); - //typename ArrayTypes::t_efloat_1d d_dpdTheta = atomKK->k_dpdTheta.view(); d_x = atomKK->k_x.view(); d_type = atomKK->k_type.view(); d_dpdTheta = atomKK->k_dpdTheta.view(); @@ -1993,21 +1818,9 @@ void FixRxKokkos::computeLocalTemperature() nlocal = atom->nlocal; const int nghost = atom->nghost; - //printf("Inside FixRxKokkos::computeLocalTemperature: %d %d %d %d %d %d %d\n", WT_FLAG, LOCAL_TEMP_FLAG, NEWTON_PAIR, (int)lmp->kokkos->neighflag, NEIGHFLAG, nlocal, nghost); - - // Pull from pairDPDE. The pairDPDEKK objects are protected so recreate here for now. - //pairDPDEKK->k_cutsq.template sync(); - //typename ArrayTypes::t_ffloat_2d d_cutsq = pairDPDEKK->k_cutsq.template view::tdual_ffloat_2d k_cutsq; - //typename ArrayTypes::t_ffloat_2d d_cutsq; - //double **h_cutsq; - { const int ntypes = atom->ntypes; - //memoryKK->create_kokkos (k_cutsq, h_cutsq, ntypes+1, ntypes+1, "pair:cutsq"); if (ntypes+1 > (int) k_cutsq.extent(0)) { memoryKK->destroy_kokkos (k_cutsq); memoryKK->create_kokkos (k_cutsq, ntypes+1, ntypes+1, "FixRxKokkos::k_cutsq"); @@ -2028,7 +1841,6 @@ void FixRxKokkos::computeLocalTemperature() // Initialize the local temperature weight array int sumWeightsCt = nlocal + (NEWTON_PAIR ? nghost : 0); - //memoryKK->create_kokkos (k_sumWeights, sumWeights, sumWeightsCt, "FixRxKokkos::sumWeights"); if (sumWeightsCt > (int)k_sumWeights.template view().extent(0)) { memoryKK->destroy_kokkos(k_sumWeights, sumWeights); memoryKK->create_kokkos (k_sumWeights, sumWeightsCt, "FixRxKokkos::sumWeights"); @@ -2036,14 +1848,6 @@ void FixRxKokkos::computeLocalTemperature() h_sumWeights = k_sumWeights.h_view; } - // Initialize the accumulator to zero ... - //Kokkos::parallel_for (sumWeightsCt, - // LAMMPS_LAMBDA(const int i) - // { - // d_sumWeights(i) = 0.0; - // } - // ); - Kokkos::parallel_for (Kokkos::RangePolicy(0, sumWeightsCt), *this); // Local list views. (This isn't working!) @@ -2051,9 +1855,6 @@ void FixRxKokkos::computeLocalTemperature() if (!list->kokkos) error->one(FLERR,"list is not a Kokkos list\n"); - //typename ArrayTypes::t_neighbors_2d d_neighbors = k_list->d_neighbors; - //typename ArrayTypes::t_int_1d d_ilist = k_list->d_ilist; - //typename ArrayTypes::t_int_1d d_numneigh = k_list->d_numneigh; d_neighbors = k_list->d_neighbors; d_ilist = k_list->d_ilist; d_numneigh = k_list->d_numneigh; @@ -2067,7 +1868,6 @@ void FixRxKokkos::computeLocalTemperature() { // Create an atomic view of sumWeights and dpdThetaLocal. Only needed // for Half/thread scenarios. - //typedef Kokkos::View< E_FLOAT*, typename DAT::t_efloat_1d::array_layout, typename KKDevice::value, Kokkos::MemoryTraits< AtomicF< NEIGHFLAG >::value> > AtomicViewType; typedef Kokkos::View< E_FLOAT*, typename DAT::t_efloat_1d::array_layout, typename KKDevice::value, Kokkos::MemoryTraits< AtomicF< NEIGHFLAG >::value> > AtomicViewType; AtomicViewType a_dpdThetaLocal = d_dpdThetaLocal; @@ -2174,8 +1974,6 @@ void FixRxKokkos::computeLocalTemperature() template int FixRxKokkos::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) { - //printf("inside FixRxKokkos::pack_forward_comm %d\n", comm->me); - HAT::t_float_2d h_dvector = atomKK->k_dvector.h_view; int m = 0; @@ -2186,9 +1984,6 @@ int FixRxKokkos::pack_forward_comm(int n, int *list, double *buf, in buf[m++] = h_dvector(ispecies+nspecies,jj); } } - - //printf("done with FixRxKokkos::pack_forward_comm %d\n", comm->me); - return m; } @@ -2197,8 +1992,6 @@ int FixRxKokkos::pack_forward_comm(int n, int *list, double *buf, in template void FixRxKokkos::unpack_forward_comm(int n, int first, double *buf) { - //printf("inside FixRxKokkos::unpack_forward_comm %d\n", comm->me); - HAT::t_float_2d h_dvector = atomKK->k_dvector.h_view; const int last = first + n ; @@ -2209,8 +2002,6 @@ void FixRxKokkos::unpack_forward_comm(int n, int first, double *buf) h_dvector(ispecies+nspecies,ii) = buf[m++]; } } - - //printf("done with FixRxKokkos::unpack_forward_comm %d\n", comm->me); } /* ---------------------------------------------------------------------- */ @@ -2218,7 +2009,6 @@ void FixRxKokkos::unpack_forward_comm(int n, int first, double *buf) template int FixRxKokkos::pack_reverse_comm(int n, int first, double *buf) { - //printf("inside FixRxKokkos::pack_reverse_comm %d %d %d\n", comm->me, first, n); // Sync the host view. k_dpdThetaLocal.template sync(); k_sumWeights. template sync(); @@ -2230,8 +2020,6 @@ int FixRxKokkos::pack_reverse_comm(int n, int first, double *buf) buf[m++] = h_dpdThetaLocal(i); buf[m++] = h_sumWeights(i); } - //printf("done with FixRxKokkos::pack_reverse_comm %d\n", comm->me); - return m; } @@ -2240,7 +2028,6 @@ int FixRxKokkos::pack_reverse_comm(int n, int first, double *buf) template void FixRxKokkos::unpack_reverse_comm(int n, int *list, double *buf) { - // printf("inside FixRxKokkos::unpack_reverse_comm %d\n", comm->me); int m = 0; for (int i = 0; i < n; i++) { const int j = list[i]; @@ -2252,8 +2039,6 @@ void FixRxKokkos::unpack_reverse_comm(int n, int *list, double *buf) // Signal that the host view has been modified. k_dpdThetaLocal.template modify(); k_sumWeights. template modify(); - - // printf("done with FixRxKokkos::unpack_reverse_comm %d\n", comm->me); } namespace LAMMPS_NS { diff --git a/src/fix_property_atom.h b/src/fix_property_atom.h index f2f26d511b..92497d6188 100644 --- a/src/fix_property_atom.h +++ b/src/fix_property_atom.h @@ -56,10 +56,10 @@ class FixPropertyAtom : public Fix { int nvalue, border; int molecule_flag, q_flag, rmass_flag; // flags for specific fields int temperature_flag, heatflow_flag; - int *styles; // style of each value, see enum - int *index; // indices into atom custom data structs - int *cols; // columns per value, for arrays - char *astyle; // atom style at instantiation + int *styles; // style of each value, see enum + int *index; // indices into atom custom data structs + int *cols; // columns per value, for arrays + char *astyle; // atom style at instantiation int values_peratom; // # of values per atom, including multiple for arrays int nmax_old; // length of peratom arrays the last time they grew