diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index 453067130d..5f79015aaf 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -172,12 +172,14 @@ OPT. * :doc:`lj/cut/dipole/long (g) ` * :doc:`lj/cut/dipole/sf (go) ` * :doc:`lj/cut/soft (o) ` + * :doc:`lj/cut/sphere (o) ` * :doc:`lj/cut/thole/long (o) ` * :doc:`lj/cut/tip4p/cut (o) ` * :doc:`lj/cut/tip4p/long (got) ` * :doc:`lj/cut/tip4p/long/soft (o) ` * :doc:`lj/expand (gko) ` * :doc:`lj/expand/coul/long (gk) ` + * :doc:`lj/expand/sphere (o) ` * :doc:`lj/gromacs (gko) ` * :doc:`lj/gromacs/coul/gromacs (ko) ` * :doc:`lj/long/coul/long (iot) ` @@ -186,7 +188,6 @@ OPT. * :doc:`lj/long/tip4p/long (o) ` * :doc:`lj/mdf ` * :doc:`lj/relres (o) ` - * :doc:`lj/sphere (o) ` * :doc:`lj/spica (gko) ` * :doc:`lj/spica/coul/long (go) ` * :doc:`lj/spica/coul/msm (o) ` diff --git a/doc/src/pair_lj_sphere.rst b/doc/src/pair_lj_cut_sphere.rst similarity index 84% rename from doc/src/pair_lj_sphere.rst rename to doc/src/pair_lj_cut_sphere.rst index 6f030e506a..1987ac895f 100644 --- a/doc/src/pair_lj_sphere.rst +++ b/doc/src/pair_lj_cut_sphere.rst @@ -1,10 +1,10 @@ -.. index:: pair_style lj/sphere -.. index:: pair_style lj/sphere/omp +.. index:: pair_style lj/cut/sphere +.. index:: pair_style lj/cut/sphere/omp -pair_style lj/sphere command -============================ +pair_style lj/cut/sphere command +================================ -Accelerator Variant: *lj/sphere/omp* +Accelerator Variant: *lj/cut/sphere/omp* Syntax """""" @@ -13,12 +13,12 @@ Syntax pair_style style args -* style = *lj/sphere* +* style = *lj/cut/sphere* * args = list of arguments for a particular style .. parsed-literal:: - *lj/sphere* args = cutoff ratio + *lj/cut/sphere* args = cutoff ratio cutoff = global cutoff ratio for Lennard Jones interactions (unitless) Examples @@ -26,14 +26,14 @@ Examples .. code-block:: LAMMPS - pair_style lj/sphere 2.5 + pair_style lj/cut/sphere 2.5 pair_coeff * * 1.0 pair_coeff 1 1 1.1 2.8 Description """"""""""" -The *lj/sphere* style compute the standard 12/6 Lennard-Jones potential, +The *lj/cut/sphere* style compute the standard 12/6 Lennard-Jones potential, given by .. math:: @@ -106,7 +106,7 @@ is at :math:`2^{\frac{1}{6}} \sigma_{ij}`. group large variable large set group largea type 2 - pair_style lj/sphere 2.5 + pair_style lj/cut/sphere 2.5 pair_coeff * * 1.0 neighbor 0.3 multi @@ -139,18 +139,18 @@ Mixing, shift, table, tail correction, restart, rRESPA info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" For atom type pairs I,J and I != J, the epsilon coefficients and cutoff -ratio for the *lj/sphere* pair style can be mixed. The default mixing +ratio for the *lj/cut/sphere* pair style can be mixed. The default mixing style is *geometric*. See the :doc:`pair_modify command ` for details. -The *lj/sphere* pair style supports the :doc:`pair_modify shift ` +The *lj/cut/sphere* pair style supports the :doc:`pair_modify shift ` option for the energy of the Lennard-Jones portion of the pair interaction. -The *lj/sphere* pair style does *not* support the :doc:`pair_modify +The *lj/cut/sphere* pair style does *not* support the :doc:`pair_modify ` tail option for adding a long-range tail corrections to the energy and pressure. -The *lj/sphere* pair style writes its information to :doc:`binary +The *lj/cut/sphere* pair style writes its information to :doc:`binary restart files `, so pair_style and pair_coeff commands do not need to be specified in an input script that reads a restart file. @@ -163,11 +163,11 @@ This pair style can only be used via the *pair* keyword of the Restrictions """""""""""" -The *lj/sphere* pair style is only enabled if LAMMPS was built with the +The *lj/cut/sphere* pair style is only enabled if LAMMPS was built with the EXTRA-PAIR package. See the :doc:`Build package ` page for more info. -The *lj/sphere* pair style does not support the *sixthpower* mixing rule. +The *lj/cut/sphere* pair style does not support the *sixthpower* mixing rule. ---------- @@ -176,6 +176,7 @@ Related commands * :doc:`pair_coeff ` * :doc:`pair_style lj/cut ` +* :doc:`pair_style lj/expnd/sphere ` Default """"""" diff --git a/doc/src/pair_lj_expand_sphere.rst b/doc/src/pair_lj_expand_sphere.rst new file mode 100644 index 0000000000..0eb3d497af --- /dev/null +++ b/doc/src/pair_lj_expand_sphere.rst @@ -0,0 +1,184 @@ +.. index:: pair_style lj/expand/sphere +.. index:: pair_style lj/expand/sphere/omp + +pair_style lj/expand/sphere command +=================================== + +Accelerator Variant: *lj/expand/sphere/omp* + +Syntax +"""""" + +.. code-block:: LAMMPS + + pair_style style args + +* style = *lj/expand/sphere* +* args = list of arguments for a particular style + +.. parsed-literal:: + + *lj/expand/sphere* args = cutoff ratio + cutoff = global cutoff ratio for Lennard Jones interactions (unitless) + +Examples +"""""""" + +.. code-block:: LAMMPS + + pair_style lj/expand/sphere 2.5 + pair_coeff * * 1.0 + pair_coeff 1 1 1.1 2.8 + +Description +""""""""""" + +The *lj/expand/sphere* style compute the standard 12/6 Lennard-Jones potential, +given by + +.. math:: + + E = 4 \epsilon \left[ \left(\frac{\sigma_{ij}}{r}\right)^{12} - + \left(\frac{\sigma_{ij}}{r}\right)^6 \right] \qquad r < r_c * \sigma_{ij} + +:math:`r_c` is the cutoff ratio. + +This is the same potential function used by the :doc:`lj/cut +` pair style, but the :math:`\sigma_{ij}` parameter is not +set as a per-type parameter via the :doc:`pair_coeff command +`. Instead it is calculated individually for each pair +using the per-atom diameter attribute of :doc:`atom_style sphere +` for the two atoms as :math:`\sigma_{i}` and +:math:`\sigma_{j}`; :math:`\sigma_{ij}` is then computed by the mixing +rule for pair coefficients as set by the :doc:`pair_modify mix +` command (defaults to geometric mixing). The cutoff is +not specified as a distance, but as ratio that is internally +multiplied by :math:`\sigma_{ij}` to obtain the actual cutoff for each +pair of atoms. + +Note that :math:`\sigma_{ij}` is defined in the LJ formula above as the +zero-crossing distance for the potential, *not* as the energy minimum which +is at :math:`2^{\frac{1}{6}} \sigma_{ij}`. + +.. admonition:: Notes on cutoffs, neighbor lists, and efficiency + :class: note + + If your system is mildly polydisperse, meaning the ratio of the + diameter of the largest particle to the smallest is less than 2.0, + then the neighbor lists built by the code should be resonably + efficient. Which means they will not contain too many particle + pairs that do not interact. However, if your system is highly + polydisperse (ratio > 2.0), the neighbor list build and force + computations may be inefficient. There are two ways to try and + speed up the computations. + + The first is to assign multiple atom types so that atoms of each + type are similar in size. E.g. if particle diameters range from 1 + to 5 use 4 atom types, ensuring atoms of type 1 have diameters from + 1.0-2.0, type 2 from 2.0-3.0, etc. + + The second is to try the :doc:`neighbor multi ` command + which uses a different algorithm for buliding neighbor lists. This + will also require that you assign multiple atom types as in the + preceeding paragraph. + + Here are example input script commands using both ideas for a + highly polydisperse system: + + .. code-block:: c++ + + units lj + atom_style sphere + lattice fcc 0.8442 + region box block 0 10 0 10 0 10 + create_box 2 box + create_atoms 1 box + + # create atoms with random diameters from bimodal distribution + + variable switch atom random(0.0,1.0,345634) + variable diam atom (v_switch<0.75)*normal(0.4,0.075,325)+(v_switch>=0.7)*normal(1.2,0.2,453) + set group all diameter v_diam + + # assign type 2 to atoms with diameter > 0.5 + + variable large atom 2.0*radius>0.5 + group large variable large + set group largea type 2 + + pair_style lj/expand/sphere 2.5 + pair_coeff * * 1.0 + + neighbor 0.3 multi + +Coefficients +"""""""""""" + +The following coefficients must be defined for each pair of atoms types via the +:doc:`pair_coeff ` command as in the examples above, or in the data +file or restart files read by the :doc:`read_data ` or +:doc:`read_restart ` commands, or by mixing as described below: + +* :math:`\epsilon` (energy units) +* LJ cutoff ratio (unitless) (optional) + +The last coefficient is optional. If not specified, the global LJ +cutoff ratio specified in the :doc:`pair_style command ` is +used. + +If a repulsive only LJ interaction is desired, the coefficient for the cutoff +ratio should be set to the minimum of the LJ potential using ``$(2.0^(1.0/6.0))`` + +---------- + +.. include:: accel_styles.rst + +---------- + +Mixing, shift, table, tail correction, restart, rRESPA info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +For atom type pairs I,J and I != J, the epsilon coefficients and cutoff +ratio for the *lj/expand/sphere* pair style can be mixed. The default mixing +style is *geometric*. See the :doc:`pair_modify command ` +for details. + +The *lj/expand/sphere* pair style supports the :doc:`pair_modify shift ` +option for the energy of the Lennard-Jones portion of the pair interaction. + +The *lj/expand/sphere* pair style does *not* support the :doc:`pair_modify +` tail option for adding a long-range tail corrections to +the energy and pressure. + +The *lj/expand/sphere* pair style writes its information to :doc:`binary +restart files `, so pair_style and pair_coeff commands do not +need to be specified in an input script that reads a restart file. + +This pair style can only be used via the *pair* keyword of the +:doc:`run_style respa ` command. It does *not* support the +*inner*, *middle*, *outer* keywords. + +---------- + +Restrictions +"""""""""""" + +The *lj/expand/sphere* pair style is only enabled if LAMMPS was built with the +EXTRA-PAIR package. See the :doc:`Build package ` page +for more info. + +The *lj/expand/sphere* pair style does not support the *sixthpower* mixing rule. + +---------- + +Related commands +"""""""""""""""" + +* :doc:`pair_coeff ` +* :doc:`pair_style lj/cut ` +* :doc:`pair_style lj/cut/sphere ` + +Default +""""""" + +none diff --git a/doc/src/pair_style.rst b/doc/src/pair_style.rst index 96488ee154..49e20a93b0 100644 --- a/doc/src/pair_style.rst +++ b/doc/src/pair_style.rst @@ -250,12 +250,14 @@ accelerated styles exist. * :doc:`lj/cut/dipole/cut ` - point dipoles with cutoff * :doc:`lj/cut/dipole/long ` - point dipoles with long-range Ewald * :doc:`lj/cut/soft ` - LJ with a soft core +* :doc:`lj/cut/sphere ` - LJ where per-atom radius is used as LJ sigma * :doc:`lj/cut/thole/long ` - LJ with Coulomb with thole damping * :doc:`lj/cut/tip4p/cut ` - LJ with cutoff Coulomb for TIP4P water * :doc:`lj/cut/tip4p/long ` - LJ with long-range Coulomb for TIP4P water * :doc:`lj/cut/tip4p/long/soft ` - LJ with cutoff Coulomb for TIP4P water with a soft core * :doc:`lj/expand ` - Lennard-Jones for variable size particles * :doc:`lj/expand/coul/long ` - Lennard-Jones for variable size particles with long-range Coulomb +* :doc:`lj/expand/sphere ` - Variable size LJ where per-atom radius is used as delta (size) * :doc:`lj/gromacs ` - GROMACS-style Lennard-Jones potential * :doc:`lj/gromacs/coul/gromacs ` - GROMACS-style LJ and Coulomb potential * :doc:`lj/long/coul/long ` - long-range LJ and long-range Coulomb @@ -264,7 +266,6 @@ accelerated styles exist. * :doc:`lj/long/tip4p/long ` - long-range LJ and long-range Coulomb for TIP4P water * :doc:`lj/mdf ` - LJ potential with a taper function * :doc:`lj/relres ` - LJ using multiscale Relative Resolution (RelRes) methodology :ref:`(Chaimovich) `. -* :doc:`lj/sphere ` - LJ where per-atom radius is used as LJ sigma * :doc:`lj/spica ` - LJ for SPICA coarse-graining * :doc:`lj/spica/coul/long ` - LJ for SPICA coarse-graining with long-range Coulomb * :doc:`lj/spica/coul/msm ` - LJ for SPICA coarse-graining with long-range Coulomb via MSM diff --git a/src/.gitignore b/src/.gitignore index 9988ace999..a0d3d63de9 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -1249,6 +1249,8 @@ /pair_lj_cut_dipole_long.h /pair_lj_cut_*hars_*.cpp /pair_lj_cut_*hars_*.h +/pair_lj_cut_sphere.cpp +/pair_lj_cut_sphere.h /pair_lj_cut_soft.cpp /pair_lj_cut_soft.h /pair_lj_cut_tip4p_long.cpp @@ -1263,6 +1265,8 @@ /pair_lj_long_tip4p_long.h /pair_lj_cut_tgpu.cpp /pair_lj_cut_tgpu.h +/pair_lj_expand_sphere.cpp +/pair_lj_expand_sphere.h /pair_lj_spica.cpp /pair_lj_spica.h /pair_lj_spica_coul_long.cpp diff --git a/src/EXTRA-PAIR/pair_lj_sphere.cpp b/src/EXTRA-PAIR/pair_lj_cut_sphere.cpp similarity index 90% rename from src/EXTRA-PAIR/pair_lj_sphere.cpp rename to src/EXTRA-PAIR/pair_lj_cut_sphere.cpp index 2d54d96bb0..3baff2252f 100644 --- a/src/EXTRA-PAIR/pair_lj_sphere.cpp +++ b/src/EXTRA-PAIR/pair_lj_cut_sphere.cpp @@ -11,7 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include "pair_lj_sphere.h" +#include "pair_lj_cut_sphere.h" #include "atom.h" #include "comm.h" @@ -32,14 +32,15 @@ using MathSpecial::square; /* ---------------------------------------------------------------------- */ -PairLJSphere::PairLJSphere(LAMMPS *lmp) : Pair(lmp), rmax(nullptr), cut(nullptr), epsilon(nullptr) +PairLJCutSphere::PairLJCutSphere(LAMMPS *lmp) : + Pair(lmp), rmax(nullptr), cut(nullptr), epsilon(nullptr) { writedata = 1; } /* ---------------------------------------------------------------------- */ -PairLJSphere::~PairLJSphere() +PairLJCutSphere::~PairLJCutSphere() { if (allocated) { memory->destroy(setflag); @@ -53,7 +54,7 @@ PairLJSphere::~PairLJSphere() /* ---------------------------------------------------------------------- */ -void PairLJSphere::compute(int eflag, int vflag) +void PairLJCutSphere::compute(int eflag, int vflag) { int i, j, ii, jj, inum, jnum, itype, jtype; double xtmp, ytmp, ztmp, rtmp, delx, dely, delz, evdwl, sigma, sigma6, fpair; @@ -145,7 +146,7 @@ void PairLJSphere::compute(int eflag, int vflag) allocate all arrays ------------------------------------------------------------------------- */ -void PairLJSphere::allocate() +void PairLJCutSphere::allocate() { allocated = 1; int n = atom->ntypes + 1; @@ -165,7 +166,7 @@ void PairLJSphere::allocate() global settings ------------------------------------------------------------------------- */ -void PairLJSphere::settings(int narg, char **arg) +void PairLJCutSphere::settings(int narg, char **arg) { if (narg != 1) error->all(FLERR, "Illegal pair_style command"); @@ -185,7 +186,7 @@ void PairLJSphere::settings(int narg, char **arg) set coeffs for one or more type pairs ------------------------------------------------------------------------- */ -void PairLJSphere::coeff(int narg, char **arg) +void PairLJCutSphere::coeff(int narg, char **arg) { if (narg < 3 || narg > 4) error->all(FLERR, "Incorrect args for pair coefficients"); if (!allocated) allocate(); @@ -215,13 +216,14 @@ void PairLJSphere::coeff(int narg, char **arg) init specific to this pair style ------------------------------------------------------------------------- */ -void PairLJSphere::init_style() +void PairLJCutSphere::init_style() { Pair::init_style(); - if (!atom->radius_flag) error->all(FLERR, "Pair style lj/sphere requires atom attribute radius"); + if (!atom->radius_flag) + error->all(FLERR, "Pair style lj/cut/sphere requires atom attribute radius"); if (mix_flag == SIXTHPOWER) - error->all(FLERR, "Pair_modify mix sixthpower is not compatible with pair style lj/sphere"); + error->all(FLERR, "Pair_modify mix sixthpower is not compatible with pair style lj/cut/sphere"); // determine max radius per type @@ -241,7 +243,7 @@ void PairLJSphere::init_style() init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ -double PairLJSphere::init_one(int i, int j) +double PairLJCutSphere::init_one(int i, int j) { if (setflag[i][j] == 0) { epsilon[i][j] = mix_energy(epsilon[i][i], epsilon[j][j], 0.0, 0.0); @@ -260,7 +262,7 @@ double PairLJSphere::init_one(int i, int j) proc 0 writes to restart file ------------------------------------------------------------------------- */ -void PairLJSphere::write_restart(FILE *fp) +void PairLJCutSphere::write_restart(FILE *fp) { write_restart_settings(fp); @@ -279,7 +281,7 @@ void PairLJSphere::write_restart(FILE *fp) proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ -void PairLJSphere::read_restart(FILE *fp) +void PairLJCutSphere::read_restart(FILE *fp) { read_restart_settings(fp); allocate(); @@ -305,7 +307,7 @@ void PairLJSphere::read_restart(FILE *fp) proc 0 writes to restart file ------------------------------------------------------------------------- */ -void PairLJSphere::write_restart_settings(FILE *fp) +void PairLJCutSphere::write_restart_settings(FILE *fp) { fwrite(&cut_global, sizeof(double), 1, fp); fwrite(&offset_flag, sizeof(int), 1, fp); @@ -316,7 +318,7 @@ void PairLJSphere::write_restart_settings(FILE *fp) proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ -void PairLJSphere::read_restart_settings(FILE *fp) +void PairLJCutSphere::read_restart_settings(FILE *fp) { int me = comm->me; if (me == 0) { @@ -333,7 +335,7 @@ void PairLJSphere::read_restart_settings(FILE *fp) proc 0 writes to data file ------------------------------------------------------------------------- */ -void PairLJSphere::write_data(FILE *fp) +void PairLJCutSphere::write_data(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g\n", i, epsilon[i][i]); } @@ -342,7 +344,7 @@ void PairLJSphere::write_data(FILE *fp) proc 0 writes all pairs to data file ------------------------------------------------------------------------- */ -void PairLJSphere::write_data_all(FILE *fp) +void PairLJCutSphere::write_data_all(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) @@ -351,8 +353,8 @@ void PairLJSphere::write_data_all(FILE *fp) /* ---------------------------------------------------------------------- */ -double PairLJSphere::single(int i, int j, int itype, int jtype, double rsq, double /*factor_coul*/, - double factor_lj, double &fforce) +double PairLJCutSphere::single(int i, int j, int itype, int jtype, double rsq, + double /*factor_coul*/, double factor_lj, double &fforce) { double r2inv, r6inv, rcutsq, sigma, sigma6, forcelj, philj; @@ -374,7 +376,7 @@ double PairLJSphere::single(int i, int j, int itype, int jtype, double rsq, doub /* ---------------------------------------------------------------------- */ -void *PairLJSphere::extract(const char *str, int &dim) +void *PairLJCutSphere::extract(const char *str, int &dim) { dim = 2; if (strcmp(str, "epsilon") == 0) return (void *) epsilon; diff --git a/src/EXTRA-PAIR/pair_lj_sphere.h b/src/EXTRA-PAIR/pair_lj_cut_sphere.h similarity index 87% rename from src/EXTRA-PAIR/pair_lj_sphere.h rename to src/EXTRA-PAIR/pair_lj_cut_sphere.h index 55932e5db2..4f2fc9785b 100644 --- a/src/EXTRA-PAIR/pair_lj_sphere.h +++ b/src/EXTRA-PAIR/pair_lj_cut_sphere.h @@ -13,21 +13,21 @@ #ifdef PAIR_CLASS // clang-format off -PairStyle(lj/sphere,PairLJSphere); +PairStyle(lj/cut/sphere,PairLJCutSphere); // clang-format on #else -#ifndef LMP_PAIR_LJ_SPHERE_H -#define LMP_PAIR_LJ_SPHERE_H +#ifndef LMP_PAIR_LJ_CUT_SPHERE_H +#define LMP_PAIR_LJ_CUT_SPHERE_H #include "pair.h" namespace LAMMPS_NS { -class PairLJSphere : public Pair { +class PairLJCutSphere : public Pair { public: - PairLJSphere(class LAMMPS *); - ~PairLJSphere() override; + PairLJCutSphere(class LAMMPS *); + ~PairLJCutSphere() override; void compute(int, int) override; void settings(int, char **) override; void coeff(int, char **) override; diff --git a/src/EXTRA-PAIR/pair_lj_expand_sphere.cpp b/src/EXTRA-PAIR/pair_lj_expand_sphere.cpp new file mode 100644 index 0000000000..d0b6bc0258 --- /dev/null +++ b/src/EXTRA-PAIR/pair_lj_expand_sphere.cpp @@ -0,0 +1,385 @@ +/* ---------------------------------------------------------------------- + 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 "pair_lj_expand_sphere.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "math_special.h" +#include "memory.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; +using MathSpecial::powint; +using MathSpecial::square; + +/* ---------------------------------------------------------------------- */ + +PairLJExpandSphere::PairLJExpandSphere(LAMMPS *lmp) : + Pair(lmp), rmax(nullptr), cut(nullptr), epsilon(nullptr) +{ + writedata = 1; +} + +/* ---------------------------------------------------------------------- */ + +PairLJExpandSphere::~PairLJExpandSphere() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(rmax); + memory->destroy(cut); + memory->destroy(epsilon); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLJExpandSphere::compute(int eflag, int vflag) +{ + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, rtmp, delx, dely, delz, evdwl, sigma, sigma6, fpair; + double rcutsq, rsq, r2inv, r6inv, forcelj, factor_lj; + int *ilist, *jlist, *numneigh, **firstneigh; + + evdwl = 0.0; + ev_init(eflag, vflag); + + double **x = atom->x; + double **f = atom->f; + double *radius = atom->radius; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + rtmp = radius[i]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + + // cutsq is maximum cutoff per type. Now compute and apply real cutoff + + sigma = 2.0 * mix_distance(rtmp, radius[j]); + rcutsq = square(cut[itype][jtype] * sigma); + if (rsq < rcutsq) { + + r2inv = 1.0 / rsq; + r6inv = r2inv * r2inv * r2inv; + sigma6 = powint(sigma, 6); + forcelj = r6inv * 24.0 * epsilon[itype][jtype] * (2.0 * sigma6 * sigma6 * r6inv - sigma6); + fpair = factor_lj * forcelj * r2inv; + + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; + } + + if (eflag) { + evdwl = r6inv * 4.0 * epsilon[itype][jtype]; + evdwl *= sigma6 * sigma6 * r6inv - sigma6; + if (offset_flag && (rcutsq > 0.0)) { + double ratio6 = sigma6 / powint(rcutsq, 3); + evdwl -= 4.0 * epsilon[itype][jtype] * (ratio6 * ratio6 - ratio6); + } + evdwl *= factor_lj; + } + + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz); + } + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairLJExpandSphere::allocate() +{ + allocated = 1; + int n = atom->ntypes + 1; + + memory->create(setflag, n, n, "pair:setflag"); + for (int i = 1; i < n; i++) + for (int j = i; j < n; j++) setflag[i][j] = 0; + + memory->create(cutsq, n, n, "pair:cutsq"); + + memory->create(rmax, n, "pair:rmax"); + memory->create(cut, n, n, "pair:cut"); + memory->create(epsilon, n, n, "pair:epsilon"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairLJExpandSphere::settings(int narg, char **arg) +{ + if (narg != 1) error->all(FLERR, "Illegal pair_style command"); + + cut_global = utils::numeric(FLERR, arg[0], false, lmp); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i, j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) + if (setflag[i][j]) cut[i][j] = cut_global; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairLJExpandSphere::coeff(int narg, char **arg) +{ + if (narg < 3 || narg > 4) error->all(FLERR, "Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo, ihi, jlo, jhi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); + + double epsilon_one = utils::numeric(FLERR, arg[2], false, lmp); + double cut_one = cut_global; + if (narg == 5) cut_one = utils::numeric(FLERR, arg[4], false, lmp); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { + epsilon[i][j] = epsilon_one; + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairLJExpandSphere::init_style() +{ + Pair::init_style(); + + if (!atom->radius_flag) + error->all(FLERR, "Pair style lj/expand/sphere requires atom attribute radius"); + if (mix_flag == SIXTHPOWER) + error->all(FLERR, + "Pair_modify mix sixthpower is not compatible with pair style lj/expand/sphere"); + + // determine max radius per type + + int *type = atom->type; + double *radius = atom->radius; + rmax[0] = 0.0; + for (int itype = 1; itype <= atom->ntypes; ++itype) { + double rmax_one = 0.0; + for (int i = 0; i < atom->nlocal; ++i) { + if (type[i] == itype) rmax_one = MAX(rmax_one, radius[i]); + } + MPI_Allreduce(&rmax_one, &rmax[itype], 1, MPI_DOUBLE, MPI_MAX, world); + } +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairLJExpandSphere::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + epsilon[i][j] = mix_energy(epsilon[i][i], epsilon[j][j], 0.0, 0.0); + cut[i][j] = mix_distance(cut[i][i], cut[j][j]); + } + + epsilon[j][i] = epsilon[i][j]; + cut[j][i] = cut[i][j]; + + // since cut is a scaled by the mixed diameter, report maximum possible cutoff. + + return cut[i][j] * 2.0 * mix_distance(rmax[i], rmax[j]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJExpandSphere::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i, j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j], sizeof(int), 1, fp); + if (setflag[i][j]) { + fwrite(&epsilon[i][j], sizeof(double), 1, fp); + fwrite(&cut[i][j], sizeof(double), 1, fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJExpandSphere::read_restart(FILE *fp) +{ + read_restart_settings(fp); + allocate(); + + int i, j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); + if (setflag[i][j]) { + if (me == 0) { + utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error); + } + MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJExpandSphere::write_restart_settings(FILE *fp) +{ + fwrite(&cut_global, sizeof(double), 1, fp); + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJExpandSphere::read_restart_settings(FILE *fp) +{ + int me = comm->me; + if (me == 0) { + utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); + } + MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void PairLJExpandSphere::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g\n", i, epsilon[i][i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairLJExpandSphere::write_data_all(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i; j <= atom->ntypes; j++) + fprintf(fp, "%d %d %g %g\n", i, j, epsilon[i][j], cut[i][j]); +} + +/* ---------------------------------------------------------------------- */ + +double PairLJExpandSphere::single(int i, int j, int itype, int jtype, double rsq, + double /*factor_coul*/, double factor_lj, double &fforce) +{ + double r2inv, r6inv, rcutsq, sigma, sigma6, forcelj, philj; + + sigma = 2.0 * mix_distance(atom->radius[i], atom->radius[j]); + rcutsq = square(cut[itype][jtype] * sigma); + sigma6 = powint(sigma, 6); + r2inv = 1.0 / rsq; + r6inv = r2inv * r2inv * r2inv; + forcelj = r6inv * 24.0 * epsilon[itype][jtype] * (sigma6 * sigma6 * r6inv - sigma6); + fforce = factor_lj * forcelj * r2inv; + + philj = r6inv * 4.0 * epsilon[itype][jtype] * (sigma6 * sigma6 * r6inv - sigma6); + if (offset_flag && (rcutsq > 0.0)) { + double ratio6 = sigma6 / powint(rcutsq, 3); + philj -= 4.0 * epsilon[itype][jtype] * (ratio6 * ratio6 - ratio6); + } + return factor_lj * philj; +} + +/* ---------------------------------------------------------------------- */ + +void *PairLJExpandSphere::extract(const char *str, int &dim) +{ + dim = 2; + if (strcmp(str, "epsilon") == 0) return (void *) epsilon; + return nullptr; +} diff --git a/src/EXTRA-PAIR/pair_lj_expand_sphere.h b/src/EXTRA-PAIR/pair_lj_expand_sphere.h new file mode 100644 index 0000000000..48d4e9c754 --- /dev/null +++ b/src/EXTRA-PAIR/pair_lj_expand_sphere.h @@ -0,0 +1,58 @@ +/* -*- 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 PAIR_CLASS +// clang-format off +PairStyle(lj/expand/sphere,PairLJExpandSphere); +// clang-format on +#else + +#ifndef LMP_PAIR_LJ_EXPAND_SPHERE_H +#define LMP_PAIR_LJ_EXPAND_SPHERE_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairLJExpandSphere : public Pair { + public: + PairLJExpandSphere(class LAMMPS *); + ~PairLJExpandSphere() override; + void compute(int, int) override; + void settings(int, char **) override; + void coeff(int, char **) override; + void init_style() override; + double init_one(int, int) override; + void write_restart(FILE *) override; + void read_restart(FILE *) override; + void write_restart_settings(FILE *) override; + void read_restart_settings(FILE *) override; + void write_data(FILE *) override; + void write_data_all(FILE *) override; + double single(int, int, int, int, double, double, double, double &) override; + void *extract(const char *, int &) override; + + protected: + double cut_global; + double *rmax; + double **cut; + double **epsilon; + double **sigma; + + virtual void allocate(); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/OPENMP/pair_lj_sphere_omp.cpp b/src/OPENMP/pair_lj_cut_sphere_omp.cpp similarity index 93% rename from src/OPENMP/pair_lj_sphere_omp.cpp rename to src/OPENMP/pair_lj_cut_sphere_omp.cpp index bb6d4620cc..e71f5add59 100644 --- a/src/OPENMP/pair_lj_sphere_omp.cpp +++ b/src/OPENMP/pair_lj_cut_sphere_omp.cpp @@ -12,7 +12,7 @@ Contributing author: Axel Kohlmeyer (Temple U) ------------------------------------------------------------------------- */ -#include "pair_lj_sphere_omp.h" +#include "pair_lj_cut_sphere_omp.h" #include "atom.h" #include "comm.h" @@ -28,14 +28,14 @@ using MathSpecial::square; /* ---------------------------------------------------------------------- */ -PairLJSphereOMP::PairLJSphereOMP(LAMMPS *lmp) : PairLJSphere(lmp), ThrOMP(lmp, THR_PAIR) +PairLJCutSphereOMP::PairLJCutSphereOMP(LAMMPS *lmp) : PairLJCutSphere(lmp), ThrOMP(lmp, THR_PAIR) { suffix_flag |= Suffix::OMP; } /* ---------------------------------------------------------------------- */ -void PairLJSphereOMP::compute(int eflag, int vflag) +void PairLJCutSphereOMP::compute(int eflag, int vflag) { ev_init(eflag, vflag); @@ -78,7 +78,7 @@ void PairLJSphereOMP::compute(int eflag, int vflag) } template -void PairLJSphereOMP::eval(int iifrom, int iito, ThrData *const thr) +void PairLJCutSphereOMP::eval(int iifrom, int iito, ThrData *const thr) { const auto *_noalias const x = (dbl3_t *) atom->x[0]; auto *_noalias const f = (dbl3_t *) thr->get_f()[0]; @@ -172,10 +172,10 @@ void PairLJSphereOMP::eval(int iifrom, int iito, ThrData *const thr) /* ---------------------------------------------------------------------- */ -double PairLJSphereOMP::memory_usage() +double PairLJCutSphereOMP::memory_usage() { double bytes = memory_usage_thr(); - bytes += PairLJSphere::memory_usage(); + bytes += PairLJCutSphere::memory_usage(); return bytes; } diff --git a/src/OPENMP/pair_lj_sphere_omp.h b/src/OPENMP/pair_lj_cut_sphere_omp.h similarity index 82% rename from src/OPENMP/pair_lj_sphere_omp.h rename to src/OPENMP/pair_lj_cut_sphere_omp.h index 1f6ac5a0f1..7fe7044229 100644 --- a/src/OPENMP/pair_lj_sphere_omp.h +++ b/src/OPENMP/pair_lj_cut_sphere_omp.h @@ -17,22 +17,22 @@ #ifdef PAIR_CLASS // clang-format off -PairStyle(lj/sphere/omp,PairLJSphereOMP); +PairStyle(lj/cut/sphere/omp,PairLJCutSphereOMP); // clang-format on #else -#ifndef LMP_PAIR_LJ_SPHERE_OMP_H -#define LMP_PAIR_LJ_SPHERE_OMP_H +#ifndef LMP_PAIR_LJ_CUT_SPHERE_OMP_H +#define LMP_PAIR_LJ_CUT_SPHERE_OMP_H -#include "pair_lj_sphere.h" +#include "pair_lj_cut_sphere.h" #include "thr_omp.h" namespace LAMMPS_NS { -class PairLJSphereOMP : public PairLJSphere, public ThrOMP { +class PairLJCutSphereOMP : public PairLJCutSphere, public ThrOMP { public: - PairLJSphereOMP(class LAMMPS *); + PairLJCutSphereOMP(class LAMMPS *); void compute(int, int) override; double memory_usage() override; diff --git a/src/OPENMP/pair_lj_expand_sphere_omp.cpp b/src/OPENMP/pair_lj_expand_sphere_omp.cpp new file mode 100644 index 0000000000..309e4b9df3 --- /dev/null +++ b/src/OPENMP/pair_lj_expand_sphere_omp.cpp @@ -0,0 +1,181 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + This software is distributed under the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#include "pair_lj_expand_sphere_omp.h" + +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "math_special.h" +#include "neigh_list.h" +#include "suffix.h" + +#include "omp_compat.h" +using namespace LAMMPS_NS; +using MathSpecial::powint; +using MathSpecial::square; + +/* ---------------------------------------------------------------------- */ + +PairLJExpandSphereOMP::PairLJExpandSphereOMP(LAMMPS *lmp) : PairLJExpandSphere(lmp), ThrOMP(lmp, THR_PAIR) +{ + suffix_flag |= Suffix::OMP; +} + +/* ---------------------------------------------------------------------- */ + +void PairLJExpandSphereOMP::compute(int eflag, int vflag) +{ + ev_init(eflag, vflag); + + const int nall = atom->nlocal + atom->nghost; + const int nthreads = comm->nthreads; + const int inum = list->inum; + +#if defined(_OPENMP) +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag, vflag) +#endif + { + int ifrom, ito, tid; + + loop_setup_thr(ifrom, ito, tid, inum, nthreads); + ThrData *thr = fix->get_thr(tid); + thr->timer(Timer::START); + ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr); + + if (evflag) { + if (eflag) { + if (force->newton_pair) + eval<1, 1, 1>(ifrom, ito, thr); + else + eval<1, 1, 0>(ifrom, ito, thr); + } else { + if (force->newton_pair) + eval<1, 0, 1>(ifrom, ito, thr); + else + eval<1, 0, 0>(ifrom, ito, thr); + } + } else { + if (force->newton_pair) + eval<0, 0, 1>(ifrom, ito, thr); + else + eval<0, 0, 0>(ifrom, ito, thr); + } + thr->timer(Timer::PAIR); + reduce_thr(this, eflag, vflag, thr); + } // end of omp parallel region +} + +template +void PairLJExpandSphereOMP::eval(int iifrom, int iito, ThrData *const thr) +{ + const auto *_noalias const x = (dbl3_t *) atom->x[0]; + auto *_noalias const f = (dbl3_t *) thr->get_f()[0]; + const double *_noalias const radius = atom->radius; + const int *_noalias const type = atom->type; + const double *_noalias const special_lj = force->special_lj; + const int *_noalias const ilist = list->ilist; + const int *_noalias const numneigh = list->numneigh; + const int *const *const firstneigh = list->firstneigh; + + double xtmp, ytmp, ztmp, rtmp, delx, dely, delz, fxtmp, fytmp, fztmp; + double rcutsq, rsq, r2inv, r6inv, forcelj, factor_lj, evdwl, sigma, sigma6, fpair; + + const int nlocal = atom->nlocal; + int j, jj, jnum, jtype; + + evdwl = 0.0; + + // loop over neighbors of my atoms + + for (int ii = iifrom; ii < iito; ++ii) { + const int i = ilist[ii]; + const int itype = type[i]; + const int *_noalias const jlist = firstneigh[i]; + const double *_noalias const epsiloni = epsilon[itype]; + const double *_noalias const cutsqi = cutsq[itype]; + + xtmp = x[i].x; + ytmp = x[i].y; + ztmp = x[i].z; + rtmp = radius[i]; + jnum = numneigh[i]; + fxtmp = fytmp = fztmp = 0.0; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j].x; + dely = ytmp - x[j].y; + delz = ztmp - x[j].z; + rsq = delx * delx + dely * dely + delz * delz; + jtype = type[j]; + + if (rsq < cutsqi[jtype]) { + + // cutsq is maximum cutoff per type. Now compute and apply real cutoff + + sigma = 2.0 * mix_distance(rtmp, radius[j]); + rcutsq = square(cut[itype][jtype] * sigma); + + if (rsq < rcutsq) { + + r2inv = 1.0 / rsq; + r6inv = r2inv * r2inv * r2inv; + + sigma6 = powint(sigma, 6); + forcelj = r6inv * 24.0 * epsiloni[jtype] * (2.0 * sigma6 * sigma6 * r6inv - sigma6); + fpair = factor_lj * forcelj * r2inv; + + fxtmp += delx * fpair; + fytmp += dely * fpair; + fztmp += delz * fpair; + if (NEWTON_PAIR || j < nlocal) { + f[j].x -= delx * fpair; + f[j].y -= dely * fpair; + f[j].z -= delz * fpair; + } + + if (EFLAG) { + evdwl = r6inv * 4.0 * epsiloni[jtype]; + evdwl *= sigma6 * sigma6 * r6inv - sigma6; + if (offset_flag && (cutsqi[jtype] > 0.0)) { + const double ratio6 = sigma6 / powint(rcutsq, 3); + evdwl -= 4.0 * epsiloni[jtype] * (ratio6 * ratio6 - ratio6); + } + evdwl *= factor_lj; + } + + if (EVFLAG) + ev_tally_thr(this, i, j, nlocal, NEWTON_PAIR, evdwl, 0.0, fpair, delx, dely, delz, thr); + } + } + } + f[i].x += fxtmp; + f[i].y += fytmp; + f[i].z += fztmp; + } +} + +/* ---------------------------------------------------------------------- */ + +double PairLJExpandSphereOMP::memory_usage() +{ + double bytes = memory_usage_thr(); + bytes += PairLJExpandSphere::memory_usage(); + + return bytes; +} diff --git a/src/OPENMP/pair_lj_expand_sphere_omp.h b/src/OPENMP/pair_lj_expand_sphere_omp.h new file mode 100644 index 0000000000..d160e49f1f --- /dev/null +++ b/src/OPENMP/pair_lj_expand_sphere_omp.h @@ -0,0 +1,46 @@ +/* -*- 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(lj/expand/sphere/omp,PairLJExpandSphereOMP); +// clang-format on +#else + +#ifndef LMP_PAIR_LJ_EXPAND_SPHERE_OMP_H +#define LMP_PAIR_LJ_EXPAND_SPHERE_OMP_H + +#include "pair_lj_expand_sphere.h" +#include "thr_omp.h" + +namespace LAMMPS_NS { + +class PairLJExpandSphereOMP : public PairLJExpandSphere, public ThrOMP { + + public: + PairLJExpandSphereOMP(class LAMMPS *); + + void compute(int, int) override; + double memory_usage() override; + + private: + template + void eval(int ifrom, int ito, ThrData *const thr); +}; +} // namespace LAMMPS_NS +#endif +#endif diff --git a/unittest/force-styles/tests/atomic-pair-lj_sphere.yaml b/unittest/force-styles/tests/atomic-pair-lj_cut_sphere.yaml similarity index 99% rename from unittest/force-styles/tests/atomic-pair-lj_sphere.yaml rename to unittest/force-styles/tests/atomic-pair-lj_cut_sphere.yaml index 4205868be4..3a5122a896 100644 --- a/unittest/force-styles/tests/atomic-pair-lj_sphere.yaml +++ b/unittest/force-styles/tests/atomic-pair-lj_cut_sphere.yaml @@ -4,7 +4,7 @@ date_generated: Thu Mar 30 14:38:22 2023 epsilon: 7.5e-13 skip_tests: single prerequisites: ! | - pair lj/sphere + pair lj/cut/sphere atom sphere pre_commands: ! | echo screen @@ -21,7 +21,7 @@ pre_commands: ! | velocity all create 3.0 4534624 loop geom post_commands: ! "" input_file: in.empty -pair_style: lj/sphere 2.5 +pair_style: lj/cut/sphere 2.5 pair_coeff: ! | 1 1 1.0 extract: ! | diff --git a/unittest/force-styles/tests/atomic-pair-lj_expand_sphere.yaml b/unittest/force-styles/tests/atomic-pair-lj_expand_sphere.yaml new file mode 100644 index 0000000000..027ff82306 --- /dev/null +++ b/unittest/force-styles/tests/atomic-pair-lj_expand_sphere.yaml @@ -0,0 +1,104 @@ +--- +lammps_version: 28 Mar 2023 +date_generated: Thu Mar 30 14:38:22 2023 +epsilon: 7.5e-13 +skip_tests: single +prerequisites: ! | + pair lj/expand/sphere + atom sphere +pre_commands: ! | + echo screen + atom_modify map array + units lj + atom_style sphere + lattice fcc 0.8442 + region box block 0 2 0 2 0 2 + create_box 1 box + create_atoms 1 box + displace_atoms all random 0.1 0.1 0.1 623426 + set group all diameter 1.0 + set group all mass 1.0 + velocity all create 3.0 4534624 loop geom +post_commands: ! "" +input_file: in.empty +pair_style: lj/expand/sphere 2.5 +pair_coeff: ! | + 1 1 1.0 +extract: ! | + epsilon 2 +natoms: 32 +init_vdwl: -116.44208632277189 +init_coul: 0 +init_stress: ! |2- + 2.8031747808711788e+02 3.3753332096206151e+02 3.6161491443553058e+02 9.1207933030798330e+01 1.5243342905404037e+02 -9.9498392667113833e+01 +init_forces: ! |2 + 1 -1.4586568952535859e+02 -2.1931945697523867e+02 8.3668114906961620e+01 + 2 1.5110060734448868e+02 1.6194983298501683e+02 5.1680305209387697e+01 + 3 -4.1775824251599438e+01 1.2867107093830445e+00 -3.6313865233596999e+01 + 4 2.9472326459801721e+01 8.8486535373136732e+01 -4.2617580722552610e+01 + 5 5.2397837003192160e+00 2.0459375918115779e+01 -2.3426042702551680e+01 + 6 9.6723821245494683e+00 -4.1806082480411575e+01 1.2511584321172965e+01 + 7 -5.1675096093601567e+01 -4.0551926756618926e+01 -1.3294327178693157e+01 + 8 -8.2578663162446531e+01 -2.0603180343277316e+02 1.1080781077571790e+02 + 9 6.6563366218520699e+00 -1.3576270514928492e+01 -6.2990879039987373e+00 + 10 1.2562028508263918e+01 7.1383899379197500e+00 -7.4200177576219133e+00 + 11 -4.0616225569260138e+01 2.9426535484086656e+01 -1.5194870287132140e+01 + 12 -1.2185920630572705e+01 2.0903918930070052e+01 2.3935223952337630e+01 + 13 -1.0168289428829605e+02 1.2181428476742342e+02 -2.4970463716793029e+02 + 14 1.1276607360598639e+01 2.0932524903784913e+00 -8.4445936145127263e+00 + 15 2.1147814932644917e+02 7.2455411745058896e+01 1.3994698736930633e+02 + 16 2.1307381939437682e+00 -9.7110437674312617e+00 1.6772426774878742e+01 + 17 1.3917131309134271e+01 -4.3979124892489672e+01 5.9765040597753625e+01 + 18 2.0607405082480078e+01 -3.0878340639707460e+01 1.4167531514432428e+01 + 19 4.6190203936710283e+01 -1.1202702997067846e+01 -3.9476561053715500e+01 + 20 -4.0392151835430383e+01 7.6413322558301815e+01 -9.9315439072912369e+01 + 21 3.9812181220755775e+01 -8.0229922011073427e+00 4.8494181799577049e+01 + 22 3.5117168233481337e+01 -5.0880287326141598e+01 8.1590667313257761e+00 + 23 -3.3183438691899802e+01 5.8452017613468186e+01 1.3689689859305707e+01 + 24 -7.0085612135333859e+00 -1.4882951194443576e+01 -1.6978104605753799e+01 + 25 -2.6945478276353324e+01 1.4989629752987154e+01 2.2837312716809768e+00 + 26 1.5934827401819398e+00 1.0517784981514938e+01 1.1915938970456459e+01 + 27 2.8342468546108526e+00 3.3183558717337318e+00 -9.9597072230404908e+00 + 28 1.0973021677369854e+02 -1.8959241094027011e+01 1.1838770808181975e+02 + 29 -3.7698386363756697e+01 4.4468578018428147e+01 1.1140943290101237e+01 + 30 -8.2712093475021490e+01 -4.7423589066212351e+01 -1.3330169888378890e+02 + 31 -2.8868167685733068e+00 2.1403398282492340e+01 -3.1026283179188514e+01 + 32 -2.1837556456164062e+00 1.6484779190829122e+00 5.4465311607738585e+00 +run_vdwl: -117.1001111727232 +run_coul: 0 +run_stress: ! |2- + 2.7695385711446124e+02 3.3278546187873974e+02 3.6006476859072728e+02 8.6715925372879269e+01 1.5264294056452928e+02 -9.6728681863809229e+01 +run_forces: ! |2 + 1 -1.4371784794329676e+02 -2.1498566778457018e+02 8.3095590725041319e+01 + 2 1.4875978675032550e+02 1.5837418598684059e+02 5.1936572842869772e+01 + 3 -4.1902045381131956e+01 1.3746278031477663e+00 -3.7012679382131516e+01 + 4 2.7963282998007497e+01 8.7029018636752809e+01 -4.2734730369263040e+01 + 5 5.0038082469556810e+00 2.0394292459977670e+01 -2.3158066316063287e+01 + 6 9.8078855865137822e+00 -4.1221556143424770e+01 1.2111606884041064e+01 + 7 -4.9916678949523160e+01 -3.9032632064977967e+01 -1.3057662213781221e+01 + 8 -7.9404298765990461e+01 -2.0228919740312625e+02 1.0842014530532373e+02 + 9 6.5305851467989928e+00 -1.3836534088482027e+01 -6.3753564084228440e+00 + 10 1.2948519162479368e+01 7.1663216684672859e+00 -7.9026502950718713e+00 + 11 -4.0485527318703632e+01 2.9658129253905681e+01 -1.5413396670143056e+01 + 12 -1.2384564654464873e+01 2.0906797996653175e+01 2.4583702026262241e+01 + 13 -1.0117689275037893e+02 1.1973355069764031e+02 -2.4623920962306576e+02 + 14 1.1129073511061154e+01 1.9554924379279703e+00 -8.3457798988948593e+00 + 15 2.0816516562230623e+02 7.0803370473935900e+01 1.3842506215462157e+02 + 16 2.1310828504091881e+00 -9.8135507177243841e+00 1.6536268967032719e+01 + 17 1.3641628242862716e+01 -4.4119911040585336e+01 5.9614054383878042e+01 + 18 2.1006853310301640e+01 -3.1211293326535454e+01 1.4241007813718358e+01 + 19 4.6009040579610449e+01 -1.1345944054890328e+01 -3.9459071459686321e+01 + 20 -4.0397130717929159e+01 7.5279993851872220e+01 -9.8511341807527316e+01 + 21 4.0039315283046079e+01 -8.0137748724521511e+00 4.8660642711420813e+01 + 22 3.5708168073593527e+01 -5.1821322839622454e+01 8.2819277512096683e+00 + 23 -3.3113759560430822e+01 5.9072079885315674e+01 1.4264033179943622e+01 + 24 -6.9593081581622727e+00 -1.4886965489878571e+01 -1.6947529042442767e+01 + 25 -2.7405937853181872e+01 1.5364157291023506e+01 2.2982650066647481e+00 + 26 1.5890405119239828e+00 1.0786544007531514e+01 1.2070193510601317e+01 + 27 2.8643954209448790e+00 3.4431707338380022e+00 -1.0003156259622884e+01 + 28 1.1184386324630830e+02 -1.8203684205576618e+01 1.2042827913401095e+02 + 29 -3.8541921364351666e+01 4.5525124123081802e+01 1.1427605135155794e+01 + 30 -8.4505523479161951e+01 -4.8548646933958288e+01 -1.3594721117648746e+02 + 31 -2.9919458949399576e+00 2.0910279091646661e+01 -3.0672195920243109e+01 + 32 -2.2381117518014042e+00 1.5535445662463561e+00 5.3850793110514772e+00 +...