Merge remote-tracking branch 'github/develop' into collected-small-fixes

This commit is contained in:
Axel Kohlmeyer
2023-10-26 20:35:36 -04:00
14 changed files with 304 additions and 38 deletions

View File

@ -24,6 +24,7 @@ Syntax
c1,c2 = coords of cone axis in other 2 dimensions (distance units)
radlo,radhi = cone radii at lo and hi end (distance units)
lo,hi = bounds of cone in dim (distance units)
c1,c2,radlo,radhi,lo,hi can be a variable (see below)
*cylinder* args = dim c1 c2 radius lo hi
dim = *x* or *y* or *z* = axis of cylinder
c1,c2 = coords of cylinder axis in other 2 dimensions (distance units)
@ -206,7 +207,8 @@ equal-style :doc:`variable <variable>`. Likewise, for style *sphere*
and *ellipsoid* the x-, y-, and z- coordinates of the center of the
sphere/ellipsoid can be specified as an equal-style variable. And for
style *cylinder* the two center positions c1 and c2 for the location
of the cylinder axes can be specified as a equal-style variable.
of the cylinder axes can be specified as a equal-style variable. For style *cone*
all properties can be defined via equal-style variables.
If the value is a variable, it should be specified as v_name, where
name is the variable name. In this case, the variable will be

View File

