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

rename lj/sphere to lj/cut/sphere and add (yet non-functional) lj/expand/sphere

Browse Source
This commit is contained in:
Axel Kohlmeyer
2023-04-06 08:17:17 -04:00
parent 8e3ec4d567
commit e82fd31bd4
15 changed files with 1025 additions and 58 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
doc/src/Commands_pair.rst
View File

@ -172,12 +172,14 @@ OPT.
* :doc:`lj/cut/dipole/long (g) <pair_dipole>` * :doc:`lj/cut/dipole/long (g) <pair_dipole>`
* :doc:`lj/cut/dipole/sf (go) <pair_dipole>` * :doc:`lj/cut/dipole/sf (go) <pair_dipole>`
* :doc:`lj/cut/soft (o) <pair_fep_soft>` * :doc:`lj/cut/soft (o) <pair_fep_soft>`
* :doc:`lj/cut/sphere (o) <pair_lj_cut_sphere>`
* :doc:`lj/cut/thole/long (o) <pair_thole>` * :doc:`lj/cut/thole/long (o) <pair_thole>`
* :doc:`lj/cut/tip4p/cut (o) <pair_lj_cut_tip4p>` * :doc:`lj/cut/tip4p/cut (o) <pair_lj_cut_tip4p>`
* :doc:`lj/cut/tip4p/long (got) <pair_lj_cut_tip4p>` * :doc:`lj/cut/tip4p/long (got) <pair_lj_cut_tip4p>`
* :doc:`lj/cut/tip4p/long/soft (o) <pair_fep_soft>` * :doc:`lj/cut/tip4p/long/soft (o) <pair_fep_soft>`
* :doc:`lj/expand (gko) <pair_lj_expand>` * :doc:`lj/expand (gko) <pair_lj_expand>`
* :doc:`lj/expand/coul/long (gk) <pair_lj_expand>` * :doc:`lj/expand/coul/long (gk) <pair_lj_expand>`
* :doc:`lj/expand/sphere (o) <pair_lj_expand_sphere>`
* :doc:`lj/gromacs (gko) <pair_gromacs>` * :doc:`lj/gromacs (gko) <pair_gromacs>`
* :doc:`lj/gromacs/coul/gromacs (ko) <pair_gromacs>` * :doc:`lj/gromacs/coul/gromacs (ko) <pair_gromacs>`
* :doc:`lj/long/coul/long (iot) <pair_lj_long>` * :doc:`lj/long/coul/long (iot) <pair_lj_long>`
@ -186,7 +188,6 @@ OPT.
* :doc:`lj/long/tip4p/long (o) <pair_lj_long>` * :doc:`lj/long/tip4p/long (o) <pair_lj_long>`
* :doc:`lj/mdf <pair_mdf>` * :doc:`lj/mdf <pair_mdf>`
* :doc:`lj/relres (o) <pair_lj_relres>` * :doc:`lj/relres (o) <pair_lj_relres>`
* :doc:`lj/sphere (o) <pair_lj_sphere>`
* :doc:`lj/spica (gko) <pair_spica>` * :doc:`lj/spica (gko) <pair_spica>`
* :doc:`lj/spica/coul/long (go) <pair_spica>` * :doc:`lj/spica/coul/long (go) <pair_spica>`
* :doc:`lj/spica/coul/msm (o) <pair_spica>` * :doc:`lj/spica/coul/msm (o) <pair_spica>`

33
doc/src/pair_lj_sphere.rst → doc/src/pair_lj_cut_sphere.rst
View File

