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

Merge branch 'develop' into add_set_time

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-05-04 15:47:28 -04:00
parent f8742d599b ad6cebb049
commit 059c84bf41
17 changed files with 1752 additions and 154 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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();