From c2f4c8d23a623afb85bbcf66319f352186659b99 Mon Sep 17 00:00:00 2001 From: mkanski <20713012+mkanski@users.noreply.github.com> Date: Thu, 10 Nov 2022 14:13:44 +0100 Subject: [PATCH 01/12] First version of KOKKOS fix viscous --- src/KOKKOS/fix_viscous_kokkos.cpp | 104 ++++++++++++++++++++++++++++++ src/KOKKOS/fix_viscous_kokkos.h | 61 ++++++++++++++++++ src/fix_viscous.cpp | 2 + 3 files changed, 167 insertions(+) create mode 100644 src/KOKKOS/fix_viscous_kokkos.cpp create mode 100644 src/KOKKOS/fix_viscous_kokkos.h diff --git a/src/KOKKOS/fix_viscous_kokkos.cpp b/src/KOKKOS/fix_viscous_kokkos.cpp new file mode 100644 index 0000000000..c66553fb13 --- /dev/null +++ b/src/KOKKOS/fix_viscous_kokkos.cpp @@ -0,0 +1,104 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "fix_viscous_kokkos.h" + +#include "atom_kokkos.h" +#include "update.h" +#include "modify.h" +#include "domain.h" +#include "input.h" +#include "variable.h" +#include "memory_kokkos.h" +#include "error.h" +#include "atom_masks.h" +#include "kokkos_base.h" + +#include + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +template +FixViscousKokkos::FixViscousKokkos(LAMMPS *lmp, int narg, char **arg) : + FixViscous(lmp, narg, arg) +{ + kokkosable = 1; + atomKK = (AtomKokkos *) atom; + execution_space = ExecutionSpaceFromDevice::space; + datamask_read = EMPTY_MASK; + datamask_modify = EMPTY_MASK; +} + +/* ---------------------------------------------------------------------- */ + +template +FixViscousKokkos::~FixViscousKokkos() +{ + if (copymode) return; +} + +/* ---------------------------------------------------------------------- */ + +template +void FixViscousKokkos::init() +{ + FixViscous::init(); + + if (utils::strmatch(update->integrate_style,"^respa")) + error->all(FLERR,"Cannot (yet) use respa with Kokkos"); +} + +/* ---------------------------------------------------------------------- */ + +template +void FixViscousKokkos::post_force(int /*vflag*/) +{ + atomKK->sync(execution_space, V_MASK | F_MASK | MASK_MASK); + + v = atomKK->k_v.view(); + f = atomKK->k_f.view(); + mask = atomKK->k_mask.view(); + type = atomKK->k_type.view(); + + int nlocal = atom->nlocal; + + copymode = 1; + Kokkos::parallel_for(Kokkos::RangePolicy(0,nlocal),*this); + copymode = 0; + + atomKK->modified(execution_space, F_MASK); + +} + +template +KOKKOS_INLINE_FUNCTION +void FixViscousKokkos::operator()(TagFixViscous, const int &i) const { + if (mask[i] & groupbit) { + double drag = gamma[type[i]]; + f(i,0) -= drag*v(i,0); + f(i,1) -= drag*v(i,1); + f(i,2) -= drag*v(i,2); + } +} + +namespace LAMMPS_NS { +template class FixViscousKokkos; +#ifdef LMP_KOKKOS_GPU +template class FixViscousKokkos; +#endif +} + diff --git a/src/KOKKOS/fix_viscous_kokkos.h b/src/KOKKOS/fix_viscous_kokkos.h new file mode 100644 index 0000000000..f08c2a5fd6 --- /dev/null +++ b/src/KOKKOS/fix_viscous_kokkos.h @@ -0,0 +1,61 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS +// clang-format off +FixStyle(viscous/kk,FixViscousKokkos); +FixStyle(viscous/kk/device,FixViscousKokkos); +FixStyle(viscous/kk/host,FixViscousKokkos); +// clang-format on +#else + +// clang-format off +#ifndef LMP_FIX_VISCOUS_KOKKOS_H +#define LMP_FIX_VISCOUS_KOKKOS_H + +#include "fix_viscous.h" +#include "kokkos_type.h" + +namespace LAMMPS_NS { + +struct TagFixViscous{}; + +struct TagFixViscousNonConstant{}; + +template +class FixViscousKokkos : public FixViscous { + public: + typedef DeviceType device_type; +// typedef double_3 value_type; + typedef ArrayTypes AT; + + FixViscousKokkos(class LAMMPS *, int, char **); + ~FixViscousKokkos() override; + void init() override; + void post_force(int) override; + + KOKKOS_INLINE_FUNCTION + void operator()(TagFixViscous, const int&) const; + + private: + typename AT::t_v_array_randomread v; + typename AT::t_f_array f; + typename AT::t_int_1d_randomread mask; + typename AT::t_int_1d_randomread type; +}; + +} + +#endif +#endif + diff --git a/src/fix_viscous.cpp b/src/fix_viscous.cpp index 740dbd91fb..fb0b38f264 100644 --- a/src/fix_viscous.cpp +++ b/src/fix_viscous.cpp @@ -61,6 +61,8 @@ FixViscous::FixViscous(LAMMPS *lmp, int narg, char **arg) : FixViscous::~FixViscous() { + if (copymode) return; + delete [] gamma; } From ee5d40984fc3316de8f89c441eec9d17110502e2 Mon Sep 17 00:00:00 2001 From: mkanski <20713012+mkanski@users.noreply.github.com> Date: Thu, 10 Nov 2022 15:43:08 +0100 Subject: [PATCH 02/12] Change the way gamma is accessed --- src/KOKKOS/fix_viscous_kokkos.cpp | 7 ++++++- src/KOKKOS/fix_viscous_kokkos.h | 4 +++- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/src/KOKKOS/fix_viscous_kokkos.cpp b/src/KOKKOS/fix_viscous_kokkos.cpp index c66553fb13..7d3d5f042a 100644 --- a/src/KOKKOS/fix_viscous_kokkos.cpp +++ b/src/KOKKOS/fix_viscous_kokkos.cpp @@ -49,6 +49,7 @@ template FixViscousKokkos::~FixViscousKokkos() { if (copymode) return; + } /* ---------------------------------------------------------------------- */ @@ -58,6 +59,10 @@ void FixViscousKokkos::init() { FixViscous::init(); + k_gamma = Kokkos::DualView("FixViscousKokkos:gamma",atom->ntypes+1); + + for (int i = 1; i <= atom->ntypes; i++) k_gamma.h_view(i) = gamma[i]; + if (utils::strmatch(update->integrate_style,"^respa")) error->all(FLERR,"Cannot (yet) use respa with Kokkos"); } @@ -88,7 +93,7 @@ template KOKKOS_INLINE_FUNCTION void FixViscousKokkos::operator()(TagFixViscous, const int &i) const { if (mask[i] & groupbit) { - double drag = gamma[type[i]]; + double drag = k_gamma.h_view(type[i]); f(i,0) -= drag*v(i,0); f(i,1) -= drag*v(i,1); f(i,2) -= drag*v(i,2); diff --git a/src/KOKKOS/fix_viscous_kokkos.h b/src/KOKKOS/fix_viscous_kokkos.h index f08c2a5fd6..904e3bd6c4 100644 --- a/src/KOKKOS/fix_viscous_kokkos.h +++ b/src/KOKKOS/fix_viscous_kokkos.h @@ -48,10 +48,12 @@ class FixViscousKokkos : public FixViscous { void operator()(TagFixViscous, const int&) const; private: - typename AT::t_v_array_randomread v; + typename AT::t_v_array v; typename AT::t_f_array f; typename AT::t_int_1d_randomread mask; typename AT::t_int_1d_randomread type; + + Kokkos::DualView k_gamma; }; } From 9a78f45c09ec083797e7ae3898de6d531065ca2a Mon Sep 17 00:00:00 2001 From: mkanski <20713012+mkanski@users.noreply.github.com> Date: Thu, 10 Nov 2022 17:03:55 +0100 Subject: [PATCH 03/12] Second try for CUDA --- src/KOKKOS/fix_viscous_kokkos.cpp | 6 +++++- src/KOKKOS/fix_viscous_kokkos.h | 1 + 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/KOKKOS/fix_viscous_kokkos.cpp b/src/KOKKOS/fix_viscous_kokkos.cpp index 7d3d5f042a..6fb2bd9d6a 100644 --- a/src/KOKKOS/fix_viscous_kokkos.cpp +++ b/src/KOKKOS/fix_viscous_kokkos.cpp @@ -61,8 +61,12 @@ void FixViscousKokkos::init() k_gamma = Kokkos::DualView("FixViscousKokkos:gamma",atom->ntypes+1); + gamma2 = k_gamma.template view(); + for (int i = 1; i <= atom->ntypes; i++) k_gamma.h_view(i) = gamma[i]; + k_gamma.template sync(); + if (utils::strmatch(update->integrate_style,"^respa")) error->all(FLERR,"Cannot (yet) use respa with Kokkos"); } @@ -93,7 +97,7 @@ template KOKKOS_INLINE_FUNCTION void FixViscousKokkos::operator()(TagFixViscous, const int &i) const { if (mask[i] & groupbit) { - double drag = k_gamma.h_view(type[i]); + double drag = gamma2(type[i]); f(i,0) -= drag*v(i,0); f(i,1) -= drag*v(i,1); f(i,2) -= drag*v(i,2); diff --git a/src/KOKKOS/fix_viscous_kokkos.h b/src/KOKKOS/fix_viscous_kokkos.h index 904e3bd6c4..0498b3ab77 100644 --- a/src/KOKKOS/fix_viscous_kokkos.h +++ b/src/KOKKOS/fix_viscous_kokkos.h @@ -54,6 +54,7 @@ class FixViscousKokkos : public FixViscous { typename AT::t_int_1d_randomread type; Kokkos::DualView k_gamma; + typename Kokkos::DualView::t_dev_const_um gamma2; }; } From 27bd28bf34842ff506c19d8c0668768ec50960a7 Mon Sep 17 00:00:00 2001 From: mkanski <20713012+mkanski@users.noreply.github.com> Date: Thu, 10 Nov 2022 18:52:28 +0100 Subject: [PATCH 04/12] CUDA should work now --- src/KOKKOS/fix_viscous_kokkos.cpp | 15 ++++----------- src/KOKKOS/fix_viscous_kokkos.h | 4 ---- 2 files changed, 4 insertions(+), 15 deletions(-) diff --git a/src/KOKKOS/fix_viscous_kokkos.cpp b/src/KOKKOS/fix_viscous_kokkos.cpp index 6fb2bd9d6a..80ddff2fce 100644 --- a/src/KOKKOS/fix_viscous_kokkos.cpp +++ b/src/KOKKOS/fix_viscous_kokkos.cpp @@ -17,16 +17,12 @@ #include "atom_kokkos.h" #include "update.h" #include "modify.h" -#include "domain.h" #include "input.h" -#include "variable.h" #include "memory_kokkos.h" #include "error.h" #include "atom_masks.h" #include "kokkos_base.h" -#include - using namespace LAMMPS_NS; using namespace FixConst; @@ -49,7 +45,6 @@ template FixViscousKokkos::~FixViscousKokkos() { if (copymode) return; - } /* ---------------------------------------------------------------------- */ @@ -61,10 +56,9 @@ void FixViscousKokkos::init() k_gamma = Kokkos::DualView("FixViscousKokkos:gamma",atom->ntypes+1); - gamma2 = k_gamma.template view(); - - for (int i = 1; i <= atom->ntypes; i++) k_gamma.h_view(i) = gamma[i]; + for (int i = 1; i <= atom->ntypes; i++) k_gamma.h_view(i) = gamma[i]; + k_gamma.template modify(); k_gamma.template sync(); if (utils::strmatch(update->integrate_style,"^respa")) @@ -76,7 +70,7 @@ void FixViscousKokkos::init() template void FixViscousKokkos::post_force(int /*vflag*/) { - atomKK->sync(execution_space, V_MASK | F_MASK | MASK_MASK); + atomKK->sync(execution_space, V_MASK | F_MASK | MASK_MASK | TYPE_MASK); v = atomKK->k_v.view(); f = atomKK->k_f.view(); @@ -90,14 +84,13 @@ void FixViscousKokkos::post_force(int /*vflag*/) copymode = 0; atomKK->modified(execution_space, F_MASK); - } template KOKKOS_INLINE_FUNCTION void FixViscousKokkos::operator()(TagFixViscous, const int &i) const { if (mask[i] & groupbit) { - double drag = gamma2(type[i]); + double drag = k_gamma.d_view(type[i]); f(i,0) -= drag*v(i,0); f(i,1) -= drag*v(i,1); f(i,2) -= drag*v(i,2); diff --git a/src/KOKKOS/fix_viscous_kokkos.h b/src/KOKKOS/fix_viscous_kokkos.h index 0498b3ab77..2594c0cf85 100644 --- a/src/KOKKOS/fix_viscous_kokkos.h +++ b/src/KOKKOS/fix_viscous_kokkos.h @@ -30,13 +30,10 @@ namespace LAMMPS_NS { struct TagFixViscous{}; -struct TagFixViscousNonConstant{}; - template class FixViscousKokkos : public FixViscous { public: typedef DeviceType device_type; -// typedef double_3 value_type; typedef ArrayTypes AT; FixViscousKokkos(class LAMMPS *, int, char **); @@ -54,7 +51,6 @@ class FixViscousKokkos : public FixViscous { typename AT::t_int_1d_randomread type; Kokkos::DualView k_gamma; - typename Kokkos::DualView::t_dev_const_um gamma2; }; } From fabfd863388a73b8b246709328bbdec7af29f799 Mon Sep 17 00:00:00 2001 From: mkanski <20713012+mkanski@users.noreply.github.com> Date: Thu, 10 Nov 2022 20:23:55 +0100 Subject: [PATCH 05/12] Base for KOKKOS dt/reset --- src/KOKKOS/fix_dt_reset_kokkos.cpp | 136 +++++++++++++++++++++++++++++ src/KOKKOS/fix_dt_reset_kokkos.h | 63 +++++++++++++ src/fix_dt_reset.h | 2 +- 3 files changed, 200 insertions(+), 1 deletion(-) create mode 100644 src/KOKKOS/fix_dt_reset_kokkos.cpp create mode 100644 src/KOKKOS/fix_dt_reset_kokkos.h diff --git a/src/KOKKOS/fix_dt_reset_kokkos.cpp b/src/KOKKOS/fix_dt_reset_kokkos.cpp new file mode 100644 index 0000000000..c8d37e7341 --- /dev/null +++ b/src/KOKKOS/fix_dt_reset_kokkos.cpp @@ -0,0 +1,136 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "fix_dt_reset_kokkos.h" + +#include "atom_kokkos.h" +#include "update.h" +#include "modify.h" +#include "input.h" +#include "memory_kokkos.h" +#include "error.h" +#include "atom_masks.h" +#include "kokkos_base.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +#define BIG 1.0e20 + +/* ---------------------------------------------------------------------- */ + +template +FixDtResetKokkos::FixDtResetKokkos(LAMMPS *lmp, int narg, char **arg) : + FixDtReset(lmp, narg, arg) +{ + kokkosable = 1; + atomKK = (AtomKokkos *) atom; + execution_space = ExecutionSpaceFromDevice::space; + datamask_read = EMPTY_MASK; + datamask_modify = EMPTY_MASK; +} + +/* ---------------------------------------------------------------------- */ + +template +void FixDtResetKokkos::init() +{ + FixDtReset::init(); + + k_params = Kokkos::DualView("FixDtResetKokkos:gamma", 1); + + k_params.h_view(0) = emax; + + + k_params.template modify(); + k_params.template sync(); + + if (utils::strmatch(update->integrate_style,"^respa")) + error->all(FLERR,"Cannot (yet) use respa with Kokkos"); +} + +/* ---------------------------------------------------------------------- */ + +template +void FixDtResetKokkos::end_of_step() +{ + atomKK->sync(execution_space, V_MASK | F_MASK | MASK_MASK | TYPE_MASK | RMASS_MASK); + + v = atomKK->k_v.view(); + f = atomKK->k_f.view(); + mask = atomKK->k_mask.view(); + type = atomKK->k_type.view(); + if (atomKK->rmass) + rmass = atomKK->k_rmass.view(); + else + mass = atomKK->k_mass.view(); + + int nlocal = atom->nlocal; + + double dtmin = BIG; + +// Kokkos::DualView dt = +// Kokkos::DualView("FixViscousKokkos:gamma", 1); + + double dt; + + copymode = 1; + Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal), *this, Kokkos::Min(dt)); +// Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal), *this, dt); + copymode = 0; + + printf ("Dt: %f \n", dt); + atomKK->modified(execution_space, F_MASK); +} + +template +KOKKOS_INLINE_FUNCTION +void FixDtResetKokkos::operator()(TagFixDtResetMass, const int &i, double &k_dt) const { + + double dtv, dtf, dte, dtsq; + double vsq, fsq, massinv; + double delx, dely, delz, delr; + + double k_emax = k_params.d_view(0); + + if (mask[i] & groupbit) { + + massinv = 1.0 / mass[type[i]]; + vsq = v(i,0) * v(i,0) + v(i,1) * v(i,1) + v(i,2) * v(i,2); + fsq = f(i,0) * f(i,0) + f(i,1) * f(i,1) + f(i,2) * f(i,2); + dtv = dtf = dte = BIG; + if (vsq > 0.0) dtv = xmax / sqrt(vsq); + if (fsq > 0.0) dtf = sqrt(2.0 * xmax / (ftm2v * sqrt(fsq) * massinv)); + k_dt = MIN(dtv, dtf); + if ((k_emax > 0.0) && (fsq * vsq > 0.0)) { + dte = k_emax / sqrt(fsq * vsq) / sqrt(ftm2v * mvv2e); + k_dt = MIN(dt, dte); + } + dtsq = k_dt * k_dt; + delx = k_dt * v(i,0) + 0.5 * dtsq * massinv * f(i,0) * ftm2v; + dely = k_dt * v(i,1) + 0.5 * dtsq * massinv * f(i,1) * ftm2v; + delz = k_dt * v(i,2) + 0.5 * dtsq * massinv * f(i,2) * ftm2v; + delr = sqrt(delx * delx + dely * dely + delz * delz); + if (delr > xmax) k_dt *= xmax / delr; + } + + } + +namespace LAMMPS_NS { +template class FixDtResetKokkos; +#ifdef LMP_KOKKOS_GPU +template class FixDtResetKokkos; +#endif +} + diff --git a/src/KOKKOS/fix_dt_reset_kokkos.h b/src/KOKKOS/fix_dt_reset_kokkos.h new file mode 100644 index 0000000000..cef6ac85e9 --- /dev/null +++ b/src/KOKKOS/fix_dt_reset_kokkos.h @@ -0,0 +1,63 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS +// clang-format off +FixStyle(dt/reset/kk,FixDtResetKokkos); +FixStyle(dt/reset/kk/device,FixDtResetKokkos); +FixStyle(dt/reset/kk/host,FixDtResetKokkos); +// clang-format on +#else + +// clang-format off +#ifndef LMP_FIX_DT_RESET_KOKKOS_H +#define LMP_FIX_DT_RESET_KOKKOS_H + +#include "fix_dt_reset.h" +#include "kokkos_type.h" + +namespace LAMMPS_NS { + +struct TagFixDtResetMass{}; + +template +class FixDtResetKokkos : public FixDtReset { + public: + typedef DeviceType device_type; + typedef ArrayTypes AT; + + FixDtResetKokkos(class LAMMPS *, int, char **); +// ~FixDtResetKokkos() override; + void init() override; + void end_of_step() override; + + KOKKOS_INLINE_FUNCTION + void operator()(TagFixDtResetMass, const int&, double&) const; + + private: + typename AT::t_v_array v; + typename AT::t_f_array f; + typename AT::t_int_1d_randomread mask; + typename AT::t_int_1d_randomread type; + typename ArrayTypes::t_float_1d_randomread rmass; + typename ArrayTypes::t_float_1d_randomread mass; + + + Kokkos::DualView k_params; +}; + +} + +#endif +#endif + diff --git a/src/fix_dt_reset.h b/src/fix_dt_reset.h index f0a0564685..00070dac4b 100644 --- a/src/fix_dt_reset.h +++ b/src/fix_dt_reset.h @@ -34,7 +34,7 @@ class FixDtReset : public Fix { void end_of_step() override; double compute_scalar() override; - private: + protected: bigint laststep; int minbound, maxbound; double tmin, tmax, xmax, emax; From 5dbc41c16837ed3522a59865ebddebca04304aef Mon Sep 17 00:00:00 2001 From: mkanski <20713012+mkanski@users.noreply.github.com> Date: Thu, 10 Nov 2022 22:33:19 +0100 Subject: [PATCH 06/12] First working version (rmass not tested yet) --- src/KOKKOS/fix_dt_reset_kokkos.cpp | 95 +++++++++++++++++++++++------- src/KOKKOS/fix_dt_reset_kokkos.h | 5 +- 2 files changed, 79 insertions(+), 21 deletions(-) diff --git a/src/KOKKOS/fix_dt_reset_kokkos.cpp b/src/KOKKOS/fix_dt_reset_kokkos.cpp index c8d37e7341..d27fa676bc 100644 --- a/src/KOKKOS/fix_dt_reset_kokkos.cpp +++ b/src/KOKKOS/fix_dt_reset_kokkos.cpp @@ -15,13 +15,17 @@ #include "fix_dt_reset_kokkos.h" #include "atom_kokkos.h" -#include "update.h" -#include "modify.h" -#include "input.h" -#include "memory_kokkos.h" -#include "error.h" #include "atom_masks.h" +#include "error.h" +#include "force.h" +#include "input.h" +#include "integrate.h" #include "kokkos_base.h" +#include "memory_kokkos.h" +#include "modify.h" +#include "output.h" +#include "pair.h" +#include "update.h" using namespace LAMMPS_NS; using namespace FixConst; @@ -48,13 +52,13 @@ void FixDtResetKokkos::init() { FixDtReset::init(); - k_params = Kokkos::DualView("FixDtResetKokkos:gamma", 1); + k_emax = Kokkos::DualView("FixDtResetKokkos:gamma", 1); - k_params.h_view(0) = emax; + k_emax.h_view(0) = emax; - k_params.template modify(); - k_params.template sync(); + k_emax.template modify(); + k_emax.template sync(); if (utils::strmatch(update->integrate_style,"^respa")) error->all(FLERR,"Cannot (yet) use respa with Kokkos"); @@ -78,20 +82,38 @@ void FixDtResetKokkos::end_of_step() int nlocal = atom->nlocal; - double dtmin = BIG; - -// Kokkos::DualView dt = -// Kokkos::DualView("FixViscousKokkos:gamma", 1); - double dt; copymode = 1; - Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal), *this, Kokkos::Min(dt)); -// Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal), *this, dt); + if (atomKK->rmass) + Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal), *this, Kokkos::Min(dt)); + else + Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal), *this, Kokkos::Min(dt)); copymode = 0; - printf ("Dt: %f \n", dt); atomKK->modified(execution_space, F_MASK); + + if (minbound) dt = MAX(dt, tmin); + if (maxbound) dt = MIN(dt, tmax); + + // if timestep didn't change, just return + // else reset update->dt and other classes that depend on it + // rRESPA, pair style, fixes + + if (dt == update->dt) return; + + laststep = update->ntimestep; + + // calls to other classes that need to know timestep size changed + // similar logic is in Input::timestep() + + update->update_time(); + update->dt = dt; + update->dt_default = 0; + if (force->pair) force->pair->reset_dt(); + for (int i = 0; i < modify->nfix; i++) modify->fix[i]->reset_dt(); + output->reset_dt(); + } template @@ -102,7 +124,7 @@ void FixDtResetKokkos::operator()(TagFixDtResetMass, const int &i, d double vsq, fsq, massinv; double delx, dely, delz, delr; - double k_emax = k_params.d_view(0); + double emax = k_emax.d_view(0); if (mask[i] & groupbit) { @@ -113,8 +135,41 @@ void FixDtResetKokkos::operator()(TagFixDtResetMass, const int &i, d if (vsq > 0.0) dtv = xmax / sqrt(vsq); if (fsq > 0.0) dtf = sqrt(2.0 * xmax / (ftm2v * sqrt(fsq) * massinv)); k_dt = MIN(dtv, dtf); - if ((k_emax > 0.0) && (fsq * vsq > 0.0)) { - dte = k_emax / sqrt(fsq * vsq) / sqrt(ftm2v * mvv2e); + if ((emax > 0.0) && (fsq * vsq > 0.0)) { + dte = emax / sqrt(fsq * vsq) / sqrt(ftm2v * mvv2e); + k_dt = MIN(dt, dte); + } + dtsq = k_dt * k_dt; + delx = k_dt * v(i,0) + 0.5 * dtsq * massinv * f(i,0) * ftm2v; + dely = k_dt * v(i,1) + 0.5 * dtsq * massinv * f(i,1) * ftm2v; + delz = k_dt * v(i,2) + 0.5 * dtsq * massinv * f(i,2) * ftm2v; + delr = sqrt(delx * delx + dely * dely + delz * delz); + if (delr > xmax) k_dt *= xmax / delr; + } + + } + +template +KOKKOS_INLINE_FUNCTION +void FixDtResetKokkos::operator()(TagFixDtResetRMass, const int &i, double &k_dt) const { + + double dtv, dtf, dte, dtsq; + double vsq, fsq, massinv; + double delx, dely, delz, delr; + + double emax = k_emax.d_view(0); + + if (mask[i] & groupbit) { + + massinv = 1.0 / rmass[i]; + vsq = v(i,0) * v(i,0) + v(i,1) * v(i,1) + v(i,2) * v(i,2); + fsq = f(i,0) * f(i,0) + f(i,1) * f(i,1) + f(i,2) * f(i,2); + dtv = dtf = dte = BIG; + if (vsq > 0.0) dtv = xmax / sqrt(vsq); + if (fsq > 0.0) dtf = sqrt(2.0 * xmax / (ftm2v * sqrt(fsq) * massinv)); + k_dt = MIN(dtv, dtf); + if ((emax > 0.0) && (fsq * vsq > 0.0)) { + dte = emax / sqrt(fsq * vsq) / sqrt(ftm2v * mvv2e); k_dt = MIN(dt, dte); } dtsq = k_dt * k_dt; diff --git a/src/KOKKOS/fix_dt_reset_kokkos.h b/src/KOKKOS/fix_dt_reset_kokkos.h index cef6ac85e9..a269b58137 100644 --- a/src/KOKKOS/fix_dt_reset_kokkos.h +++ b/src/KOKKOS/fix_dt_reset_kokkos.h @@ -29,6 +29,7 @@ FixStyle(dt/reset/kk/host,FixDtResetKokkos); namespace LAMMPS_NS { struct TagFixDtResetMass{}; +struct TagFixDtResetRMass{}; template class FixDtResetKokkos : public FixDtReset { @@ -43,6 +44,8 @@ class FixDtResetKokkos : public FixDtReset { KOKKOS_INLINE_FUNCTION void operator()(TagFixDtResetMass, const int&, double&) const; + KOKKOS_INLINE_FUNCTION + void operator()(TagFixDtResetRMass, const int&, double&) const; private: typename AT::t_v_array v; @@ -53,7 +56,7 @@ class FixDtResetKokkos : public FixDtReset { typename ArrayTypes::t_float_1d_randomread mass; - Kokkos::DualView k_params; + Kokkos::DualView k_emax; }; } From 4647096a97b7e3ff694ec4115efcdb60035b5ec1 Mon Sep 17 00:00:00 2001 From: mkanski <20713012+mkanski@users.noreply.github.com> Date: Mon, 14 Nov 2022 13:30:41 +0100 Subject: [PATCH 07/12] Add docs --- doc/src/fix_dt_reset.rst | 6 ++++++ doc/src/fix_viscous.rst | 6 ++++++ 2 files changed, 12 insertions(+) diff --git a/doc/src/fix_dt_reset.rst b/doc/src/fix_dt_reset.rst index 368a3dcd70..e99bf37b75 100644 --- a/doc/src/fix_dt_reset.rst +++ b/doc/src/fix_dt_reset.rst @@ -3,6 +3,8 @@ fix dt/reset command ==================== +Accelerator Variants: *dt/reset/kk* + Syntax """""" @@ -86,6 +88,10 @@ allows dump files to be written at intervals specified by simulation time, rather than by timesteps. Simulation time is in time units; see the :doc:`units ` doc page for details. +---------- + +.. include:: accel_styles.rst + Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/doc/src/fix_viscous.rst b/doc/src/fix_viscous.rst index 6c48cb5f12..5e68ec0a8e 100644 --- a/doc/src/fix_viscous.rst +++ b/doc/src/fix_viscous.rst @@ -3,6 +3,8 @@ fix viscous command =================== +Accelerator Variants: *setforce/kk* + Syntax """""" @@ -84,6 +86,10 @@ more easily be used as a thermostat. ---------- +.. include:: accel_styles.rst + +---------- + Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" From a7de83d2899c3c12bd0dd5af21e933282a8ef985 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 14 Nov 2022 08:24:03 -0700 Subject: [PATCH 08/12] Update Kokkos Install.sh and fix typo in docs --- doc/src/fix_viscous.rst | 2 +- src/KOKKOS/Install.sh | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/doc/src/fix_viscous.rst b/doc/src/fix_viscous.rst index 5e68ec0a8e..723bc15085 100644 --- a/doc/src/fix_viscous.rst +++ b/doc/src/fix_viscous.rst @@ -3,7 +3,7 @@ fix viscous command =================== -Accelerator Variants: *setforce/kk* +Accelerator Variants: *viscous/kk* Syntax """""" diff --git a/src/KOKKOS/Install.sh b/src/KOKKOS/Install.sh index 063ea44720..53d2093976 100755 --- a/src/KOKKOS/Install.sh +++ b/src/KOKKOS/Install.sh @@ -117,6 +117,8 @@ action fix_acks2_reaxff_kokkos.cpp fix_acks2_reaxff.cpp action fix_acks2_reaxff_kokkos.h fix_acks2_reaxff.h action fix_deform_kokkos.cpp action fix_deform_kokkos.h +action fix_dt_reset_kokkos.cpp +action fix_dt_reset_kokkos.h action fix_enforce2d_kokkos.cpp action fix_enforce2d_kokkos.h action fix_eos_table_rx_kokkos.cpp fix_eos_table_rx.cpp @@ -165,6 +167,8 @@ action fix_wall_lj93_kokkos.cpp action fix_wall_lj93_kokkos.h action fix_wall_reflect_kokkos.cpp action fix_wall_reflect_kokkos.h +action fix_viscous_kokkos.cpp +action fix_viscous_kokkos.h action fix_dpd_energy_kokkos.cpp fix_dpd_energy.cpp action fix_dpd_energy_kokkos.h fix_dpd_energy.h action fix_rx_kokkos.cpp fix_rx.cpp From 859403e2af55aeb08c5ada361817c1c292732b9f Mon Sep 17 00:00:00 2001 From: mkanski <20713012+mkanski@users.noreply.github.com> Date: Tue, 15 Nov 2022 12:01:29 +0100 Subject: [PATCH 09/12] dt/reset now works in parallel --- src/KOKKOS/fix_dt_reset_kokkos.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/KOKKOS/fix_dt_reset_kokkos.cpp b/src/KOKKOS/fix_dt_reset_kokkos.cpp index d27fa676bc..6ffa7fcf56 100644 --- a/src/KOKKOS/fix_dt_reset_kokkos.cpp +++ b/src/KOKKOS/fix_dt_reset_kokkos.cpp @@ -91,6 +91,8 @@ void FixDtResetKokkos::end_of_step() Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal), *this, Kokkos::Min(dt)); copymode = 0; + MPI_Allreduce(MPI_IN_PLACE, &dt, 1, MPI_DOUBLE, MPI_MIN, world); + atomKK->modified(execution_space, F_MASK); if (minbound) dt = MAX(dt, tmin); From 2de997b52dacb5af2b888ea5fa4c6121c5c97e70 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 16 Dec 2022 01:11:52 -0500 Subject: [PATCH 10/12] import tabulate scripts for table files --- tools/tabulate/README.md | 39 +++ tools/tabulate/angle_harmonic_tabulate.py | 21 ++ tools/tabulate/bond_morse_tabulate.py | 28 ++ tools/tabulate/dihedral_harmonic_tabulate.py | 19 ++ tools/tabulate/pair_hybrid_tabulate.py | 31 +++ tools/tabulate/pair_lj_tabulate.py | 26 ++ tools/tabulate/tabulate.py | 279 +++++++++++++++++++ 7 files changed, 443 insertions(+) create mode 100644 tools/tabulate/README.md create mode 100755 tools/tabulate/angle_harmonic_tabulate.py create mode 100755 tools/tabulate/bond_morse_tabulate.py create mode 100755 tools/tabulate/dihedral_harmonic_tabulate.py create mode 100755 tools/tabulate/pair_hybrid_tabulate.py create mode 100755 tools/tabulate/pair_lj_tabulate.py create mode 100755 tools/tabulate/tabulate.py diff --git a/tools/tabulate/README.md b/tools/tabulate/README.md new file mode 100644 index 0000000000..cf6a92d166 --- /dev/null +++ b/tools/tabulate/README.md @@ -0,0 +1,39 @@ +# Python scripts to generate tabulated potential files for LAMMPS + +This directory contains a Python module 'tabulate' that can be used to +generate tabulated potential files for pair style table, bond style +table, and angle style table + +To create tables, you need to define your energy and - optionally - +force functions and then an instance of either the PairTabulate(), +BondTabulate(), AngleTabulate(), or DihedralTabulate() class from the +tabulate module and call its run() method to generate the table. Most +of the settings (number of points, minimum, maximum etc.) are provided +as command line flags. The run() method may be called multiple times to +generate multiple tables, for instance after changing parameters of the +energy/force functions. If the force function is not provided, the +force will be determined through numerical differentiation. + +Please see the individual tabulation scripts in this folder for examples: + +| ------------------------------|----------------------------------------------------------------------------| +| File | Description | +| ------------------------------|----------------------------------------------------------------------------| +| pair_lj_tabulate.py | creates two Lennard-Jones pair potential tables with different parameters | +| bond_morse_tabulate.py | creates a table for a Morse bond potential table | +| angle_harmonic_tabulate.py | creates a table for a harmonic angle potential table | +| dihedral_harmonic_tabulate.py | creates a table for a harmonic dihedral potential table | +| pair_hybrid_tabulate.py | creates a Morse/Lennard-Jones hybrid potential table with smooth switching | +| ------------------------------|----------------------------------------------------------------------------| + +Common command line flags: + +``` +options: + -h, --help show this help message and exit + --num-points NUM, -n NUM Number of tabulated points (default: 1000) + --filename FILENAME, -f FILENAME Name of output file (default: -) + --diff-num, -d Differentiate energy function numerically + --inner XMIN, -i XMIN Inner cutoff of table (required for pair) + --outer XMAX, -o XMAX Outer cutoff of table (required) +``` diff --git a/tools/tabulate/angle_harmonic_tabulate.py b/tools/tabulate/angle_harmonic_tabulate.py new file mode 100755 index 0000000000..665b392198 --- /dev/null +++ b/tools/tabulate/angle_harmonic_tabulate.py @@ -0,0 +1,21 @@ +#!/usr/bin/env python3 + +from tabulate import AngleTabulate + +################################################################################ +import math + +def harmonic_energy(theta): + k = 100.0 + thetazero = 120.0 + # the force constant in LAMMPS is in energy per radians^2 so convert from degrees to radians + deg2rad = math.pi / 180.0 + t = (theta - thetazero) * deg2rad + f = k * t * t + return f + +################################################################################ + +if __name__ == "__main__": + atable = AngleTabulate(harmonic_energy, units='real') + atable.run('HARM') diff --git a/tools/tabulate/bond_morse_tabulate.py b/tools/tabulate/bond_morse_tabulate.py new file mode 100755 index 0000000000..2301343d4e --- /dev/null +++ b/tools/tabulate/bond_morse_tabulate.py @@ -0,0 +1,28 @@ +#!/usr/bin/env python3 + +from tabulate import BondTabulate + +################################################################################ +import math + +def morse_energy(r): + depth = 1.0 + width = 2.0 + rzero = 1.2 + ralpha = math.exp(-width*(r-rzero)) + f = depth * (-1.0 + (1.0 - ralpha) * (1.0 - ralpha)) + return f + +def morse_force(r): + depth = 1.0 + width = 2.0 + rzero = 1.2 + ralpha = math.exp(-width*(r-rzero)) + f = -2.0 * depth * width * (1.0 -ralpha) * ralpha + return f + +################################################################################ + +if __name__ == "__main__": + btable = BondTabulate(morse_energy, morse_force, units='lj') + btable.run('MORSE') diff --git a/tools/tabulate/dihedral_harmonic_tabulate.py b/tools/tabulate/dihedral_harmonic_tabulate.py new file mode 100755 index 0000000000..085cc73c3b --- /dev/null +++ b/tools/tabulate/dihedral_harmonic_tabulate.py @@ -0,0 +1,19 @@ +#!/usr/bin/env python3 + +from tabulate import DihedralTabulate + +################################################################################ +import math + +def harmonic_energy(theta): + k = 100.0 + # the force constant in LAMMPS is in energy per radians^2 so convert from degrees to radians + deg2rad = math.pi / 180.0 + f = k * (1 - math.cos(2.0 * deg2rad * theta)) + return f + +################################################################################ + +if __name__ == "__main__": + dtable = DihedralTabulate(harmonic_energy, units='metal') + dtable.run('HARM') diff --git a/tools/tabulate/pair_hybrid_tabulate.py b/tools/tabulate/pair_hybrid_tabulate.py new file mode 100755 index 0000000000..dabc0c7e85 --- /dev/null +++ b/tools/tabulate/pair_hybrid_tabulate.py @@ -0,0 +1,31 @@ +#!/usr/bin/env python3 + +from tabulate import PairTabulate + +################################################################################ +import math + +# hybrid potential using Morse for repulsive and LJ for attractive. +# using tanh() as smooth switching function between the two. + +def hybrid_energy(r): + depth = 1.0 + width = 5.0 + epsilon = 1.0 + sigma = 1.0 + rzero = 1.2 + + f1 = 4.0*epsilon*(math.pow(sigma/r,12.0) - math.pow(sigma/r,6.0)) + + ralpha = math.exp(-width*(r-rzero)) + f2 = depth * (-1.0 + (1.0 - ralpha) * (1.0 - ralpha)) + + switch = 0.5*math.tanh(10.0*(r - rzero)) + 0.5 + f = switch*f1 + (1.0 - switch)*f2 + return f + +################################################################################ + +if __name__ == "__main__": + ptable = PairTabulate(hybrid_energy, units='lj', comment='Morse repulsion + LJ attraction') + ptable.run('MORSE_LJ') diff --git a/tools/tabulate/pair_lj_tabulate.py b/tools/tabulate/pair_lj_tabulate.py new file mode 100755 index 0000000000..713b10a761 --- /dev/null +++ b/tools/tabulate/pair_lj_tabulate.py @@ -0,0 +1,26 @@ +#!/usr/bin/env python + +from tabulate import PairTabulate + +################################################################################ +import math +epsilon = 1.0 +sigma = 1.0 + +def lj_energy(r): + f = 4.0*epsilon*(math.pow(sigma/r,12.0) - math.pow(sigma/r,6.0)) + return f + +def lj_force(r): + epsilon = 1.0 + sigma = 1.0 + f = -4.0*epsilon*(-12.0*math.pow(sigma/r,12.0)/r + 6.0*math.pow(sigma/r,6.0)/r) + return f +################################################################################ + +if __name__ == "__main__": + ptable = PairTabulate(lj_energy, lj_force) + ptable.run('LJ_11') + epsilon = 1.0 + sigma = 1.5 + ptable.run('LJ_12') diff --git a/tools/tabulate/tabulate.py b/tools/tabulate/tabulate.py new file mode 100755 index 0000000000..58b7055791 --- /dev/null +++ b/tools/tabulate/tabulate.py @@ -0,0 +1,279 @@ +# ---------------------------------------------------------------------- +# LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator +# https://www.lammps.org/ Sandia National Laboratories +# LAMMPS Development team: developers@lammps.org +# +# Copyright (2003) Sandia Corporation. Under the terms of Contract +# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains +# certain rights in this software. This software is distributed under +# the GNU General Public License. +# +# See the README file in the top-level LAMMPS directory. +# ------------------------------------------------------------------------- +""" +Python classes and support functions to generate tabulated potential files for LAMMPS. +""" + +# for python2/3 compatibility +from __future__ import print_function + +import sys +import argparse +from os import path +from datetime import datetime + +######################################################################## + +def numdiff(x, func): + """ Get the value of the first derivative of a function 'func(x)' at 'x' + from numerical differentiation.""" + + # optimal delta x value for 2-point numerical differentiation of two floating point numbers + epsilon = 6.05504e-6 + fval1 = func(x - epsilon) + fval2 = func(x + epsilon) + return 0.5 * (fval2-fval1) / epsilon + +######################################################################## + +def mktable(tstyle, label, num, xmin, xmax, efunc, diff=False, ffunc=None): + """ Do the tabulation of the provided energy function. Compute force from + numerical differentiation if no force function is provided. Also detect + minimum for use to determine potential shifting in bonded potentials.""" + + # must use numerical differentiation if no force function provided + if not ffunc: + diff = True + + print("# Creating %s table %s with %d points from %g to %g" % (tstyle, label, num, xmin, xmax)) + + table = [] + delx = (xmax - xmin) / (float(num) - 1.0) + emin = 999999999.0 + xzero = 0.0 + + for i in range(0,num): + x = xmin + i*delx + energy = efunc(x) + if energy < emin: + emin = energy + xzero = x + if diff: + force = -numdiff(x, efunc) + else: + force = ffunc(x) + table.append([i+1, x, energy, force]) + + return table, xzero + + +######################################################################## +# base class with shared functionality +class Tabulate(object): + """Base tabulation class. Contains all shared functionality: common argument parsing, + output file handling, table writing""" + + def __init__(self, style, efunc, ffunc=None, units=None, comment=None): + self.fp = sys.stdout + self.tstyle = style + self.energyfunc = efunc + self.forcefunc = ffunc + self.units = units + self.comment = comment + self.eshift = 0.0 + self.args = None + self.diff = True + + self.parser = argparse.ArgumentParser(description='Tool to generate tabulated ' + + self.tstyle + ' potential files for LAMMPS') + self.parser.add_argument('--num-points', '-n', dest='num', default=1000, type=int, + help="Number of tabulated points") + self.parser.add_argument('--filename', '-f', dest='filename', default='-', + help="Name of output file") + self.parser.add_argument('--diff-num', '-d', dest='diff',default=False, + action='store_true', + help="Differentiate energy function numerically") + self.parser.add_argument('--inner', '-i', dest='xmin', required=True, type=float, + help="Inner cutoff of table") + self.parser.add_argument('--outer', '-o', dest='xmax', required=True, type=float, + help="Outer cutoff of table") + + def openfile(self, label): + """Open table file, if needed and print label for new table entry""" + if self.args and self.args.filename != '-': + try: + if path.isfile(self.args.filename): + self.fp = open(self.args.filename, 'a') + print("# Appending table to file " + self.args.filename) + else: + self.fp = open(self.args.filename, 'w') + print("# Writing table to new file " + self.args.filename) + self.fp.write('# DATE: ' + datetime.now().date().isoformat()) + if self.units: + self.fp.write(' UNITS: ' + str(self.units)) + if self.comment: + self.fp.write(' COMMENT: ' + str(self.comment) + '\n') + except IOError: + sys.exit("Cannot open file %s for writing table data" % self.args.filename) + self.fp.write('\n' + label + '\n') + + def writetable(self, table, offset): + """ Formatted output tabulated data with 4 columns""" + for i,r,energy,force in table: + self.fp.write("%8d %- 22.15g %- 22.15g %- 22.15g\n" % (i, r, energy - offset, force)) + + def helpexit(self, text): + """ Convenience function to exit program with error and help message""" + sys.exit('\n' + text + '\n\n' + self.parser.format_help()) + +################################################################################ +# create tabulation for pair styles +class PairTabulate(Tabulate): + def __init__(self, efunc, ffunc=None, units=None, comment=None): + super(PairTabulate, self).__init__('pair', efunc, ffunc, units, comment) + self.parser.add_argument('--eshift', '-e', dest='eshift', default=False, + action='store_true', + help="Shift potential energy to be zero at outer cutoff") + try: + self.args = self.parser.parse_args() + except argparse.ArgumentError: + sys.exit() + + def run(self, label): + # sanity checks + if self.args.num < 2: + self.helpexit('Expect 2 or more points in table for tabulation') + if self.args.xmin <= 0.0: + self.helpexit('Inner tabulation cutoff must be > 0 for pair style table') + if self.args.xmax <= self.args.xmin: + self.helpexit('Outer cutoff must be larger than inner cutoff') + + self.diff = self.args.diff + if not self.forcefunc: + self.diff = True + offset = 0.0 + if self.args.eshift: + offset=self.energyfunc(self.args.xmax) + + table, dummy = mktable(self.tstyle, label, self.args.num, self.args.xmin, self.args.xmax, + self.energyfunc, self.args.diff, self.forcefunc) + + # open table file + self.openfile(label) + + # write pair style specific header + if self.forcefunc: + diffmin = -numdiff(self.args.xmin, self.forcefunc) + diffmax = -numdiff(self.args.xmax, self.forcefunc) + self.fp.write("N %d R %g %g FPRIME %- 22.15g %- 22.15g\n\n" + % (self.args.num, self.args.xmin, self.args.xmax, diffmin, diffmax)) + else: + self.fp.write("N %d R %g %g\n\n" % (self.args.num, self.args.xmin, self.args.xmax)) + + self.writetable(table, offset) + if self.args.filename != '-': + self.fp.close() + + +################################################################################ +# shared functionality to create tabulation for bond or angle styles +class BondAngleTabulate(Tabulate): + def __init__(self, style, efunc, ffunc=None, units=None, comment=None): + super(BondAngleTabulate, self).__init__(style, efunc, ffunc, units, comment) + self.parser.add_argument('--eshift', '-e', dest='eshift', default=False, + action='store_true', + help="Shift potential energy to be zero at minimum") + idx = [a.dest for a in self.parser._actions].index('xmin') + self.parser._actions[idx].required=False + self.parser._actions[idx].default=0.0 + if style == 'angle': + idx = [a.dest for a in self.parser._actions].index('xmax') + self.parser._actions[idx].required=False + self.parser._actions[idx].default=180.0 + try: + self.args = self.parser.parse_args() + except argparse.ArgumentError: + sys.exit() + + def run(self, label): + # sanity checks + if self.args.num < 2: + self.helpexit('Expect 2 or more points in table for tabulation') + if self.args.xmin < 0.0: + self.helpexit('Inner cutoff must not be negative') + if self.tstyle == 'angle' and self.args.xmax > 180.0: + self.helpexit('Outer cutoff must not be larger than 180.0 degrees') + + self.diff = self.args.diff + if not self.forcefunc: + self.diff = True + + table, xzero = mktable(self.tstyle, label, self.args.num, self.args.xmin, self.args.xmax, + self.energyfunc, self.args.diff, self.forcefunc) + print("# Minimum energy of tabulated potential is at %g" % xzero) + offset = 0.0 + if self.args.eshift: + offset=self.energyfunc(xzero) + + self.openfile(label) + + if self.forcefunc: + diffmin = -numdiff(self.args.xmin, self.forcefunc) + diffmax = -numdiff(self.args.xmax, self.forcefunc) + self.fp.write("N %d FP %- 22.15g %- 22.15g EQ %g\n\n" % + (self.args.num, diffmin, diffmax, xzero)) + else: + self.fp.write("N %d EQ %g\n\n" % (self.args.num, xzero)) + + self.writetable(table, offset) + if self.args.filename != '-': + self.fp.close() + +################################################################################ +class BondTabulate(BondAngleTabulate): + def __init__(self, efunc, ffunc=None, units=None, comment=None): + super(BondTabulate, self).__init__('bond', efunc, ffunc, units, comment) + +################################################################################ +class AngleTabulate(BondAngleTabulate): + def __init__(self, efunc, ffunc=None, units=None, comment=None): + super(AngleTabulate, self).__init__('angle', efunc, ffunc, units, comment) + +################################################################################ +# create tabulation for dihdedral +class DihedralTabulate(Tabulate): + def __init__(self, efunc, ffunc=None, units=None, comment=None): + super(DihedralTabulate, self).__init__('dihedral', efunc, ffunc, units, comment) + idx = [a.dest for a in self.parser._actions].index('xmin') + self.parser._actions[idx].required=False + self.parser._actions[idx].default=-180.0 + idx = [a.dest for a in self.parser._actions].index('xmax') + self.parser._actions[idx].required=False + self.parser._actions[idx].default=179.999999 + try: + self.args = self.parser.parse_args() + except argparse.ArgumentError: + sys.exit() + + def run(self, label): + # sanity checks + if self.args.num < 2: + self.helpexit('Expect 2 or more points in table for tabulation') + if self.args.xmin < -180 or self.args.xmin > 0.0: + self.helpexit('Inner cutoff must be within -180.0 and 0.0 degrees') + if self.args.xmax < 0.0 or self.args.xmin > 360.0: + self.helpexit('Outer cutoff must be within 0.0 and 360.0 degrees') + if (self.args.xmax - self.args.xmin) >= 360.0: + self.helpexit('Inner and outer cutoff range must be less than 360.0 degrees') + + self.diff = self.args.diff + if not self.forcefunc: + self.diff = True + + table, dummy = mktable(self.tstyle, label, self.args.num, self.args.xmin, self.args.xmax, + self.energyfunc, self.args.diff, self.forcefunc) + self.openfile(label) + self.fp.write("N %d DEGREES \n\n" % (self.args.num)) + self.writetable(table, 0.0) + if self.args.filename != '-': + self.fp.close() From 8998ef23de538c1d835ec839287aae0a01fbc159 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 16 Dec 2022 02:49:39 -0500 Subject: [PATCH 11/12] update docs for tools/tabulate scripts --- doc/src/Tools.rst | 20 +++++++++++++++++--- doc/src/angle_table.rst | 4 ++++ doc/src/bond_table.rst | 10 +++++++++- doc/src/dihedral_table.rst | 4 ++++ doc/src/pair_table.rst | 7 +++++++ tools/README | 1 + 6 files changed, 42 insertions(+), 4 deletions(-) diff --git a/doc/src/Tools.rst b/doc/src/Tools.rst index b70c7677e0..73ca11e689 100644 --- a/doc/src/Tools.rst +++ b/doc/src/Tools.rst @@ -57,7 +57,7 @@ Pre-processing tools * :ref:`msi2lmp ` * :ref:`polybond ` * :ref:`stl_bin2txt ` - + * :ref:`tabulate ` Post-processing tools ===================== @@ -1159,13 +1159,27 @@ For illustration purposes below is a part of the Tcl example script. ---------- +.. _tabulate: + +tabulate tool +-------------- + +.. versionadded:: TBD + +The ``tabulate`` folder contains Python scripts scripts to generate tabulated +potential files for LAMMPS. The bulk of the code is in the ``tabulate`` module +in the ``tabulate.py`` file. Some example files demonstrating its use are +included. See the README file for more information. + +---------- + .. _vim: vim tool ------------------ -The files in the tools/vim directory are add-ons to the VIM editor -that allow easier editing of LAMMPS input scripts. See the README.txt +The files in the ``tools/vim`` directory are add-ons to the VIM editor +that allow easier editing of LAMMPS input scripts. See the ``README.txt`` file for details. These files were provided by Gerolf Ziegenhain (gerolf at diff --git a/doc/src/angle_table.rst b/doc/src/angle_table.rst index b8aeae780b..c0faa3f046 100644 --- a/doc/src/angle_table.rst +++ b/doc/src/angle_table.rst @@ -59,6 +59,10 @@ format of this file is described below. ---------- +Suitable tables for use with this angle style can be created using the +Python code in the ``tools/tabulate`` folder of the LAMMPS source code +distribution. + The format of a tabulated file is as follows (without the parenthesized comments): diff --git a/doc/src/bond_table.rst b/doc/src/bond_table.rst index 898503b951..51e677d74a 100644 --- a/doc/src/bond_table.rst +++ b/doc/src/bond_table.rst @@ -59,6 +59,13 @@ this file is described below. ---------- +Suitable tables for use with this bond style can be created by LAMMPS +itself from existing bond styles using the :doc:`bond_write +` command. This can be useful to have a template file for +testing the bond style settings and to build a compatible custom file. +Another option to generate tables is the Python code in the +``tools/tabulate`` folder of the LAMMPS source code distribution. + The format of a tabulated file is as follows (without the parenthesized comments): @@ -149,7 +156,8 @@ info. Related commands """""""""""""""" -:doc:`bond_coeff `, :doc:`delete_bonds ` +:doc:`bond_coeff `, :doc:`delete_bonds `, +:doc:`bond_write ` Default """"""" diff --git a/doc/src/dihedral_table.rst b/doc/src/dihedral_table.rst index e0830b523c..5e8d589c55 100644 --- a/doc/src/dihedral_table.rst +++ b/doc/src/dihedral_table.rst @@ -114,6 +114,10 @@ below. ---------- +Suitable tables for use with this dihedral style can be created using +the Python code in the ``tools/tabulate`` folder of the LAMMPS source +code distribution. + The format of a tabulated file is as follows (without the parenthesized comments). It can begin with one or more comment or blank lines. diff --git a/doc/src/pair_table.rst b/doc/src/pair_table.rst index 79eb8f5b42..e5b05a3a23 100644 --- a/doc/src/pair_table.rst +++ b/doc/src/pair_table.rst @@ -120,6 +120,13 @@ best effect: ---------- +Suitable tables in the correct format for use with these pair styles can +be created by LAMMPS itself using the :doc:`pair_write ` +command. In combination with the :doc:`pair style python ` +this can be a powerful mechanism to implement and test tables for use +with LAMMPS. Another option to generate tables is the Python code in +the ``tools/tabulate`` folder of the LAMMPS source code distribution. + The format of a tabulated file has an (optional) header followed by a series of one or more sections, defined as follows (without the parenthesized comments). The header must start with a `#` character diff --git a/tools/README b/tools/README index 22d9c8ea18..82af6f3f66 100644 --- a/tools/README +++ b/tools/README @@ -48,6 +48,7 @@ smd convert Smooth Mach Dynamics triangles to VTK spin perform a cubic polynomial interpolation of a GNEB MEP stl_bin2txt convert binary STL files to ASCII swig Interface file and demo scripts for SWIG wrappers for the LAMMPS C library interface +tabulate Python scripts to generate tabulated potential files for LAMMPS tinker2lmp.py convert Tinker input files to LAMMPS input files valgrind suppression files for use with valgrind's memcheck tool vim add-ons to VIM editor for editing LAMMPS input scripts From c421a445bda06f59d2bc11ec7ebf3caf50d97833 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Fri, 16 Dec 2022 15:18:57 -0700 Subject: [PATCH 12/12] Small cleanup, docs --- doc/src/Commands_fix.rst | 4 ++-- src/KOKKOS/fix_dt_reset_kokkos.cpp | 12 +++++++----- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index ff6e402a2e..0d3f69602e 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -65,7 +65,7 @@ OPT. * :doc:`drude ` * :doc:`drude/transform/direct ` * :doc:`drude/transform/inverse ` - * :doc:`dt/reset ` + * :doc:`dt/reset (k) ` * :doc:`edpd/source ` * :doc:`efield ` * :doc:`ehex ` @@ -250,7 +250,7 @@ OPT. * :doc:`tune/kspace ` * :doc:`vector ` * :doc:`viscosity ` - * :doc:`viscous ` + * :doc:`viscous (k) ` * :doc:`viscous/sphere ` * :doc:`wall/body/polygon ` * :doc:`wall/body/polyhedron ` diff --git a/src/KOKKOS/fix_dt_reset_kokkos.cpp b/src/KOKKOS/fix_dt_reset_kokkos.cpp index 6ffa7fcf56..f3435e711e 100644 --- a/src/KOKKOS/fix_dt_reset_kokkos.cpp +++ b/src/KOKKOS/fix_dt_reset_kokkos.cpp @@ -118,6 +118,8 @@ void FixDtResetKokkos::end_of_step() } +/* ---------------------------------------------------------------------- */ + template KOKKOS_INLINE_FUNCTION void FixDtResetKokkos::operator()(TagFixDtResetMass, const int &i, double &k_dt) const { @@ -127,7 +129,7 @@ void FixDtResetKokkos::operator()(TagFixDtResetMass, const int &i, d double delx, dely, delz, delr; double emax = k_emax.d_view(0); - + if (mask[i] & groupbit) { massinv = 1.0 / mass[type[i]]; @@ -148,9 +150,10 @@ void FixDtResetKokkos::operator()(TagFixDtResetMass, const int &i, d delr = sqrt(delx * delx + dely * dely + delz * delz); if (delr > xmax) k_dt *= xmax / delr; } - } +/* ---------------------------------------------------------------------- */ + template KOKKOS_INLINE_FUNCTION void FixDtResetKokkos::operator()(TagFixDtResetRMass, const int &i, double &k_dt) const { @@ -160,7 +163,7 @@ void FixDtResetKokkos::operator()(TagFixDtResetRMass, const int &i, double delx, dely, delz, delr; double emax = k_emax.d_view(0); - + if (mask[i] & groupbit) { massinv = 1.0 / rmass[i]; @@ -181,8 +184,7 @@ void FixDtResetKokkos::operator()(TagFixDtResetRMass, const int &i, delr = sqrt(delx * delx + dely * dely + delz * delz); if (delr > xmax) k_dt *= xmax / delr; } - - } +} namespace LAMMPS_NS { template class FixDtResetKokkos;