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

Merge branch 'lammps:develop' into mliap_null_fix

Browse Source
This commit is contained in:
bnebgen-LANL
2023-10-21 04:08:19 -06:00
committed by GitHub
parent 1064c82c00 58d299d08e
commit 710cd1c001
62 changed files with 1295 additions and 568 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
.github/CODEOWNERS vendored
View File

@ -135,6 +135,7 @@ src/timer.* @akohlmey
src/utils.* @akohlmey @rbberger
src/verlet.* @sjplimp @stanmoore1
src/math_eigen_impl.h @jewettaij
src/fix_press_langevin.* @Bibobu
# tools
tools/coding_standard/* @akohlmey @rbberger

3
doc/src/Howto_triclinic.rst
View File

@ -12,7 +12,8 @@ is created, e.g. by the :doc:`create_box <create_box>` or
:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
commands. Additionally, LAMMPS defines box size parameters lx,ly,lz
where lx = xhi-xlo, and similarly in the y and z dimensions. The 6
parameters, as well as lx,ly,lz, can be output via the :doc:`thermo_style custom <thermo_style>` command.
parameters, as well as lx,ly,lz, can be output via the
:doc:`thermo_style custom <thermo_style>` command.
LAMMPS also allows simulations to be performed in triclinic
(non-orthogonal) simulation boxes shaped as a parallelepiped with

5
doc/src/atom_modify.rst
View File

@ -65,6 +65,11 @@ switch. This is described on the :doc:`Build_settings <Build_settings>`
doc page. If atom IDs are not used, they must be specified as 0 for
all atoms, e.g. in a data or restart file.
.. note::
If a :doc:`triclinic simulation box <Howto_triclinic>` is used,
atom IDs are required, due to how neighbor lists are built.
The *map* keyword determines how atoms with specific IDs are found
when required. An example are the bond (angle, etc) methods which
need to find the local index of an atom with a specific global ID

2
doc/src/compute.rst
View File

@ -74,7 +74,7 @@ Global, per-atom, local, and per-grid quantities can also be of three
for each atom, each local entity, or each grid cell.
Note that a single compute can produce any combination of global,
per-atom, local, or per-grid values. Likewise it can prouduce any
per-atom, local, or per-grid values. Likewise it can produce any
combination of scalar, vector, or array output for each style. The
exception is that for per-atom, local, and per-grid output, either a
vector or array can be produced, but not both. The doc page for each

2
doc/src/compute_voronoi_atom.rst
View File

@ -232,4 +232,4 @@ Related commands
Default
"""""""
The default for the neighobrs keyword is no.
The default for the neighbors keyword is no.

2
doc/src/fix.rst
View File

@ -106,7 +106,7 @@ Global, per-atom, local, and per-grid quantities can also be of three
for each atom, each local entity, or each grid cell.
Note that a single fix can produce any combination of global,
per-atom, local, or per-grid values. Likewise it can prouduce any
per-atom, local, or per-grid values. Likewise it can produce any
combination of scalar, vector, or array output for each style. The
exception is that for per-atom, local, and per-grid output, either a
vector or array can be produced, but not both. The doc page for each

4
doc/src/fix_ave_histo.rst
View File

@ -106,7 +106,7 @@ attributes are per-atom vector values. See the page for individual
generate.
Note that a compute or fix can produce multiple kinds of data (global,
per-atom, local). If LAMMPS cannot unambiguosly determine which kind
per-atom, local). If LAMMPS cannot unambiguously determine which kind
of data to use, the optional *kind* keyword discussed below can force
the desired disambiguation.
@ -263,7 +263,7 @@ keyword is set to *vector*, then all input values must be global or
per-atom or local vectors, or columns of global or per-atom or local
arrays.
The *kind* keyword only needs to be used if any of the specfied input
The *kind* keyword only needs to be used if any of the specified input
computes or fixes produce more than one kind of output (global,
per-atom, local). If not, LAMMPS will determine the kind of data all
the inputs produce and verify it is all the same kind. If not, an

51
doc/src/fix_srd.rst
View File

@ -61,25 +61,30 @@ Description
Treat a group of particles as stochastic rotation dynamics (SRD)
particles that serve as a background solvent when interacting with big
(colloidal) particles in groupbig-ID. The SRD formalism is described
in :ref:`(Hecht) <Hecht>`. The key idea behind using SRD particles as a
cheap coarse-grained solvent is that SRD particles do not interact
with each other, but only with the solute particles, which in LAMMPS
can be spheroids, ellipsoids, or line segments, or triangles, or rigid
bodies containing multiple spheroids or ellipsoids or line segments
or triangles. The collision and rotation properties of the model
imbue the SRD particles with fluid-like properties, including an
effective viscosity. Thus simulations with large solute particles can
be run more quickly, to measure solute properties like diffusivity
and viscosity in a background fluid. The usual LAMMPS fixes for such
simulations, such as :doc:`fix deform <fix_deform>`,
:doc:`fix viscosity <fix_viscosity>`, and :doc:`fix nvt/sllod <fix_nvt_sllod>`,
can be used in conjunction with the SRD model.
in :ref:`(Hecht) <Hecht>`. The same methodology is also called
multi-particle collision dynamics (MPCD) in the literature.
For more details on how the SRD model is implemented in LAMMPS,
:ref:`(Petersen) <Petersen1>` describes the implementation and usage of
pure SRD fluids. See the ``examples/srd`` directory for sample input
scripts using SRD particles for that and for mixture systems (solute
particles in an SRD fluid).
The key idea behind using SRD particles as a cheap coarse-grained
solvent is that SRD particles do not interact with each other, but
only with the solute particles, which in LAMMPS can be spheroids,
ellipsoids, or line segments, or triangles, or rigid bodies containing
multiple spheroids or ellipsoids or line segments or triangles. The
collision and rotation properties of the model imbue the SRD particles
with fluid-like properties, including an effective viscosity. Thus
simulations with large solute particles can be run more quickly, to
measure solute properties like diffusivity and viscosity in a
background fluid. The usual LAMMPS fixes for such simulations, such
as :doc:`fix deform <fix_deform>`, :doc:`fix viscosity
<fix_viscosity>`, and :doc:`fix nvt/sllod <fix_nvt_sllod>`, can be
used in conjunction with the SRD model.
These 3 papers give more details on how the SRD model is implemented
in LAMMPS. :ref:`(Petersen) <Petersen1>` describes pure SRD fluid
systems. :ref:`(Bolintineanu1) <Bolintineanu1>` describes models
where pure SRD fluids :ref:interact with boundary walls.
:ref:`(Bolintineanu2) <Bolintineanu2>` describes mixture models where
large colloidal particles are solvated by an SRD fluid. See the
``examples/srd`` :ref:directory for sample input scripts.
This fix does two things:
@ -404,3 +409,13 @@ no, and rescale = yes.
**(Petersen)** Petersen, Lechman, Plimpton, Grest, in' t Veld, Schunk, J
Chem Phys, 132, 174106 (2010).
.. _Bolintineanu1:
**(Bolintineanu1)**
Bolintineanu, Lechman, Plimpton, Grest, Phys Rev E, 86, 066703 (2012).
.. _Bolintineanu2:
**(Bolintineanu2)** Bolintineanu, Grest, Lechman, Pierce, Plimpton,
Schunk, Comp Particle Mechanics, 1, 321-356 (2014).

2
doc/src/thermo_style.rst
View File

@ -442,7 +442,7 @@ equal-style and vector-style variables can be referenced; the latter
requires a bracketed term to specify the Ith element of the vector
calculated by the variable. However, an equal-style variable can use
an atom-style variable in its formula indexed by the ID of an
individual atom. This is a way to output a speciic atom's per-atom
individual atom. This is a way to output a specific atom's per-atom
coordinates or other per-atom properties in thermo output. See the
:doc:`variable <variable>` command for details. Note that variables
of style *equal* and *vector* and *atom* define a formula which can

8
doc/src/variable.rst
View File

@ -1167,7 +1167,7 @@ variables), or global vectors of values. The latter can also be a
column of a global array.
Atom-style variables can use scalar values (same as for equal-style
varaibles), or per-atom vectors of values. The latter can also be a
variables), or per-atom vectors of values. The latter can also be a
column of a per-atom array.
The various allowed compute references in the variable formulas for
@ -1183,7 +1183,7 @@ table:
+--------+------------+------------------------------------------+
| vector | c_ID | global vector |
| vector | c_ID[I] | column of global array |
---------+------------+------------------------------------------+
+--------+------------+------------------------------------------+
| atom | c_ID | per-atom vector |
| atom | c_ID[I] | column of per-atom array |
+--------+------------+------------------------------------------+
@ -1232,7 +1232,7 @@ variables), or global vectors of values. The latter can also be a
column of a global array.
Atom-style variables can use scalar values (same as for equal-style
varaibles), or per-atom vectors of values. The latter can also be a
variables), or per-atom vectors of values. The latter can also be a
column of a per-atom array.
The allowed fix references in variable formulas for equal-, vector-,
@ -1247,7 +1247,7 @@ and atom-style variables are listed in the following table:
+--------+------------+------------------------------------------+
| vector | f_ID | global vector |
| vector | f_ID[I] | column of global array |
---------+------------+------------------------------------------+
+--------+------------+------------------------------------------+
| atom | f_ID | per-atom vector |
| atom | f_ID[I] | column of per-atom array |
+--------+------------+------------------------------------------+

9
src/BOCS/fix_bocs.cpp
View File

@ -1024,7 +1024,10 @@ void FixBocs::final_integrate()
if (pstat_flag) {
if (pstyle == ISO) pressure->compute_scalar();
else pressure->compute_vector();
else {
temperature->compute_vector();
pressure->compute_vector();
}
couple();
pressure->addstep(update->ntimestep+1);
}
@ -1961,6 +1964,7 @@ void FixBocs::nhc_press_integrate()
int ich,i,pdof;
double expfac,factor_etap,kecurrent;
double kt = boltz * t_target;
double lkt_press;
// Update masses, to preserve initial freq, if flag set
@ -2006,7 +2010,8 @@ void FixBocs::nhc_press_integrate()
}
}
double lkt_press = pdof * kt;
if (pstyle == ISO) lkt_press = kt;
else lkt_press = pdof * kt;
etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];
double ncfac = 1.0/nc_pchain;

16
src/EXTRA-DUMP/dump_yaml.cpp
View File

@ -24,6 +24,8 @@
using namespace LAMMPS_NS;
static constexpr char special_chars[] = "{}[],&:*#?|-<>=!%@\\";
/* ---------------------------------------------------------------------- */
DumpYAML::DumpYAML(class LAMMPS *_lmp, int narg, char **args) :
DumpCustom(_lmp, narg, args), thermo(false)
@ -67,7 +69,12 @@ void DumpYAML::write_header(bigint ndump)
const auto &fields = th->get_fields();
thermo_data += "thermo:\n - keywords: [ ";
for (int i = 0; i < nfield; ++i) thermo_data += fmt::format("{}, ", keywords[i]);
for (int i = 0; i < nfield; ++i) {
if (keywords[i].find_first_of(special_chars) == std::string::npos)
thermo_data += fmt::format("{}, ", keywords[i]);
else
thermo_data += fmt::format("'{}', ", keywords[i]);
}
thermo_data += "]\n - data: [ ";
for (int i = 0; i < nfield; ++i) {
@ -107,7 +114,12 @@ void DumpYAML::write_header(bigint ndump)
if (domain->triclinic) fmt::print(fp, " - [ {}, {}, {} ]\n", boxxy, boxxz, boxyz);
fmt::print(fp, "keywords: [ ");
for (const auto &item : utils::split_words(columns)) fmt::print(fp, "{}, ", item);
for (const auto &item : utils::split_words(columns)) {
if (item.find_first_of(special_chars) == std::string::npos)
fmt::print(fp, "{}, ", item);
else
fmt::print(fp, "'{}', ", item);
}
fputs(" ]\ndata:\n", fp);
} else // reset so that the remainder of the output is not multi-proc
filewriter = 0;

2
src/INTEL/fix_intel.cpp
View File

@ -20,6 +20,7 @@
#include "fix_intel.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "neighbor.h"
@ -470,6 +471,7 @@ void FixIntel::pair_init_check(const bool cdmessage)
int need_tag = 0;
if (atom->molecular != Atom::ATOMIC || three_body_neighbor()) need_tag = 1;
if (domain->triclinic && force->newton_pair) need_tag = 1;
// Clear buffers used for pair style
char kmode[80];

66
src/INTEL/npair_halffull_newton_intel.cpp
View File

