tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'develop' into collected-small-changes

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-05-04 15:45:38 -04:00
parent 87b99306ba ad6cebb049
commit eec86c8038
17 changed files with 1752 additions and 154 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Commands_pair.rst
View File

@ -34,7 +34,7 @@ OPT.
*
*
*
* :doc:`adp (o) <pair_adp>`
* :doc:`adp (ko) <pair_adp>`
* :doc:`agni (o) <pair_agni>`
* :doc:`airebo (io) <pair_airebo>`
* :doc:`airebo/morse (io) <pair_airebo>`

11
doc/src/compute_reduce.rst
View File

@ -23,7 +23,7 @@ Syntax
*reduce/region* arg = region-ID
region-ID = ID of region to use for choosing atoms
* mode = *sum* or *min* or *max* or *ave* or *sumsq* or *avesq*
* mode = *sum* or *min* or *max* or *ave* or *sumsq* or *avesq* or *sumabs* or *aveabs*
* one or more inputs can be listed
* input = x, y, z, vx, vy, vz, fx, fy, fz, c_ID, c_ID[N], f_ID, f_ID[N], v_name
@ -77,7 +77,10 @@ option sums the square of the values in the vector into a global
total. The *avesq* setting does the same as *sumsq*, then divides the
sum of squares by the number of values. The last two options can be
useful for calculating the variance of some quantity, e.g. variance =
sumsq - ave\^2.
sumsq - ave\^2. The *sumabs* option sums the absolute values in the
vector into a global total. The *aveabs* setting does the same as
*sumabs*, then divides the sum of absolute values by the number of
values.
Each listed input is operated on independently. For per-atom inputs,
the group specified with this command means only atoms within the
@ -189,7 +192,7 @@ value. If multiple inputs are specified, this compute produces a
global vector of values, the length of which is equal to the number of
inputs specified.
As discussed below, for the *sum* and *sumsq* modes, the value(s)
As discussed below, for the *sum*, *sumabs* and *sumsq* modes, the value(s)
produced by this compute are all "extensive", meaning their value
scales linearly with the number of atoms involved. If normalized
values are desired, this compute can be accessed by the :doc:`thermo_style custom <thermo_style>` command with :doc:`thermo_modify norm yes <thermo_modify>` set as an option. Or it can be accessed by a
@ -208,7 +211,7 @@ compute as input. See the :doc:`Howto output <Howto_output>` doc page
for an overview of LAMMPS output options.
All the scalar or vector values calculated by this compute are
"intensive", except when the *sum* or *sumsq* modes are used on
"intensive", except when the *sum*, *sumabs* or *sumsq* modes are used on
per-atom or local vectors, in which case the calculated values are
"extensive".

11
doc/src/delete_atoms.rst
View File

