diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 1fa84fe370..0169826751 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -60,6 +60,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`erotate/sphere ` * :doc:`erotate/sphere/atom ` * :doc:`event/displace ` + * :doc:`fabric ` * :doc:`fep ` * :doc:`force/tally ` * :doc:`fragment/atom ` diff --git a/doc/src/Howto_granular.rst b/doc/src/Howto_granular.rst index 70c2b2fbb7..9492e5755e 100644 --- a/doc/src/Howto_granular.rst +++ b/doc/src/Howto_granular.rst @@ -17,6 +17,11 @@ This compute * :doc:`compute erotate/sphere ` calculates rotational kinetic energy which can be :doc:`output with thermodynamic info `. +The compute + +* :doc:`compute fabric ` + +calculates various versions of the fabric tensor for granular and non-granular pair styles. Use one of these 4 pair potentials, which compute forces and torques between interacting pairs of particles: diff --git a/doc/src/Library_create.rst b/doc/src/Library_create.rst index 350569e54e..b3fea4b89e 100644 --- a/doc/src/Library_create.rst +++ b/doc/src/Library_create.rst @@ -9,6 +9,7 @@ This section documents the following functions: - :cpp:func:`lammps_close` - :cpp:func:`lammps_mpi_init` - :cpp:func:`lammps_mpi_finalize` +- :cpp:func:`lammps_kokkos_finalize` -------------------- diff --git a/doc/src/Python_create.rst b/doc/src/Python_create.rst index 00b1c08814..ec4241f36a 100644 --- a/doc/src/Python_create.rst +++ b/doc/src/Python_create.rst @@ -134,7 +134,10 @@ compiled with. The :py:func:`lmp.close() ` call is optional since the LAMMPS class instance will also be deleted automatically during the :py:class:`lammps ` class -destructor. +destructor. Instead of :py:func:`lmp.close() ` +it is also possible to call :py:func:`lmp.finalize() `; +this will destruct the LAMMPS instance, but also finalized and release +the MPI and/or Kokkos environment if enabled and active. Note that you can create multiple LAMMPS objects in your Python script, and coordinate and run multiple simulations, e.g. diff --git a/doc/src/compute.rst b/doc/src/compute.rst index afe54a69be..c2b2b92aa1 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -206,6 +206,7 @@ The individual style names on the :doc:`Commands compute ` doc * :doc:`erotate/sphere ` - rotational energy of spherical particles * :doc:`erotate/sphere/atom ` - rotational energy for each spherical particle * :doc:`event/displace ` - detect event on atom displacement +* :doc:`fabric ` - calculates fabric tensors from pair interactions * :doc:`fep ` - * :doc:`force/tally ` - * :doc:`fragment/atom ` - fragment ID for each atom diff --git a/doc/src/compute_fabric.rst b/doc/src/compute_fabric.rst new file mode 100644 index 0000000000..3ded3d7adc --- /dev/null +++ b/doc/src/compute_fabric.rst @@ -0,0 +1,185 @@ +.. index:: compute fabric + +compute fabric command +====================== + +Syntax +"""""" + +.. parsed-literal:: + + compute ID group-ID fabric cutoff attribute1 attribute2 ... keyword values ... + +* ID, group-ID are documented in :doc:`compute ` command +* fabric = style name of this compute command +* cutoff = *type* or *radius* + + .. parsed-literal:: + + *type* = cutoffs determined based on atom types + *radius* = cutoffs determined based on atom diameters (atom style sphere) + +* one or more attributes may be appended + + .. parsed-literal:: + + *contact* = contact tensor + *branch* = branch tensor + *force/normal* = normal force tensor + *force/tangential* = tangential force tensor + +* zero or more keyword/value pairs may be appended +* keyword = *type/include* + + .. parsed-literal:: + + *type/include* value = arg1 arg2 + arg = separate lists of types (see below) + +Examples +"""""""" + +.. code-block:: LAMMPS + + compute 1 all fabric type contact force/normal type/include 1,2 3*4 + compute 1 all fabric radius force/normal force/tangential + +Description +""""""""""" + +Define a compute that calculates various fabric tensors for pairwise +interaction :ref:`(Ouadfel) `. Fabric tensors are commonly used +to quantify the anisotropy or orientation of granular contacts but can also +be used to characterize the direction of pairwise interactions in general +systems. The *type* and *radius* settings are used to select whether interactions +cutoffs are determined by atom types or by the sum of atomic radii (atom +style sphere), respectively. Calling this compute is roughly the cost of a +pair style invocation as it involves a loop over the neighbor list. If the +normal or tangential force tensors are requested, it will be more expensive +than a pair style invocation as it will also recalculate all pair forces. + +Four fabric tensors are available: the contact, branch, normal force, or +tangential force tensor. The contact tensor is calculated as + +.. math:: + + C_{ab} = \frac{15}{2} (\phi_{ab} - \mathrm{tr}(\phi) \delta_{ab}) + +where :math:`a` and :math:`b` are the :math:`x`, :math:`y`, :math:`z` +directions, :math:`\delta_{ab}` is the Kronecker delta function, and +the tensor :math:`\phi` is defined as + +.. math:: + + \phi_{ab} = \sum_{n = 1}^{N_p} \frac{r_{a} r_{b}}{r^2} + +where :math:`n` loops over the :math:`N_p` pair interactions in the simulation, +:math:`r_{a}` is the :math:`a` component of the radial vector between the +two pairwise interacting particles, and :math:`r` is the magnitude of the radial vector. + +The branch tensor is calculated as + +.. math:: + + B_{ab} = \frac{15}{6 \mathrm{tr}(D)} (D_{ab} - \mathrm{tr}(D) \delta_{ab}) + +where the tensor :math:`D` is defined as + +.. math:: + + D_{ab} = \sum_{n = 1}^{N_p} + \frac{1}{N_c (r^2 + C_{cd} r_c r_d)} + \frac{r_{a} r_{b}}{r} + +where :math:`N_c` is the total number of contacts in the system and the subscripts +:math:`c` and :math:`d` indices are summed according to Einstein notation. + +The normal force fabric tensor is calculated as + +.. math:: + + F^n_{ab} = \frac{15}{6 \mathrm{tr}(N)} (N_{ab} - \mathrm{tr}(N) \delta_{ab}) + +where the tensor :math:`N` is defined as + +.. math:: + + N_{ab} = \sum_{n = 1}^{N_p} + \frac{1}{N_c (r^2 + C_{cd} r_c r_d)} + \frac{r_{a} r_{b}}{r^2} f_n + +and :math:`f_n` is the magnitude of the normal, central-body force between the two atoms. + +Finally, the tangential force fabric tensor is only defined for pair styles that +apply tangential forces to particles, namely granular pair styles. It is calculated +as + +.. math:: + + F^t_{ab} = \frac{15}{9 \mathrm{tr}(N)} (T_{ab} - \mathrm{tr}(T) \delta_{ab}) + +where the tensor :math:`T` is defined as + +.. math:: + + T_{ab} = \sum_{n = 1}^{N_p} + \frac{1}{N_c (r^2 + C_{cd} r_c r_d)} + \frac{r_{a} r_{b}}{r^2} f_t + +and :math:`f_t` is the magnitude of the tangential force between the two atoms. + +The *type/include* keyword filters interactions based on the types of the two atoms. +Interactions between two atoms are only included in calculations if the atom types +are in the two lists. Each list consists of a series of type +ranges separated by commas. The range can be specified as a +single numeric value, or a wildcard asterisk can be used to specify a range +of values. This takes the form "\*" or "\*n" or "n\*" or "m\*n". For +example, if M = the number of atom types, then an asterisk with no numeric +values means all types from 1 to M. A leading asterisk means all types +from 1 to n (inclusive). A trailing asterisk means all types from n to M +(inclusive). A middle asterisk means all types from m to n (inclusive). +Multiple *type/include* keywords may be added. + +Output info +""""""""""" + +This compute calculates a local vector of doubles and a scalar. The vector stores the +unique components of the first requested tensor in the order xx, yy, zz, xy, xz, yz +followed by the same components for all subsequent tensors. The length of the vector +is therefore six times the number of requested tensors. The scalar output is the +number of pairwise interactions included in the calculation of the fabric tensor. + +Restrictions +"""""""""""" + +This fix is part of the GRANULAR package. It is only enabled if LAMMPS +was built with that package. See the :doc:`Build package ` +doc page for more info. + +Currently, compute *fabric* does not support pair styles +with many-body interactions. It also does not +support models with long-range Coulombic or dispersion forces, +i.e. the kspace_style command in LAMMPS. It also does not support the +following fixes which add rigid-body constraints: :doc:`fix shake +`, :doc:`fix rattle `, :doc:`fix rigid +`, :doc:`fix rigid/small `. It does not support +granular pair styles that extend beyond the contact of atomic radii +(e.g. JKR and DMT). + +Related commands +"""""""""""""""" + +none + +Default +""""""" + +none + +---------- + +.. _Ouadfel: + +**(Ouadfel)** Ouadfel and Rothenburg +"Stress-force-fabric relationship for assemblies of ellipsoids", +Mechanics of Materials (2001). (`link to paper `_) diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 12e6fc88f0..e57f81dbc8 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -2680,6 +2680,7 @@ qoverride qqr qqrd qtb +Quadfel quadratically quadrupolar Quant @@ -2861,6 +2862,7 @@ Rossky rosybrown rotationally Rotenberg +Rothenburg Rovigatti royalblue rozero diff --git a/fortran/lammps.f90 b/fortran/lammps.f90 index 21909e1288..2e78cdd00b 100644 --- a/fortran/lammps.f90 +++ b/fortran/lammps.f90 @@ -76,17 +76,15 @@ MODULE LIBLAMMPS TYPE(c_ptr), VALUE :: handle END SUBROUTINE lammps_close - SUBROUTINE lammps_mpi_init(handle) BIND(C, name='lammps_mpi_init') - IMPORT :: c_ptr - TYPE(c_ptr), VALUE :: handle + SUBROUTINE lammps_mpi_init() BIND(C, name='lammps_mpi_init') END SUBROUTINE lammps_mpi_init - SUBROUTINE lammps_mpi_finalize(handle) & - BIND(C, name='lammps_mpi_finalize') - IMPORT :: c_ptr - TYPE(c_ptr), VALUE :: handle + SUBROUTINE lammps_mpi_finalize() BIND(C, name='lammps_mpi_finalize') END SUBROUTINE lammps_mpi_finalize + SUBROUTINE lammps_kokkos_finalize() BIND(C, name='lammps_kokkos_finalize') + END SUBROUTINE lammps_kokkos_finalize + SUBROUTINE lammps_file(handle,filename) BIND(C, name='lammps_file') IMPORT :: c_ptr TYPE(c_ptr), VALUE :: handle @@ -188,7 +186,8 @@ CONTAINS IF (PRESENT(finalize)) THEN IF (finalize) THEN - CALL lammps_mpi_finalize(self%handle) + CALL lammps_kokkos_finalize() + CALL lammps_mpi_finalize() END IF END IF END SUBROUTINE lmp_close diff --git a/python/lammps/core.py b/python/lammps/core.py index df5e578683..2f101f4eab 100644 --- a/python/lammps/core.py +++ b/python/lammps/core.py @@ -460,10 +460,16 @@ class lammps(object): # ------------------------------------------------------------------------- def finalize(self): - """Shut down the MPI communication through the library interface by calling :cpp:func:`lammps_finalize`. + """Shut down the MPI communication and Kokkos environment (if active) through the + library interface by calling :cpp:func:`lammps_mpi_finalize` and + :cpp:func:`lammps_kokkos_finalize`. + + You cannot create or use any LAMMPS instances after this function is called + unless LAMMPS was compiled without MPI and without Kokkos support. """ self.close() - self.lib.lammps_finalize() + self.lib.lammps_kokkos_finalize() + self.lib.lammps_mpi_finalize() # ------------------------------------------------------------------------- diff --git a/src/.gitignore b/src/.gitignore index 67072d5137..c562f89e61 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -386,6 +386,8 @@ /compute_erotate_rigid.h /compute_event_displace.cpp /compute_event_displace.h +/compute_fabric.cpp +/compute_fabric.h /compute_fep.cpp /compute_fep.h /compute_force_tally.cpp diff --git a/src/GRANULAR/compute_fabric.cpp b/src/GRANULAR/compute_fabric.cpp new file mode 100644 index 0000000000..61d3b632a5 --- /dev/null +++ b/src/GRANULAR/compute_fabric.cpp @@ -0,0 +1,595 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 "compute_fabric.h" + +#include "atom.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "neighbor.h" +#include "pair.h" +#include "tokenizer.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; + +enum { OTHER, GRANULAR }; +enum { TYPE, RADIUS }; +enum { CN, BR, FN, FT }; + +/* ---------------------------------------------------------------------- */ + +ComputeFabric::ComputeFabric(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), tensor_style(NULL) +{ + if (narg < 3) error->all(FLERR, "Illegal compute fabric command"); + + if (strcmp(arg[3], "type") == 0) + cutstyle = TYPE; + else if (strcmp(arg[3], "radius") == 0) + cutstyle = RADIUS; + else + error->all(FLERR, "Illegal compute fabric command"); + + if (cutstyle == RADIUS && !atom->radius_flag) + error->all(FLERR, "Compute fabric radius style requires atom attribute radius"); + + // If optional arguments included, this will be oversized + ntensors = narg - 4; + tensor_style = new int[ntensors]; + + cn_flag = 0; + br_flag = 0; + fn_flag = 0; + ft_flag = 0; + type_filter = nullptr; + + ntensors = 0; + int iarg = 4; + while (iarg < narg) { + if (strcmp(arg[iarg], "contact") == 0) { + cn_flag = 1; + tensor_style[ntensors++] = CN; + } else if (strcmp(arg[iarg], "branch") == 0) { + br_flag = 1; + tensor_style[ntensors++] = BR; + } else if (strcmp(arg[iarg], "force/normal") == 0) { + fn_flag = 1; + tensor_style[ntensors++] = FN; + } else if (strcmp(arg[iarg], "force/tangential") == 0) { + ft_flag = 1; + tensor_style[ntensors++] = FT; + } else if (strcmp(arg[iarg], "type/include") == 0) { + if (iarg + 1 >= narg) error->all(FLERR, "Invalid keyword in compute fabric command"); + int ntypes = atom->ntypes; + + int i, j, itype, jtype, in, jn, infield, jnfield; + int inlo, inhi, jnlo, jnhi; + char *istr, *jstr; + if (!type_filter) { + memory->create(type_filter, ntypes + 1, ntypes + 1, "compute/fabric:type_filter"); + + for (i = 0; i <= ntypes; i++) { + for (j = 0; j <= ntypes; j++) { type_filter[i][j] = 0; } + } + } + + in = strlen(arg[iarg + 1]) + 1; + istr = new char[in]; + strcpy(istr, arg[iarg + 1]); + std::vector iwords = Tokenizer(istr, ",").as_vector(); + infield = iwords.size(); + + jn = strlen(arg[iarg + 2]) + 1; + jstr = new char[jn]; + strcpy(jstr, arg[iarg + 2]); + std::vector jwords = Tokenizer(jstr, ",").as_vector(); + jnfield = jwords.size(); + + for (i = 0; i < infield; i++) { + const char *ifield = iwords[i].c_str(); + utils::bounds(FLERR, ifield, 1, ntypes, inlo, inhi, error); + + for (j = 0; j < jnfield; j++) { + const char *jfield = jwords[j].c_str(); + utils::bounds(FLERR, jfield, 1, ntypes, jnlo, jnhi, error); + + for (itype = inlo; itype <= inhi; itype++) { + for (jtype = jnlo; jtype <= jnhi; jtype++) { + type_filter[itype][jtype] = 1; + type_filter[jtype][itype] = 1; + } + } + } + } + + delete[] istr; + delete[] jstr; + + iarg += 2; + } else + error->all(FLERR, "Illegal compute fabric command"); + iarg++; + } + + vector_flag = 1; + size_vector = ntensors * 6; + extvector = 0; + + scalar_flag = 1; + extscalar = 1; + + vector = new double[size_vector]; +} + +/* ---------------------------------------------------------------------- */ + +ComputeFabric::~ComputeFabric() +{ + delete[] vector; + delete[] tensor_style; + memory->destroy(type_filter); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeFabric::init() +{ + if (force->pair == NULL) error->all(FLERR, "No pair style is defined for compute fabric"); + if (force->pair->single_enable == 0 && (fn_flag || ft_flag)) + error->all(FLERR, "Pair style does not support compute fabric normal or tangential force"); + + // Find if granular or gran + pstyle = OTHER; + if (force->pair_match("^granular", 0) || force->pair_match("^gran/", 0)) pstyle = GRANULAR; + + if (pstyle != GRANULAR && ft_flag) + error->all(FLERR, "Pair style does not calculate tangential forces for compute fabric"); + + if (force->pair->beyond_contact) + error->all(FLERR, "Compute fabric does not support pair styles that extend beyond contact"); + + // need an occasional half neighbor list + // set size to same value as request made by force->pair + // this should enable it to always be a copy list (e.g. for granular pstyle) + + int irequest = neighbor->request(this, instance_me); + neighbor->requests[irequest]->pair = 0; + neighbor->requests[irequest]->compute = 1; + neighbor->requests[irequest]->occasional = 1; + NeighRequest *pairrequest = neighbor->find_request((void *) force->pair); + if (pairrequest) neighbor->requests[irequest]->size = pairrequest->size; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeFabric::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeFabric::compute_vector() +{ + invoked_vector = update->ntimestep; + + int i, j, ii, jj, inum, jnum, itype, jtype; + tagint itag, jtag; + double xtmp, ytmp, ztmp, delx, dely, delz; + double r, rinv, rsq, radsum, eng, fpair; + + double nx, ny, nz; + double ncinv, denom, fn, ft, prefactor; + double br_tensor[6], ft_tensor[6], fn_tensor[6]; + double trace_phi, trace_D, trace_Xfn, trace_Xft; + double phi_ij[6] = {0.0}; + double Ac_ij[6] = {0.0}; + double D_ij[6] = {0.0}; + double Xfn_ij[6] = {0.0}; + double Xft_ij[6] = {0.0}; + double temp_dbl[6]; + + int *ilist, *jlist, *numneigh, **firstneigh; + + double **x = atom->x; + double *radius = atom->radius; + tagint *tag = atom->tag; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + // invoke half neighbor list (will copy or build if necessary) + neighbor->build_one(list); + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + Pair *pair = force->pair; + double **cutsq = force->pair->cutsq; + + // invoke compute_scalar() to update the number of contacts, if needed + nc = compute_scalar(); + + // If no contacts, everything will be zero + if (nc == 0) { + for (i = 0; i < size_vector; i++) vector[i] = 0.0; + return; + } + ncinv = 1.0 / nc; + + // First loop through and calculate contact fabric tensor + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(mask[i] & groupbit)) continue; + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itag = tag[i]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + if (!(mask[j] & groupbit)) continue; + + // itag = jtag is possible for long cutoffs that include images of self + + if (newton_pair == 0 && j >= nlocal) { + jtag = tag[j]; + if (itag > jtag) { + if ((itag + jtag) % 2 == 0) continue; + } else if (itag < jtag) { + if ((itag + jtag) % 2 == 1) continue; + } else { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + } + } + + jtype = type[j]; + + if (type_filter) + if (type_filter[itype][jtype] == 0) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + + if (cutstyle == TYPE) { + if (rsq >= cutsq[itype][jtype]) continue; + } else { + radsum = radius[i] + radius[j]; + if (rsq >= radsum * radsum) continue; + } + + r = sqrt(rsq); + rinv = 1.0 / r; + nx = delx * rinv; + ny = dely * rinv; + nz = delz * rinv; + + phi_ij[0] += nx * nx; + phi_ij[1] += ny * ny; + phi_ij[2] += nz * nz; + phi_ij[3] += nx * ny; + phi_ij[4] += nx * nz; + phi_ij[5] += ny * nz; + } + } + + //Sum phi across processors + MPI_Allreduce(phi_ij, temp_dbl, 6, MPI_DOUBLE, MPI_SUM, world); + for (i = 0; i < 6; i++) phi_ij[i] = temp_dbl[i] * ncinv; + + trace_phi = (1.0 / 3.0) * (phi_ij[0] + phi_ij[1] + phi_ij[2]); + + Ac_ij[0] = (15.0 / 2.0) * (phi_ij[0] - trace_phi); + Ac_ij[1] = (15.0 / 2.0) * (phi_ij[1] - trace_phi); + Ac_ij[2] = (15.0 / 2.0) * (phi_ij[2] - trace_phi); + Ac_ij[3] = (15.0 / 2.0) * (phi_ij[3]); + Ac_ij[4] = (15.0 / 2.0) * (phi_ij[4]); + Ac_ij[5] = (15.0 / 2.0) * (phi_ij[5]); + + // If needed, loop through and calculate other fabric tensors + if (br_flag || fn_flag || ft_flag) { + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(mask[i] & groupbit)) continue; + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itag = tag[i]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + if (!(mask[j] & groupbit)) continue; + + // itag = jtag is possible for long cutoffs that include images of self + + if (newton_pair == 0 && j >= nlocal) { + jtag = tag[j]; + if (itag > jtag) { + if ((itag + jtag) % 2 == 0) continue; + } else if (itag < jtag) { + if ((itag + jtag) % 2 == 1) continue; + } else { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + } + } + + jtype = type[j]; + + if (type_filter) + if (type_filter[itype][jtype] == 0) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + + if (cutstyle == TYPE) { + if (rsq >= cutsq[itype][jtype]) continue; + } else { + radsum = radius[i] + radius[j]; + if (rsq >= radsum * radsum) continue; + } + + if (fn_flag || ft_flag) eng = pair->single(i, j, itype, jtype, rsq, 1.0, 1.0, fpair); + + r = sqrt(rsq); + rinv = 1.0 / r; + nx = delx * rinv; + ny = dely * rinv; + nz = delz * rinv; + + denom = 1 + Ac_ij[0] * nx * nx + Ac_ij[1] * ny * ny + Ac_ij[2] * nz * nz; + denom += 2 * Ac_ij[3] * nx * ny + 2 * Ac_ij[4] * nx * nz + 2 * Ac_ij[5] * ny * nz; + prefactor = ncinv / denom; + + if (br_flag) { + D_ij[0] += prefactor * nx * nx * r; + D_ij[1] += prefactor * ny * ny * r; + D_ij[2] += prefactor * nz * nz * r; + D_ij[3] += prefactor * nx * ny * r; + D_ij[4] += prefactor * nx * nz * r; + D_ij[5] += prefactor * ny * nz * r; + } + + if (fn_flag || ft_flag) { + fn = r * fpair; + + Xfn_ij[0] += prefactor * nx * nx * fn; + Xfn_ij[1] += prefactor * ny * ny * fn; + Xfn_ij[2] += prefactor * nz * nz * fn; + Xfn_ij[3] += prefactor * nx * ny * fn; + Xfn_ij[4] += prefactor * nx * nz * fn; + Xfn_ij[5] += prefactor * ny * nz * fn; + + if (ft_flag) { + ft = force->pair->svector[3]; + + Xft_ij[0] += prefactor * nx * nx * ft; + Xft_ij[1] += prefactor * ny * ny * ft; + Xft_ij[2] += prefactor * nz * nz * ft; + Xft_ij[3] += prefactor * nx * ny * ft; + Xft_ij[4] += prefactor * nx * nz * ft; + Xft_ij[5] += prefactor * ny * nz * ft; + } + } + } + } + } + + // Output results + + if (cn_flag) { + for (i = 0; i < ntensors; i++) { + if (tensor_style[i] == CN) { + for (j = 0; j < 6; j++) vector[6 * i + j] = Ac_ij[j]; + } + } + } + + if (br_flag) { + MPI_Allreduce(D_ij, temp_dbl, 6, MPI_DOUBLE, MPI_SUM, world); + for (i = 0; i < 6; i++) D_ij[i] = temp_dbl[i]; + + trace_D = (1.0 / 3.0) * (D_ij[0] + D_ij[1] + D_ij[2]); + + br_tensor[0] = (15.0 / (6.0 * trace_D)) * (D_ij[0] - trace_D); + br_tensor[1] = (15.0 / (6.0 * trace_D)) * (D_ij[1] - trace_D); + br_tensor[2] = (15.0 / (6.0 * trace_D)) * (D_ij[2] - trace_D); + br_tensor[3] = (15.0 / (6.0 * trace_D)) * (D_ij[3]); + br_tensor[4] = (15.0 / (6.0 * trace_D)) * (D_ij[4]); + br_tensor[5] = (15.0 / (6.0 * trace_D)) * (D_ij[5]); + + for (i = 0; i < ntensors; i++) { + if (tensor_style[i] == BR) { + for (j = 0; j < 6; j++) vector[6 * i + j] = br_tensor[j]; + } + } + } + + if (fn_flag || ft_flag) { + MPI_Allreduce(Xfn_ij, temp_dbl, 6, MPI_DOUBLE, MPI_SUM, world); + for (i = 0; i < 6; i++) Xfn_ij[i] = temp_dbl[i]; + + trace_Xfn = (1.0 / 3.0) * (Xfn_ij[0] + Xfn_ij[1] + Xfn_ij[2]); + } + + if (fn_flag) { + + fn_tensor[0] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[0] - trace_Xfn); + fn_tensor[1] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[1] - trace_Xfn); + fn_tensor[2] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[2] - trace_Xfn); + fn_tensor[3] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[3]); + fn_tensor[4] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[4]); + fn_tensor[5] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[5]); + + for (i = 0; i < ntensors; i++) { + if (tensor_style[i] == FN) { + for (j = 0; j < 6; j++) vector[6 * i + j] = fn_tensor[j]; + } + } + } + + if (ft_flag) { + MPI_Allreduce(Xft_ij, temp_dbl, 6, MPI_DOUBLE, MPI_SUM, world); + for (i = 0; i < 6; i++) Xft_ij[i] = temp_dbl[i]; + + trace_Xft = (1.0 / 3.0) * (Xft_ij[0] + Xft_ij[1] + Xft_ij[2]); + + ft_tensor[0] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[0] - trace_Xft); + ft_tensor[1] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[1] - trace_Xft); + ft_tensor[2] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[2] - trace_Xft); + ft_tensor[3] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[3]); + ft_tensor[4] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[4]); + ft_tensor[5] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[5]); + + for (i = 0; i < ntensors; i++) { + if (tensor_style[i] == FT) { + for (j = 0; j < 6; j++) vector[6 * i + j] = ft_tensor[j]; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +double ComputeFabric::compute_scalar() +{ + // Skip if already calculated on this timestep + if (invoked_scalar == update->ntimestep) return nc; + + invoked_scalar = update->ntimestep; + + int i, j, ii, jj, inum, jnum, itype, jtype; + tagint itag, jtag; + double xtmp, ytmp, ztmp, delx, dely, delz; + double rsq, radsum, temp_dbl; + + int *ilist, *jlist, *numneigh, **firstneigh; + + double **x = atom->x; + double *radius = atom->radius; + tagint *tag = atom->tag; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + // invoke half neighbor list (will copy or build if necessary) + neighbor->build_one(list); + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + double **cutsq = force->pair->cutsq; + + // First loop through and calculate contact fabric tensor + nc = 0; + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(mask[i] & groupbit)) continue; + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itag = tag[i]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + if (!(mask[j] & groupbit)) continue; + + // itag = jtag is possible for long cutoffs that include images of self + + if (newton_pair == 0 && j >= nlocal) { + jtag = tag[j]; + if (itag > jtag) { + if ((itag + jtag) % 2 == 0) continue; + } else if (itag < jtag) { + if ((itag + jtag) % 2 == 1) continue; + } else { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + } + } + + jtype = type[j]; + + if (type_filter) + if (type_filter[itype][jtype] == 0) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + + if (cutstyle == TYPE) { + if (rsq >= cutsq[itype][jtype]) continue; + } else { + radsum = radius[i] + radius[j]; + if (rsq >= radsum * radsum) continue; + } + + nc += 1.0; + } + } + //Count total contacts across processors + MPI_Allreduce(&nc, &temp_dbl, 1, MPI_DOUBLE, MPI_SUM, world); + nc = temp_dbl; + + scalar = nc; + return nc; +} + diff --git a/src/GRANULAR/compute_fabric.h b/src/GRANULAR/compute_fabric.h new file mode 100644 index 0000000000..5e20ea6562 --- /dev/null +++ b/src/GRANULAR/compute_fabric.h @@ -0,0 +1,79 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 COMPUTE_CLASS +// clang-format off +ComputeStyle(fabric,ComputeFabric); +// clang-format on +#else + +#ifndef LMP_COMPUTE_FABRIC_H +#define LMP_COMPUTE_FABRIC_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeFabric : public Compute { + public: + ComputeFabric(class LAMMPS *, int, char **); + ~ComputeFabric(); + void init(); + void init_list(int, class NeighList *); + void compute_vector(); + double compute_scalar(); + + private: + int ntensors, pstyle, cutstyle; + double nc; + int *tensor_style; + int **type_filter; + class NeighList *list; + + int cn_flag, br_flag, fn_flag, ft_flag; +}; + +} // namespace LAMMPS_NS + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Compute fabric radius style requires atom attribute radius + +Self-explanatory. + +E: No pair style is defined for compute fabric + +Self-explanatory. + +E: Pair style does not support compute fabric normal or tangential force + +Pair style must be single enabled to calculate the normal or tangential force tensors + +E: Pair style does not calculate tangential forces for compute fabric + +The tangential force tensor can only be calculated for granular pair styles with tangential forces + +E: Compute fabric does not support pair styles that extend beyond contact + +Granular pair styles that extend beyond contact such as JKR or DMT are not supported + +*/ diff --git a/src/KOKKOS/kokkos.cpp b/src/KOKKOS/kokkos.cpp index fd171b61af..2588bf5d49 100644 --- a/src/KOKKOS/kokkos.cpp +++ b/src/KOKKOS/kokkos.cpp @@ -69,6 +69,10 @@ GPU_AWARE_UNKNOWN using namespace LAMMPS_NS; +Kokkos::InitArguments KokkosLMP::args{-1, -1, -1, false}; +int KokkosLMP::is_finalized = 0; +int KokkosLMP::init_ngpus = 0; + /* ---------------------------------------------------------------------- */ KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) @@ -155,6 +159,10 @@ KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) } else if (strcmp(arg[iarg],"t") == 0 || strcmp(arg[iarg],"threads") == 0) { nthreads = atoi(arg[iarg+1]); + + if (nthreads <= 0) + error->all(FLERR,"Invalid number of threads requested for Kokkos: must be 1 or greater"); + iarg += 2; } else if (strcmp(arg[iarg],"n") == 0 || @@ -165,13 +173,27 @@ KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) } else error->all(FLERR,"Invalid Kokkos command-line args"); } - // initialize Kokkos + // Initialize Kokkos. However, we cannot change any + // Kokkos library parameters after the first initalization - if (me == 0) { - if (screen) fprintf(screen," will use up to %d GPU(s) per node\n",ngpus); - if (logfile) fprintf(logfile," will use up to %d GPU(s) per node\n",ngpus); + if (args.num_threads != -1) { + if (args.num_threads != nthreads || args.num_numa != numa || args.device_id != device) + if (me == 0) + error->warning(FLERR,"Kokkos package already initalized, cannot reinitialize with different parameters"); + nthreads = args.num_threads; + numa = args.num_numa; + device = args.device_id; + ngpus = init_ngpus; + } else { + args.num_threads = nthreads; + args.num_numa = numa; + args.device_id = device; + init_ngpus = ngpus; } + if (me == 0) + utils::logmesg(lmp, " will use up to {} GPU(s) per node\n",ngpus); + #ifdef LMP_KOKKOS_GPU if (ngpus <= 0) error->all(FLERR,"Kokkos has been compiled for CUDA, HIP, or SYCL but no GPUs are requested"); @@ -184,12 +206,7 @@ KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) "than the OpenMP backend"); #endif - Kokkos::InitArguments args; - args.num_threads = nthreads; - args.num_numa = numa; - args.device_id = device; - - Kokkos::initialize(args); + KokkosLMP::initialize(args,error); // default settings for package kokkos command @@ -299,9 +316,27 @@ KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) KokkosLMP::~KokkosLMP() { - // finalize Kokkos - Kokkos::finalize(); +} + +/* ---------------------------------------------------------------------- */ + +void KokkosLMP::initialize(Kokkos::InitArguments args, Error *error) +{ + if (!Kokkos::is_initialized()) { + if (is_finalized) + error->all(FLERR,"Kokkos package already finalized, cannot re-initialize"); + Kokkos::initialize(args); + } +} + +/* ---------------------------------------------------------------------- */ + +void KokkosLMP::finalize() +{ + if (Kokkos::is_initialized() && !is_finalized) + Kokkos::finalize(); + is_finalized = 1; } /* ---------------------------------------------------------------------- diff --git a/src/KOKKOS/kokkos.h b/src/KOKKOS/kokkos.h index 22060ecc09..65990544ad 100644 --- a/src/KOKKOS/kokkos.h +++ b/src/KOKKOS/kokkos.h @@ -49,8 +49,14 @@ class KokkosLMP : protected Pointers { int newtonflag; double binsize; + static int is_finalized; + static Kokkos::InitArguments args; + static int init_ngpus; + KokkosLMP(class LAMMPS *, int, char **); ~KokkosLMP(); + static void initialize(Kokkos::InitArguments, Error *); + static void finalize(); void accelerator(int, char **); int neigh_count(int); @@ -84,13 +90,21 @@ because MPI library not recognized The local MPI rank was not found in one of four supported environment variables. +E: Invalid number of threads requested for Kokkos: must be 1 or greater + +Self-explanatory. + E: GPUs are requested but Kokkos has not been compiled for CUDA Recompile Kokkos with CUDA support to use GPUs. -E: Kokkos has been compiled for CUDA but no GPUs are requested +E: Kokkos has been compiled for CUDA, HIP, or SYCL but no GPUs are requested -One or more GPUs must be used when Kokkos is compiled for CUDA. +One or more GPUs must be used when Kokkos is compiled for CUDA/HIP/SYCL. + +W: Kokkos package already initalized, cannot reinitialize with different parameters + +Self-explanatory. E: Illegal ... command diff --git a/src/SPIN/fix_nve_spin.cpp b/src/SPIN/fix_nve_spin.cpp index e6929d2d72..c64609cfe3 100644 --- a/src/SPIN/fix_nve_spin.cpp +++ b/src/SPIN/fix_nve_spin.cpp @@ -176,6 +176,7 @@ void FixNVESpin::init() // loop 1: obtain # of Pairs, and # of Pair/Spin styles + npairspin = 0; PairHybrid *hybrid = (PairHybrid *)force->pair_match("^hybrid",0); if (force->pair_match("^spin",0,0)) { // only one Pair/Spin style pair = force->pair_match("^spin",0,0); @@ -232,6 +233,7 @@ void FixNVESpin::init() // loop 1: obtain # of fix precession/spin styles int iforce; + nprecspin = 0; for (iforce = 0; iforce < modify->nfix; iforce++) { if (utils::strmatch(modify->fix[iforce]->style,"^precession/spin")) { nprecspin++; @@ -264,6 +266,7 @@ void FixNVESpin::init() // loop 1: obtain # of fix langevin/spin styles + nlangspin = 0; for (iforce = 0; iforce < modify->nfix; iforce++) { if (utils::strmatch(modify->fix[iforce]->style,"^langevin/spin")) { nlangspin++; diff --git a/src/accelerator_kokkos.h b/src/accelerator_kokkos.h index 5601e8b855..0cea855f1e 100644 --- a/src/accelerator_kokkos.h +++ b/src/accelerator_kokkos.h @@ -56,16 +56,12 @@ class KokkosLMP { KokkosLMP(class LAMMPS *, int, char **) { kokkos_exists = 0; } ~KokkosLMP() {} + static void finalize() {} void accelerator(int, char **) {} int neigh_list_kokkos(int) { return 0; } int neigh_count(int) { return 0; } }; -class Kokkos { - public: - static void finalize() {} -}; - class AtomKokkos : public Atom { public: tagint **k_special; diff --git a/src/error.cpp b/src/error.cpp index 491c72bdbb..9811a1d3eb 100644 --- a/src/error.cpp +++ b/src/error.cpp @@ -81,7 +81,7 @@ void Error::universe_all(const std::string &file, int line, const std::string &s throw LAMMPSException(mesg); #else - if (lmp->kokkos) Kokkos::finalize(); + KokkosLMP::finalize(); MPI_Finalize(); exit(1); #endif @@ -107,6 +107,7 @@ void Error::universe_one(const std::string &file, int line, const std::string &s throw LAMMPSAbortException(mesg, universe->uworld); #else + KokkosLMP::finalize(); MPI_Abort(universe->uworld,1); exit(1); // to trick "smart" compilers into believing this does not return #endif @@ -173,8 +174,8 @@ void Error::all(const std::string &file, int line, const std::string &str) if (screen && screen != stdout) fclose(screen); if (logfile) fclose(logfile); + KokkosLMP::finalize(); if (universe->nworlds > 1) MPI_Abort(universe->uworld,1); - if (lmp->kokkos) Kokkos::finalize(); MPI_Finalize(); exit(1); #endif @@ -213,6 +214,7 @@ void Error::one(const std::string &file, int line, const std::string &str) #else if (screen) fflush(screen); if (logfile) fflush(logfile); + KokkosLMP::finalize(); MPI_Abort(world,1); exit(1); // to trick "smart" compilers into believing this does not return #endif @@ -315,7 +317,7 @@ void Error::done(int status) if (screen && screen != stdout) fclose(screen); if (logfile) fclose(logfile); - if (lmp->kokkos) Kokkos::finalize(); + KokkosLMP::finalize(); MPI_Finalize(); exit(status); } diff --git a/src/library.cpp b/src/library.cpp index 436892683b..c0dcb4f328 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -19,6 +19,7 @@ #include "library.h" #include +#include "accelerator_kokkos.h" #include "atom.h" #include "atom_vec.h" #include "comm.h" @@ -333,8 +334,8 @@ The MPI standard requires that any MPI application calls do any MPI calls, MPI is still initialized internally to avoid errors accessing any MPI functions. This function should then be called right before exiting the program to wait until all (parallel) tasks are -completed and then MPI is cleanly shut down. After this function no -more MPI calls may be made. +completed and then MPI is cleanly shut down. After calling this +function no more MPI calls may be made. .. versionadded:: 18Sep2020 @@ -353,6 +354,28 @@ void lammps_mpi_finalize() } } +/* ---------------------------------------------------------------------- */ + +/** Shut down the Kokkos library environment. + * +\verbatim embed:rst + +The Kokkos library may only be initialized once during the execution of +a process. This is done automatically the first time Kokkos +functionality is used. This requires that the Kokkos environment +must be explicitly shut down after any LAMMPS instance using it is +closed (to release associated resources). +After calling this function no Kokkos functionality may be used. + +.. versionadded:: TBD + +\endverbatim */ + +void lammps_kokkos_finalize() +{ + KokkosLMP::finalize(); +} + // ---------------------------------------------------------------------- // Library functions to process commands // ---------------------------------------------------------------------- diff --git a/src/library.h b/src/library.h index 70ece96716..a337d7b510 100644 --- a/src/library.h +++ b/src/library.h @@ -94,6 +94,7 @@ void lammps_close(void *handle); void lammps_mpi_init(); void lammps_mpi_finalize(); +void lammps_kokkos_finalize(); /* ---------------------------------------------------------------------- * Library functions to process commands diff --git a/src/main.cpp b/src/main.cpp index 8399ecab45..76ae4fde09 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -14,6 +14,7 @@ #include "lammps.h" #include "input.h" +#include "accelerator_kokkos.h" #if defined(LAMMPS_EXCEPTIONS) #include "exceptions.h" #endif @@ -77,13 +78,16 @@ int main(int argc, char **argv) lammps->input->file(); delete lammps; } catch (LAMMPSAbortException &ae) { + KokkosLMP::finalize(); MPI_Abort(ae.universe, 1); } catch (LAMMPSException &e) { + KokkosLMP::finalize(); MPI_Barrier(lammps_comm); MPI_Finalize(); exit(1); } catch (fmt::format_error &fe) { fprintf(stderr, "fmt::format_error: %s\n", fe.what()); + KokkosLMP::finalize(); MPI_Abort(MPI_COMM_WORLD, 1); exit(1); } @@ -94,10 +98,12 @@ int main(int argc, char **argv) delete lammps; } catch (fmt::format_error &fe) { fprintf(stderr, "fmt::format_error: %s\n", fe.what()); + KokkosLMP::finalize(); MPI_Abort(MPI_COMM_WORLD, 1); exit(1); } #endif + KokkosLMP::finalize(); MPI_Barrier(lammps_comm); MPI_Finalize(); } diff --git a/tools/swig/lammps.i b/tools/swig/lammps.i index ec37120c07..56547dda53 100644 --- a/tools/swig/lammps.i +++ b/tools/swig/lammps.i @@ -63,6 +63,7 @@ extern void *lammps_open_fortran(int argc, char **argv, int f_comm); extern void lammps_close(void *handle); extern void lammps_mpi_init(); extern void lammps_mpi_finalize(); +extern void lammps_kokkos_finalize(); extern void lammps_file(void *handle, const char *file); extern char *lammps_command(void *handle, const char *cmd); extern void lammps_commands_list(void *handle, int ncmd, const char **cmds); @@ -185,6 +186,7 @@ extern void *lammps_open_fortran(int argc, char **argv, int f_comm); extern void lammps_close(void *handle); extern void lammps_mpi_init(); extern void lammps_mpi_finalize(); +extern void lammps_kokkos_finalize(); extern void lammps_file(void *handle, const char *file); extern char *lammps_command(void *handle, const char *cmd); extern void lammps_commands_list(void *handle, int ncmd, const char **cmds); diff --git a/unittest/python/python-capabilities.py b/unittest/python/python-capabilities.py index a0b8bd2618..8f72c39670 100644 --- a/unittest/python/python-capabilities.py +++ b/unittest/python/python-capabilities.py @@ -165,6 +165,13 @@ class PythonCapabilities(unittest.TestCase): if self.cmake_cache['GPU_PREC'].lower() == 'single': self.assertIn('single',settings['GPU']['precision']) + if self.cmake_cache['PKG_KOKKOS']: + if self.cmake_cache['Kokkos_ENABLE_OPENMP']: + self.assertIn('openmp',settings['KOKKOS']['api']) + if self.cmake_cache['Kokkos_ENABLE_SERIAL']: + self.assertIn('serial',settings['KOKKOS']['api']) + self.assertIn('double',settings['KOKKOS']['precision']) + def test_gpu_device(self): info = self.lmp.get_gpu_device_info()