@ -85,6 +85,7 @@ void MLIAPDataKokkos<DeviceType>::generate_neighdata(class NeighList *list_in, i
// clear gradforce and elems arrays
int nall = atom->nlocal + atom->nghost;
nlocal = atom->nlocal;
ntotal = nall;
if (gradgradflag > -1){
auto d_gradforce = k_gradforce.template view<DeviceType>();

View File

@ -118,6 +118,7 @@ public:
egradient(nullptr),
ntotal(base.ntotal),
nlistatoms(base.nlistatoms),
nlocal(base.nlocal),
natomneigh(base.natomneigh),
numneighs(base.numneighs),
iatoms(base.k_iatoms.d_view.data()),
@ -171,6 +172,7 @@ public:
// Neighborlist stuff
const int ntotal;
const int nlistatoms;
const int nlocal;
const int natomneigh;
int *numneighs;
int *iatoms;
@ -191,7 +193,7 @@ public:
int dev;
#ifdef LMP_KOKKOS_GPU
MLIAPDataKokkosDevice(MLIAPDataKokkos<LMPHostType> &base) : ndescriptors(-1),nparams(-1),nelements(-1),ntotal(-1),nlistatoms(-1),natomneigh(-1),
MLIAPDataKokkosDevice(MLIAPDataKokkos<LMPHostType> &base) : ndescriptors(-1),nparams(-1),nelements(-1),ntotal(-1),nlistatoms(-1),nlocal(-1),natomneigh(-1),
nneigh_max(-1),npairs(-1)
{
// It cannot get here, but needed for compilation

View File

@ -25,6 +25,7 @@ cdef extern from "mliap_data_kokkos.h" namespace "LAMMPS_NS":
cdef cppclass MLIAPDataKokkosDevice:
# Array shapes
int nlistatoms
int nlocal
int ndescriptors
# Input data

View File

@ -59,6 +59,7 @@ cdef extern from "mliap_data_kokkos.h" namespace "LAMMPS_NS":
int ntotal # total number of owned and ghost atoms on this proc
int nlistatoms # current number of atoms in local atom lists
int nlocal
int natomneigh # current number of atoms and ghosts in atom neighbor arrays
int * numneighs # neighbors count for each atom
int * iatoms # index of each atom
@ -290,6 +291,10 @@ cdef class MLIAPDataPy:
@property
def nlistatoms(self):
return self.data.nlistatoms
@property
def nlocal(self):
return self.data.nlocal
@property
def natomneigh(self):

View File

@ -271,8 +271,8 @@ void LAMMPS_NS::update_pair_energy(MLIAPDataKokkosDevice *data, double *eij)
{
auto d_eatoms = data->eatoms;
auto d_pair_i= data->pair_i;
const auto nlistatoms = data->nlistatoms;
Kokkos::parallel_for(nlistatoms, KOKKOS_LAMBDA(int ii){
const auto nlocal = data->nlocal;
Kokkos::parallel_for(nlocal, KOKKOS_LAMBDA(int ii){
d_eatoms[ii] = 0;
});
@ -281,7 +281,7 @@ void LAMMPS_NS::update_pair_energy(MLIAPDataKokkosDevice *data, double *eij)
double e = 0.5 * eij[ii];
// must not count any contribution where i is not a local atom
if (i < nlistatoms) {
if (i < nlocal) {
Kokkos::atomic_add(&d_eatoms[i], e);
local_sum += e;
}
@ -294,7 +294,7 @@ void LAMMPS_NS::update_pair_energy(MLIAPDataKokkosDevice *data, double *eij)
void LAMMPS_NS::update_pair_forces(MLIAPDataKokkosDevice *data, double *fij)
{
const auto nlistatoms = data->nlistatoms;
const auto nlocal = data->nlocal;
auto *f = data->f;
auto pair_i = data->pair_i;
auto j_atoms = data->jatoms;
@ -315,7 +315,7 @@ void LAMMPS_NS::update_pair_forces(MLIAPDataKokkosDevice *data, double *fij)
int i = pair_i[ii];
int j = j_atoms[ii];
// must not count any contribution where i is not a local atom
if (i < nlistatoms) {
if (i < nlocal) {
Kokkos::atomic_add(&f[i*3+0], fij[ii3+0]);
Kokkos::atomic_add(&f[i*3+1], fij[ii3+1]);
Kokkos::atomic_add(&f[i*3+2], fij[ii3+2]);
@ -378,12 +378,12 @@ void LAMMPS_NS::update_pair_forces(MLIAPDataKokkosDevice *data, double *fij)
void LAMMPS_NS::update_atom_energy(MLIAPDataKokkosDevice *data, double *ei)
{
auto d_eatoms = data->eatoms;
const auto nlistatoms = data->nlistatoms;
const auto nlocal = data->nlocal;
Kokkos::parallel_reduce(nlistatoms, KOKKOS_LAMBDA(int i, double &local_sum){
Kokkos::parallel_reduce(nlocal, KOKKOS_LAMBDA(int i, double &local_sum){
double e = ei[i];
// must not count any contribution where i is not a local atom
if (i < nlistatoms) {
if (i < nlocal) {
d_eatoms[i] = e;
local_sum += e;
}

View File

@ -138,6 +138,7 @@ template<class DeviceType>
void PairMLIAPKokkos<DeviceType>::allocate()
{
int n = atom->ntypes;
memoryKK->destroy_kokkos(k_map, map);
memoryKK->destroy_kokkos(k_cutsq, cutsq);
memoryKK->destroy_kokkos(k_setflag, setflag);
@ -275,7 +276,10 @@ void PairMLIAPKokkos<DeviceType>::coeff(int narg, char **arg) {
auto h_cutsq=k_cutsq.template view<LMPHostType>();
for (int itype=1; itype <= atom->ntypes; ++itype)
for (int jtype=1; jtype <= atom->ntypes; ++jtype)
h_cutsq(itype,jtype) = descriptor->cutsq[map[itype]][map[jtype]];
// do not set cuts for NULL atoms
if (map[itype] >= 0 && map[jtype] >= 0) {
h_cutsq(itype,jtype) = descriptor->cutsq[map[itype]][map[jtype]];
}
k_cutsq.modify<LMPHostType>();
k_cutsq.sync<DeviceType>();
constexpr int gradgradflag = -1;

View File

@ -115,6 +115,7 @@ void MLIAPData::generate_neighdata(NeighList *list_in, int eflag_in, int vflag_i
int **firstneigh = list->firstneigh;
int nall = atom->nlocal + atom->nghost;
nlocal = atom->nlocal;
ntotal = nall;
// grow nmax gradforce, elems arrays if necessary

View File

@ -60,6 +60,7 @@ class MLIAPData : protected Pointers {
int ntotal; // total number of owned and ghost atoms on this proc
int nlistatoms; // current number of atoms in local atom lists
int nlocal;
int nlistatoms_max; // allocated size of descriptor array
int natomneigh; // current number of atoms and ghosts in atom neighbor arrays
int natomneigh_max; // allocated size of atom neighbor arrays

View File

@ -18,6 +18,7 @@ cdef extern from "mliap_data.h" namespace "LAMMPS_NS":
cdef cppclass MLIAPData:
# Array shapes
int nlistatoms
int nlocal
int ndescriptors
# Input data

View File

@ -246,6 +246,7 @@ void LAMMPS_NS::update_pair_energy(MLIAPData *data, double *eij)
{
double e_total = 0.0;
const auto nlistatoms = data->nlistatoms;
const auto nlocal = data->nlocal;
for (int ii = 0; ii < nlistatoms; ii++) data->eatoms[ii] = 0;
for (int ii = 0; ii < data->npairs; ii++) {
@ -253,7 +254,7 @@ void LAMMPS_NS::update_pair_energy(MLIAPData *data, double *eij)
double e = 0.5 * eij[ii];
// must not count any contribution where i is not a local atom
if (i < nlistatoms) {
if (i < nlocal) {
data->eatoms[i] += e;
e_total += e;
}
@ -267,7 +268,9 @@ void LAMMPS_NS::update_pair_energy(MLIAPData *data, double *eij)
void LAMMPS_NS::update_pair_forces(MLIAPData *data, double *fij)
{
const auto nlistatoms = data->nlistatoms;
//Bugfix: need to account for Null atoms in local atoms
//const auto nlistatoms = data->nlistatoms;
const auto nlocal = data->nlocal;
double **f = data->f;
for (int ii = 0; ii < data->npairs; ii++) {
int ii3 = ii * 3;
@ -275,7 +278,7 @@ void LAMMPS_NS::update_pair_forces(MLIAPData *data, double *fij)
int j = data->jatoms[ii];
// must not count any contribution where i is not a local atom
if (i < nlistatoms) {
if (i < nlocal) {
f[i][0] += fij[ii3];
f[i][1] += fij[ii3 + 1];
f[i][2] += fij[ii3 + 2];

View File

@ -53,7 +53,8 @@ cdef extern from "mliap_data.h" namespace "LAMMPS_NS":
# only neighbors strictly inside descriptor cutoff
int ntotal # total number of owned and ghost atoms on this proc
int nlistatoms # current number of atoms in local atom lists
int nlistatoms # current number of non-NULL atoms in local atom lists
int nlocal # current number of NULL and normal atoms in local atom lists
int natomneigh # current number of atoms and ghosts in atom neighbor arrays
int * numneighs # neighbors count for each atom
int * iatoms # index of each atom
@ -226,6 +227,10 @@ cdef class MLIAPDataPy:
@property
def nlistatoms(self):
return self.data.nlistatoms
@property
def nlocal(self):
return self.data.nlocal
@property
def natomneigh(self):

View File

@ -19,17 +19,22 @@
#include "domain.h"
#include "error.h"
#include "input.h"
#include "variable.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
enum { CONSTANT, VARIABLE };
static constexpr double BIG = 1.0e20;
/* ---------------------------------------------------------------------- */
RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg), lo(0.0), hi(0.0)
RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg), lo(0.0), hi(0.0),
c1str(nullptr), c2str(nullptr), rlostr(nullptr), rhistr(nullptr), lostr(nullptr), histr(nullptr)
{
options(narg - 9, &arg[9]);
@ -44,22 +49,133 @@ RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg), lo
axis = arg[2][0];
if (axis == 'x') {
c1 = yscale * utils::numeric(FLERR, arg[3], false, lmp);
c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp);
radiuslo = yscale * utils::numeric(FLERR, arg[5], false, lmp);
radiushi = yscale * utils::numeric(FLERR, arg[6], false, lmp);
if (utils::strmatch(arg[3], "^v_")) {
c1str = utils::strdup(arg[3] + 2);
c1 = 0.0;
c1style = VARIABLE;
varshape = 1;
} else {
c1 = yscale * utils::numeric(FLERR, arg[3], false, lmp);
c1style = CONSTANT;
}
if (utils::strmatch(arg[4], "^v_")) {
c2str = utils::strdup(arg[4] + 2);
c2 = 0.0;
c2style = VARIABLE;
varshape = 1;
} else {
c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp);
c2style = CONSTANT;
}
if (utils::strmatch(arg[5], "^v_")) {
rlostr = utils::strdup(arg[5] + 2);
radiuslo = 0.0;
rlostyle = VARIABLE;
varshape = 1;
} else {
radiuslo = yscale * utils::numeric(FLERR, arg[5], false, lmp);
rlostyle = CONSTANT;
}
if (utils::strmatch(arg[6], "^v_")) {
rhistr = utils::strdup(arg[6] + 2);
radiushi = 0.0;
rhistyle = VARIABLE;
varshape = 1;
} else {
radiushi = yscale * utils::numeric(FLERR, arg[6], false, lmp);
rhistyle = CONSTANT;
}
} else if (axis == 'y') {
c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp);
c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp);
radiuslo = xscale * utils::numeric(FLERR, arg[5], false, lmp);
radiushi = xscale * utils::numeric(FLERR, arg[6], false, lmp);
if (utils::strmatch(arg[3], "^v_")) {
c1str = utils::strdup(arg[3] + 2);
c1 = 0.0;
c1style = VARIABLE;
varshape = 1;
} else {
c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp);
c1style = CONSTANT;
}
if (utils::strmatch(arg[4], "^v_")) {
c2str = utils::strdup(arg[4] + 2);
c2 = 0.0;
c2style = VARIABLE;
varshape = 1;
} else {
c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp);
c2style = CONSTANT;
}
if (utils::strmatch(arg[5], "^v_")) {
rlostr = utils::strdup(arg[5] + 2);
radiuslo = 0.0;
rlostyle = VARIABLE;
varshape = 1;
} else {
radiuslo = xscale * utils::numeric(FLERR, arg[5], false, lmp);
rlostyle = CONSTANT;
}
if (utils::strmatch(arg[6], "^v_")) {
rhistr = utils::strdup(arg[6] + 2);
radiushi = 0.0;
rhistyle = VARIABLE;
varshape = 1;
} else {
radiushi = xscale * utils::numeric(FLERR, arg[6], false, lmp);
rhistyle = CONSTANT;
}
} else if (axis == 'z') {
c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp);
c2 = yscale * utils::numeric(FLERR, arg[4], false, lmp);
radiuslo = xscale * utils::numeric(FLERR, arg[5], false, lmp);
radiushi = xscale * utils::numeric(FLERR, arg[6], false, lmp);
if (utils::strmatch(arg[3], "^v_")) {
c1str = utils::strdup(arg[3] + 2);
c1 = 0.0;
c1style = VARIABLE;
varshape = 1;
} else {
c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp);
c1style = CONSTANT;
}
if (utils::strmatch(arg[4], "^v_")) {
c2str = utils::strdup(arg[4] + 2);
c2 = 0.0;
c2style = VARIABLE;
varshape = 1;
} else {
c2 = yscale * utils::numeric(FLERR, arg[4], false, lmp);
c2style = CONSTANT;
}
if (utils::strmatch(arg[5], "^v_")) {
rlostr = utils::strdup(arg[5] + 2);
radiuslo = 0.0;
rlostyle = VARIABLE;
varshape = 1;
} else {
radiuslo = xscale * utils::numeric(FLERR, arg[5], false, lmp);
rlostyle = CONSTANT;
}
if (utils::strmatch(arg[6], "^v_")) {
rhistr = utils::strdup(arg[6] + 2);
radiushi = 0.0;
rhistyle = VARIABLE;
varshape = 1;
} else {
radiushi = xscale * utils::numeric(FLERR, arg[6], false, lmp);
rhistyle = CONSTANT;
}
}
lostyle = CONSTANT;
if (strcmp(arg[7], "INF") == 0 || strcmp(arg[7], "EDGE") == 0) {
if (domain->box_exist == 0)
error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
@ -87,6 +203,11 @@ RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg), lo
else
lo = domain->boxlo_bound[2];
}
} else if (utils::strmatch(arg[7], "^v_")) {
lostr = utils::strdup(arg[7] + 2);
lo = 0.0;
lostyle = VARIABLE;
varshape = 1;
} else {
if (axis == 'x') lo = xscale * utils::numeric(FLERR, arg[7], false, lmp);
if (axis == 'y') lo = yscale * utils::numeric(FLERR, arg[7], false, lmp);
@ -119,26 +240,32 @@ RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg), lo
else
hi = domain->boxhi_bound[2];
}
} else if (utils::strmatch(arg[8], "^v_")) {
histr = utils::strdup(arg[8] + 2);
hi = 0.0;
histyle = VARIABLE;
varshape = 1;
} else {
if (axis == 'x') hi = xscale * utils::numeric(FLERR, arg[8], false, lmp);
if (axis == 'y') hi = yscale * utils::numeric(FLERR, arg[8], false, lmp);
if (axis == 'z') hi = zscale * utils::numeric(FLERR, arg[8], false, lmp);
}
if (varshape) {
variable_check();
RegCone::shape_update();
}
// error check
if (radiuslo < 0.0) error->all(FLERR, "Illegal radius in region cone command");
if (radiushi < 0.0) error->all(FLERR, "Illegal radius in region cone command");
if (radiuslo == 0.0 && radiushi == 0.0)
error->all(FLERR, "Illegal radius in region cone command");
if (hi == lo) error->all(FLERR, "Illegal cone length in region cone command");
if (hi <= lo) error->all(FLERR, "Illegal cone length in region cone command");
// extent of cone
if (radiuslo > radiushi)
maxradius = radiuslo;
else
maxradius = radiushi;
maxradius = ( (radiuslo > radiushi) ? radiuslo : radiushi);
if (interior) {
bboxflag = 1;
@ -174,19 +301,30 @@ RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg), lo
cmax = 3;
contact = new Contact[cmax];
if (interior)
tmax = 2;
else
tmax = 1;
tmax = (interior ? 2 : 1);
}
/* ---------------------------------------------------------------------- */
RegCone::~RegCone()
{
delete[] c1str;
delete[] c2str;
delete[] rlostr;
delete[] rhistr;
delete[] lostr;
delete[] histr;
delete[] contact;
}
/* ---------------------------------------------------------------------- */
void RegCone::init()
{
Region::init();
if (varshape) variable_check();
}
/* ----------------------------------------------------------------------
inside = 1 if x,y,z is inside or on surface
inside = 0 if x,y,z is outside and not on surface
@ -622,6 +760,101 @@ int RegCone::surface_exterior(double *x, double cutoff)
}
}
/* ----------------------------------------------------------------------
change region shape via variable evaluation
------------------------------------------------------------------------- */
void RegCone::shape_update()
{
if (c1style == VARIABLE) c1 = input->variable->compute_equal(c1var);
if (c2style == VARIABLE) c2 = input->variable->compute_equal(c2var);
if (rlostyle == VARIABLE) {
radiuslo = input->variable->compute_equal(rlovar);
if (radiuslo < 0.0) error->one(FLERR, "Variable evaluation in region gave bad value");
}
if (rhistyle == VARIABLE) {
radiushi = input->variable->compute_equal(rhivar);
if (radiushi < 0.0) error->one(FLERR, "Variable evaluation in region gave bad value");
}
if (lostyle == VARIABLE) lo = input->variable->compute_equal(lovar);
if (histyle == VARIABLE) hi = input->variable->compute_equal(hivar);
if (radiuslo == 0.0 && radiushi == 0.0)
error->all(FLERR, "dtion in region gave bad value");
if (axis == 'x') {
if (c1style == VARIABLE) c1 *= yscale;
if (c2style == VARIABLE) c2 *= zscale;
if (rlostyle == VARIABLE) radiuslo *= yscale;
if (rhistyle == VARIABLE) radiushi *= yscale;
if (lostyle == VARIABLE) lo *= xscale;
if (histyle == VARIABLE) hi *= xscale;
} else if (axis == 'y') {
if (c1style == VARIABLE) c1 *= xscale;
if (c2style == VARIABLE) c2 *= zscale;
if (rlostyle == VARIABLE) radiuslo *= xscale;
if (rhistyle == VARIABLE) radiushi *= xscale;
if (lostyle == VARIABLE) lo *= yscale;
if (histyle == VARIABLE) hi *= yscale;
} else {
if (c1style == VARIABLE) c1 *= xscale;
if (c2style == VARIABLE) c2 *= yscale;
if (rlostyle == VARIABLE) radiuslo *= xscale;
if (rhistyle == VARIABLE) radiushi *= xscale;
if (lostyle == VARIABLE) lo *= zscale;
if (histyle == VARIABLE) hi *= zscale;
}
}
/* ----------------------------------------------------------------------
error check on existence of variable
------------------------------------------------------------------------- */
void RegCone::variable_check()
{
if (c1style == VARIABLE) {
c1var = input->variable->find(c1str);
if (c1var < 0) error->all(FLERR, "Variable {} for region cone does not exist", c1str);
if (!input->variable->equalstyle(c1var))
error->all(FLERR, "Variable {} for region cone is invalid style", c1str);
}
if (c2style == VARIABLE) {
c2var = input->variable->find(c2str);
if (c2var < 0) error->all(FLERR, "Variable {} for region cone does not exist", c2str);
if (!input->variable->equalstyle(c2var))
error->all(FLERR, "Variable {} for region cone is invalid style", c2str);
}
if (rlostyle == VARIABLE) {
rlovar = input->variable->find(rlostr);
if (rlovar < 0) error->all(FLERR, "Variable {} for region cone does not exist", rlostr);
if (!input->variable->equalstyle(rlovar))
error->all(FLERR, "Variable {} for region cone is invalid style", rlostr);
}
if (rhistyle == VARIABLE) {
rhivar = input->variable->find(rhistr);
if (rhivar < 0) error->all(FLERR, "Variable {} for region cone does not exist", rhistr);
if (!input->variable->equalstyle(rhivar))
error->all(FLERR, "Variable {} for region cone is invalid style", rhistr);
}
if (lostyle == VARIABLE) {
lovar = input->variable->find(lostr);
if (lovar < 0) error->all(FLERR, "Variable {} for region cone does not exist", lostr);
if (!input->variable->equalstyle(lovar))
error->all(FLERR, "Variable {} for region cone is invalid style", lostr);
}
if (histyle == VARIABLE) {
hivar = input->variable->find(histr);
if (hivar < 0) error->all(FLERR, "Variable {} for region cone does not exist", histr);
if (!input->variable->equalstyle(hivar))
error->all(FLERR, "Variable {} for region cone is invalid style", histr);
}
}
/* ---------------------------------------------------------------------- */
double RegCone::closest(double *x, double *near, double *nearest, double dsq)

View File

@ -28,9 +28,11 @@ class RegCone : public Region {
public:
RegCone(class LAMMPS *, int, char **);
~RegCone() override;
void init() override;
int inside(double, double, double) override;
int surface_interior(double *, double) override;
int surface_exterior(double *, double) override;
void shape_update() override;
private:
char axis;
@ -39,7 +41,12 @@ class RegCone : public Region {
double lo, hi;
double maxradius;
int c1style, c2style, rlostyle, rhistyle, lostyle, histyle;
int c1var, c2var, rlovar, rhivar, lovar, hivar;
char *c1str, *c2str, *rlostr, *rhistr, *lostr, *histr;
double closest(double *, double *, double *, double);
void variable_check();
};
} // namespace LAMMPS_NS