@ -10,7 +10,7 @@ Syntax
delete_atoms style args keyword value ...
* style = *group* or *region* or *overlap* or *porosity*
* style = *group* or *region* or *overlap* or *porosity* or *variable*
.. parsed-literal::
@ -26,6 +26,7 @@ Syntax
or NULL to only impose the group criterion
fraction = delete this fraction of atoms
seed = random number seed (positive integer)
*variable* args = variable-name
* zero or more keyword/value pairs may be appended
* keyword = *compress* or *bond* or *mol*
@ -47,6 +48,7 @@ Examples
delete_atoms overlap 0.5 solvent colloid
delete_atoms porosity all cube 0.1 482793 bond yes
delete_atoms porosity polymer cube 0.1 482793 bond yes
detele_atoms variable checkers
Description
"""""""""""
@ -91,6 +93,13 @@ guarantee that the exact fraction of atoms will be deleted, or that
the same atoms will be deleted when running on different numbers of
processors.
For style *variable*, all atoms for which the atom-style variable with
the given name evaluates to non-zero will be deleted. Additional atoms
can be deleted if they are in a molecule for which one or more atoms
were deleted within the region; see the *mol* keyword discussion below.
This option allows complex selections of atoms not covered by the
other options listed above.
If the *compress* keyword is set to *yes*, then after atoms are
deleted, then atom IDs are re-assigned so that they run from 1 to the
number of atoms in the system. Note that this is not done for

3
doc/src/pair_adp.rst
View File

@ -1,10 +1,11 @@
.. index:: pair_style adp
.. index:: pair_style adp/kk
.. index:: pair_style adp/omp
pair_style adp command
======================
Accelerator Variants: *adp/omp*
Accelerator Variants: *adp/kk*, *adp/omp*
Syntax
""""""

2
src/KOKKOS/Install.sh
View File

@ -214,6 +214,8 @@ action min_kokkos.h
action min_linesearch_kokkos.cpp
action min_linesearch_kokkos.h
action pack_kokkos.h pack.h
action pair_adp_kokkos.cpp pair_adp.cpp
action pair_adp_kokkos.h pair_adp.h
action pair_buck_coul_cut_kokkos.cpp
action pair_buck_coul_cut_kokkos.h
action pair_buck_coul_long_kokkos.cpp pair_buck_coul_long.cpp

1182
src/KOKKOS/pair_adp_kokkos.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

203
src/KOKKOS/pair_adp_kokkos.h Normal file
View File

@ -0,0 +1,203 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Vladislav Galigerov (HSE), Vsevolod Nikolskiy (HSE)
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
// clang-format off
PairStyle(adp/kk,PairADPKokkos<LMPDeviceType>);
PairStyle(adp/kk/device,PairADPKokkos<LMPDeviceType>);
PairStyle(adp/kk/host,PairADPKokkos<LMPHostType>);
// clang-format on
#else
// clang-format off
#ifndef LMP_PAIR_ADP_KOKKOS_H
#define LMP_PAIR_ADP_KOKKOS_H
#include "kokkos_base.h"
#include "pair_kokkos.h"
#include "pair_adp.h"
#include "neigh_list_kokkos.h"
namespace LAMMPS_NS {
struct TagPairADPPackForwardComm{};
struct TagPairADPUnpackForwardComm{};
struct TagPairADPInitialize{};
template<int NEIGHFLAG, int NEWTON_PAIR>
struct TagPairADPKernelA{};
template<int EFLAG>
struct TagPairADPKernelB{};
template<int EFLAG>
struct TagPairADPKernelAB{};
template<int NEIGHFLAG, int NEWTON_PAIR, int EVFLAG>
struct TagPairADPKernelC{};
template<class DeviceType>
class PairADPKokkos : public PairADP, public KokkosBase
{
public:
enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF};
enum {COUL_FLAG=0};
typedef DeviceType device_type;
typedef ArrayTypes<DeviceType> AT;
typedef EV_FLOAT value_type;
PairADPKokkos(class LAMMPS *);
~PairADPKokkos() override;
void compute(int, int) override;
void init_style() override;
KOKKOS_INLINE_FUNCTION
void operator()(TagPairADPPackForwardComm, const int&) const;
KOKKOS_INLINE_FUNCTION
void operator()(TagPairADPUnpackForwardComm, const int&) const;
KOKKOS_INLINE_FUNCTION
void operator()(TagPairADPInitialize, const int&) const;
template<int NEIGHFLAG, int NEWTON_PAIR>
KOKKOS_INLINE_FUNCTION
void operator()(TagPairADPKernelA<NEIGHFLAG,NEWTON_PAIR>, const int&) const;
template<int EFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagPairADPKernelB<EFLAG>, const int&, EV_FLOAT&) const;
template<int EFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagPairADPKernelB<EFLAG>, const int&) const;
template<int EFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagPairADPKernelAB<EFLAG>, const int&, EV_FLOAT&) const;
template<int EFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagPairADPKernelAB<EFLAG>, const int&) const;
template<int NEIGHFLAG, int NEWTON_PAIR, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagPairADPKernelC<NEIGHFLAG,NEWTON_PAIR,EVFLAG>, const int&, EV_FLOAT&) const;
template<int NEIGHFLAG, int NEWTON_PAIR, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagPairADPKernelC<NEIGHFLAG,NEWTON_PAIR,EVFLAG>, const int&) const;
template<int NEIGHFLAG, int NEWTON_PAIR>
KOKKOS_INLINE_FUNCTION
void ev_tally_xyz(EV_FLOAT &ev, const int &i, const int &j,
const F_FLOAT &epair, const F_FLOAT &fx, const F_FLOAT &fy, const F_FLOAT &fz, const F_FLOAT &delx,
const F_FLOAT &dely, const F_FLOAT &delz) const;
int pack_forward_comm_kokkos(int, DAT::tdual_int_2d, int, DAT::tdual_xfloat_1d&,
int, int *) override;
void unpack_forward_comm_kokkos(int, int, DAT::tdual_xfloat_1d&) override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
int pack_reverse_comm(int, int, double *) override;
void unpack_reverse_comm(int, int *, double *) override;
protected:
typename AT::t_x_array x;
typename AT::t_f_array f;
typename AT::t_int_1d type;
typename AT::t_tagint_1d tag;
DAT::tdual_efloat_1d k_eatom;
DAT::tdual_virial_array k_vatom;
typename AT::t_efloat_1d d_eatom;
typename AT::t_virial_array d_vatom;
int need_dup;
using KKDeviceType = typename KKDevice<DeviceType>::value;
template<typename DataType, typename Layout>
using DupScatterView = KKScatterView<DataType, Layout, KKDeviceType, KKScatterSum, KKScatterDuplicated>;
template<typename DataType, typename Layout>
using NonDupScatterView = KKScatterView<DataType, Layout, KKDeviceType, KKScatterSum, KKScatterNonDuplicated>;
DupScatterView<F_FLOAT*, typename DAT::t_ffloat_1d::array_layout> dup_rho;
DupScatterView<F_FLOAT*[3], typename DAT::t_f_array::array_layout> dup_mu;
DupScatterView<F_FLOAT*[6], typename DAT::t_virial_array::array_layout> dup_lambda;
DupScatterView<F_FLOAT*[3], typename DAT::t_f_array::array_layout> dup_f;
DupScatterView<E_FLOAT*, typename DAT::t_efloat_1d::array_layout> dup_eatom;
DupScatterView<F_FLOAT*[6], typename DAT::t_virial_array::array_layout> dup_vatom;
NonDupScatterView<F_FLOAT*, typename DAT::t_ffloat_1d::array_layout> ndup_rho;
NonDupScatterView<F_FLOAT*[3], typename DAT::t_f_array::array_layout> ndup_mu;
NonDupScatterView<F_FLOAT*[6], typename DAT::t_virial_array::array_layout> ndup_lambda;
NonDupScatterView<F_FLOAT*[3], typename DAT::t_f_array::array_layout> ndup_f;
NonDupScatterView<E_FLOAT*, typename DAT::t_efloat_1d::array_layout> ndup_eatom;
NonDupScatterView<F_FLOAT*[6], typename DAT::t_virial_array::array_layout> ndup_vatom;
DAT::tdual_ffloat_1d k_rho;
DAT::tdual_ffloat_1d k_fp;
DAT::tdual_f_array k_mu;
DAT::tdual_virial_array k_lambda;
typename AT::t_ffloat_1d d_rho;
typename AT::t_ffloat_1d d_fp;
typename AT::t_f_array d_mu;
typename AT::t_virial_array d_lambda;
HAT::t_ffloat_1d h_rho;
HAT::t_ffloat_1d h_fp;
HAT::t_f_array h_mu;
HAT::t_virial_array h_lambda;
typename AT::t_int_1d d_type2frho;
typename AT::t_int_2d_dl d_type2rhor;
typename AT::t_int_2d_dl d_type2z2r;
typename AT::t_int_2d_dl d_type2u2r;
typename AT::t_int_2d_dl d_type2w2r;
typedef Kokkos::DualView<F_FLOAT**[7],DeviceType> tdual_ffloat_2d_n7;
typedef typename tdual_ffloat_2d_n7::t_dev_const t_ffloat_2d_n7;
typedef typename tdual_ffloat_2d_n7::t_host t_host_ffloat_2d_n7;
t_ffloat_2d_n7 d_frho_spline;
t_ffloat_2d_n7 d_rhor_spline;
t_ffloat_2d_n7 d_z2r_spline;
t_ffloat_2d_n7 d_u2r_spline, d_w2r_spline;
void interpolate(int, double, double *, t_host_ffloat_2d_n7, int);
void file2array() override;
void array2spline() override;
typename AT::t_neighbors_2d d_neighbors;
typename AT::t_int_1d d_ilist;
typename AT::t_int_1d d_numneigh;
int iswap;
int first;
typename AT::t_int_2d d_sendlist;
typename AT::t_xfloat_1d_um v_buf;
int neighflag,newton_pair;
int nlocal,nall,eflag,vflag;
friend void pair_virial_fdotr_compute<PairADPKokkos>(PairADPKokkos*);
};
}
#endif
#endif

9
src/KOKKOS/pair_eam_kokkos.cpp
View File

@ -54,10 +54,10 @@ PairEAMKokkos<DeviceType>::PairEAMKokkos(LAMMPS *lmp) : PairEAM(lmp)
template<class DeviceType>
PairEAMKokkos<DeviceType>::~PairEAMKokkos()
{
if (!copymode) {
memoryKK->destroy_kokkos(k_eatom,eatom);
memoryKK->destroy_kokkos(k_vatom,vatom);
}
if (copymode) return;
memoryKK->destroy_kokkos(k_eatom,eatom);
memoryKK->destroy_kokkos(k_vatom,vatom);
}
/* ---------------------------------------------------------------------- */
@ -640,6 +640,7 @@ void PairEAMKokkos<DeviceType>::operator()(TagPairEAMKernelAB<EFLAG>, const int
for (int jj = 0; jj < jnum; jj++) {
int j = d_neighbors(i,jj);
j &= NEIGHMASK;
const X_FLOAT delx = xtmp - x(j,0);
const X_FLOAT dely = ytmp - x(j,1);
const X_FLOAT delz = ztmp - x(j,2);

5
src/KOKKOS/pair_eam_kokkos.h
View File

@ -59,7 +59,6 @@ class PairEAMKokkos : public PairEAM, public KokkosBase {
~PairEAMKokkos() override;
void compute(int, int) override;
void init_style() override;
void *extract(const char *, int &) override { return nullptr; }
KOKKOS_INLINE_FUNCTION
void operator()(TagPairEAMPackForwardComm, const int&) const;
@ -161,21 +160,18 @@ class PairEAMKokkos : public PairEAM, public KokkosBase {
t_ffloat_2d_n7 d_rhor_spline;
t_ffloat_2d_n7 d_z2r_spline;
void interpolate(int, double, double *, t_host_ffloat_2d_n7, int);
void file2array() override;
void array2spline() override;
typename AT::t_neighbors_2d d_neighbors;
typename AT::t_int_1d d_ilist;
typename AT::t_int_1d d_numneigh;
//NeighListKokkos<DeviceType> k_list;
int iswap;
int first;
typename AT::t_int_2d d_sendlist;
typename AT::t_xfloat_1d_um v_buf;
int neighflag,newton_pair;
int nlocal,nall,eflag,vflag;
@ -183,7 +179,6 @@ class PairEAMKokkos : public PairEAM, public KokkosBase {
};
}
#endif
#endif

2
src/MANYBODY/pair_adp.cpp
View File

@ -76,6 +76,8 @@ PairADP::PairADP(LAMMPS *lmp) : Pair(lmp)
PairADP::~PairADP()
{
if (copymode) return;
memory->destroy(rho);
memory->destroy(fp);
memory->destroy(mu);

4
src/MANYBODY/pair_adp.h
View File

@ -78,11 +78,11 @@ class PairADP : public Pair {
Setfl *setfl;
void allocate();
void array2spline();
virtual void array2spline();
void interpolate(int, double, double *, double **);
void read_file(char *);
void file2array();
virtual void file2array();
};
} // namespace LAMMPS_NS

18
src/compute_reduce.cpp
View File

@ -54,6 +54,8 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
mode = SUM;
else if (strcmp(arg[iarg], "sumsq") == 0)
mode = SUMSQ;
else if (strcmp(arg[iarg], "sumabs") == 0)
mode = SUMABS;
else if (strcmp(arg[iarg], "min") == 0)
mode = MINN;
else if (strcmp(arg[iarg], "max") == 0)
@ -62,6 +64,8 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
mode = AVE;
else if (strcmp(arg[iarg], "avesq") == 0)
mode = AVESQ;
else if (strcmp(arg[iarg], "aveabs") == 0)
mode = AVEABS;
else
error->all(FLERR, "Illegal compute {} operation {}", style, arg[iarg]);
iarg++;
@ -253,7 +257,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
if (nvalues == 1) {
scalar_flag = 1;
if (mode == SUM || mode == SUMSQ)
if (mode == SUM || mode == SUMSQ || mode == SUMABS)
extscalar = 1;
else
extscalar = 0;
@ -262,7 +266,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
} else {
vector_flag = 1;
size_vector = nvalues;
if (mode == SUM || mode == SUMSQ)
if (mode == SUM || mode == SUMSQ || mode == SUMABS)
extvector = 1;
else
extvector = 0;
@ -339,13 +343,13 @@ double ComputeReduce::compute_scalar()
double one = compute_one(0, -1);
if (mode == SUM || mode == SUMSQ) {
if (mode == SUM || mode == SUMSQ || mode == SUMABS) {
MPI_Allreduce(&one, &scalar, 1, MPI_DOUBLE, MPI_SUM, world);
} else if (mode == MINN) {
MPI_Allreduce(&one, &scalar, 1, MPI_DOUBLE, MPI_MIN, world);
} else if (mode == MAXX) {
MPI_Allreduce(&one, &scalar, 1, MPI_DOUBLE, MPI_MAX, world);
} else if (mode == AVE || mode == AVESQ) {
} else if (mode == AVE || mode == AVESQ || mode == AVEABS) {
MPI_Allreduce(&one, &scalar, 1, MPI_DOUBLE, MPI_SUM, world);
bigint n = count(0);
if (n) scalar /= n;
@ -366,7 +370,7 @@ void ComputeReduce::compute_vector()
indices[m] = index;
}
if (mode == SUM || mode == SUMSQ) {
if (mode == SUM || mode == SUMSQ || mode == AVEABS) {
for (int m = 0; m < nvalues; m++)
MPI_Allreduce(&onevec[m], &vector[m], 1, MPI_DOUBLE, MPI_SUM, world);
@ -412,7 +416,7 @@ void ComputeReduce::compute_vector()
}
}
} else if (mode == AVE || mode == AVESQ) {
} else if (mode == AVE || mode == AVESQ || mode == AVEABS) {
for (int m = 0; m < nvalues; m++) {
MPI_Allreduce(&onevec[m], &vector[m], 1, MPI_DOUBLE, MPI_SUM, world);
bigint n = count(m);
@ -646,6 +650,8 @@ void ComputeReduce::combine(double &one, double two, int i)
one += two;
else if (mode == SUMSQ || mode == AVESQ)
one += two * two;
else if (mode == SUMABS || mode == AVEABS)
one += std::fabs(two);
else if (mode == MINN) {
if (two < one) {
one = two;

2
src/compute_reduce.h
View File

@ -26,7 +26,7 @@ namespace LAMMPS_NS {
class ComputeReduce : public Compute {
public:
enum { SUM, SUMSQ, MINN, MAXX, AVE, AVESQ };
enum { SUM, SUMSQ, SUMABS, MINN, MAXX, AVE, AVESQ, AVEABS };
enum { PERATOM, LOCAL };
ComputeReduce(class LAMMPS *, int, char **);

296
src/delete_atoms.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -25,6 +24,7 @@
#include "error.h"
#include "force.h"
#include "group.h"
#include "input.h"
#include "memory.h"
#include "modify.h"
#include "molecule.h"
@ -32,6 +32,7 @@
#include "neighbor.h"
#include "random_mars.h"
#include "region.h"
#include "variable.h"
#include <cstring>
#include <map>
@ -48,10 +49,9 @@ DeleteAtoms::DeleteAtoms(LAMMPS *lmp) : Command(lmp) {}
void DeleteAtoms::command(int narg, char **arg)
{
if (domain->box_exist == 0)
error->all(FLERR,"Delete_atoms command before simulation box is defined");
if (narg < 1) error->all(FLERR,"Illegal delete_atoms command");
if (atom->tag_enable == 0)
error->all(FLERR,"Cannot use delete_atoms unless atoms have IDs");
error->all(FLERR, "Delete_atoms command before simulation box is defined");
if (narg < 1) utils::missing_cmd_args(FLERR, "delete_atoms", error);
if (atom->tag_enable == 0) error->all(FLERR, "Cannot use delete_atoms unless atoms have IDs");
// store state before delete
@ -65,26 +65,33 @@ void DeleteAtoms::command(int narg, char **arg)
allflag = 0;
if (strcmp(arg[0],"group") == 0) delete_group(narg,arg);
else if (strcmp(arg[0],"region") == 0) delete_region(narg,arg);
else if (strcmp(arg[0],"overlap") == 0) delete_overlap(narg,arg);
else if (strcmp(arg[0],"porosity") == 0) delete_porosity(narg,arg);
else error->all(FLERR,"Illegal delete_atoms command");
if (strcmp(arg[0], "group") == 0)
delete_group(narg, arg);
else if (strcmp(arg[0], "region") == 0)
delete_region(narg, arg);
else if (strcmp(arg[0], "overlap") == 0)
delete_overlap(narg, arg);
else if (strcmp(arg[0], "porosity") == 0)
delete_porosity(narg, arg);
else if (strcmp(arg[0], "variable") == 0)
delete_variable(narg, arg);
else
error->all(FLERR, "Unknown delete_atoms sub-command: {}", arg[0]);
if (allflag) {
int igroup = group->find("all");
if ((igroup >= 0) &&
modify->check_rigid_group_overlap(group->bitmask[igroup]))
error->warning(FLERR,"Attempting to delete atoms in rigid bodies");
if ((igroup >= 0) && modify->check_rigid_group_overlap(group->bitmask[igroup]))
error->warning(FLERR, "Attempting to delete atoms in rigid bodies");
} else {
if (modify->check_rigid_list_overlap(dlist))
error->warning(FLERR,"Attempting to delete atoms in rigid bodies");
error->warning(FLERR, "Attempting to delete atoms in rigid bodies");
}
// if allflag = 1, just reset atom->nlocal
// else delete atoms one by one
if (allflag) atom->nlocal = 0;
if (allflag)
atom->nlocal = 0;
else {
// optionally delete additional bonds or atoms in molecules
@ -101,10 +108,11 @@ void DeleteAtoms::command(int narg, char **arg)
int i = 0;
while (i < nlocal) {
if (dlist[i]) {
avec->copy(nlocal-1,i,1);
dlist[i] = dlist[nlocal-1];
avec->copy(nlocal - 1, i, 1);
dlist[i] = dlist[nlocal - 1];
nlocal--;
} else i++;
} else
i++;
}
atom->nlocal = nlocal;
@ -122,38 +130,37 @@ void DeleteAtoms::command(int narg, char **arg)
for (int i = 0; i < nlocal; i++) tag[i] = 0;
atom->tag_extend();
} else if (comm->me == 0)
error->warning(FLERR,"Ignoring 'compress yes' for molecular system");
error->warning(FLERR, "Ignoring 'compress yes' for molecular system");
}
// reset atom->natoms and also topology counts
bigint nblocal = atom->nlocal;
MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
MPI_Allreduce(&nblocal, &atom->natoms, 1, MPI_LMP_BIGINT, MPI_SUM, world);
// reset bonus data counts
auto avec_ellipsoid =
dynamic_cast<AtomVecEllipsoid *>( atom->style_match("ellipsoid"));
auto avec_line = dynamic_cast<AtomVecLine *>( atom->style_match("line"));
auto avec_tri = dynamic_cast<AtomVecTri *>( atom->style_match("tri"));
auto avec_body = dynamic_cast<AtomVecBody *>( atom->style_match("body"));
auto avec_ellipsoid = dynamic_cast<AtomVecEllipsoid *>(atom->style_match("ellipsoid"));
auto avec_line = dynamic_cast<AtomVecLine *>(atom->style_match("line"));
auto avec_tri = dynamic_cast<AtomVecTri *>(atom->style_match("tri"));
auto avec_body = dynamic_cast<AtomVecBody *>(atom->style_match("body"));
bigint nlocal_bonus;
if (atom->nellipsoids > 0) {
nlocal_bonus = avec_ellipsoid->nlocal_bonus;
MPI_Allreduce(&nlocal_bonus,&atom->nellipsoids,1,MPI_LMP_BIGINT,MPI_SUM,world);
MPI_Allreduce(&nlocal_bonus, &atom->nellipsoids, 1, MPI_LMP_BIGINT, MPI_SUM, world);
}
if (atom->nlines > 0) {
nlocal_bonus = avec_line->nlocal_bonus;
MPI_Allreduce(&nlocal_bonus,&atom->nlines,1,MPI_LMP_BIGINT,MPI_SUM,world);
MPI_Allreduce(&nlocal_bonus, &atom->nlines, 1, MPI_LMP_BIGINT, MPI_SUM, world);
}
if (atom->ntris > 0) {
nlocal_bonus = avec_tri->nlocal_bonus;
MPI_Allreduce(&nlocal_bonus,&atom->ntris,1,MPI_LMP_BIGINT,MPI_SUM,world);
MPI_Allreduce(&nlocal_bonus, &atom->ntris, 1, MPI_LMP_BIGINT, MPI_SUM, world);
}
if (atom->nbodies > 0) {
nlocal_bonus = avec_body->nlocal_bonus;
MPI_Allreduce(&nlocal_bonus,&atom->nbodies,1,MPI_LMP_BIGINT,MPI_SUM,world);
MPI_Allreduce(&nlocal_bonus, &atom->nbodies, 1, MPI_LMP_BIGINT, MPI_SUM, world);
}
// reset atom->map if it exists
@ -176,23 +183,20 @@ void DeleteAtoms::command(int narg, char **arg)
bigint ndelete_impropers = nimpropers_previous - atom->nimpropers;
if (comm->me == 0) {
std::string mesg = fmt::format("Deleted {} atoms, new total = {}\n",
ndelete,atom->natoms);
std::string mesg = fmt::format("Deleted {} atoms, new total = {}\n", ndelete, atom->natoms);
if (bond_flag || mol_flag) {
if (nbonds_previous)
mesg += fmt::format("Deleted {} bonds, new total = {}\n",
ndelete_bonds,atom->nbonds);
mesg += fmt::format("Deleted {} bonds, new total = {}\n", ndelete_bonds, atom->nbonds);
if (nangles_previous)
mesg += fmt::format("Deleted {} angles, new total = {}\n",
ndelete_angles,atom->nangles);
mesg += fmt::format("Deleted {} angles, new total = {}\n", ndelete_angles, atom->nangles);
if (ndihedrals_previous)
mesg += fmt::format("Deleted {} dihedrals, new total = {}\n",
ndelete_dihedrals,atom->ndihedrals);
mesg += fmt::format("Deleted {} dihedrals, new total = {}\n", ndelete_dihedrals,
atom->ndihedrals);
if (nimpropers_previous)
mesg += fmt::format("Deleted {} impropers, new total = {}\n",
ndelete_impropers,atom->nimpropers);
mesg += fmt::format("Deleted {} impropers, new total = {}\n", ndelete_impropers,
atom->nimpropers);
}
utils::logmesg(lmp,mesg);
utils::logmesg(lmp, mesg);
}
}
@ -202,15 +206,15 @@ void DeleteAtoms::command(int narg, char **arg)
void DeleteAtoms::delete_group(int narg, char **arg)
{
if (narg < 2) error->all(FLERR,"Illegal delete_atoms command");
if (narg < 2) utils::missing_cmd_args(FLERR, "delete_atoms group", error);
int igroup = group->find(arg[1]);
if (igroup == -1) error->all(FLERR,"Could not find delete_atoms group ID");
options(narg-2,&arg[2]);
if (igroup == -1) error->all(FLERR, "Could not find delete_atoms group ID {}", arg[1]);
options(narg - 2, &arg[2]);
// check for special case of group = all
if (strcmp(arg[1],"all") == 0) {
if (strcmp(arg[1], "all") == 0) {
allflag = 1;
return;
}
@ -218,7 +222,7 @@ void DeleteAtoms::delete_group(int narg, char **arg)
// allocate and initialize deletion list
int nlocal = atom->nlocal;
memory->create(dlist,nlocal,"delete_atoms:dlist");
memory->create(dlist, nlocal, "delete_atoms:dlist");
for (int i = 0; i < nlocal; i++) dlist[i] = 0;
int *mask = atom->mask;
@ -234,24 +238,24 @@ void DeleteAtoms::delete_group(int narg, char **arg)
void DeleteAtoms::delete_region(int narg, char **arg)
{
if (narg < 2) error->all(FLERR,"Illegal delete_atoms command");
if (narg < 2) utils::missing_cmd_args(FLERR, "delete_atoms region", error);
auto iregion = domain->get_region_by_id(arg[1]);
if (!iregion) error->all(FLERR,"Could not find delete_atoms region ID");
if (!iregion) error->all(FLERR, "Could not find delete_atoms region ID {}", arg[1]);
iregion->prematch();
options(narg-2,&arg[2]);
options(narg - 2, &arg[2]);
// allocate and initialize deletion list
int nlocal = atom->nlocal;
memory->create(dlist,nlocal,"delete_atoms:dlist");
memory->create(dlist, nlocal, "delete_atoms:dlist");
for (int i = 0; i < nlocal; i++) dlist[i] = 0;
double **x = atom->x;
for (int i = 0; i < nlocal; i++)
if (iregion->match(x[i][0],x[i][1],x[i][2])) dlist[i] = 1;
if (iregion->match(x[i][0], x[i][1], x[i][2])) dlist[i] = 1;
}
/* ----------------------------------------------------------------------
@ -263,23 +267,23 @@ void DeleteAtoms::delete_region(int narg, char **arg)
void DeleteAtoms::delete_overlap(int narg, char **arg)
{
if (narg < 4) error->all(FLERR,"Illegal delete_atoms command");
if (narg < 4) utils::missing_cmd_args(FLERR, "delete_atoms overlap", error);
// read args
double cut = utils::numeric(FLERR,arg[1],false,lmp);
double cutsq = cut*cut;
double cut = utils::numeric(FLERR, arg[1], false, lmp);
double cutsq = cut * cut;
int igroup1 = group->find(arg[2]);
int igroup2 = group->find(arg[3]);
if (igroup1 < 0 || igroup2 < 0)
error->all(FLERR,"Could not find delete_atoms group ID");
options(narg-4,&arg[4]);
error->all(FLERR, "Could not find delete_atoms group ID {}", arg[1]);
options(narg - 4, &arg[4]);
int group1bit = group->bitmask[igroup1];
int group2bit = group->bitmask[igroup2];
if (comm->me == 0) utils::logmesg(lmp,"System init for delete_atoms ...\n");
if (comm->me == 0) utils::logmesg(lmp, "System init for delete_atoms ...\n");
// request a full neighbor list for use by this command
@ -293,12 +297,10 @@ void DeleteAtoms::delete_overlap(int narg, char **arg)
// error check on cutoff
// if no pair style, neighbor list will be empty
if (force->pair == nullptr)
error->all(FLERR,"Delete_atoms requires a pair style be defined");
if (cut > neighbor->cutneighmax)
error->all(FLERR,"Delete_atoms cutoff > max neighbor cutoff");
if (force->pair == nullptr) error->all(FLERR, "Delete_atoms requires a pair style be defined");
if (cut > neighbor->cutneighmax) error->all(FLERR, "Delete_atoms cutoff > max neighbor cutoff");
if (cut > neighbor->cutneighmin && comm->me == 0)
error->warning(FLERR,"Delete_atoms cutoff > minimum neighbor cutoff");
error->warning(FLERR, "Delete_atoms cutoff > minimum neighbor cutoff");
// setup domain, communication and neighboring
// acquire ghosts and build standard neighbor lists
@ -310,7 +312,7 @@ void DeleteAtoms::delete_overlap(int narg, char **arg)
if (neighbor->style) neighbor->setup_bins();
comm->exchange();
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
if (domain->triclinic) domain->lamda2x(atom->nlocal + atom->nghost);
neighbor->build(1);
// build neighbor list this command needs based on the earlier request
@ -322,7 +324,7 @@ void DeleteAtoms::delete_overlap(int narg, char **arg)
// must be after exchange potentially changes nlocal
int nlocal = atom->nlocal;
memory->create(dlist,nlocal,"delete_atoms:dlist");
memory->create(dlist, nlocal, "delete_atoms:dlist");
for (int i = 0; i < nlocal; i++) dlist[i] = 0;
// double loop over owned atoms and their full neighbor list
@ -335,10 +337,10 @@ void DeleteAtoms::delete_overlap(int narg, char **arg)
double *special_coul = force->special_coul;
double *special_lj = force->special_lj;
int i,j,ii,jj,inum,jnum;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *ilist,*jlist,*numneigh,**firstneigh;
double factor_lj,factor_coul;
int i, j, ii, jj, inum, jnum;
double xtmp, ytmp, ztmp, delx, dely, delz, rsq;
int *ilist, *jlist, *numneigh, **firstneigh;
double factor_lj, factor_coul;
inum = list->inum;
ilist = list->ilist;
@ -378,7 +380,7 @@ void DeleteAtoms::delete_overlap(int narg, char **arg)
dely = x[j][1] - ytmp;
delz = x[j][2] - ztmp;
}
rsq = delx*delx + dely*dely + delz*delz;
rsq = delx * delx + dely * dely + delz * delz;
if (rsq >= cutsq) continue;
// only consider deletion if I,J are in groups 1,2 respectively
@ -415,25 +417,25 @@ void DeleteAtoms::delete_overlap(int narg, char **arg)
void DeleteAtoms::delete_porosity(int narg, char **arg)
{
if (narg < 5) error->all(FLERR,"Illegal delete_atoms command");
if (narg < 5) utils::missing_cmd_args(FLERR, "delete_atoms porosity", error);
int igroup = group->find(arg[1]);
if (igroup == -1) error->all(FLERR,"Could not find delete_atoms group ID");
if (igroup == -1) error->all(FLERR, "Could not find delete_atoms porosity group ID {}", arg[1]);
auto iregion = domain->get_region_by_id(arg[2]);
if (!iregion && (strcmp(arg[2],"NULL") != 0))
error->all(FLERR,"Could not find delete_atoms region ID");
if (!iregion && (strcmp(arg[2], "NULL") != 0))
error->all(FLERR, "Could not find delete_atoms porosity region ID {}", arg[2]);
double porosity_fraction = utils::numeric(FLERR,arg[3],false,lmp);
int seed = utils::inumeric(FLERR,arg[4],false,lmp);
options(narg-5,&arg[5]);
double porosity_fraction = utils::numeric(FLERR, arg[3], false, lmp);
int seed = utils::inumeric(FLERR, arg[4], false, lmp);
options(narg - 5, &arg[5]);
auto random = new RanMars(lmp,seed + comm->me);
auto random = new RanMars(lmp, seed + comm->me);
// allocate and initialize deletion list
int nlocal = atom->nlocal;
memory->create(dlist,nlocal,"delete_atoms:dlist");
memory->create(dlist, nlocal, "delete_atoms:dlist");
for (int i = 0; i < nlocal; i++) dlist[i] = 0;
// delete fraction of atoms which are in both group and region
@ -446,13 +448,45 @@ void DeleteAtoms::delete_porosity(int narg, char **arg)
for (int i = 0; i < nlocal; i++) {
if (!(mask[i] & groupbit)) continue;
if (iregion && !iregion->match(x[i][0],x[i][1],x[i][2])) continue;
if (iregion && !iregion->match(x[i][0], x[i][1], x[i][2])) continue;
if (random->uniform() <= porosity_fraction) dlist[i] = 1;
}
delete random;
}
/* ----------------------------------------------------------------------
delete all as flagged by non-zero atom style variable
------------------------------------------------------------------------- */
void DeleteAtoms::delete_variable(int narg, char **arg)
{
if (narg < 2) utils::missing_cmd_args(FLERR, "delete_atoms variable", error);
int ivar = input->variable->find(arg[1]);
if (ivar < 0) error->all(FLERR, "Variable name {} for delete_atoms does not exist", arg[1]);
if (!input->variable->atomstyle(ivar))
error->all(FLERR, "Variable {} for delete_atoms is invalid style", arg[1]);
// consume remaining options
options(narg - 2, &arg[2]);
// aflag = evaluation of per-atom variable
const int nlocal = atom->nlocal;
double *aflag;
memory->create(dlist, nlocal, "delete_atoms:dlist");
memory->create(aflag, nlocal, "group:aflag");
input->variable->compute_atom(ivar, 0, aflag, 1, 0);
// delete if per-atom variable evaluated to non-zero
for (int i = 0; i < nlocal; i++) dlist[i] = (aflag[i] == 0.0) ? 0 : 1;
memory->destroy(aflag);
}
/* ----------------------------------------------------------------------
delete all topology interactions that include deleted atoms
------------------------------------------------------------------------- */
@ -463,7 +497,7 @@ void DeleteAtoms::delete_bond()
// list of these IDs is sent around ring
// at each stage of ring pass, hash is re-populated with received IDs
hash = new std::map<tagint,int>();
hash = new std::map<tagint, int>();
// list = set of unique molecule IDs from which I deleted atoms
// pass list to all other procs via comm->ring()
@ -475,13 +509,13 @@ void DeleteAtoms::delete_bond()
for (int i = 0; i < nlocal; i++)
if (dlist[i]) n++;
tagint *list;
memory->create(list,n,"delete_atoms:list");
memory->create(list, n, "delete_atoms:list");
n = 0;
for (int i = 0; i < nlocal; i++)
if (dlist[i]) list[n++] = tag[i];
comm->ring(n,sizeof(tagint),list,1,bondring,nullptr,(void *)this);
comm->ring(n, sizeof(tagint), list, 1, bondring, nullptr, (void *) this);
delete hash;
memory->destroy(list);
@ -497,15 +531,14 @@ void DeleteAtoms::delete_molecule()
{
// hash = unique molecule IDs from which I deleted atoms
hash = new std::map<tagint,int>();
hash = new std::map<tagint, int>();
tagint *molecule = atom->molecule;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (molecule[i] == 0) continue;
if (dlist[i] && hash->find(molecule[i]) == hash->end())
(*hash)[molecule[i]] = 1;
if (dlist[i] && hash->find(molecule[i]) == hash->end()) (*hash)[molecule[i]] = 1;
}
// list = set of unique molecule IDs from which I deleted atoms
@ -513,13 +546,13 @@ void DeleteAtoms::delete_molecule()
int n = hash->size();
tagint *list;
memory->create(list,n,"delete_atoms:list");
memory->create(list, n, "delete_atoms:list");
n = 0;
std::map<tagint,int>::iterator pos;
std::map<tagint, int>::iterator pos;
for (pos = hash->begin(); pos != hash->end(); ++pos) list[n++] = pos->first;
comm->ring(n,sizeof(tagint),list,1,molring,nullptr,(void *)this);
comm->ring(n, sizeof(tagint), list, 1, molring, nullptr, (void *) this);
delete hash;
memory->destroy(list);
@ -557,7 +590,7 @@ void DeleteAtoms::recount_topology()
int *molatom = atom->molatom;
int nlocal = atom->nlocal;
int imol,iatom;
int imol, iatom;
for (int i = 0; i < nlocal; i++) {
imol = molindex[i];
@ -571,19 +604,19 @@ void DeleteAtoms::recount_topology()
}
if (atom->avec->bonds_allow) {
MPI_Allreduce(&nbonds,&atom->nbonds,1,MPI_LMP_BIGINT,MPI_SUM,world);
MPI_Allreduce(&nbonds, &atom->nbonds, 1, MPI_LMP_BIGINT, MPI_SUM, world);
if (!force->newton_bond) atom->nbonds /= 2;
}
if (atom->avec->angles_allow) {
MPI_Allreduce(&nangles,&atom->nangles,1,MPI_LMP_BIGINT,MPI_SUM,world);
MPI_Allreduce(&nangles, &atom->nangles, 1, MPI_LMP_BIGINT, MPI_SUM, world);
if (!force->newton_bond) atom->nangles /= 3;
}
if (atom->avec->dihedrals_allow) {
MPI_Allreduce(&ndihedrals,&atom->ndihedrals,1,MPI_LMP_BIGINT,MPI_SUM,world);
MPI_Allreduce(&ndihedrals, &atom->ndihedrals, 1, MPI_LMP_BIGINT, MPI_SUM, world);
if (!force->newton_bond) atom->ndihedrals /= 4;
}
if (atom->avec->impropers_allow) {
MPI_Allreduce(&nimpropers,&atom->nimpropers,1,MPI_LMP_BIGINT,MPI_SUM,world);
MPI_Allreduce(&nimpropers, &atom->nimpropers, 1, MPI_LMP_BIGINT, MPI_SUM, world);
if (!force->newton_bond) atom->nimpropers /= 4;
}
}
@ -596,7 +629,7 @@ void DeleteAtoms::bondring(int nbuf, char *cbuf, void *ptr)
{
auto daptr = (DeleteAtoms *) ptr;
auto list = (tagint *) cbuf;
std::map<tagint,int> *hash = daptr->hash;
std::map<tagint, int> *hash = daptr->hash;
int *num_bond = daptr->atom->num_bond;
int *num_angle = daptr->atom->num_angle;
@ -633,17 +666,18 @@ void DeleteAtoms::bondring(int nbuf, char *cbuf, void *ptr)
// loop over my atoms and their bond topology lists
// if any atom in an interaction matches atom ID in hash, delete interaction
int m,n;
int m, n;
for (int i = 0; i < nlocal; i++) {
if (num_bond) {
m = 0;
n = num_bond[i];
while (m < n) {
if (hash->find(bond_atom[i][m]) != hash->end()) {
bond_type[i][m] = bond_type[i][n-1];
bond_atom[i][m] = bond_atom[i][n-1];
bond_type[i][m] = bond_type[i][n - 1];
bond_atom[i][m] = bond_atom[i][n - 1];
n--;
} else m++;
} else
m++;
}
num_bond[i] = n;
}
@ -655,12 +689,13 @@ void DeleteAtoms::bondring(int nbuf, char *cbuf, void *ptr)
if (hash->find(angle_atom1[i][m]) != hash->end() ||
hash->find(angle_atom2[i][m]) != hash->end() ||
hash->find(angle_atom3[i][m]) != hash->end()) {
angle_type[i][m] = angle_type[i][n-1];
angle_atom1[i][m] = angle_atom1[i][n-1];
angle_atom2[i][m] = angle_atom2[i][n-1];
angle_atom3[i][m] = angle_atom3[i][n-1];
angle_type[i][m] = angle_type[i][n - 1];
angle_atom1[i][m] = angle_atom1[i][n - 1];
angle_atom2[i][m] = angle_atom2[i][n - 1];
angle_atom3[i][m] = angle_atom3[i][n - 1];
n--;
} else m++;
} else
m++;
}
num_angle[i] = n;
}
@ -673,13 +708,14 @@ void DeleteAtoms::bondring(int nbuf, char *cbuf, void *ptr)
hash->find(dihedral_atom2[i][m]) != hash->end() ||
hash->find(dihedral_atom3[i][m]) != hash->end() ||
hash->find(dihedral_atom4[i][m]) != hash->end()) {
dihedral_type[i][m] = dihedral_type[i][n-1];
dihedral_atom1[i][m] = dihedral_atom1[i][n-1];
dihedral_atom2[i][m] = dihedral_atom2[i][n-1];
dihedral_atom3[i][m] = dihedral_atom3[i][n-1];
dihedral_atom4[i][m] = dihedral_atom4[i][n-1];
dihedral_type[i][m] = dihedral_type[i][n - 1];
dihedral_atom1[i][m] = dihedral_atom1[i][n - 1];
dihedral_atom2[i][m] = dihedral_atom2[i][n - 1];
dihedral_atom3[i][m] = dihedral_atom3[i][n - 1];
dihedral_atom4[i][m] = dihedral_atom4[i][n - 1];
n--;
} else m++;
} else
m++;
}
num_dihedral[i] = n;
}
@ -692,13 +728,14 @@ void DeleteAtoms::bondring(int nbuf, char *cbuf, void *ptr)
hash->find(improper_atom2[i][m]) != hash->end() ||
hash->find(improper_atom3[i][m]) != hash->end() ||
hash->find(improper_atom4[i][m]) != hash->end()) {
improper_type[i][m] = improper_type[i][n-1];
improper_atom1[i][m] = improper_atom1[i][n-1];
improper_atom2[i][m] = improper_atom2[i][n-1];
improper_atom3[i][m] = improper_atom3[i][n-1];
improper_atom4[i][m] = improper_atom4[i][n-1];
improper_type[i][m] = improper_type[i][n - 1];
improper_atom1[i][m] = improper_atom1[i][n - 1];
improper_atom2[i][m] = improper_atom2[i][n - 1];
improper_atom3[i][m] = improper_atom3[i][n - 1];
improper_atom4[i][m] = improper_atom4[i][n - 1];
n--;
} else m++;
} else
m++;
}
num_improper[i] = n;
}
@ -711,10 +748,10 @@ void DeleteAtoms::bondring(int nbuf, char *cbuf, void *ptr)
void DeleteAtoms::molring(int n, char *cbuf, void *ptr)
{
auto daptr = (DeleteAtoms *)ptr;
auto daptr = (DeleteAtoms *) ptr;
auto list = (tagint *) cbuf;
int *dlist = daptr->dlist;
std::map<tagint,int> *hash = daptr->hash;
std::map<tagint, int> *hash = daptr->hash;
int nlocal = daptr->atom->nlocal;
tagint *molecule = daptr->atom->molecule;
@ -740,24 +777,25 @@ void DeleteAtoms::options(int narg, char **arg)
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"compress") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal delete_atoms command");
compress_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
if (strcmp(arg[iarg], "compress") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "delete_atoms compress", error);
compress_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"bond") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal delete_atoms command");
} else if (strcmp(arg[iarg], "bond") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "delete_atoms bond", error);
if (atom->molecular == Atom::ATOMIC)
error->all(FLERR,"Cannot delete_atoms bond yes for non-molecular systems");
error->all(FLERR, "Cannot use delete_atoms bond yes for non-molecular systems");
if (atom->molecular == Atom::TEMPLATE)
error->all(FLERR,"Cannot use delete_atoms bond yes with atom_style template");
bond_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
error->all(FLERR, "Cannot use delete_atoms bond yes with atom_style template");
bond_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"mol") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal delete_atoms command");
} else if (strcmp(arg[iarg], "mol") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "delete_atoms mol", error);
if (atom->molecule_flag == 0)
error->all(FLERR,"Delete_atoms mol yes requires atom attribute molecule");
mol_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
error->all(FLERR, "Delete_atoms mol yes requires atom attribute molecule");
mol_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else error->all(FLERR,"Illegal delete_atoms command");
} else
error->all(FLERR, "Unknown delete_atoms option: {}", arg[iarg]);
}
}

1
src/delete_atoms.h
View File

@ -39,6 +39,7 @@ class DeleteAtoms : public Command {
void delete_region(int, char **);
void delete_overlap(int, char **);
void delete_porosity(int, char **);
void delete_variable(int, char **);
void delete_bond();
void delete_molecule();

4
unittest/commands/CMakeLists.txt
View File

@ -25,6 +25,10 @@ add_executable(test_regions test_regions.cpp)
target_link_libraries(test_regions PRIVATE lammps GTest::GMock)
add_test(NAME Regions COMMAND test_regions)
add_executable(test_delete_atoms test_delete_atoms.cpp)
target_link_libraries(test_delete_atoms PRIVATE lammps GTest::GMock)
add_test(NAME DeleteAtoms COMMAND test_delete_atoms)
add_executable(test_variables test_variables.cpp)
target_link_libraries(test_variables PRIVATE lammps GTest::GMock)
add_test(NAME Variables COMMAND test_variables)

151
unittest/commands/test_delete_atoms.cpp Normal file
View File

@ -0,0 +1,151 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "lammps.h"
#include "atom.h"
#include "group.h"
#include "info.h"
#include "input.h"
#include "math_const.h"
#include "region.h"
#include "variable.h"
#include "../testing/core.h"
#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include <cstring>
#include <vector>
// whether to print verbose output (i.e. not capturing LAMMPS screen output).
bool verbose = false;
using LAMMPS_NS::MathConst::MY_PI;
using LAMMPS_NS::utils::split_words;
namespace LAMMPS_NS {
using ::testing::ContainsRegex;
using ::testing::ExitedWithCode;
using ::testing::StrEq;
class DeleteAtomsTest : public LAMMPSTest {
protected:
Atom *atom;
void SetUp() override
{
testbinary = "DeleteAtomsTest";
args = {"-log", "none", "-echo", "screen", "-nocite", "-v", "num", "1"};
LAMMPSTest::SetUp();
atom = lmp->atom;
}
void TearDown() override
{
LAMMPSTest::TearDown();
platform::unlink("test_variable.file");
platform::unlink("test_variable.atomfile");
}
void atomic_system()
{
BEGIN_HIDE_OUTPUT();
command("units real");
command("lattice sc 1.0 origin 0.125 0.125 0.125");
command("region box block -4 4 -4 4 -4 4");
command("create_box 8 box");
command("create_atoms 1 box");
command("mass * 1.0");
command("region left block -2.0 -1.0 INF INF INF INF");
command("region right block 0.5 2.0 INF INF INF INF");
command("region top block INF INF -2.0 -1.0 INF INF");
command("set region left type 2");
command("set region right type 3");
command("group top region top");
END_HIDE_OUTPUT();
}
void molecular_system()
{
HIDE_OUTPUT([&] {
command("fix props all property/atom mol rmass q");
});
atomic_system();
BEGIN_HIDE_OUTPUT();
command("variable molid atom floor(id/4)+1");
command("variable charge atom 2.0*sin(PI/32*id)");
command("set atom * mol v_molid");
command("set atom * charge v_charge");
command("set type 1 mass 0.5");
command("set type 2*4 mass 2.0");
END_HIDE_OUTPUT();
}
};
TEST_F(DeleteAtomsTest, Simple)
{
atomic_system();
ASSERT_EQ(atom->natoms, 512);
HIDE_OUTPUT([&] {
command("delete_atoms group top");
});
ASSERT_EQ(atom->natoms, 448);
HIDE_OUTPUT([&] {
command("delete_atoms region left");
});
ASSERT_EQ(atom->natoms, 392);
HIDE_OUTPUT([&] {
command("delete_atoms porosity all right 0.5 43252");
});
ASSERT_EQ(atom->natoms, 362);
HIDE_OUTPUT([&] {
command("variable checker atom sin(4*PI*x/lx)*sin(4*PI*y/ly)*sin(4*PI*z/lz)>0");
command("delete_atoms variable checker");
});
ASSERT_EQ(atom->natoms, 177);
TEST_FAILURE(".*ERROR: Illegal delete_atoms command: missing argument.*",
command("delete_atoms"););
TEST_FAILURE(".*ERROR: Unknown delete_atoms sub-command: xxx.*", command("delete_atoms xxx"););
}
} // namespace LAMMPS_NS
int main(int argc, char **argv)
{
MPI_Init(&argc, &argv);
::testing::InitGoogleMock(&argc, argv);
if (platform::mpi_vendor() == "Open MPI" && !LAMMPS_NS::Info::has_exceptions())
std::cout << "Warning: using OpenMPI without exceptions. Death tests will be skipped\n";
// handle arguments passed via environment variable
if (const char *var = getenv("TEST_ARGS")) {
std::vector<std::string> env = split_words(var);
for (auto arg : env) {
if (arg == "-v") {
verbose = true;
}
}
}
if ((argc > 1) && (strcmp(argv[1], "-v") == 0)) verbose = true;
int rv = RUN_ALL_TESTS();
MPI_Finalize();
return rv;
}