@ -20,7 +20,9 @@
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "modify.h"
#include "my_page.h"
#include "neigh_list.h"
@ -56,6 +58,9 @@ void NPairHalffullNewtonIntel::build_t(NeighList *list,
const int * _noalias const numneigh_full = list->listfull->numneigh;
const int ** _noalias const firstneigh_full = (const int ** const)list->listfull->firstneigh; // NOLINT
const double delta = 0.01 * force->angstrom;
const int triclinic = domain->triclinic;
#if defined(_OPENMP)
#pragma omp parallel
#endif
@ -82,25 +87,50 @@ void NPairHalffullNewtonIntel::build_t(NeighList *list,
const int * _noalias const jlist = firstneigh_full[i];
const int jnum = numneigh_full[i];
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma ivdep
#endif
for (int jj = 0; jj < jnum; jj++) {
const int joriginal = jlist[jj];
const int j = joriginal & NEIGHMASK;
int addme = 1;
if (j < nlocal) {
if (i > j) addme = 0;
} else {
if (x[j].z < ztmp) addme = 0;
if (x[j].z == ztmp) {
if (x[j].y < ytmp) addme = 0;
if (x[j].y == ytmp && x[j].x < xtmp) addme = 0;
if (!triclinic) {
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma ivdep
#endif
for (int jj = 0; jj < jnum; jj++) {
const int joriginal = jlist[jj];
const int j = joriginal & NEIGHMASK;
int addme = 1;
if (j < nlocal) {
if (i > j) addme = 0;
} else {
if (x[j].z < ztmp) addme = 0;
if (x[j].z == ztmp) {
if (x[j].y < ytmp) addme = 0;
if (x[j].y == ytmp && x[j].x < xtmp) addme = 0;
}
}
if (addme)
neighptr[n++] = joriginal;
}
} else {
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma ivdep
#endif
for (int jj = 0; jj < jnum; jj++) {
const int joriginal = jlist[jj];
const int j = joriginal & NEIGHMASK;
int addme = 1;
if (j < nlocal) {
if (i > j) addme = 0;
} else {
if (fabs(x[j].z-ztmp) > delta) {
if (x[j].z < ztmp) addme = 0;
} else if (fabs(x[j].y-ytmp) > delta) {
if (x[j].y < ytmp) addme = 0;
} else {
if (x[j].x < xtmp) addme = 0;
}
}
if (addme)
neighptr[n++] = joriginal;
}
if (addme)
neighptr[n++] = joriginal;
}
ilist[ii] = i;
@ -203,7 +233,7 @@ void NPairHalffullNewtonIntel::build_t3(NeighList *list, int *numhalf)
void NPairHalffullNewtonIntel::build(NeighList *list)
{
if (_fix->three_body_neighbor() == 0) {
if (_fix->three_body_neighbor() == 0 || domain->triclinic) {
if (_fix->precision() == FixIntel::PREC_MODE_MIXED)
build_t(list, _fix->get_mixed_buffers());
else if (_fix->precision() == FixIntel::PREC_MODE_DOUBLE)

87
src/INTEL/npair_halffull_newton_trim_intel.cpp
View File

@ -20,7 +20,9 @@
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "modify.h"
#include "my_page.h"
#include "neigh_list.h"
@ -57,6 +59,8 @@ void NPairHalffullNewtonTrimIntel::build_t(NeighList *list,
const int ** _noalias const firstneigh_full = (const int ** const)list->listfull->firstneigh; // NOLINT
const flt_t cutsq_custom = cutoff_custom * cutoff_custom;
const double delta = 0.01 * force->angstrom;
const int triclinic = domain->triclinic;
#if defined(_OPENMP)
#pragma omp parallel
@ -84,35 +88,70 @@ void NPairHalffullNewtonTrimIntel::build_t(NeighList *list,
const int * _noalias const jlist = firstneigh_full[i];
const int jnum = numneigh_full[i];
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma ivdep
#endif
for (int jj = 0; jj < jnum; jj++) {
const int joriginal = jlist[jj];
const int j = joriginal & NEIGHMASK;
int addme = 1;
if (j < nlocal) {
if (i > j) addme = 0;
} else {
if (x[j].z < ztmp) addme = 0;
if (x[j].z == ztmp) {
if (x[j].y < ytmp) addme = 0;
if (x[j].y == ytmp && x[j].x < xtmp) addme = 0;
if (!triclinic) {
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma ivdep
#endif
for (int jj = 0; jj < jnum; jj++) {
const int joriginal = jlist[jj];
const int j = joriginal & NEIGHMASK;
int addme = 1;
if (j < nlocal) {
if (i > j) addme = 0;
} else {
if (x[j].z < ztmp) addme = 0;
if (x[j].z == ztmp) {
if (x[j].y < ytmp) addme = 0;
if (x[j].y == ytmp && x[j].x < xtmp) addme = 0;
}
}
// trim to shorter cutoff
const flt_t delx = xtmp - x[j].x;
const flt_t dely = ytmp - x[j].y;
const flt_t delz = ztmp - x[j].z;
const flt_t rsq = delx * delx + dely * dely + delz * delz;
if (rsq > cutsq_custom) addme = 0;
if (addme)
neighptr[n++] = joriginal;
}
} else {
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma ivdep
#endif
for (int jj = 0; jj < jnum; jj++) {
const int joriginal = jlist[jj];
const int j = joriginal & NEIGHMASK;
int addme = 1;
if (j < nlocal) {
if (i > j) addme = 0;
} else {
if (fabs(x[j].z-ztmp) > delta) {
if (x[j].z < ztmp) addme = 0;
} else if (fabs(x[j].y-ytmp) > delta) {
if (x[j].y < ytmp) addme = 0;
} else {
if (x[j].x < xtmp) addme = 0;
}
}
// trim to shorter cutoff
// trim to shorter cutoff
const flt_t delx = xtmp - x[j].x;
const flt_t dely = ytmp - x[j].y;
const flt_t delz = ztmp - x[j].z;
const flt_t rsq = delx * delx + dely * dely + delz * delz;
const flt_t delx = xtmp - x[j].x;
const flt_t dely = ytmp - x[j].y;
const flt_t delz = ztmp - x[j].z;
const flt_t rsq = delx * delx + dely * dely + delz * delz;
if (rsq > cutsq_custom) addme = 0;
if (rsq > cutsq_custom) addme = 0;
if (addme)
neighptr[n++] = joriginal;
if (addme)
neighptr[n++] = joriginal;
}
}
ilist[ii] = i;
@ -235,7 +274,7 @@ void NPairHalffullNewtonTrimIntel::build_t3(NeighList *list, int *numhalf,
void NPairHalffullNewtonTrimIntel::build(NeighList *list)
{
if (_fix->three_body_neighbor() == 0) {
if (_fix->three_body_neighbor() == 0 || domain->triclinic) {
if (_fix->precision() == FixIntel::PREC_MODE_MIXED)
build_t(list, _fix->get_mixed_buffers());
else if (_fix->precision() == FixIntel::PREC_MODE_DOUBLE)

40
src/INTEL/npair_intel.cpp
View File

@ -204,6 +204,8 @@ void NPairIntel::bin_newton(const int offload, NeighList *list,
}
const int special_bound = sb;
const double delta = 0.01 * force->angstrom;
#ifdef _LMP_INTEL_OFFLOAD
const int * _noalias const binhead = this->binhead;
const int * _noalias const bins = this->bins;
@ -229,7 +231,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list,
in(ncache_stride,maxnbors,nthreads,maxspecial,nstencil,e_nall,offload) \
in(offload_end,separate_buffers,astart,aend,nlocal,molecular) \
in(ntypes,xperiodic,yperiodic,zperiodic,xprd_half,yprd_half,zprd_half) \
in(pack_width,special_bound) \
in(pack_width,special_bound,delta) \
out(overflow:length(5) alloc_if(0) free_if(0)) \
out(timer_compute:length(1) alloc_if(0) free_if(0)) \
signal(tag)
@ -331,7 +333,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list,
const flt_t ztmp = x[i].z;
const int itype = x[i].w;
tagint itag;
if (THREE) itag = tag[i];
if (THREE || (TRI && !FULL)) itag = tag[i];
const int ioffset = ntypes * itype;
const int ibin = atombin[i];
@ -365,7 +367,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list,
ty[u] = x[j].y;
tz[u] = x[j].z;
tjtype[u] = x[j].w;
if (THREE) ttag[u] = tag[j];
if (THREE || (TRI && !FULL)) ttag[u] = tag[j];
}
if (FULL == 0 && TRI != 1) {
@ -486,12 +488,32 @@ void NPairIntel::bin_newton(const int offload, NeighList *list,
// Triclinic
if (TRI) {
if (tz[u] < ztmp) addme = 0;
if (tz[u] == ztmp) {
if (ty[u] < ytmp) addme = 0;
if (ty[u] == ytmp) {
if (tx[u] < xtmp) addme = 0;
if (tx[u] == xtmp && j <= i) addme = 0;
if (FULL) {
if (tz[u] < ztmp) addme = 0;
if (tz[u] == ztmp) {
if (ty[u] < ytmp) addme = 0;
if (ty[u] == ytmp) {
if (tx[u] < xtmp) addme = 0;
if (tx[u] == xtmp && j <= i) addme = 0;
}
}
} else {
if (j <= i) addme = 0;
if (j >= nlocal) {
const tagint jtag = ttag[u];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) addme = 0;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) addme = 0;
} else {
if (fabs(tz[u]-ztmp) > delta) {
if (tz[u] < ztmp) addme = 0;
} else if (fabs(ty[u]-ytmp) > delta) {
if (ty[u] < ytmp) addme = 0;
} else {
if (tx[u] < xtmp) addme = 0;
}
}
}
}
}

3
src/KOKKOS/min_kokkos.cpp
View File

@ -59,6 +59,9 @@ void MinKokkos::init()
{
Min::init();
if (!fix_minimize->kokkosable)
error->all(FLERR,"KOKKOS package requires fix minimize/kk");
fix_minimize_kk = (FixMinimizeKokkos*) fix_minimize;
}

49
src/KOKKOS/npair_halffull_kokkos.cpp
View File

@ -18,6 +18,7 @@
#include "atom_masks.h"
#include "atom_vec.h"
#include "domain.h"
#include "force.h"
#include "neigh_list_kokkos.h"
#include <cmath>
@ -26,8 +27,8 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
template<class DeviceType, int NEWTON, int TRIM>
NPairHalffullKokkos<DeviceType,NEWTON,TRIM>::NPairHalffullKokkos(LAMMPS *lmp) : NPair(lmp) {
template<class DeviceType, int NEWTON, int TRI, int TRIM>
NPairHalffullKokkos<DeviceType,NEWTON,TRI,TRIM>::NPairHalffullKokkos(LAMMPS *lmp) : NPair(lmp) {
atomKK = (AtomKokkos *) atom;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
}
@ -41,13 +42,14 @@ NPairHalffullKokkos<DeviceType,NEWTON,TRIM>::NPairHalffullKokkos(LAMMPS *lmp) :
if ghost, also store neighbors of ghost atoms & set inum,gnum correctly
------------------------------------------------------------------------- */
template<class DeviceType, int NEWTON, int TRIM>
void NPairHalffullKokkos<DeviceType,NEWTON,TRIM>::build(NeighList *list)
template<class DeviceType, int NEWTON, int TRI, int TRIM>
void NPairHalffullKokkos<DeviceType,NEWTON,TRI,TRIM>::build(NeighList *list)
{
if (NEWTON || TRIM) {
x = atomKK->k_x.view<DeviceType>();
atomKK->sync(execution_space,X_MASK);
}
nlocal = atom->nlocal;
cutsq_custom = cutoff_custom*cutoff_custom;
@ -66,6 +68,8 @@ void NPairHalffullKokkos<DeviceType,NEWTON,TRIM>::build(NeighList *list)
d_numneigh = k_list->d_numneigh;
d_neighbors = k_list->d_neighbors;
delta = 0.01 * force->angstrom;
// loop over parent full list
copymode = 1;
@ -78,9 +82,9 @@ void NPairHalffullKokkos<DeviceType,NEWTON,TRIM>::build(NeighList *list)
k_list->k_ilist.template modify<DeviceType>();
}
template<class DeviceType, int NEWTON, int TRIM>
template<class DeviceType, int NEWTON, int TRI, int TRIM>
KOKKOS_INLINE_FUNCTION
void NPairHalffullKokkos<DeviceType,NEWTON,TRIM>::operator()(TagNPairHalffullCompute, const int &ii) const {
void NPairHalffullKokkos<DeviceType,NEWTON,TRI,TRIM>::operator()(TagNPairHalffullCompute, const int &ii) const {
int n = 0;
const int i = d_ilist_full(ii);
@ -92,6 +96,11 @@ void NPairHalffullKokkos<DeviceType,NEWTON,TRIM>::operator()(TagNPairHalffullCom
}
// loop over full neighbor list
// use i < j < nlocal to eliminate half the local/local interactions
// for triclinic, must use delta to eliminate half the local/ghost interactions
// cannot use I/J exact coord comparision as for orthog
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
const int jnum = d_numneigh_full(i);
const AtomNeighbors neighbors_i = AtomNeighbors(&d_neighbors(i,0),d_numneigh(i),
@ -103,6 +112,14 @@ void NPairHalffullKokkos<DeviceType,NEWTON,TRIM>::operator()(TagNPairHalffullCom
if (NEWTON) {
if (j < nlocal) {
if (i > j) continue;
} else if (TRI) {
if (fabs(x(j,2)-ztmp) > delta) {
if (x(j,2) < ztmp) continue;
} else if (fabs(x(j,1)-ytmp) > delta) {
if (x(j,1) < ytmp) continue;
} else {
if (x(j,0) < xtmp) continue;
}
} else {
if (x(j,2) < ztmp) continue;
if (x(j,2) == ztmp) {
@ -141,14 +158,18 @@ void NPairHalffullKokkos<DeviceType,NEWTON,TRIM>::operator()(TagNPairHalffullCom
}
namespace LAMMPS_NS {
template class NPairHalffullKokkos<LMPDeviceType,0,0>;
template class NPairHalffullKokkos<LMPDeviceType,0,1>;
template class NPairHalffullKokkos<LMPDeviceType,1,0>;
template class NPairHalffullKokkos<LMPDeviceType,1,1>;
template class NPairHalffullKokkos<LMPDeviceType,0,0,0>;
template class NPairHalffullKokkos<LMPDeviceType,0,0,1>;
template class NPairHalffullKokkos<LMPDeviceType,1,0,0>;
template class NPairHalffullKokkos<LMPDeviceType,1,0,1>;
template class NPairHalffullKokkos<LMPDeviceType,1,1,0>;
template class NPairHalffullKokkos<LMPDeviceType,1,1,1>;
#ifdef LMP_KOKKOS_GPU
template class NPairHalffullKokkos<LMPHostType,0,0>;
template class NPairHalffullKokkos<LMPHostType,0,1>;
template class NPairHalffullKokkos<LMPHostType,1,0>;
template class NPairHalffullKokkos<LMPHostType,1,1>;
template class NPairHalffullKokkos<LMPHostType,0,0,0>;
template class NPairHalffullKokkos<LMPHostType,0,0,1>;
template class NPairHalffullKokkos<LMPHostType,1,0,0>;
template class NPairHalffullKokkos<LMPHostType,1,0,1>;
template class NPairHalffullKokkos<LMPHostType,1,1,0>;
template class NPairHalffullKokkos<LMPHostType,1,1,1>;
#endif
}

252
src/KOKKOS/npair_halffull_kokkos.h
View File

@ -16,53 +16,79 @@
// Trim off
// Newton
// Newton, no triclinic
typedef NPairHalffullKokkos<LMPDeviceType,1,0> NPairKokkosHalffullNewtonDevice;
typedef NPairHalffullKokkos<LMPDeviceType,1,0,0> NPairKokkosHalffullNewtonDevice;
NPairStyle(halffull/newton/kk/device,
NPairKokkosHalffullNewtonDevice,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_KOKKOS_DEVICE);
NP_ORTHO | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,1,0> NPairKokkosHalffullNewtonHost;
typedef NPairHalffullKokkos<LMPHostType,1,0,0> NPairKokkosHalffullNewtonHost;
NPairStyle(halffull/newton/kk/host,
NPairKokkosHalffullNewtonHost,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_KOKKOS_HOST);
NP_ORTHO | NP_KOKKOS_HOST);
typedef NPairHalffullKokkos<LMPDeviceType,1,0> NPairKokkosHalffullNewtonDevice;
typedef NPairHalffullKokkos<LMPDeviceType,1,0,0> NPairKokkosHalffullNewtonDevice;
NPairStyle(halffull/newton/skip/kk/device,
NPairKokkosHalffullNewtonDevice,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_SKIP | NP_KOKKOS_DEVICE);
NP_ORTHO | NP_SKIP | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,1,0> NPairKokkosHalffullNewtonHost;
typedef NPairHalffullKokkos<LMPHostType,1,0,0> NPairKokkosHalffullNewtonHost;
NPairStyle(halffull/newton/skip/kk/host,
NPairKokkosHalffullNewtonHost,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_SKIP | NP_KOKKOS_HOST);
// Newton, triclinic
typedef NPairHalffullKokkos<LMPDeviceType,1,1,0> NPairKokkosHalffullNewtonTriDevice;
NPairStyle(halffull/newton/tri/kk/device,
NPairKokkosHalffullNewtonTriDevice,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,1,1,0> NPairKokkosHalffullNewtonTriHost;
NPairStyle(halffull/newton/tri/kk/host,
NPairKokkosHalffullNewtonTriHost,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_KOKKOS_HOST);
typedef NPairHalffullKokkos<LMPDeviceType,1,1,0> NPairKokkosHalffullNewtonTriDevice;
NPairStyle(halffull/newton/tri/skip/kk/device,
NPairKokkosHalffullNewtonTriDevice,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_SKIP | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,1,1,0> NPairKokkosHalffullNewtonTriHost;
NPairStyle(halffull/newton/tri/skip/kk/host,
NPairKokkosHalffullNewtonTriHost,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_SKIP | NP_KOKKOS_HOST);
// Newtoff
// Newtoff (can be triclinic but template param always set to 0)
typedef NPairHalffullKokkos<LMPDeviceType,0,0> NPairKokkosHalffullNewtoffDevice;
typedef NPairHalffullKokkos<LMPDeviceType,0,0,0> NPairKokkosHalffullNewtoffDevice;
NPairStyle(halffull/newtoff/kk/device,
NPairKokkosHalffullNewtoffDevice,
NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,0,0> NPairKokkosHalffullNewtoffHost;
typedef NPairHalffullKokkos<LMPHostType,0,0,0> NPairKokkosHalffullNewtoffHost;
NPairStyle(halffull/newtoff/kk/host,
NPairKokkosHalffullNewtoffHost,
NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_KOKKOS_HOST);
typedef NPairHalffullKokkos<LMPDeviceType,0,0> NPairKokkosHalffullNewtoffDevice;
typedef NPairHalffullKokkos<LMPDeviceType,0,0,0> NPairKokkosHalffullNewtoffDevice;
NPairStyle(halffull/newtoff/skip/kk/device,
NPairKokkosHalffullNewtoffDevice,
NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_SKIP | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,0,0> NPairKokkosHalffullNewtoffHost;
typedef NPairHalffullKokkos<LMPHostType,0,0,0> NPairKokkosHalffullNewtoffHost;
NPairStyle(halffull/newtoff/skip/kk/host,
NPairKokkosHalffullNewtoffHost,
NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
@ -70,166 +96,244 @@ NPairStyle(halffull/newtoff/skip/kk/host,
//************ Ghost **************
// Newton
// Newton, no triclinic
typedef NPairHalffullKokkos<LMPDeviceType,1,0> NPairKokkosHalffullNewtonGhostDevice;
typedef NPairHalffullKokkos<LMPDeviceType,1,0,0> NPairKokkosHalffullNewtonDevice;
NPairStyle(halffull/newton/ghost/kk/device,
NPairKokkosHalffullNewtonGhostDevice,
NPairKokkosHalffullNewtonDevice,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_KOKKOS_DEVICE);
NP_ORTHO | NP_GHOST | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,1,0> NPairKokkosHalffullNewtonHost;
typedef NPairHalffullKokkos<LMPHostType,1,0,0> NPairKokkosHalffullNewtonHost;
NPairStyle(halffull/newton/ghost/kk/host,
NPairKokkosHalffullNewtonHost,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_KOKKOS_HOST);
NP_ORTHO | NP_GHOST | NP_KOKKOS_HOST);
typedef NPairHalffullKokkos<LMPDeviceType,1,0> NPairKokkosHalffullNewtonGhostDevice;
typedef NPairHalffullKokkos<LMPDeviceType,1,0,0> NPairKokkosHalffullNewtonDevice;
NPairStyle(halffull/newton/skip/ghost/kk/device,
NPairKokkosHalffullNewtonGhostDevice,
NPairKokkosHalffullNewtonDevice,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_KOKKOS_DEVICE);
NP_ORTHO | NP_GHOST | NP_SKIP | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,1,0> NPairKokkosHalffullNewtonHost;
typedef NPairHalffullKokkos<LMPHostType,1,0,0> NPairKokkosHalffullNewtonHost;
NPairStyle(halffull/newton/skip/ghost/kk/host,
NPairKokkosHalffullNewtonHost,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_GHOST | NP_SKIP | NP_KOKKOS_HOST);
// Newton, triclinic
typedef NPairHalffullKokkos<LMPDeviceType,1,1,0> NPairKokkosHalffullNewtonTriDevice;
NPairStyle(halffull/newton/tri/ghost/kk/device,
NPairKokkosHalffullNewtonTriDevice,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,1,1,0> NPairKokkosHalffullNewtonTriHost;
NPairStyle(halffull/newton/tri/ghost/kk/host,
NPairKokkosHalffullNewtonTriHost,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_KOKKOS_HOST);
typedef NPairHalffullKokkos<LMPDeviceType,1,1,0> NPairKokkosHalffullNewtonTriDevice;
NPairStyle(halffull/newton/tri/skip/ghost/kk/device,
NPairKokkosHalffullNewtonTriDevice,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,1,1,0> NPairKokkosHalffullNewtonTriHost;
NPairStyle(halffull/newton/tri/skip/ghost/kk/host,
NPairKokkosHalffullNewtonTriHost,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_KOKKOS_HOST);
// Newtoff
// Newtoff (can be triclinic but template param always set to 0)
typedef NPairHalffullKokkos<LMPDeviceType,0,0> NPairKokkosHalffullNewtoffGhostDevice;
typedef NPairHalffullKokkos<LMPDeviceType,0,0,0> NPairKokkosHalffullNewtoffDevice;
NPairStyle(halffull/newtoff/ghost/kk/device,
NPairKokkosHalffullNewtoffGhostDevice,
NPairKokkosHalffullNewtoffDevice,
NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,0,0> NPairKokkosHalffullNewtoffHost;
typedef NPairHalffullKokkos<LMPHostType,0,0,0> NPairKokkosHalffullNewtoffHost;
NPairStyle(halffull/newtoff/ghost/kk/host,
NPairKokkosHalffullNewtoffHost,
NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_KOKKOS_HOST);
typedef NPairHalffullKokkos<LMPDeviceType,0,0> NPairKokkosHalffullNewtoffGhostDevice;
typedef NPairHalffullKokkos<LMPDeviceType,0,0,0> NPairKokkosHalffullNewtoffDevice;
NPairStyle(halffull/newtoff/skip/ghost/kk/device,
NPairKokkosHalffullNewtoffGhostDevice,
NPairKokkosHalffullNewtoffDevice,
NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,0,0> NPairKokkosHalffullNewtoffHost;
typedef NPairHalffullKokkos<LMPHostType,0,0,0> NPairKokkosHalffullNewtoffHost;
NPairStyle(halffull/newtoff/skip/ghost/kk/host,
NPairKokkosHalffullNewtoffHost,
NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_KOKKOS_HOST);
//************ Trim **************
// Newton
// Newton, no triclinic
typedef NPairHalffullKokkos<LMPDeviceType,1,1> NPairKokkosHalffullNewtonTrimDevice;
typedef NPairHalffullKokkos<LMPDeviceType,1,0,1> NPairKokkosHalffullNewtonTrimDevice;
NPairStyle(halffull/newton/trim/kk/device,
NPairKokkosHalffullNewtonTrimDevice,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_TRIM | NP_KOKKOS_DEVICE);
NP_ORTHO | NP_TRIM | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,1,1> NPairKokkosHalffullNewtonTrimHost;
typedef NPairHalffullKokkos<LMPHostType,1,0,1> NPairKokkosHalffullNewtonTrimHost;
NPairStyle(halffull/newton/trim/kk/host,
NPairKokkosHalffullNewtonTrimHost,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRIM | NP_KOKKOS_HOST);
typedef NPairHalffullKokkos<LMPDeviceType,1,0,1> NPairKokkosHalffullNewtonTrimDevice;
NPairStyle(halffull/newton/trim/skip/kk/device,
NPairKokkosHalffullNewtonTrimDevice,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_SKIP | NP_TRIM | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,1,0,1> NPairKokkosHalffullNewtonTrimHost;
NPairStyle(halffull/newton/trim/skip/kk/host,
NPairKokkosHalffullNewtonTrimHost,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_SKIP | NP_TRIM | NP_KOKKOS_HOST);
// Newton, triclinic
typedef NPairHalffullKokkos<LMPDeviceType,1,1,1> NPairKokkosHalffullNewtonTriTrimDevice;
NPairStyle(halffull/newton/tri/trim/kk/device,
NPairKokkosHalffullNewtonTriTrimDevice,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_TRIM | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,1,1,1> NPairKokkosHalffullNewtonTriTrimHost;
NPairStyle(halffull/newton/tri/trim/kk/host,
NPairKokkosHalffullNewtonTriTrimHost,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_TRIM | NP_KOKKOS_HOST);
typedef NPairHalffullKokkos<LMPDeviceType,1,1> NPairKokkosHalffullNewtonTrimDevice;
NPairStyle(halffull/newton/skip/trim/kk/device,
typedef NPairHalffullKokkos<LMPDeviceType,1,1,1> NPairKokkosHalffullNewtonTriTrimDevice;
NPairStyle(halffull/newton/tri/trim/skip/kk/device,
NPairKokkosHalffullNewtonTrimDevice,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_SKIP | NP_TRIM | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,1,1> NPairKokkosHalffullNewtonTrimHost;
NPairStyle(halffull/newton/skip/trim/kk/host,
NPairKokkosHalffullNewtonTrimHost,
typedef NPairHalffullKokkos<LMPHostType,1,1,1> NPairKokkosHalffullNewtonTriTrimHost;
NPairStyle(halffull/newton/tri/trim/skip/kk/host,
NPairKokkosHalffullNewtonTriTrimHost,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_SKIP | NP_TRIM | NP_KOKKOS_HOST);
// Newtoff
// Newtoff (can be triclinic but template param always set to 0)
typedef NPairHalffullKokkos<LMPDeviceType,0,1> NPairKokkosHalffullNewtoffTrimDevice;
typedef NPairHalffullKokkos<LMPDeviceType,0,0,1> NPairKokkosHalffullNewtoffTrimDevice;
NPairStyle(halffull/newtoff/trim/kk/device,
NPairKokkosHalffullNewtoffTrimDevice,
NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_TRIM | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,0,1> NPairKokkosHalffullNewtoffTrimHost;
typedef NPairHalffullKokkos<LMPHostType,0,0,1> NPairKokkosHalffullNewtoffTrimHost;
NPairStyle(halffull/newtoff/trim/kk/host,
NPairKokkosHalffullNewtoffTrimHost,
NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_TRIM | NP_KOKKOS_HOST);
typedef NPairHalffullKokkos<LMPDeviceType,0,1> NPairKokkosHalffullNewtoffTrimDevice;
NPairStyle(halffull/newtoff/skip/trim/kk/device,
typedef NPairHalffullKokkos<LMPDeviceType,0,0,1> NPairKokkosHalffullNewtoffTrimDevice;
NPairStyle(halffull/newtoff/trim/skip/kk/device,
NPairKokkosHalffullNewtoffTrimDevice,
NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_SKIP | NP_TRIM | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,0,1> NPairKokkosHalffullNewtoffTrimHost;
NPairStyle(halffull/newtoff/skip/trim/kk/host,
typedef NPairHalffullKokkos<LMPHostType,0,0,1> NPairKokkosHalffullNewtoffTrimHost;
NPairStyle(halffull/newtoff/trim/skip/kk/host,
NPairKokkosHalffullNewtoffTrimHost,
NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_SKIP | NP_TRIM | NP_KOKKOS_HOST);
//************ Ghost **************
// Newton
// Newton, no triclinic
typedef NPairHalffullKokkos<LMPDeviceType,1,1> NPairKokkosHalffullNewtonGhostTrimDevice;
NPairStyle(halffull/newton/ghost/trim/kk/device,
NPairKokkosHalffullNewtonGhostTrimDevice,
typedef NPairHalffullKokkos<LMPDeviceType,1,0,1> NPairKokkosHalffullNewtonTrimDevice;
NPairStyle(halffull/newton/tri/trim/ghost/kk/device,
NPairKokkosHalffullNewtonTrimDevice,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_GHOST | NP_TRIM | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,1,0,1> NPairKokkosHalffullNewtonTrimHost;
NPairStyle(halffull/newton/trim/ghost/kk/host,
NPairKokkosHalffullNewtonTrimHost,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_GHOST | NP_TRIM | NP_KOKKOS_HOST);
typedef NPairHalffullKokkos<LMPDeviceType,1,0,1> NPairKokkosHalffullNewtonTrimDevice;
NPairStyle(halffull/newton/trim/skip/ghost/kk/device,
NPairKokkosHalffullNewtonTrimDevice,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_GHOST | NP_SKIP | NP_TRIM | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,1,0,1> NPairKokkosHalffullNewtonTrimHost;
NPairStyle(halffull/newton/trim/skip/ghost/kk/host,
NPairKokkosHalffullNewtonTrimHost,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_GHOST | NP_SKIP | NP_TRIM | NP_KOKKOS_HOST);
// Newton, triclinic
typedef NPairHalffullKokkos<LMPDeviceType,1,1,1> NPairKokkosHalffullNewtonTriTrimDevice;
NPairStyle(halffull/newton/tri/trim/ghost/kk/device,
NPairKokkosHalffullNewtonTriTrimDevice,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_TRIM | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,1,1> NPairKokkosHalffullNewtonTrimHost;
NPairStyle(halffull/newton/ghost/trim/kk/host,
NPairKokkosHalffullNewtonTrimHost,
typedef NPairHalffullKokkos<LMPHostType,1,1,1> NPairKokkosHalffullNewtonTriTrimHost;
NPairStyle(halffull/newton/tri/trim/ghost/kk/host,
NPairKokkosHalffullNewtonTriTrimHost,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_TRIM | NP_KOKKOS_HOST);
typedef NPairHalffullKokkos<LMPDeviceType,1,1> NPairKokkosHalffullNewtonGhostTrimDevice;
NPairStyle(halffull/newton/skip/ghost/trim/kk/device,
NPairKokkosHalffullNewtonGhostTrimDevice,
typedef NPairHalffullKokkos<LMPDeviceType,1,1,1> NPairKokkosHalffullNewtonTriTrimDevice;
NPairStyle(halffull/newton/tri/trim/skip/ghost/kk/device,
NPairKokkosHalffullNewtonTriTrimDevice,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_TRIM | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,1,1> NPairKokkosHalffullNewtonTrimHost;
NPairStyle(halffull/newton/skip/ghost/trim/kk/host,
NPairKokkosHalffullNewtonTrimHost,
typedef NPairHalffullKokkos<LMPHostType,1,1,1> NPairKokkosHalffullNewtonTriTrimHost;
NPairStyle(halffull/newton/tri/trim/skip/ghost/kk/host,
NPairKokkosHalffullNewtonTriTrimHost,
NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_TRIM | NP_KOKKOS_HOST);
// Newtoff
// Newtoff (can be triclinic but template param always set to 0)
typedef NPairHalffullKokkos<LMPDeviceType,0,1> NPairKokkosHalffullNewtoffGhostTrimDevice;
NPairStyle(halffull/newtoff/ghost/trim/kk/device,
NPairKokkosHalffullNewtoffGhostTrimDevice,
typedef NPairHalffullKokkos<LMPDeviceType,0,0,1> NPairKokkosHalffullNewtoffTrimDevice;
NPairStyle(halffull/newtoff/trim/ghost/kk/device,
NPairKokkosHalffullNewtoffTrimDevice,
NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_TRIM | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,0,1> NPairKokkosHalffullNewtoffTrimHost;
NPairStyle(halffull/newtoff/ghost/trim/kk/host,
typedef NPairHalffullKokkos<LMPHostType,0,0,1> NPairKokkosHalffullNewtoffTrimHost;
NPairStyle(halffull/newtoff/trim/ghost/kk/host,
NPairKokkosHalffullNewtoffTrimHost,
NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_TRIM | NP_KOKKOS_HOST);
typedef NPairHalffullKokkos<LMPDeviceType,0,1> NPairKokkosHalffullNewtoffGhostTrimDevice;
NPairStyle(halffull/newtoff/skip/ghost/trim/kk/device,
NPairKokkosHalffullNewtoffGhostTrimDevice,
typedef NPairHalffullKokkos<LMPDeviceType,0,0,1> NPairKokkosHalffullNewtoffTrimDevice;
NPairStyle(halffull/newtoff/trim/skip/ghost/kk/device,
NPairKokkosHalffullNewtoffTrimDevice,
NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_TRIM | NP_KOKKOS_DEVICE);
typedef NPairHalffullKokkos<LMPHostType,0,1> NPairKokkosHalffullNewtoffTrimHost;
NPairStyle(halffull/newtoff/skip/ghost/trim/kk/host,
typedef NPairHalffullKokkos<LMPHostType,0,0,1> NPairKokkosHalffullNewtoffTrimHost;
NPairStyle(halffull/newtoff/trim/skip/ghost/kk/host,
NPairKokkosHalffullNewtoffTrimHost,
NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_TRIM | NP_KOKKOS_HOST);
// clang-format on
#else
@ -244,7 +348,7 @@ namespace LAMMPS_NS {
struct TagNPairHalffullCompute{};
template<class DeviceType, int NEWTON, int TRIM>
template<class DeviceType, int NEWTON, int TRI, int TRIM>
class NPairHalffullKokkos : public NPair {
public:
typedef DeviceType device_type;
@ -257,8 +361,8 @@ class NPairHalffullKokkos : public NPair {
void operator()(TagNPairHalffullCompute, const int&) const;
private:
int nlocal;
double cutsq_custom;
int nlocal,triclinic;
double cutsq_custom,delta;
typename AT::t_x_array_randomread x;

161
src/KOKKOS/npair_kokkos.cpp
View File

@ -155,6 +155,8 @@ void NPairKokkos<DeviceType,HALF,NEWTON,GHOST,TRI,SIZE>::build(NeighList *list_)
list->grow(nall);
const double delta = 0.01 * force->angstrom;
NeighborKokkosExecute<DeviceType>
data(*list,
k_cutneighsq.view<DeviceType>(),
@ -176,7 +178,7 @@ void NPairKokkos<DeviceType,HALF,NEWTON,GHOST,TRI,SIZE>::build(NeighList *list_)
atomKK->molecular,
nbinx,nbiny,nbinz,mbinx,mbiny,mbinz,mbinxlo,mbinylo,mbinzlo,
bininvx,bininvy,bininvz,
exclude, nex_type,
delta, exclude, nex_type,
k_ex1_type.view<DeviceType>(),
k_ex2_type.view<DeviceType>(),
k_ex_type.view<DeviceType>(),
@ -217,6 +219,8 @@ void NPairKokkos<DeviceType,HALF,NEWTON,GHOST,TRI,SIZE>::build(NeighList *list_)
atomKK->sync(Device,X_MASK|RADIUS_MASK|TYPE_MASK);
}
if (HALF && NEWTON && TRI) atomKK->sync(Device,TAG_MASK);
data.special_flag[0] = special_flag[0];
data.special_flag[1] = special_flag[1];
data.special_flag[2] = special_flag[2];
@ -261,7 +265,7 @@ void NPairKokkos<DeviceType,HALF,NEWTON,GHOST,TRI,SIZE>::build(NeighList *list_)
//#endif
} else {
if (SIZE) {
NPairKokkosBuildFunctorSize<DeviceType,HALF,NEWTON,TRI> f(data,atoms_per_bin * 6 * sizeof(X_FLOAT) * factor);
NPairKokkosBuildFunctorSize<DeviceType,HALF,NEWTON,TRI> f(data,atoms_per_bin * 7 * sizeof(X_FLOAT) * factor);
#ifdef LMP_KOKKOS_GPU
if (ExecutionSpaceFromDevice<DeviceType>::space == Device) {
int team_size = atoms_per_bin*factor;
@ -279,7 +283,7 @@ void NPairKokkos<DeviceType,HALF,NEWTON,GHOST,TRI,SIZE>::build(NeighList *list_)
Kokkos::parallel_for(nall, f);
#endif
} else {
NPairKokkosBuildFunctor<DeviceType,HALF,NEWTON,TRI> f(data,atoms_per_bin * 5 * sizeof(X_FLOAT) * factor);
NPairKokkosBuildFunctor<DeviceType,HALF,NEWTON,TRI> f(data,atoms_per_bin * 6 * sizeof(X_FLOAT) * factor);
#ifdef LMP_KOKKOS_GPU
if (ExecutionSpaceFromDevice<DeviceType>::space == Device) {
int team_size = atoms_per_bin*factor;
@ -414,6 +418,8 @@ void NeighborKokkosExecute<DeviceType>::
const X_FLOAT ytmp = x(i, 1);
const X_FLOAT ztmp = x(i, 2);
const int itype = type(i);
tagint itag;
if (HalfNeigh && Newton && Tri) itag = tag(i);
const int ibin = c_atom2bin(i);
@ -484,13 +490,29 @@ void NeighborKokkosExecute<DeviceType>::
if (HalfNeigh && !Newton && j <= i) continue;
if (!HalfNeigh && j == i) continue;
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
if (HalfNeigh && Newton && Tri) {
if (x(j,2) < ztmp) continue;
if (x(j,2) == ztmp) {
if (x(j,1) < ytmp) continue;
if (x(j,1) == ytmp) {
if (x(j,0) < xtmp) continue;
if (x(j,0) == xtmp && j <= i) continue;
if (j <= i) continue;
if (j >= nlocal) {
const tagint jtag = tag(j);
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x(j,2)-ztmp) > delta) {
if (x(j,2) < ztmp) continue;
} else if (fabs(x(j,1)-ytmp) > delta) {
if (x(j,1) < ytmp) continue;
} else {
if (x(j,0) < xtmp) continue;
}
}
}
}
@ -568,8 +590,9 @@ void NeighborKokkosExecute<DeviceType>::build_ItemGPU(typename Kokkos::TeamPolic
size_t sharedsize) const
{
auto* sharedmem = static_cast<X_FLOAT *>(dev.team_shmem().get_shmem(sharedsize));
/* loop over atoms in i's bin,
*/
// loop over atoms in i's bin
const int atoms_per_bin = c_bins.extent(1);
const int BINS_PER_TEAM = dev.team_size()/atoms_per_bin <1?1:dev.team_size()/atoms_per_bin;
const int TEAMS_PER_BIN = atoms_per_bin/dev.team_size()<1?1:atoms_per_bin/dev.team_size();
@ -579,15 +602,14 @@ void NeighborKokkosExecute<DeviceType>::build_ItemGPU(typename Kokkos::TeamPolic
if (ibin >= mbins) return;
X_FLOAT* other_x = sharedmem + 5*atoms_per_bin*MY_BIN;
int* other_id = (int*) &other_x[4 * atoms_per_bin];
X_FLOAT* other_x = sharedmem + 6*atoms_per_bin*MY_BIN;
int* other_id = (int*) &other_x[5 * atoms_per_bin];
int bincount_current = c_bincount[ibin];
for (int kk = 0; kk < TEAMS_PER_BIN; kk++) {
const int MY_II = dev.team_rank()%atoms_per_bin+kk*dev.team_size();
const int i = MY_II < bincount_current ? c_bins(ibin, MY_II) : -1;
/* if necessary, goto next page and add pages */
int n = 0;
@ -595,6 +617,7 @@ void NeighborKokkosExecute<DeviceType>::build_ItemGPU(typename Kokkos::TeamPolic
X_FLOAT ytmp;
X_FLOAT ztmp;
int itype;
tagint itag;
const int index = (i >= 0 && i < nlocal) ? i : 0;
const AtomNeighbors neighbors_i = neigh_transpose ?
neigh_list.get_neighbors_transpose(index) : neigh_list.get_neighbors(index);
@ -608,6 +631,10 @@ void NeighborKokkosExecute<DeviceType>::build_ItemGPU(typename Kokkos::TeamPolic
other_x[MY_II + atoms_per_bin] = ytmp;
other_x[MY_II + 2 * atoms_per_bin] = ztmp;
other_x[MY_II + 3 * atoms_per_bin] = itype;
if (HalfNeigh && Newton && Tri) {
itag = tag(i);
other_x[MY_II + 4 * atoms_per_bin] = itag;
}
}
other_id[MY_II] = i;
@ -695,6 +722,8 @@ void NeighborKokkosExecute<DeviceType>::build_ItemGPU(typename Kokkos::TeamPolic
other_x[MY_II + atoms_per_bin] = x(j, 1);
other_x[MY_II + 2 * atoms_per_bin] = x(j, 2);
other_x[MY_II + 3 * atoms_per_bin] = type(j);
if (HalfNeigh && Newton && Tri)
other_x[MY_II + 4 * atoms_per_bin] = tag(j);
}
other_id[MY_II] = j;
@ -708,13 +737,29 @@ void NeighborKokkosExecute<DeviceType>::build_ItemGPU(typename Kokkos::TeamPolic
if (HalfNeigh && !Newton && j <= i) continue;
if (!HalfNeigh && j == i) continue;
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
if (HalfNeigh && Newton && Tri) {
if (x(j,2) < ztmp) continue;
if (x(j,2) == ztmp) {
if (x(j,1) < ytmp) continue;
if (x(j,1) == ytmp) {
if (x(j,0) < xtmp) continue;
if (x(j,0) == xtmp && j <= i) continue;
if (j <= i) continue;
if (j >= nlocal) {
const tagint jtag = other_x[m + 4 * atoms_per_bin];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x(j,2)-ztmp) > delta) {
if (x(j,2) < ztmp) continue;
} else if (fabs(x(j,1)-ytmp) > delta) {
if (x(j,1) < ytmp) continue;
} else {
if (x(j,0) < xtmp) continue;
}
}
}
}
@ -905,6 +950,7 @@ void NeighborKokkosExecute<DeviceType>::build_ItemGhostGPU(typename Kokkos::Team
size_t sharedsize) const
{
auto* sharedmem = static_cast<X_FLOAT *>(dev.team_shmem().get_shmem(sharedsize));
// loop over atoms in i's bin
const int atoms_per_bin = c_bins.extent(1);
@ -1084,6 +1130,8 @@ void NeighborKokkosExecute<DeviceType>::
const X_FLOAT ztmp = x(i, 2);
const X_FLOAT radi = radius(i);
const int itype = type(i);
tagint itag;
if (HalfNeigh && Newton && Tri) itag = tag(i);
const int ibin = c_atom2bin(i);
@ -1167,13 +1215,29 @@ void NeighborKokkosExecute<DeviceType>::
if (HalfNeigh && !Newton && j <= i) continue;
if (!HalfNeigh && j == i) continue;
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
if (HalfNeigh && Newton && Tri) {
if (x(j,2) < ztmp) continue;
if (x(j,2) == ztmp) {
if (x(j,1) < ytmp) continue;
if (x(j,1) == ytmp) {
if (x(j,0) < xtmp) continue;
if (x(j,0) == xtmp && j <= i) continue;
if (j <= i) continue;
if (j >= nlocal) {
const tagint jtag = tag(j);
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x(j,2)-ztmp) > delta) {
if (x(j,2) < ztmp) continue;
} else if (fabs(x(j,1)-ytmp) > delta) {
if (x(j,1) < ytmp) continue;
} else {
if (x(j,0) < xtmp) continue;
}
}
}
}
@ -1245,8 +1309,9 @@ void NeighborKokkosExecute<DeviceType>::build_ItemSizeGPU(typename Kokkos::TeamP
size_t sharedsize) const
{
auto* sharedmem = static_cast<X_FLOAT *>(dev.team_shmem().get_shmem(sharedsize));
/* loop over atoms in i's bin,
*/
// loop over atoms in i's bin
const int atoms_per_bin = c_bins.extent(1);
const int BINS_PER_TEAM = dev.team_size()/atoms_per_bin <1?1:dev.team_size()/atoms_per_bin;
const int TEAMS_PER_BIN = atoms_per_bin/dev.team_size()<1?1:atoms_per_bin/dev.team_size();
@ -1256,15 +1321,14 @@ void NeighborKokkosExecute<DeviceType>::build_ItemSizeGPU(typename Kokkos::TeamP
if (ibin >= mbins) return;
X_FLOAT* other_x = sharedmem + 6*atoms_per_bin*MY_BIN;
int* other_id = (int*) &other_x[5 * atoms_per_bin];
X_FLOAT* other_x = sharedmem + 7*atoms_per_bin*MY_BIN;
int* other_id = (int*) &other_x[6 * atoms_per_bin];
int bincount_current = c_bincount[ibin];
for (int kk = 0; kk < TEAMS_PER_BIN; kk++) {
const int MY_II = dev.team_rank()%atoms_per_bin+kk*dev.team_size();
const int i = MY_II < bincount_current ? c_bins(ibin, MY_II) : -1;
/* if necessary, goto next page and add pages */
int n = 0;
@ -1273,6 +1337,7 @@ void NeighborKokkosExecute<DeviceType>::build_ItemSizeGPU(typename Kokkos::TeamP
X_FLOAT ztmp;
X_FLOAT radi;
int itype;
tagint itag;
const int index = (i >= 0 && i < nlocal) ? i : 0;
const AtomNeighbors neighbors_i = neigh_transpose ?
neigh_list.get_neighbors_transpose(index) : neigh_list.get_neighbors(index);
@ -1289,6 +1354,10 @@ void NeighborKokkosExecute<DeviceType>::build_ItemSizeGPU(typename Kokkos::TeamP
other_x[MY_II + 2 * atoms_per_bin] = ztmp;
other_x[MY_II + 3 * atoms_per_bin] = itype;
other_x[MY_II + 4 * atoms_per_bin] = radi;
if (HalfNeigh && Newton && Tri) {
itag = tag(i);
other_x[MY_II + 5 * atoms_per_bin] = itag;
}
}
other_id[MY_II] = i;
#if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP)
@ -1381,6 +1450,8 @@ void NeighborKokkosExecute<DeviceType>::build_ItemSizeGPU(typename Kokkos::TeamP
other_x[MY_II + 2 * atoms_per_bin] = x(j, 2);
other_x[MY_II + 3 * atoms_per_bin] = type(j);
other_x[MY_II + 4 * atoms_per_bin] = radius(j);
if (HalfNeigh && Newton && Tri)
other_x[MY_II + 5 * atoms_per_bin] = tag(j);
}
other_id[MY_II] = j;
@ -1394,13 +1465,29 @@ void NeighborKokkosExecute<DeviceType>::build_ItemSizeGPU(typename Kokkos::TeamP
if (HalfNeigh && !Newton && j <= i) continue;
if (!HalfNeigh && j == i) continue;
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
if (HalfNeigh && Newton && Tri) {
if (x(j,2) < ztmp) continue;
if (x(j,2) == ztmp) {
if (x(j,1) < ytmp) continue;
if (x(j,1) == ytmp) {
if (x(j,0) < xtmp) continue;
if (x(j,0) == xtmp && j <= i) continue;
if (j <= i) continue;
if (j >= nlocal) {
const tagint jtag = other_x[m + 5 * atoms_per_bin];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x(j,2)-ztmp) > delta) {
if (x(j,2) < ztmp) continue;
} else if (fabs(x(j,1)-ytmp) > delta) {
if (x(j,1) < ytmp) continue;
} else {
if (x(j,0) < xtmp) continue;
}
}
}
}

6
src/KOKKOS/npair_kokkos.h
View File

@ -189,6 +189,8 @@ class NeighborKokkosExecute
public:
NeighListKokkos<DeviceType> neigh_list;
const double delta;
// data from Neighbor class
const typename AT::t_xfloat_2d_randomread cutneighsq;
@ -282,7 +284,7 @@ class NeighborKokkosExecute
const int & _mbinx,const int & _mbiny,const int & _mbinz,
const int & _mbinxlo,const int & _mbinylo,const int & _mbinzlo,
const X_FLOAT &_bininvx,const X_FLOAT &_bininvy,const X_FLOAT &_bininvz,
const int & _exclude,const int & _nex_type,
const double &_delta,const int & _exclude,const int & _nex_type,
const typename AT::t_int_1d_const & _ex1_type,
const typename AT::t_int_1d_const & _ex2_type,
const typename AT::t_int_2d_const & _ex_type,
@ -301,7 +303,7 @@ class NeighborKokkosExecute
const typename ArrayTypes<LMPHostType>::t_int_scalar _h_resize,
const typename AT::t_int_scalar _new_maxneighs,
const typename ArrayTypes<LMPHostType>::t_int_scalar _h_new_maxneighs):
neigh_list(_neigh_list), cutneighsq(_cutneighsq),exclude(_exclude),
neigh_list(_neigh_list), cutneighsq(_cutneighsq),delta(_delta),exclude(_exclude),
nex_type(_nex_type),ex1_type(_ex1_type),ex2_type(_ex2_type),
ex_type(_ex_type),nex_group(_nex_group),
ex1_bit(_ex1_bit),ex2_bit(_ex2_bit),

7
src/KOKKOS/pair_snap_kokkos_impl.h
View File

@ -63,10 +63,6 @@ PairSNAPKokkos<DeviceType, real_type, vector_length>::PairSNAPKokkos(LAMMPS *lmp
datamask_read = EMPTY_MASK;
datamask_modify = EMPTY_MASK;
k_cutsq = tdual_fparams("PairSNAPKokkos::cutsq",atom->ntypes+1,atom->ntypes+1);
auto d_cutsq = k_cutsq.template view<DeviceType>();
rnd_cutsq = d_cutsq;
host_flag = (execution_space == Host);
}
@ -546,6 +542,9 @@ void PairSNAPKokkos<DeviceType, real_type, vector_length>::allocate()
int n = atom->ntypes;
MemKK::realloc_kokkos(d_map,"PairSNAPKokkos::map",n+1);
MemKK::realloc_kokkos(k_cutsq,"PairSNAPKokkos::cutsq",n+1,n+1);
rnd_cutsq = k_cutsq.template view<DeviceType>();
}

47
src/OPENMP/npair_half_bin_newton_tri_omp.cpp
View File

@ -12,16 +12,18 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "omp_compat.h"
#include "npair_half_bin_newton_tri_omp.h"
#include "npair_omp.h"
#include "neigh_list.h"
#include "omp_compat.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "domain.h"
#include "my_page.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neigh_list.h"
using namespace LAMMPS_NS;
@ -40,6 +42,7 @@ void NPairHalfBinNewtonTriOmp::build(NeighList *list)
const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
const int molecular = atom->molecular;
const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
const double delta = 0.01 * force->angstrom;
NPAIR_OMP_INIT;
#if defined(_OPENMP)
@ -48,12 +51,10 @@ void NPairHalfBinNewtonTriOmp::build(NeighList *list)
NPAIR_OMP_SETUP(nlocal);
int i,j,k,n,itype,jtype,ibin,which,imol,iatom;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr;
// loop over each atom, storing neighbors
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
@ -79,6 +80,7 @@ void NPairHalfBinNewtonTriOmp::build(NeighList *list)
n = 0;
neighptr = ipage.vget();
itag = tag[i];
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
@ -90,20 +92,31 @@ void NPairHalfBinNewtonTriOmp::build(NeighList *list)
}
// loop over all atoms in bins in stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
if (j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}

119
src/OPENMP/npair_half_multi_newton_tri_omp.cpp
View File

@ -12,17 +12,19 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "omp_compat.h"
#include "npair_half_multi_newton_tri_omp.h"
#include "npair_omp.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "domain.h"
#include "my_page.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neigh_list.h"
#include "neighbor.h"
#include "npair_omp.h"
#include "omp_compat.h"
using namespace LAMMPS_NS;
@ -43,6 +45,7 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list)
const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
const int molecular = atom->molecular;
const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
const double delta = 0.01 * force->angstrom;
NPAIR_OMP_INIT;
#if defined(_OPENMP)
@ -51,13 +54,11 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list)
NPAIR_OMP_SETUP(nlocal);
int i,j,k,n,itype,jtype,ibin,jbin,icollection,jcollection,which,ns,imol,iatom;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*s;
int js;
// loop over each atom, storing neighbors
int *collection = neighbor->collection;
double **x = atom->x;
int *type = atom->type;
@ -84,6 +85,7 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list)
n = 0;
neighptr = ipage.vget();
itag = tag[i];
itype = type[i];
icollection = collection[i];
xtmp = x[i][0];
@ -98,65 +100,80 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list)
ibin = atom2bin[i];
// loop through stencils for all collections
for (jcollection = 0; jcollection < ncollections; jcollection++) {
// if same collection use own bin
if (icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// loop over all atoms in bins in stencil
// stencil is empty if i larger than j
// stencil is half if i same size as j
// stencil is full if i smaller than j
// if half: pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// for triclinic:
// stencil is empty if i larger than j
// stencil is full if i smaller than j
// stencil is full if i same size as j
// for i smaller than j:
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
// if same size (same collection), use half stencil
if (cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]){
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
// if same size (same collection), exclude half of interactions
if (cutcollectionsq[icollection][icollection] ==
cutcollectionsq[jcollection][jcollection]) {
if (j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
}
ilist[i] = i;

48
src/OPENMP/npair_half_multi_old_newton_tri_omp.cpp
View File

@ -15,13 +15,15 @@
#include "omp_compat.h"
#include "npair_half_multi_old_newton_tri_omp.h"
#include "npair_omp.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "domain.h"
#include "my_page.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neigh_list.h"
using namespace LAMMPS_NS;
@ -42,6 +44,7 @@ void NPairHalfMultiOldNewtonTriOmp::build(NeighList *list)
const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
const int molecular = atom->molecular;
const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
const double delta = 0.01 * force->angstrom;
NPAIR_OMP_INIT;
#if defined(_OPENMP)
@ -50,13 +53,11 @@ void NPairHalfMultiOldNewtonTriOmp::build(NeighList *list)
NPAIR_OMP_SETUP(nlocal);
int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*s;
double *cutsq,*distsq;
// loop over each atom, storing neighbors
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
@ -82,6 +83,7 @@ void NPairHalfMultiOldNewtonTriOmp::build(NeighList *list)
n = 0;
neighptr = ipage.vget();
itag = tag[i];
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
@ -92,13 +94,12 @@ void NPairHalfMultiOldNewtonTriOmp::build(NeighList *list)
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms in bins, including self, in stencil
// skip if i,j neighbor cutoff is less than bin distance
// bins below self are excluded from stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// loop over all atoms in bins in stencil
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
ibin = atom2bin[i];
s = stencil_multi_old[itype];
@ -109,12 +110,21 @@ void NPairHalfMultiOldNewtonTriOmp::build(NeighList *list)
for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
jtype = type[j];
if (cutsq[jtype] < distsq[k]) continue;
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}

25
src/OPENMP/npair_half_nsq_newton_omp.cpp
View File

@ -15,14 +15,16 @@
#include "omp_compat.h"
#include "npair_half_nsq_newton_omp.h"
#include "npair_omp.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "molecule.h"
#include "domain.h"
#include "my_page.h"
#include "error.h"
#include "neigh_list.h"
using namespace LAMMPS_NS;
@ -42,6 +44,8 @@ void NPairHalfNsqNewtonOmp::build(NeighList *list)
const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0;
const int molecular = atom->molecular;
const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
const double delta = 0.01 * force->angstrom;
const int triclinic = domain->triclinic;
NPAIR_OMP_INIT;
#if defined(_OPENMP)
@ -54,8 +58,6 @@ void NPairHalfNsqNewtonOmp::build(NeighList *list)
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr;
// loop over each atom, storing neighbors
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
@ -95,7 +97,12 @@ void NPairHalfNsqNewtonOmp::build(NeighList *list)
}
// loop over remaining atoms, owned and ghost
// use itag/jtap comparision to eliminate half the interactions
// itag = jtag is possible for long cutoffs that include images of self
// for triclinic, must use delta to eliminate half the I/J interactions
// cannot use I/J exact coord comparision as for orthog
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
for (j = i+1; j < nall; j++) {
if (includegroup && !(mask[j] & bitmask)) continue;
@ -106,6 +113,14 @@ void NPairHalfNsqNewtonOmp::build(NeighList *list)
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else if (triclinic) {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {

45
src/OPENMP/npair_half_respa_bin_newton_tri_omp.cpp
View File

@ -15,13 +15,15 @@
#include "omp_compat.h"
#include "npair_half_respa_bin_newton_tri_omp.h"
#include "npair_omp.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "domain.h"
#include "my_page.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neigh_list.h"
using namespace LAMMPS_NS;
@ -42,6 +44,7 @@ void NPairHalfRespaBinNewtonTriOmp::build(NeighList *list)
const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
const int molecular = atom->molecular;
const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
const double delta = 0.01 * force->angstrom;
NPAIR_OMP_INIT;
@ -53,12 +56,10 @@ void NPairHalfRespaBinNewtonTriOmp::build(NeighList *list)
NPAIR_OMP_SETUP(nlocal);
int i,j,k,n,itype,jtype,ibin,n_inner,n_middle,imol,iatom;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*neighptr_inner,*neighptr_middle;
// loop over each atom, storing neighbors
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
@ -111,6 +112,7 @@ void NPairHalfRespaBinNewtonTriOmp::build(NeighList *list)
neighptr_middle = ipage_middle->vget();
}
itag = tag[i];
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
@ -122,20 +124,31 @@ void NPairHalfRespaBinNewtonTriOmp::build(NeighList *list)
}
// loop over all atoms in bins in stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
if (j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}

29
src/OPENMP/npair_half_respa_nsq_newton_omp.cpp
View File

@ -15,21 +15,22 @@
#include "omp_compat.h"
#include "npair_half_respa_nsq_newton_omp.h"
#include "npair_omp.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "molecule.h"
#include "domain.h"
#include "my_page.h"
#include "error.h"
#include "neigh_list.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
NPairHalfRespaNsqNewtonOmp::NPairHalfRespaNsqNewtonOmp(LAMMPS *lmp) :
NPair(lmp) {}
NPairHalfRespaNsqNewtonOmp::NPairHalfRespaNsqNewtonOmp(LAMMPS *lmp) : NPair(lmp) {}
/* ----------------------------------------------------------------------
multiple respa lists
@ -45,6 +46,8 @@ void NPairHalfRespaNsqNewtonOmp::build(NeighList *list)
const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0;
const int molecular = atom->molecular;
const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
const double delta = 0.01 * force->angstrom;
const int triclinic = domain->triclinic;
NPAIR_OMP_INIT;
@ -60,8 +63,6 @@ void NPairHalfRespaNsqNewtonOmp::build(NeighList *list)
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*neighptr_inner,*neighptr_middle;
// loop over each atom, storing neighbors
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
@ -128,6 +129,12 @@ void NPairHalfRespaNsqNewtonOmp::build(NeighList *list)
}
// loop over remaining atoms, owned and ghost
// use itag/jtap comparision to eliminate half the interactions
// itag = jtag is possible for long cutoffs that include images of self
// for triclinic, must use delta to eliminate half the I/J interactions
// cannot use I/J exact coord comparision as for orthog
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
for (j = i+1; j < nall; j++) {
if (includegroup && !(mask[j] & bitmask)) continue;
@ -138,6 +145,14 @@ void NPairHalfRespaNsqNewtonOmp::build(NeighList *list)
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else if (triclinic) {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {

38
src/OPENMP/npair_half_size_bin_newton_tri_omp.cpp
View File

@ -18,6 +18,7 @@
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neigh_list.h"
@ -46,6 +47,7 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list)
const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
const int history = list->history;
const int mask_history = 1 << HISTBITS;
const double delta = 0.01 * force->angstrom;
NPAIR_OMP_INIT;
#if defined(_OPENMP)
@ -54,13 +56,11 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list)
NPAIR_OMP_SETUP(nlocal);
int i,j,jh,k,n,ibin,which,imol,iatom;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
double radi,radsum,cutsq;
int *neighptr;
// loop over each atom, storing neighbors
double **x = atom->x;
double *radius = atom->radius;
int *type = atom->type;
@ -87,6 +87,7 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list)
n = 0;
neighptr = ipage.vget();
itag = tag[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
@ -98,20 +99,31 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list)
}
// loop over all atoms in bins in stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
if (j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}

42
src/OPENMP/npair_half_size_multi_newton_tri_omp.cpp
View File

@ -18,6 +18,7 @@
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neighbor.h"
@ -48,6 +49,7 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list)
const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
const int history = list->history;
const int mask_history = 1 << HISTBITS;
const double delta = 0.01 * force->angstrom;
NPAIR_OMP_INIT;
#if defined(_OPENMP)
@ -55,15 +57,12 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list)
#endif
NPAIR_OMP_SETUP(nlocal);
int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns;
int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns,js;
int which,imol,iatom;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
double radi,radsum,cutdistsq;
int *neighptr,*s;
int js;
// loop over each atom, storing neighbors
int *collection = neighbor->collection;
double **x = atom->x;
@ -92,6 +91,7 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list)
n = 0;
neighptr = ipage.vget();
itag = tag[i];
itype = type[i];
icollection = collection[i];
xtmp = x[i][0];
@ -107,12 +107,13 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list)
ibin = atom2bin[i];
// loop through stencils for all collections
for (jcollection = 0; jcollection < ncollections; jcollection++) {
// if same collection use own bin
if(icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
if (icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
// loop over all atoms in bins in stencil
// stencil is empty if i larger than j
@ -130,14 +131,25 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list)
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
// if same size (same collection), use half stencil
if(cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]){
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
// if same size (same collection), exclude half of interactions
if (cutcollectionsq[icollection][icollection] ==
cutcollectionsq[jcollection][jcollection]) {
if (j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}
}

41
src/OPENMP/npair_half_size_multi_old_newton_tri_omp.cpp
View File

@ -18,6 +18,7 @@
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neigh_list.h"
@ -32,7 +33,6 @@ NPairHalfSizeMultiOldNewtonTriOmp::NPairHalfSizeMultiOldNewtonTriOmp(LAMMPS *lmp
NPair(lmp) {}
/* ----------------------------------------------------------------------
size particles
binned neighbor list construction with Newton's 3rd law for triclinic
each owned atom i checks its own bin and other bins in triclinic stencil
multi-type stencil is itype dependent and is distance checked
@ -46,6 +46,7 @@ void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list)
const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
const int history = list->history;
const int mask_history = 1 << HISTBITS;
const double delta = 0.01 * force->angstrom;
NPAIR_OMP_INIT;
#if defined(_OPENMP)
@ -54,7 +55,7 @@ void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list)
NPAIR_OMP_SETUP(nlocal);
int i,j,jh,k,n,itype,jtype,ibin,ns,which,imol,iatom;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
double radi,radsum,cutdistsq;
int *neighptr,*s;
@ -86,6 +87,7 @@ void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list)
n = 0;
neighptr = ipage.vget();
itag = tag[i];
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
@ -97,13 +99,12 @@ void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list)
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms in bins, including self, in stencil
// skip if i,j neighbor cutoff is less than bin distance
// bins below self are excluded from stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// loop over all atoms in bins in stencil
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
ibin = atom2bin[i];
s = stencil_multi_old[itype];
@ -114,12 +115,22 @@ void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list)
for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
jtype = type[j];
if (cutsq[jtype] < distsq[k]) continue;
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
if (j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}

21
src/OPENMP/npair_half_size_nsq_newton_omp.cpp
View File

@ -18,6 +18,7 @@
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "group.h"
#include "my_page.h"
@ -30,13 +31,11 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
NPairHalfSizeNsqNewtonOmp::NPairHalfSizeNsqNewtonOmp(LAMMPS *lmp) :
NPair(lmp) {}
NPairHalfSizeNsqNewtonOmp::NPairHalfSizeNsqNewtonOmp(LAMMPS *lmp) : NPair(lmp) {}
/* ----------------------------------------------------------------------
size particles
N^2 / 2 search for neighbor pairs with full Newton's 3rd law
shear history must be accounted for when a neighbor pair is added
pair added to list if atoms i and j are both owned and i < j
if j is ghost only me or other proc adds pair
decision based on itag,jtag tests
@ -50,6 +49,8 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list)
const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0;
const int history = list->history;
const int mask_history = 1 << HISTBITS;
const double delta = 0.01 * force->angstrom;
const int triclinic = domain->triclinic;
NPAIR_OMP_INIT;
@ -104,6 +105,12 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list)
}
// loop over remaining atoms, owned and ghost
// use itag/jtap comparision to eliminate half the interactions
// itag = jtag is possible for long cutoffs that include images of self
// for triclinic, must use delta to eliminate half the I/J interactions
// cannot use I/J exact coord comparision as for orthog
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
for (j = i+1; j < nall; j++) {
if (includegroup && !(mask[j] & bitmask)) continue;
@ -114,6 +121,14 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list)
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else if (triclinic) {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {

13
src/OPENMP/npair_halffull_newton_omp.cpp
View File

@ -15,7 +15,9 @@
#include "npair_halffull_newton_omp.h"
#include "atom.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "my_page.h"
#include "neigh_list.h"
#include "npair_omp.h"
@ -38,6 +40,8 @@ NPairHalffullNewtonOmp::NPairHalffullNewtonOmp(LAMMPS *lmp) : NPair(lmp) {}
void NPairHalffullNewtonOmp::build(NeighList *list)
{
const int inum_full = list->listfull->inum;
const double delta = 0.01 * force->angstrom;
const int triclinic = domain->triclinic;
NPAIR_OMP_INIT;
#if defined(_OPENMP)
@ -83,8 +87,17 @@ void NPairHalffullNewtonOmp::build(NeighList *list)
for (jj = 0; jj < jnum; jj++) {
joriginal = jlist[jj];
j = joriginal & NEIGHMASK;
if (j < nlocal) {
if (i > j) continue;
} else if (triclinic) {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {

2
src/balance.cpp
View File

@ -473,7 +473,7 @@ void Balance::options(int iarg, int narg, char **arg, int sortflag_default)
}
iarg += 2+nopt;
} else if (strcmp(arg[iarg+1],"sort") == 0) {
} else if (strcmp(arg[iarg],"sort") == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "balance sort", error);
sortflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;

1
src/compute_property_local.cpp
View File

@ -405,6 +405,7 @@ int ComputePropertyLocal::count_pairs(int allflag, int forceflag)
if (!(mask[j] & groupbit)) continue;
// itag = jtag is possible for long cutoffs that include images of self
// do not need triclinic logic here b/c neighbor list itself is correct
if (newton_pair == 0 && j >= nlocal) {
jtag = tag[j];

30
src/fix_press_langevin.cpp
View File

@ -46,7 +46,8 @@ enum { ISO, ANISO, TRICLINIC };
/* ---------------------------------------------------------------------- */
FixPressLangevin::FixPressLangevin(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), id_press(nullptr), pflag(0), random(nullptr), irregular(nullptr)
Fix(lmp, narg, arg), id_temp(nullptr), id_press(nullptr), temperature(nullptr),
pressure(nullptr), irregular(nullptr), random(nullptr)
{
if (narg < 5) utils::missing_cmd_args(FLERR, "fix press/langevin", error);
@ -62,6 +63,9 @@ FixPressLangevin::FixPressLangevin(LAMMPS *lmp, int narg, char **arg) :
allremap = 1;
pre_exchange_flag = 0;
flipflag = 1;
seed = 111111;
pflag = 0;
kspace_flag = 0;
p_ltime = 0.0;
@ -239,7 +243,7 @@ FixPressLangevin::FixPressLangevin(LAMMPS *lmp, int narg, char **arg) :
t_start = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
t_stop = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
seed = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
if (seed <= 0.0) error->all(FLERR, "Fix press/langevin temp seed must be > 0");
if (seed <= 0) error->all(FLERR, "Fix press/langevin temp seed must be > 0");
iarg += 4;
}
@ -349,7 +353,7 @@ FixPressLangevin::FixPressLangevin(LAMMPS *lmp, int narg, char **arg) :
// Kinetic contribution will be added by the fix style
id_press = utils::strdup(std::string(id) + "_press");
modify->add_compute(fmt::format("{} all pressure NULL virial", id_press));
pressure = modify->add_compute(fmt::format("{} all pressure NULL virial", id_press));
pflag = 1;
// p_fric is alpha coeff from GJF
@ -482,7 +486,7 @@ void FixPressLangevin::initial_integrate(int /* vflag */)
if (delta != 0.0) delta /= update->endstep - update->beginstep;
t_target = t_start + delta * (t_stop - t_start);
couple_beta(t_target);
couple_beta();
dt = update->dt;
@ -492,11 +496,12 @@ void FixPressLangevin::initial_integrate(int /* vflag */)
displacement = dt * p_deriv[i] * gjfb[i];
displacement += 0.5 * dt * dt * f_piston[i] * gjfb[i] / p_mass[i];
displacement += 0.5 * dt * fran[i] * gjfb[i] / p_mass[i];
dl = domain->boxhi[i] - domain->boxlo[i];
if (i < 3)
if (i < 3) {
dl = domain->boxhi[i] - domain->boxlo[i];
dilation[i] = (dl + displacement) / dl;
else
} else {
dilation[i] = displacement;
}
}
}
}
@ -527,7 +532,7 @@ void FixPressLangevin::post_force(int /*vflag*/)
}
couple_pressure();
couple_kinetic(t_target);
couple_kinetic();
for (int i = 0; i < 6; i++) {
if (p_flag[i]) {
@ -594,10 +599,9 @@ void FixPressLangevin::couple_pressure()
}
/* ---------------------------------------------------------------------- */
void FixPressLangevin::couple_kinetic(double t_target)
void FixPressLangevin::couple_kinetic()
{
double pk, volume;
nktv2p = force->nktv2p;
// kinetic part
@ -607,7 +611,7 @@ void FixPressLangevin::couple_kinetic(double t_target)
volume = domain->xprd * domain->yprd;
pk = atom->natoms * force->boltz * t_target / volume;
pk *= nktv2p;
pk *= force->nktv2p;
p_current[0] += pk;
p_current[1] += pk;
@ -616,7 +620,7 @@ void FixPressLangevin::couple_kinetic(double t_target)
/* ---------------------------------------------------------------------- */
void FixPressLangevin::couple_beta(double t_target)
void FixPressLangevin::couple_beta()
{
double gamma[6];
int me = comm->me;
@ -812,7 +816,7 @@ int FixPressLangevin::modify_param(int narg, char **arg)
id_press = utils::strdup(arg[1]);
pressure = modify->get_compute_by_id(arg[1]);
if (pressure) error->all(FLERR, "Could not find fix_modify pressure compute ID: {}", arg[1]);
if (!pressure) error->all(FLERR, "Could not find fix_modify pressure compute ID: {}", arg[1]);
if (pressure->pressflag == 0)
error->all(FLERR, "Fix_modify pressure compute {} does not compute pressure", arg[1]);
return 2;

7
src/fix_press_langevin.h
View File

@ -40,10 +40,9 @@ class FixPressLangevin : public Fix {
int modify_param(int, char **) override;
protected:
int dimension, which;
int dimension;
int pstyle, pcouple, allremap;
int p_flag[6]; // 1 if control P on this dim, 0 if not
double nktv2p;
double t_start, t_stop, t_target;
double p_fric[6], p_ltime; // Friction and Langevin charac. time
double p_alpha[6];
@ -68,8 +67,8 @@ class FixPressLangevin : public Fix {
int seed;
void couple_pressure();
void couple_kinetic(double);
void couple_beta(double);
void couple_kinetic();
void couple_beta();
void remap();
};

2
src/fix_property_atom.cpp
View File

@ -46,6 +46,7 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
rmass_flag = 0;
temperature_flag = 0;
heatflow_flag = 0;
nmax_old = 0;
nvalue = 0;
values_peratom = 0;
@ -212,7 +213,6 @@ void FixPropertyAtom::post_constructor()
{
// perform initial allocation of atom-based array
nmax_old = 0;
grow_arrays(atom->nmax);
}

3
src/min.cpp
View File

@ -215,6 +215,9 @@ void Min::setup(int flag)
}
update->setupflag = 1;
if (lmp->kokkos)
error->all(FLERR,"KOKKOS package requires Kokkos-enabled min_style");
// setup extra global dof due to fixes
// cannot be done in init() b/c update init() is before modify init()

5
src/neighbor.cpp
View File

@ -313,7 +313,10 @@ void Neighbor::init()
triclinic = domain->triclinic;
newton_pair = force->newton_pair;
// error check
// error checks
if (triclinic && atom->tag_enable == 0)
error->all(FLERR, "Cannot build triclinic neighbor lists unless atoms have IDs");
if (delay > 0 && (delay % every) != 0)
error->all(FLERR,"Neighbor delay must be 0 or multiple of every setting");

44
src/npair_half_bin_newton_tri.cpp
View File

@ -13,13 +13,15 @@
------------------------------------------------------------------------- */
#include "npair_half_bin_newton_tri.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "domain.h"
#include "my_page.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neigh_list.h"
using namespace LAMMPS_NS;
@ -36,10 +38,12 @@ NPairHalfBinNewtonTri::NPairHalfBinNewtonTri(LAMMPS *lmp) : NPair(lmp) {}
void NPairHalfBinNewtonTri::build(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,which,imol,iatom,moltemplate;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr;
const double delta = 0.01 * force->angstrom;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
@ -68,6 +72,7 @@ void NPairHalfBinNewtonTri::build(NeighList *list)
n = 0;
neighptr = ipage->vget();
itag = tag[i];
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
@ -79,20 +84,31 @@ void NPairHalfBinNewtonTri::build(NeighList *list)
}
// loop over all atoms in bins in stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
if (j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}

108
src/npair_half_multi_newton_tri.cpp
View File

@ -18,10 +18,11 @@
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neighbor.h"
using namespace LAMMPS_NS;
@ -38,12 +39,14 @@ NPairHalfMultiNewtonTri::NPairHalfMultiNewtonTri(LAMMPS *lmp) : NPair(lmp) {}
void NPairHalfMultiNewtonTri::build(NeighList *list)
{
int i,j,k,n,itype,jtype,icollection,jcollection,ibin,jbin,which,ns,imol,iatom,moltemplate;
tagint tagprev;
int i,j,k,n,itype,jtype,ibin,jbin,icollection,jcollection,which,ns,imol,iatom,moltemplate;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*s;
int js;
const double delta = 0.01 * force->angstrom;
int *collection = neighbor->collection;
double **x = atom->x;
int *type = atom->type;
@ -72,6 +75,8 @@ void NPairHalfMultiNewtonTri::build(NeighList *list)
for (i = 0; i < nlocal; i++) {
n = 0;
neighptr = ipage->vget();
itag = tag[i];
itype = type[i];
icollection = collection[i];
xtmp = x[i][0];
@ -86,36 +91,51 @@ void NPairHalfMultiNewtonTri::build(NeighList *list)
ibin = atom2bin[i];
// loop through stencils for all collections
for (jcollection = 0; jcollection < ncollections; jcollection++) {
// if same collection use own bin
if (icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// loop over all atoms in bins in stencil
// stencil is empty if i larger than j
// stencil is half if i same size as j
// stencil is full if i smaller than j
// if half: pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// for triclinic:
// stencil is empty if i larger than j
// stencil is full if i smaller than j
// stencil is full if i same size as j
// for i smaller than j:
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
// if same size (same collection), use half stencil
if(cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]){
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
// if same size (same collection), exclude half of interactions
if (cutcollectionsq[icollection][icollection] ==
cutcollectionsq[jcollection][jcollection]) {
if (j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}
}
@ -123,28 +143,28 @@ void NPairHalfMultiNewtonTri::build(NeighList *list)
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
}
ilist[inum++] = i;

41
src/npair_half_multi_old_newton_tri.cpp
View File

@ -18,6 +18,7 @@
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neigh_list.h"
@ -38,11 +39,13 @@ NPairHalfMultiOldNewtonTri::NPairHalfMultiOldNewtonTri(LAMMPS *lmp) : NPair(lmp)
void NPairHalfMultiOldNewtonTri::build(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom,moltemplate;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*s;
double *cutsq,*distsq;
const double delta = 0.01 * force->angstrom;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
@ -71,6 +74,7 @@ void NPairHalfMultiOldNewtonTri::build(NeighList *list)
n = 0;
neighptr = ipage->vget();
itag = tag[i];
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
@ -81,13 +85,12 @@ void NPairHalfMultiOldNewtonTri::build(NeighList *list)
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms in bins, including self, in stencil
// skip if i,j neighbor cutoff is less than bin distance
// bins below self are excluded from stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// loop over all atoms in bins in stencil
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
ibin = atom2bin[i];
s = stencil_multi_old[itype];
@ -98,12 +101,22 @@ void NPairHalfMultiOldNewtonTri::build(NeighList *list)
for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
jtype = type[j];
if (cutsq[jtype] < distsq[k]) continue;
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
if (j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}

24
src/npair_half_nsq_newton.cpp
View File

@ -13,14 +13,16 @@
------------------------------------------------------------------------- */
#include "npair_half_nsq_newton.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "molecule.h"
#include "domain.h"
#include "my_page.h"
#include "error.h"
#include "neigh_list.h"
using namespace LAMMPS_NS;
@ -41,6 +43,9 @@ void NPairHalfNsqNewton::build(NeighList *list)
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr;
const double delta = 0.01 * force->angstrom;
const int triclinic = domain->triclinic;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
@ -85,7 +90,12 @@ void NPairHalfNsqNewton::build(NeighList *list)
}
// loop over remaining atoms, owned and ghost
// use itag/jtap comparision to eliminate half the interactions
// itag = jtag is possible for long cutoffs that include images of self
// for triclinic, must use delta to eliminate half the I/J interactions
// cannot use I/J exact coord comparision as for orthog
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
for (j = i+1; j < nall; j++) {
if (includegroup && !(mask[j] & bitmask)) continue;
@ -96,6 +106,14 @@ void NPairHalfNsqNewton::build(NeighList *list)
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else if (triclinic) {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {

44
src/npair_half_respa_bin_newton_tri.cpp
View File

@ -13,13 +13,15 @@
------------------------------------------------------------------------- */
#include "npair_half_respa_bin_newton_tri.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "domain.h"
#include "my_page.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neigh_list.h"
using namespace LAMMPS_NS;
@ -38,10 +40,12 @@ NPairHalfRespaBinNewtonTri::NPairHalfRespaBinNewtonTri(LAMMPS *lmp) :
void NPairHalfRespaBinNewtonTri::build(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,n_inner,n_middle,imol,iatom,moltemplate;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*neighptr_inner,*neighptr_middle;
const double delta = 0.01 * force->angstrom;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
@ -94,6 +98,7 @@ void NPairHalfRespaBinNewtonTri::build(NeighList *list)
neighptr_middle = ipage_middle->vget();
}
itag = tag[i];
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
@ -105,20 +110,31 @@ void NPairHalfRespaBinNewtonTri::build(NeighList *list)
}
// loop over all atoms in bins in stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
if (j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}

25
src/npair_half_respa_nsq_newton.cpp
View File

@ -13,14 +13,16 @@
------------------------------------------------------------------------- */
#include "npair_half_respa_nsq_newton.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "molecule.h"
#include "domain.h"
#include "my_page.h"
#include "error.h"
#include "neigh_list.h"
using namespace LAMMPS_NS;
@ -44,6 +46,9 @@ void NPairHalfRespaNsqNewton::build(NeighList *list)
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*neighptr_inner,*neighptr_middle;
const double delta = 0.01 * force->angstrom;
const int triclinic = domain->triclinic;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
@ -112,6 +117,12 @@ void NPairHalfRespaNsqNewton::build(NeighList *list)
}
// loop over remaining atoms, owned and ghost
// use itag/jtap comparision to eliminate half the interactions
// itag = jtag is possible for long cutoffs that include images of self
// for triclinic, must use delta to eliminate half the I/J interactions
// cannot use I/J exact coord comparision as for orthog
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
for (j = i+1; j < nall; j++) {
if (includegroup && !(mask[j] & bitmask)) continue;
@ -122,6 +133,14 @@ void NPairHalfRespaNsqNewton::build(NeighList *list)
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else if (triclinic) {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {

2
src/npair_half_respa_nsq_newton.h
View File

@ -15,7 +15,7 @@
// clang-format off
NPairStyle(half/respa/nsq/newton,
NPairHalfRespaNsqNewton,
NP_HALF | NP_RESPA | NP_NSQ | NP_NEWTON | NP_ORTHO);
NP_HALF | NP_RESPA | NP_NSQ | NP_NEWTON | NP_ORTHO | NP_TRI);
// clang-format on
#else

37
src/npair_half_size_bin_newton_tri.cpp
View File

@ -18,6 +18,7 @@
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neigh_list.h"
@ -39,11 +40,13 @@ NPairHalfSizeBinNewtonTri::NPairHalfSizeBinNewtonTri(LAMMPS *lmp) :
void NPairHalfSizeBinNewtonTri::build(NeighList *list)
{
int i,j,jh,k,n,ibin,which,imol,iatom,moltemplate;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
double radi,radsum,cutsq;
int *neighptr;
const double delta = 0.01 * force->angstrom;
double **x = atom->x;
double *radius = atom->radius;
int *type = atom->type;
@ -76,6 +79,7 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list)
n = 0;
neighptr = ipage->vget();
itag = tag[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
@ -87,20 +91,31 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list)
}
// loop over all atoms in bins in stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
if (j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}

98
src/npair_half_size_multi_newton_tri.cpp
View File

@ -18,6 +18,7 @@
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neighbor.h"
@ -41,11 +42,13 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
{
int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns,js;
int which,imol,iatom,moltemplate;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
double radi,radsum,cutdistsq;
int *neighptr,*s;
const double delta = 0.01 * force->angstrom;
int *collection = neighbor->collection;
double **x = atom->x;
double *radius = atom->radius;
@ -78,6 +81,8 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
for (i = 0; i < nlocal; i++) {
n = 0;
neighptr = ipage->vget();
itag = tag[i];
itype = type[i];
icollection = collection[i];
xtmp = x[i][0];
@ -93,11 +98,13 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
ibin = atom2bin[i];
// loop through stencils for all collections
for (jcollection = 0; jcollection < ncollections; jcollection++) {
// if same collection use own bin
if (icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// loop over all atoms in bins in stencil
// stencil is empty if i larger than j
@ -108,21 +115,32 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
// if same size (same collection), use half stencil
if (cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]){
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
// if same size (same collection), exclude half of interactions
if (cutcollectionsq[icollection][icollection] ==
cutcollectionsq[jcollection][jcollection]) {
if (j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}
}
@ -130,34 +148,34 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
if (rsq <= cutdistsq) {
jh = j;
if (history && rsq < radsum*radsum)
jh = jh ^ mask_history;
if (rsq <= cutdistsq) {
jh = j;
if (history && rsq < radsum*radsum)
jh = jh ^ mask_history;
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = jh;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = jh;
else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
} else neighptr[n++] = jh;
}
}
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = jh;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = jh;
else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
} else neighptr[n++] = jh;
}
}
}
}
ilist[inum++] = i;

41
src/npair_half_size_multi_old_newton_tri.cpp
View File

@ -18,6 +18,7 @@
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neigh_list.h"
@ -38,12 +39,14 @@ NPairHalfSizeMultiOldNewtonTri::NPairHalfSizeMultiOldNewtonTri(LAMMPS *lmp) : NP
void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list)
{
int i,j,jh,k,n,itype,jtype,ibin,ns,which,imol,iatom,moltemplate;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
double radi,radsum,cutdistsq;
int *neighptr,*s;
double *cutsq,*distsq;
const double delta = 0.01 * force->angstrom;
double **x = atom->x;
double *radius = atom->radius;
int *type = atom->type;
@ -76,6 +79,7 @@ void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list)
n = 0;
neighptr = ipage->vget();
itag = tag[i];
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
@ -87,13 +91,12 @@ void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list)
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms in bins, including self, in stencil
// skip if i,j neighbor cutoff is less than bin distance
// bins below self are excluded from stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// loop over all atoms in bins in stencil
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
ibin = atom2bin[i];
s = stencil_multi_old[itype];
@ -104,12 +107,22 @@ void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list)
for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
jtype = type[j];
if (cutsq[jtype] < distsq[k]) continue;
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
if (j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}

18
src/npair_half_size_nsq_newton.cpp
View File

@ -18,6 +18,7 @@
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "group.h"
#include "my_page.h"
@ -45,6 +46,9 @@ void NPairHalfSizeNsqNewton::build(NeighList *list)
double radi,radsum,cutsq;
int *neighptr;
const double delta = 0.01 * force->angstrom;
const int triclinic = domain->triclinic;
double **x = atom->x;
double *radius = atom->radius;
tagint *tag = atom->tag;
@ -93,6 +97,12 @@ void NPairHalfSizeNsqNewton::build(NeighList *list)
}
// loop over remaining atoms, owned and ghost
// use itag/jtap comparision to eliminate half the interactions
// itag = jtag is possible for long cutoffs that include images of self
// for triclinic, must use delta to eliminate half the I/J interactions
// cannot use I/J exact coord comparision as for orthog
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
for (j = i+1; j < nall; j++) {
if (includegroup && !(mask[j] & bitmask)) continue;
@ -103,6 +113,14 @@ void NPairHalfSizeNsqNewton::build(NeighList *list)
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else if (triclinic) {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {

20
src/npair_halffull_newton.cpp
View File

@ -14,7 +14,9 @@
#include "npair_halffull_newton.h"
#include "atom.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "my_page.h"
#include "neigh_list.h"
@ -37,6 +39,9 @@ void NPairHalffullNewton::build(NeighList *list)
int *neighptr, *jlist;
double xtmp, ytmp, ztmp;
const double delta = 0.01 * force->angstrom;
const int triclinic = domain->triclinic;
double **x = atom->x;
int nlocal = atom->nlocal;
@ -65,6 +70,11 @@ void NPairHalffullNewton::build(NeighList *list)
ztmp = x[i][2];
// loop over full neighbor list
// use i < j < nlocal to eliminate half the local/local interactions
// for triclinic, must use delta to eliminate half the local/ghost interactions
// cannot use I/J exact coord comparision as for orthog
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
jlist = firstneigh_full[i];
jnum = numneigh_full[i];
@ -72,8 +82,17 @@ void NPairHalffullNewton::build(NeighList *list)
for (jj = 0; jj < jnum; jj++) {
joriginal = jlist[jj];
j = joriginal & NEIGHMASK;
if (j < nlocal) {
if (i > j) continue;
} else if (triclinic) {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
@ -81,6 +100,7 @@ void NPairHalffullNewton::build(NeighList *list)
if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
}
neighptr[n++] = joriginal;
}

19
src/npair_halffull_newton_trim.cpp
View File

@ -14,7 +14,9 @@
#include "npair_halffull_newton_trim.h"
#include "atom.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "my_page.h"
#include "neigh_list.h"
@ -38,6 +40,9 @@ void NPairHalffullNewtonTrim::build(NeighList *list)
double xtmp, ytmp, ztmp;
double delx, dely, delz, rsq;
const double delta = 0.01 * force->angstrom;
const int triclinic = domain->triclinic;
double **x = atom->x;
int nlocal = atom->nlocal;
@ -68,6 +73,11 @@ void NPairHalffullNewtonTrim::build(NeighList *list)
ztmp = x[i][2];
// loop over full neighbor list
// use i < j < nlocal to eliminate half the local/local interactions
// for triclinic, must use delta to eliminate half the local/ghost interactions
// cannot use I/J exact coord comparision as for orthog
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
jlist = firstneigh_full[i];
jnum = numneigh_full[i];
@ -75,8 +85,17 @@ void NPairHalffullNewtonTrim::build(NeighList *list)
for (jj = 0; jj < jnum; jj++) {
joriginal = jlist[jj];
j = joriginal & NEIGHMASK;
if (j < nlocal) {
if (i > j) continue;
} else if (triclinic) {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {

12
src/nstencil_half_bin_2d_tri.cpp
View File

@ -27,9 +27,17 @@ void NStencilHalfBin2dTri::create()
{
int i, j;
// for triclinic, need to use full stencil in all dims
// not a half stencil in y
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift both coords by epsilon
// thus for an I/J owned/ghost pair, the xy coords
// and bin assignments can be different on I proc vs J proc
nstencil = 0;
for (j = 0; j <= sy; j++)
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++)
if (bin_distance(i, j, 0) < cutneighmaxsq) stencil[nstencil++] = j * mbinx + i;
if (bin_distance(i, j, 0) < cutneighmaxsq)
stencil[nstencil++] = j * mbinx + i;
}

9
src/nstencil_half_bin_3d_tri.cpp
View File

@ -27,9 +27,16 @@ void NStencilHalfBin3dTri::create()
{
int i, j, k;
// for triclinic, need to use full stencil in all dims
// not a half stencil in z
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
// thus for an I/J owned/ghost pair, the xyz coords
// and bin assignments can be different on I proc vs J proc
nstencil = 0;
for (k = 0; k <= sz; k++)
for (k = -sz; k <= sz; k++)
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++)
if (bin_distance(i, j, k) < cutneighmaxsq)

2
src/nstencil_half_multi_2d_tri.cpp
View File

@ -80,7 +80,7 @@ void NStencilHalfMulti2dTri::create()
cutsq = cutcollectionsq[icollection][jcollection];
if (flag_half_multi[icollection][jcollection]) {
for (j = 0; j <= sy; j++)
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++)
if (bin_distance_multi(i, j, 0, bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] = j * mbinx + i;

2
src/nstencil_half_multi_3d_tri.cpp
View File

@ -81,7 +81,7 @@ void NStencilHalfMulti3dTri::create()
cutsq = cutcollectionsq[icollection][jcollection];
if (flag_half_multi[icollection][jcollection]) {
for (k = 0; k <= sz; k++)
for (k = -sz; k <= sz; k++)
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++)
if (bin_distance_multi(i, j, k, bin_collection) < cutsq)

2
src/nstencil_half_multi_old_2d_tri.cpp
View File

@ -37,7 +37,7 @@ void NStencilHalfMultiOld2dTri::create()
s = stencil_multi_old[itype];
distsq = distsq_multi_old[itype];
n = 0;
for (j = 0; j <= sy; j++)
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++) {
rsq = bin_distance(i, j, 0);
if (rsq < typesq) {

2
src/nstencil_half_multi_old_3d_tri.cpp
View File

@ -37,7 +37,7 @@ void NStencilHalfMultiOld3dTri::create()
s = stencil_multi_old[itype];
distsq = distsq_multi_old[itype];
n = 0;
for (k = 0; k <= sz; k++)
for (k = -sz; k <= sz; k++)
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++) {
rsq = bin_distance(i, j, k);

9
tools/lammps-gui/CMakeLists.txt
View File

@ -110,14 +110,15 @@ endif()
# we require Qt 5 and at least version 5.12 at that.
if(NOT LAMMPS_GUI_USE_QT5)
find_package(Qt6 6.2 COMPONENTS Widgets Charts)
find_package(Qt6 6.2 QUIET COMPONENTS Widgets Charts)
endif()
if(NOT Qt6_FOUND)
find_package(Qt5 5.12 REQUIRED COMPONENTS Widgets Charts)
set(QT_VERSION_MAJOR "5")
set(QT_VERSION_MAJOR 5)
else()
set(QT_VERSION_MAJOR "6")
set(QT_VERSION_MAJOR 6)
endif()
message(STATUS "Using Qt version ${Qt${QT_VERSION_MAJOR}_VERSION} for LAMMPS GUI")
set(PROJECT_SOURCES
main.cpp
@ -188,7 +189,7 @@ else()
endif()
target_include_directories(lammps-gui PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
target_compile_definitions(lammps-gui PRIVATE LAMMPS_GUI_VERSION="${PROJECT_VERSION}")
target_link_libraries(lammps-gui PRIVATE Qt${QT_VERSION_MAJOR}::Widgets Qt${VERSION_MAJOR}::Charts)
target_link_libraries(lammps-gui PRIVATE Qt${QT_VERSION_MAJOR}::Widgets Qt${QT_VERSION_MAJOR}::Charts)
if(BUILD_OMP)
find_package(OpenMP COMPONENTS CXX REQUIRED)
target_link_libraries(lammps-gui PRIVATE OpenMP::OpenMP_CXX)

6
tools/lammps-gui/lammpsgui.cpp
View File

@ -1432,12 +1432,12 @@ void LammpsGui::start_lammps()
lammps.open(narg, args);
lammpsstatus->show();
// must have at least 2 August 2023 version of LAMMPS
// must have a version newer than the 2 August 2023 release of LAMMPS
// TODO: must update this check before next feature release
if (lammps.version() < 20230802) {
if (lammps.version() <= 20230802) {
QMessageBox::critical(this, "Incompatible LAMMPS Version",
"LAMMPS-GUI version " LAMMPS_GUI_VERSION " requires\n"
"LAMMPS version 2 August 2023 or later");
"a LAMMPS version more recent than 2 August 2023");
exit(1);
}

2
unittest/formats/CMakeLists.txt
View File

@ -41,6 +41,8 @@ set_tests_properties(TextFileReader PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${
add_executable(test_file_operations test_file_operations.cpp)
target_link_libraries(test_file_operations PRIVATE lammps GTest::GMock)
add_test(NAME FileOperations COMMAND test_file_operations)
# try to mitigate possible OpenMPI bug
set_tests_properties(TextFileReader PROPERTIES ENVIRONMENT "OMPI_MCA_sharedfp=\"^sm\"")
add_executable(test_dump_atom test_dump_atom.cpp)
target_link_libraries(test_dump_atom PRIVATE lammps GTest::GMock)