@ -1,10 +1,10 @@
.. index:: pair_style lj/sphere .. index:: pair_style lj/cut/sphere
.. index:: pair_style lj/sphere/omp .. 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 Syntax
"""""" """"""
@ -13,12 +13,12 @@ Syntax
pair_style style args pair_style style args
* style = *lj/sphere* * style = *lj/cut/sphere*
* args = list of arguments for a particular style * args = list of arguments for a particular style
.. parsed-literal:: .. parsed-literal::
*lj/sphere* args = cutoff ratio *lj/cut/sphere* args = cutoff ratio
cutoff = global cutoff ratio for Lennard Jones interactions (unitless) cutoff = global cutoff ratio for Lennard Jones interactions (unitless)
Examples Examples
@ -26,14 +26,14 @@ Examples
.. code-block:: LAMMPS .. code-block:: LAMMPS
pair_style lj/sphere 2.5 pair_style lj/cut/sphere 2.5
pair_coeff * * 1.0 pair_coeff * * 1.0
pair_coeff 1 1 1.1 2.8 pair_coeff 1 1 1.1 2.8
Description 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 given by
.. math:: .. math::
@ -106,7 +106,7 @@ is at :math:`2^{\frac{1}{6}} \sigma_{ij}`.
group large variable large group large variable large
set group largea type 2 set group largea type 2
pair_style lj/sphere 2.5 pair_style lj/cut/sphere 2.5
pair_coeff * * 1.0 pair_coeff * * 1.0
neighbor 0.3 multi 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 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 <pair_modify>` style is *geometric*. See the :doc:`pair_modify command <pair_modify>`
for details. for details.
The *lj/sphere* pair style supports the :doc:`pair_modify shift <pair_modify>` The *lj/cut/sphere* pair style supports the :doc:`pair_modify shift <pair_modify>`
option for the energy of the Lennard-Jones portion of the pair interaction. 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
<pair_modify>` tail option for adding a long-range tail corrections to <pair_modify>` tail option for adding a long-range tail corrections to
the energy and pressure. 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 <restart>`, so pair_style and pair_coeff commands do not restart files <restart>`, so pair_style and pair_coeff commands do not
need to be specified in an input script that reads a restart file. 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 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 <Build_package>` page EXTRA-PAIR package. See the :doc:`Build package <Build_package>` page
for more info. 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 <pair_coeff>` * :doc:`pair_coeff <pair_coeff>`
* :doc:`pair_style lj/cut <pair_lj>` * :doc:`pair_style lj/cut <pair_lj>`
* :doc:`pair_style lj/expnd/sphere <pair_lj_expand_sphere>`
Default Default
""""""" """""""

184
doc/src/pair_lj_expand_sphere.rst Normal file
View File

@ -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_lj>` pair style, but the :math:`\sigma_{ij}` parameter is not
set as a per-type parameter via the :doc:`pair_coeff command
<pair_coeff>`. Instead it is calculated individually for each pair
using the per-atom diameter attribute of :doc:`atom_style sphere
<atom_style>` 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
<pair_modify>` 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 <neighbor>` 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 <pair_coeff>` command as in the examples above, or in the data
file or restart files read by the :doc:`read_data <read_data>` or
:doc:`read_restart <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 <pair_style>` 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 <pair_modify>`
for details.
The *lj/expand/sphere* pair style supports the :doc:`pair_modify shift <pair_modify>`
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
<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 <restart>`, 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 <run_style>` 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 <Build_package>` page
for more info.
The *lj/expand/sphere* pair style does not support the *sixthpower* mixing rule.
----------
Related commands
""""""""""""""""
* :doc:`pair_coeff <pair_coeff>`
* :doc:`pair_style lj/cut <pair_lj>`
* :doc:`pair_style lj/cut/sphere <pair_lj_cut_sphere>`
Default
"""""""
none

3
doc/src/pair_style.rst
View File

@ -250,12 +250,14 @@ accelerated styles exist.
* :doc:`lj/cut/dipole/cut <pair_dipole>` - point dipoles with cutoff * :doc:`lj/cut/dipole/cut <pair_dipole>` - point dipoles with cutoff
* :doc:`lj/cut/dipole/long <pair_dipole>` - point dipoles with long-range Ewald * :doc:`lj/cut/dipole/long <pair_dipole>` - point dipoles with long-range Ewald
* :doc:`lj/cut/soft <pair_fep_soft>` - LJ with a soft core * :doc:`lj/cut/soft <pair_fep_soft>` - LJ with a soft core
* :doc:`lj/cut/sphere <pair_lj_cut_sphere>` - LJ where per-atom radius is used as LJ sigma
* :doc:`lj/cut/thole/long <pair_thole>` - LJ with Coulomb with thole damping * :doc:`lj/cut/thole/long <pair_thole>` - LJ with Coulomb with thole damping
* :doc:`lj/cut/tip4p/cut <pair_lj_cut_tip4p>` - LJ with cutoff Coulomb for TIP4P water * :doc:`lj/cut/tip4p/cut <pair_lj_cut_tip4p>` - LJ with cutoff Coulomb for TIP4P water
* :doc:`lj/cut/tip4p/long <pair_lj_cut_tip4p>` - LJ with long-range Coulomb for TIP4P water * :doc:`lj/cut/tip4p/long <pair_lj_cut_tip4p>` - LJ with long-range Coulomb for TIP4P water
* :doc:`lj/cut/tip4p/long/soft <pair_fep_soft>` - LJ with cutoff Coulomb for TIP4P water with a soft core * :doc:`lj/cut/tip4p/long/soft <pair_fep_soft>` - LJ with cutoff Coulomb for TIP4P water with a soft core
* :doc:`lj/expand <pair_lj_expand>` - Lennard-Jones for variable size particles * :doc:`lj/expand <pair_lj_expand>` - Lennard-Jones for variable size particles
* :doc:`lj/expand/coul/long <pair_lj_expand>` - Lennard-Jones for variable size particles with long-range Coulomb * :doc:`lj/expand/coul/long <pair_lj_expand>` - Lennard-Jones for variable size particles with long-range Coulomb
* :doc:`lj/expand/sphere <pair_lj_sphere>` - Variable size LJ where per-atom radius is used as delta (size)
* :doc:`lj/gromacs <pair_gromacs>` - GROMACS-style Lennard-Jones potential * :doc:`lj/gromacs <pair_gromacs>` - GROMACS-style Lennard-Jones potential
* :doc:`lj/gromacs/coul/gromacs <pair_gromacs>` - GROMACS-style LJ and Coulomb potential * :doc:`lj/gromacs/coul/gromacs <pair_gromacs>` - GROMACS-style LJ and Coulomb potential
* :doc:`lj/long/coul/long <pair_lj_long>` - long-range LJ and long-range Coulomb * :doc:`lj/long/coul/long <pair_lj_long>` - long-range LJ and long-range Coulomb
@ -264,7 +266,6 @@ accelerated styles exist.
* :doc:`lj/long/tip4p/long <pair_lj_long>` - long-range LJ and long-range Coulomb for TIP4P water * :doc:`lj/long/tip4p/long <pair_lj_long>` - long-range LJ and long-range Coulomb for TIP4P water
* :doc:`lj/mdf <pair_mdf>` - LJ potential with a taper function * :doc:`lj/mdf <pair_mdf>` - LJ potential with a taper function
* :doc:`lj/relres <pair_lj_relres>` - LJ using multiscale Relative Resolution (RelRes) methodology :ref:`(Chaimovich) <Chaimovich2>`. * :doc:`lj/relres <pair_lj_relres>` - LJ using multiscale Relative Resolution (RelRes) methodology :ref:`(Chaimovich) <Chaimovich2>`.
* :doc:`lj/sphere <pair_lj_sphere>` - LJ where per-atom radius is used as LJ sigma
* :doc:`lj/spica <pair_spica>` - LJ for SPICA coarse-graining * :doc:`lj/spica <pair_spica>` - LJ for SPICA coarse-graining
* :doc:`lj/spica/coul/long <pair_spica>` - LJ for SPICA coarse-graining with long-range Coulomb * :doc:`lj/spica/coul/long <pair_spica>` - LJ for SPICA coarse-graining with long-range Coulomb
* :doc:`lj/spica/coul/msm <pair_spica>` - LJ for SPICA coarse-graining with long-range Coulomb via MSM * :doc:`lj/spica/coul/msm <pair_spica>` - LJ for SPICA coarse-graining with long-range Coulomb via MSM

4
src/.gitignore vendored
View File

@ -1249,6 +1249,8 @@
/pair_lj_cut_dipole_long.h /pair_lj_cut_dipole_long.h
/pair_lj_cut_*hars_*.cpp /pair_lj_cut_*hars_*.cpp
/pair_lj_cut_*hars_*.h /pair_lj_cut_*hars_*.h
/pair_lj_cut_sphere.cpp
/pair_lj_cut_sphere.h
/pair_lj_cut_soft.cpp /pair_lj_cut_soft.cpp
/pair_lj_cut_soft.h /pair_lj_cut_soft.h
/pair_lj_cut_tip4p_long.cpp /pair_lj_cut_tip4p_long.cpp
@ -1263,6 +1265,8 @@
/pair_lj_long_tip4p_long.h /pair_lj_long_tip4p_long.h
/pair_lj_cut_tgpu.cpp /pair_lj_cut_tgpu.cpp
/pair_lj_cut_tgpu.h /pair_lj_cut_tgpu.h
/pair_lj_expand_sphere.cpp
/pair_lj_expand_sphere.h
/pair_lj_spica.cpp /pair_lj_spica.cpp
/pair_lj_spica.h /pair_lj_spica.h
/pair_lj_spica_coul_long.cpp /pair_lj_spica_coul_long.cpp

42
src/EXTRA-PAIR/pair_lj_sphere.cpp → src/EXTRA-PAIR/pair_lj_cut_sphere.cpp
View File

@ -11,7 +11,7 @@
See the README file in the top-level LAMMPS directory. 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 "atom.h"
#include "comm.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; writedata = 1;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
PairLJSphere::~PairLJSphere() PairLJCutSphere::~PairLJCutSphere()
{ {
if (allocated) { if (allocated) {
memory->destroy(setflag); 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; int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, rtmp, delx, dely, delz, evdwl, sigma, sigma6, fpair; 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 allocate all arrays
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairLJSphere::allocate() void PairLJCutSphere::allocate()
{ {
allocated = 1; allocated = 1;
int n = atom->ntypes + 1; int n = atom->ntypes + 1;
@ -165,7 +166,7 @@ void PairLJSphere::allocate()
global settings 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"); 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 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 (narg < 3 || narg > 4) error->all(FLERR, "Incorrect args for pair coefficients");
if (!allocated) allocate(); if (!allocated) allocate();
@ -215,13 +216,14 @@ void PairLJSphere::coeff(int narg, char **arg)
init specific to this pair style init specific to this pair style
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairLJSphere::init_style() void PairLJCutSphere::init_style()
{ {
Pair::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) 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 // determine max radius per type
@ -241,7 +243,7 @@ void PairLJSphere::init_style()
init for one type pair i,j and corresponding j,i 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) { if (setflag[i][j] == 0) {
epsilon[i][j] = mix_energy(epsilon[i][i], epsilon[j][j], 0.0, 0.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 proc 0 writes to restart file
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairLJSphere::write_restart(FILE *fp) void PairLJCutSphere::write_restart(FILE *fp)
{ {
write_restart_settings(fp); write_restart_settings(fp);
@ -279,7 +281,7 @@ void PairLJSphere::write_restart(FILE *fp)
proc 0 reads from restart file, bcasts proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairLJSphere::read_restart(FILE *fp) void PairLJCutSphere::read_restart(FILE *fp)
{ {
read_restart_settings(fp); read_restart_settings(fp);
allocate(); allocate();
@ -305,7 +307,7 @@ void PairLJSphere::read_restart(FILE *fp)
proc 0 writes to restart file 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(&cut_global, sizeof(double), 1, fp);
fwrite(&offset_flag, sizeof(int), 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 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; int me = comm->me;
if (me == 0) { if (me == 0) {
@ -333,7 +335,7 @@ void PairLJSphere::read_restart_settings(FILE *fp)
proc 0 writes to data file 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]); 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 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 i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++) 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 PairLJCutSphere::single(int i, int j, int itype, int jtype, double rsq,
double factor_lj, double &fforce) double /*factor_coul*/, double factor_lj, double &fforce)
{ {
double r2inv, r6inv, rcutsq, sigma, sigma6, forcelj, philj; 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; dim = 2;
if (strcmp(str, "epsilon") == 0) return (void *) epsilon; if (strcmp(str, "epsilon") == 0) return (void *) epsilon;

12
src/EXTRA-PAIR/pair_lj_sphere.h → src/EXTRA-PAIR/pair_lj_cut_sphere.h
View File

@ -13,21 +13,21 @@
#ifdef PAIR_CLASS #ifdef PAIR_CLASS
// clang-format off // clang-format off
PairStyle(lj/sphere,PairLJSphere); PairStyle(lj/cut/sphere,PairLJCutSphere);
// clang-format on // clang-format on
#else #else
#ifndef LMP_PAIR_LJ_SPHERE_H #ifndef LMP_PAIR_LJ_CUT_SPHERE_H
#define LMP_PAIR_LJ_SPHERE_H #define LMP_PAIR_LJ_CUT_SPHERE_H
#include "pair.h" #include "pair.h"
namespace LAMMPS_NS { namespace LAMMPS_NS {
class PairLJSphere : public Pair { class PairLJCutSphere : public Pair {
public: public:
PairLJSphere(class LAMMPS *); PairLJCutSphere(class LAMMPS *);
~PairLJSphere() override; ~PairLJCutSphere() override;
void compute(int, int) override; void compute(int, int) override;
void settings(int, char **) override; void settings(int, char **) override;
void coeff(int, char **) override; void coeff(int, char **) override;

385
src/EXTRA-PAIR/pair_lj_expand_sphere.cpp Normal file
View File

@ -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 <cmath>
#include <cstring>
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;
}

58
src/EXTRA-PAIR/pair_lj_expand_sphere.h Normal file
View File

@ -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

12
src/OPENMP/pair_lj_sphere_omp.cpp → src/OPENMP/pair_lj_cut_sphere_omp.cpp
View File

@ -12,7 +12,7 @@
Contributing author: Axel Kohlmeyer (Temple U) Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#include "pair_lj_sphere_omp.h" #include "pair_lj_cut_sphere_omp.h"
#include "atom.h" #include "atom.h"
#include "comm.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; suffix_flag |= Suffix::OMP;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void PairLJSphereOMP::compute(int eflag, int vflag) void PairLJCutSphereOMP::compute(int eflag, int vflag)
{ {
ev_init(eflag, vflag); ev_init(eflag, vflag);
@ -78,7 +78,7 @@ void PairLJSphereOMP::compute(int eflag, int vflag)
} }
template <int EVFLAG, int EFLAG, int NEWTON_PAIR> template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
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]; const auto *_noalias const x = (dbl3_t *) atom->x[0];
auto *_noalias const f = (dbl3_t *) thr->get_f()[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(); double bytes = memory_usage_thr();
bytes += PairLJSphere::memory_usage(); bytes += PairLJCutSphere::memory_usage();
return bytes; return bytes;
} }

12
src/OPENMP/pair_lj_sphere_omp.h → src/OPENMP/pair_lj_cut_sphere_omp.h
View File

@ -17,22 +17,22 @@
#ifdef PAIR_CLASS #ifdef PAIR_CLASS
// clang-format off // clang-format off
PairStyle(lj/sphere/omp,PairLJSphereOMP); PairStyle(lj/cut/sphere/omp,PairLJCutSphereOMP);
// clang-format on // clang-format on
#else #else
#ifndef LMP_PAIR_LJ_SPHERE_OMP_H #ifndef LMP_PAIR_LJ_CUT_SPHERE_OMP_H
#define LMP_PAIR_LJ_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" #include "thr_omp.h"
namespace LAMMPS_NS { namespace LAMMPS_NS {
class PairLJSphereOMP : public PairLJSphere, public ThrOMP { class PairLJCutSphereOMP : public PairLJCutSphere, public ThrOMP {
public: public:
PairLJSphereOMP(class LAMMPS *); PairLJCutSphereOMP(class LAMMPS *);
void compute(int, int) override; void compute(int, int) override;
double memory_usage() override; double memory_usage() override;

181
src/OPENMP/pair_lj_expand_sphere_omp.cpp Normal file
View File

@ -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 <int EVFLAG, int EFLAG, int NEWTON_PAIR>
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;
}

46
src/OPENMP/pair_lj_expand_sphere_omp.h Normal file
View File

@ -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 <int EVFLAG, int EFLAG, int NEWTON_PAIR>
void eval(int ifrom, int ito, ThrData *const thr);
};
} // namespace LAMMPS_NS
#endif
#endif

4
unittest/force-styles/tests/atomic-pair-lj_sphere.yaml → unittest/force-styles/tests/atomic-pair-lj_cut_sphere.yaml
View File

@ -4,7 +4,7 @@ date_generated: Thu Mar 30 14:38:22 2023
epsilon: 7.5e-13 epsilon: 7.5e-13
skip_tests: single skip_tests: single
prerequisites: ! | prerequisites: ! |
pair lj/sphere pair lj/cut/sphere
atom sphere atom sphere
pre_commands: ! | pre_commands: ! |
echo screen echo screen
@ -21,7 +21,7 @@ pre_commands: ! |
velocity all create 3.0 4534624 loop geom velocity all create 3.0 4534624 loop geom
post_commands: ! "" post_commands: ! ""
input_file: in.empty input_file: in.empty
pair_style: lj/sphere 2.5 pair_style: lj/cut/sphere 2.5
pair_coeff: ! | pair_coeff: ! |
1 1 1.0 1 1 1.0
extract: ! | extract: ! |

104
unittest/force-styles/tests/atomic-pair-lj_expand_sphere.yaml Normal file
View File

@ -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
...