diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index 4355a12d11..9d44ee5826 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -34,7 +34,7 @@ OPT. * * * - * :doc:`adp (o) ` + * :doc:`adp (ko) ` * :doc:`agni (o) ` * :doc:`airebo (io) ` * :doc:`airebo/morse (io) ` diff --git a/doc/src/compute_reduce.rst b/doc/src/compute_reduce.rst index 1793830e0a..c2d85753db 100644 --- a/doc/src/compute_reduce.rst +++ b/doc/src/compute_reduce.rst @@ -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 ` command with :doc:`thermo_modify norm yes ` set as an option. Or it can be accessed by a @@ -208,7 +211,7 @@ compute as input. See the :doc:`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". diff --git a/doc/src/delete_atoms.rst b/doc/src/delete_atoms.rst index f78f295011..325460c8d3 100644 --- a/doc/src/delete_atoms.rst +++ b/doc/src/delete_atoms.rst @@ -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 diff --git a/doc/src/pair_adp.rst b/doc/src/pair_adp.rst index 87493724b0..63045a4792 100644 --- a/doc/src/pair_adp.rst +++ b/doc/src/pair_adp.rst @@ -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 """""" diff --git a/src/KOKKOS/Install.sh b/src/KOKKOS/Install.sh index 65ac13c061..da397cd677 100755 --- a/src/KOKKOS/Install.sh +++ b/src/KOKKOS/Install.sh @@ -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 diff --git a/src/KOKKOS/pair_adp_kokkos.cpp b/src/KOKKOS/pair_adp_kokkos.cpp new file mode 100644 index 0000000000..93806a8dc5 --- /dev/null +++ b/src/KOKKOS/pair_adp_kokkos.cpp @@ -0,0 +1,1182 @@ +// clang-format off +/* ---------------------------------------------------------------------- + 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) +------------------------------------------------------------------------- */ + +#include "pair_adp_kokkos.h" + +#include "atom_kokkos.h" +#include "atom_masks.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "kokkos.h" +#include "memory_kokkos.h" +#include "neigh_list_kokkos.h" +#include "neigh_request.h" +#include "neighbor.h" +#include "pair_kokkos.h" + +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +template +PairADPKokkos::PairADPKokkos(LAMMPS *lmp) : PairADP(lmp) +{ + respa_enable = 0; + single_enable = 0; + + kokkosable = 1; + atomKK = (AtomKokkos *) atom; + execution_space = ExecutionSpaceFromDevice::space; + datamask_read = X_MASK | F_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK; + datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK; +} + +/* ---------------------------------------------------------------------- */ + +template +PairADPKokkos::~PairADPKokkos() +{ + if (copymode) return; + + memoryKK->destroy_kokkos(k_eatom,eatom); + memoryKK->destroy_kokkos(k_vatom,vatom); +} + +/* ---------------------------------------------------------------------- */ + +template +void PairADPKokkos::compute(int eflag_in, int vflag_in) +{ + eflag = eflag_in; + vflag = vflag_in; + + if (neighflag == FULL) no_virial_fdotr_compute = 1; + + ev_init(eflag,vflag,0); + + // reallocate per-atom arrays if necessary + + if (eflag_atom) { + memoryKK->destroy_kokkos(k_eatom,eatom); + memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom"); + d_eatom = k_eatom.view(); + } + if (vflag_atom) { + memoryKK->destroy_kokkos(k_vatom,vatom); + memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom"); + d_vatom = k_vatom.view(); + } + + atomKK->sync(execution_space,datamask_read); + if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); + else atomKK->modified(execution_space,F_MASK); + + // grow energy and fp arrays if necessary + // need to be atom->nmax in length + + if (atom->nmax > nmax) { + nmax = atom->nmax; + k_rho = DAT::tdual_ffloat_1d("pair:rho",nmax); + k_fp = DAT::tdual_ffloat_1d("pair:fp",nmax); + k_mu = DAT::tdual_f_array("pair:mu", nmax); + k_lambda = DAT::tdual_virial_array("pair:lambda", nmax); + d_rho = k_rho.template view(); + d_fp = k_fp.template view(); + d_mu = k_mu.template view(); + d_lambda = k_lambda.template view(); + h_rho = k_rho.h_view; + h_fp = k_fp.h_view; + h_mu = k_mu.h_view; + h_lambda = k_lambda.h_view; + } + + x = atomKK->k_x.view(); + f = atomKK->k_f.view(); + type = atomKK->k_type.view(); + tag = atomKK->k_tag.view(); + nlocal = atom->nlocal; + nall = atom->nlocal + atom->nghost; + newton_pair = force->newton_pair; + + NeighListKokkos* k_list = static_cast*>(list); + d_numneigh = k_list->d_numneigh; + d_neighbors = k_list->d_neighbors; + d_ilist = k_list->d_ilist; + int inum = list->inum; + + need_dup = lmp->kokkos->need_dup(); + if (need_dup) { + dup_rho = Kokkos::Experimental::create_scatter_view(d_rho); + dup_mu = Kokkos::Experimental::create_scatter_view(d_mu); + dup_lambda = Kokkos::Experimental::create_scatter_view(d_lambda); + dup_f = Kokkos::Experimental::create_scatter_view(f); + dup_eatom = Kokkos::Experimental::create_scatter_view(d_eatom); + dup_vatom = Kokkos::Experimental::create_scatter_view(d_vatom); + } else { + ndup_rho = Kokkos::Experimental::create_scatter_view(d_rho); + ndup_mu = Kokkos::Experimental::create_scatter_view(d_mu); + ndup_lambda = Kokkos::Experimental::create_scatter_view(d_lambda); + ndup_f = Kokkos::Experimental::create_scatter_view(f); + ndup_eatom = Kokkos::Experimental::create_scatter_view(d_eatom); + ndup_vatom = Kokkos::Experimental::create_scatter_view(d_vatom); + } + + copymode = 1; + + // zero out density + + if (newton_pair) + Kokkos::parallel_for(Kokkos::RangePolicy(0,nall),*this); + else + Kokkos::parallel_for(Kokkos::RangePolicy(0,nlocal),*this); + + // loop over neighbors of my atoms + + EV_FLOAT ev; + + // compute kernel A + + if (neighflag == HALF || neighflag == HALFTHREAD) { + + if (neighflag == HALF) { + if (newton_pair) { + Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } else { + Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } + } else if (neighflag == HALFTHREAD) { + if (newton_pair) { + Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } else { + Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } + } + + if (need_dup) + { + Kokkos::Experimental::contribute(d_rho, dup_rho); + Kokkos::Experimental::contribute(d_mu, dup_mu); + Kokkos::Experimental::contribute(d_lambda, dup_lambda); + } + + // communicate and sum densities (on the host) + + if (newton_pair) { + k_rho.template modify(); + k_mu.template modify(); + k_lambda.template modify(); + comm->reverse_comm(this); + k_rho.template sync(); + k_mu.template sync(); + k_lambda.template sync(); + } + + // compute kernel B + + if (eflag) + Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + else + Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + + } else if (neighflag == FULL) { + + // compute kernel AB + + if (eflag) + Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + else + Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } + + if (eflag) { + eng_vdwl += ev.evdwl; + ev.evdwl = 0.0; + } + + // communicate derivative of embedding function + + k_fp.template modify(); + k_mu.template modify(); + k_lambda.template modify(); + comm->forward_comm(this); + k_fp.template sync(); + k_mu.template sync(); + k_lambda.template sync(); + + // compute kernel C + + if (evflag) { + if (neighflag == HALF) { + if (newton_pair) { + Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + } else { + Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + } + } else if (neighflag == HALFTHREAD) { + if (newton_pair) { + Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + } else { + Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + } + } else if (neighflag == FULL) { + if (newton_pair) { + Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + } else { + Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + } + } + } else { + if (neighflag == HALF) { + if (newton_pair) { + Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } else { + Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } + } else if (neighflag == HALFTHREAD) { + if (newton_pair) { + Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } else { + Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } + } else if (neighflag == FULL) { + if (newton_pair) { + Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } else { + Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } + } + } + + if (need_dup) + Kokkos::Experimental::contribute(f, dup_f); + + if (eflag_global) eng_vdwl += ev.evdwl; + if (vflag_global) { + virial[0] += ev.v[0]; + virial[1] += ev.v[1]; + virial[2] += ev.v[2]; + virial[3] += ev.v[3]; + virial[4] += ev.v[4]; + virial[5] += ev.v[5]; + } + + if (vflag_fdotr) pair_virial_fdotr_compute(this); + + if (eflag_atom) { + if (need_dup) + Kokkos::Experimental::contribute(d_eatom, dup_eatom); + k_eatom.template modify(); + k_eatom.template sync(); + } + + if (vflag_atom) { + if (need_dup) + Kokkos::Experimental::contribute(d_vatom, dup_vatom); + k_vatom.template modify(); + k_vatom.template sync(); + } + + copymode = 0; + + // free duplicated memory + if (need_dup) { + dup_rho = decltype(dup_rho)(); + dup_mu = decltype(dup_mu)(); + dup_lambda = decltype(dup_lambda)(); + dup_f = decltype(dup_f)(); + dup_eatom = decltype(dup_eatom)(); + dup_vatom = decltype(dup_vatom)(); + } +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +template +void PairADPKokkos::init_style() +{ + // convert read-in file(s) to arrays and spline them + + PairADP::init_style(); + + // adjust neighbor list request for KOKKOS + + neighflag = lmp->kokkos->neighflag; + auto request = neighbor->find_request(this); + request->set_kokkos_host(std::is_same::value && + !std::is_same::value); + request->set_kokkos_device(std::is_same::value); + if (neighflag == FULL) request->enable_full(); +} + +/* ---------------------------------------------------------------------- + convert read-in funcfl potential(s) to standard array format + interpolate all file values to a single grid and cutoff +------------------------------------------------------------------------- */ + +template +void PairADPKokkos::file2array() +{ + PairADP::file2array(); + + int i,j; + int n = atom->ntypes; + + auto k_type2frho = DAT::tdual_int_1d("pair:type2frho",n+1); + auto k_type2rhor = DAT::tdual_int_2d_dl("pair:type2rhor",n+1,n+1); + auto k_type2z2r = DAT::tdual_int_2d_dl("pair:type2z2r",n+1,n+1); + auto k_type2u2r = DAT::tdual_int_2d_dl("pair:type2u2r",n+1,n+1); + auto k_type2w2r = DAT::tdual_int_2d_dl("pair:type2w2r",n+1,n+1); + + auto h_type2frho = k_type2frho.h_view; + auto h_type2rhor = k_type2rhor.h_view; + auto h_type2z2r = k_type2z2r.h_view; + auto h_type2u2r = k_type2u2r.h_view; + auto h_type2w2r = k_type2w2r.h_view; + + for (i = 1; i <= n; i++) { + h_type2frho[i] = type2frho[i]; + for (j = 1; j <= n; j++) { + h_type2rhor(i,j) = type2rhor[i][j]; + h_type2z2r(i,j) = type2z2r[i][j]; + h_type2u2r(i,j) = type2u2r[i][j]; + h_type2w2r(i,j) = type2w2r[i][j]; + } + } + k_type2frho.template modify(); + k_type2frho.template sync(); + + k_type2rhor.template modify(); + k_type2rhor.template sync(); + + k_type2z2r.template modify(); + k_type2z2r.template sync(); + + k_type2u2r.template modify(); + k_type2u2r.template sync(); + + k_type2w2r.template modify(); + k_type2w2r.template sync(); + + + d_type2frho = k_type2frho.template view(); + d_type2rhor = k_type2rhor.template view(); + d_type2z2r = k_type2z2r.template view(); + d_type2u2r = k_type2u2r.template view(); + d_type2w2r = k_type2w2r.template view(); +} + +/* ---------------------------------------------------------------------- */ + +template +void PairADPKokkos::array2spline() +{ + rdr = 1.0/dr; + rdrho = 1.0/drho; + + tdual_ffloat_2d_n7 k_frho_spline = tdual_ffloat_2d_n7("pair:frho",nfrho,nrho+1); + tdual_ffloat_2d_n7 k_rhor_spline = tdual_ffloat_2d_n7("pair:rhor",nrhor,nr+1); + tdual_ffloat_2d_n7 k_z2r_spline = tdual_ffloat_2d_n7("pair:z2r",nz2r,nr+1); + tdual_ffloat_2d_n7 k_u2r_spline = tdual_ffloat_2d_n7("pair:z2r",nu2r,nr+1); + tdual_ffloat_2d_n7 k_w2r_spline = tdual_ffloat_2d_n7("pair:z2r",nw2r,nr+1); + + t_host_ffloat_2d_n7 h_frho_spline = k_frho_spline.h_view; + t_host_ffloat_2d_n7 h_rhor_spline = k_rhor_spline.h_view; + t_host_ffloat_2d_n7 h_z2r_spline = k_z2r_spline.h_view; + t_host_ffloat_2d_n7 h_u2r_spline = k_u2r_spline.h_view; + t_host_ffloat_2d_n7 h_w2r_spline = k_w2r_spline.h_view; + + for (int i = 0; i < nfrho; i++) + interpolate(nrho,drho,frho[i],h_frho_spline,i); + k_frho_spline.template modify(); + k_frho_spline.template sync(); + + for (int i = 0; i < nrhor; i++) + interpolate(nr,dr,rhor[i],h_rhor_spline,i); + k_rhor_spline.template modify(); + k_rhor_spline.template sync(); + + for (int i = 0; i < nz2r; i++) + interpolate(nr,dr,z2r[i],h_z2r_spline,i); + k_z2r_spline.template modify(); + k_z2r_spline.template sync(); + + for (int i = 0; i < nu2r; i++) + interpolate(nr,dr,u2r[i],h_u2r_spline,i); + k_u2r_spline.template modify(); + k_u2r_spline.template sync(); + + for (int i = 0; i < nw2r; i++) + interpolate(nr,dr,w2r[i],h_w2r_spline,i); + k_w2r_spline.template modify(); + k_w2r_spline.template sync(); + + d_frho_spline = k_frho_spline.template view(); + d_rhor_spline = k_rhor_spline.template view(); + d_z2r_spline = k_z2r_spline.template view(); + d_u2r_spline = k_u2r_spline.template view(); + d_w2r_spline = k_w2r_spline.template view(); +} + +/* ---------------------------------------------------------------------- */ + +template +void PairADPKokkos::interpolate(int n, double delta, double *f, t_host_ffloat_2d_n7 h_spline, int i) +{ + for (int m = 1; m <= n; m++) h_spline(i,m,6) = f[m]; + + h_spline(i,1,5) = h_spline(i,2,6) - h_spline(i,1,6); + h_spline(i,2,5) = 0.5 * (h_spline(i,3,6)-h_spline(i,1,6)); + h_spline(i,n-1,5) = 0.5 * (h_spline(i,n,6)-h_spline(i,n-2,6)); + h_spline(i,n,5) = h_spline(i,n,6) - h_spline(i,n-1,6); + + for (int m = 3; m <= n-2; m++) + h_spline(i,m,5) = ((h_spline(i,m-2,6)-h_spline(i,m+2,6)) + + 8.0*(h_spline(i,m+1,6)-h_spline(i,m-1,6))) / 12.0; + + for (int m = 1; m <= n-1; m++) { + h_spline(i,m,4) = 3.0*(h_spline(i,m+1,6)-h_spline(i,m,6)) - + 2.0*h_spline(i,m,5) - h_spline(i,m+1,5); + h_spline(i,m,3) = h_spline(i,m,5) + h_spline(i,m+1,5) - + 2.0*(h_spline(i,m+1,6)-h_spline(i,m,6)); + } + + h_spline(i,n,4) = 0.0; + h_spline(i,n,3) = 0.0; + + for (int m = 1; m <= n; m++) { + h_spline(i,m,2) = h_spline(i,m,5)/delta; + h_spline(i,m,1) = 2.0*h_spline(i,m,4)/delta; + h_spline(i,m,0) = 3.0*h_spline(i,m,3)/delta; + } +} + +/* ---------------------------------------------------------------------- */ + +template +int PairADPKokkos::pack_forward_comm_kokkos(int n, DAT::tdual_int_2d k_sendlist, + int iswap_in, DAT::tdual_xfloat_1d &buf, + int /*pbc_flag*/, int * /*pbc*/) +{ + d_sendlist = k_sendlist.view(); + iswap = iswap_in; + v_buf = buf.view(); + Kokkos::parallel_for(Kokkos::RangePolicy(0,n),*this); + return n*10; +} + +template +KOKKOS_INLINE_FUNCTION +void PairADPKokkos::operator()(TagPairADPPackForwardComm, const int &i) const { + int j = d_sendlist(iswap, i); + v_buf[10 * i] = d_fp(j); + v_buf[10 * i + 1] = d_mu(j, 0); + v_buf[10 * i + 2] = d_mu(j, 1); + v_buf[10 * i + 3] = d_mu(j, 2); + v_buf[10 * i + 4] = d_lambda(j, 0); + v_buf[10 * i + 5] = d_lambda(j, 1); + v_buf[10 * i + 6] = d_lambda(j, 2); + v_buf[10 * i + 7] = d_lambda(j, 3); + v_buf[10 * i + 8] = d_lambda(j, 4); + v_buf[10 * i + 9] = d_lambda(j, 5); +} + +/* ---------------------------------------------------------------------- */ + +template +void PairADPKokkos::unpack_forward_comm_kokkos(int n, int first_in, DAT::tdual_xfloat_1d &buf) +{ + first = first_in; + v_buf = buf.view(); + Kokkos::parallel_for(Kokkos::RangePolicy(0,n),*this); +} + +template +KOKKOS_INLINE_FUNCTION +void PairADPKokkos::operator()(TagPairADPUnpackForwardComm, const int &i) const { + d_fp(i + first) = v_buf[10 * i]; + d_mu(i + first, 0) = v_buf[10 * i + 1]; + d_mu(i + first, 1) = v_buf[10 * i + 2]; + d_mu(i + first, 2) = v_buf[10 * i + 3]; + d_lambda(i + first, 0) = v_buf[10 * i + 4]; + d_lambda(i + first, 1) = v_buf[10 * i + 5]; + d_lambda(i + first, 2) = v_buf[10 * i + 6]; + d_lambda(i + first, 3) = v_buf[10 * i + 7]; + d_lambda(i + first, 4) = v_buf[10 * i + 8]; + d_lambda(i + first, 5) = v_buf[10 * i + 9]; +} + +/* ---------------------------------------------------------------------- */ + +template +int PairADPKokkos::pack_forward_comm(int n, int *list, double *buf, + int /*pbc_flag*/, int * /*pbc*/) +{ + k_fp.sync_host(); + k_mu.sync_host(); + k_lambda.sync_host(); + + int i,j, m; + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + + buf[m++] = h_fp(j); + buf[m++] = h_mu(j, 0); + buf[m++] = h_mu(j, 1); + buf[m++] = h_mu(j, 2); + buf[m++] = h_lambda(j, 0); + buf[m++] = h_lambda(j, 1); + buf[m++] = h_lambda(j, 2); + buf[m++] = h_lambda(j, 3); + buf[m++] = h_lambda(j, 4); + buf[m++] = h_lambda(j, 5); + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +template +void PairADPKokkos::unpack_forward_comm(int n, int first, double *buf) +{ + k_fp.sync_host(); + k_mu.sync_host(); + k_lambda.sync_host(); + + int m, last; + m = 0; + last = n + first; + for (int i = first; i < last; i++) { + h_fp(i) = buf[m++]; + h_mu(i, 0) = buf[m++]; + h_mu(i, 1) = buf[m++]; + h_mu(i, 2) = buf[m++]; + h_lambda(i, 0) = buf[m++]; + h_lambda(i, 1) = buf[m++]; + h_lambda(i, 2) = buf[m++]; + h_lambda(i, 3) = buf[m++]; + h_lambda(i, 4) = buf[m++]; + h_lambda(i, 5) = buf[m++]; + } + + k_fp.modify_host(); + k_mu.modify_host(); + k_lambda.modify_host(); +} + +/* ---------------------------------------------------------------------- */ + +template +int PairADPKokkos::pack_reverse_comm(int n, int first, double *buf) +{ + k_rho.sync_host(); + k_mu.sync_host(); + k_lambda.sync_host(); + + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + buf[m++] = h_rho(i); + buf[m++] = h_mu(i,0); + buf[m++] = h_mu(i,1); + buf[m++] = h_mu(i,2); + buf[m++] = h_lambda(i,0); + buf[m++] = h_lambda(i,1); + buf[m++] = h_lambda(i,2); + buf[m++] = h_lambda(i,3); + buf[m++] = h_lambda(i,4); + buf[m++] = h_lambda(i,5); + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +template +void PairADPKokkos::unpack_reverse_comm(int n, int *list, double *buf) +{ + k_rho.sync_host(); + k_mu.sync_host(); + k_lambda.sync_host(); + + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + h_rho(j) += buf[m++]; + h_mu(j,0) += buf[m++]; + h_mu(j,1) += buf[m++]; + h_mu(j,2) += buf[m++]; + h_lambda(j,0) += buf[m++]; + h_lambda(j,1) += buf[m++]; + h_lambda(j,2) += buf[m++]; + h_lambda(j,3) += buf[m++]; + h_lambda(j,4) += buf[m++]; + h_lambda(j,5) += buf[m++]; + } + + k_rho.modify_host(); + k_mu.modify_host(); + k_lambda.modify_host(); +} + +/* ---------------------------------------------------------------------- */ + +template +KOKKOS_INLINE_FUNCTION +void PairADPKokkos::operator()(TagPairADPInitialize, const int &i) const { + d_rho[i] = 0.0; + d_mu(i, 0) = 0.0; + d_mu(i, 1) = 0.0; + d_mu(i, 2) = 0.0; + d_lambda(i, 0) = 0.0; + d_lambda(i, 1) = 0.0; + d_lambda(i, 2) = 0.0; + d_lambda(i, 3) = 0.0; + d_lambda(i, 4) = 0.0; + d_lambda(i, 5) = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +////Specialisation for Neighborlist types Half, HalfThread, Full +template +template +KOKKOS_INLINE_FUNCTION +void PairADPKokkos::operator()(TagPairADPKernelA, const int &ii) const { + + // rho = density at each atom + // vector mu and tensor lambda are ADP-specific + // loop over neighbors of my atoms + + // The rho array is duplicated for OpenMP, atomic for CUDA, and neither for Serial + + auto v_rho = ScatterViewHelper::value,decltype(dup_rho),decltype(ndup_rho)>::get(dup_rho,ndup_rho); + auto a_rho = v_rho.template access::value>(); + auto v_mu = ScatterViewHelper::value,decltype(dup_mu),decltype(ndup_mu)>::get(dup_mu,ndup_mu); + auto a_mu = v_mu.template access::value>(); + auto v_lambda = ScatterViewHelper::value,decltype(dup_lambda),decltype(ndup_lambda)>::get(dup_lambda,ndup_lambda); + auto a_lambda = v_lambda.template access::value>(); + + const int i = d_ilist[ii]; + const X_FLOAT xtmp = x(i,0); + const X_FLOAT ytmp = x(i,1); + const X_FLOAT ztmp = x(i,2); + const int itype = type(i); + + const int jnum = d_numneigh[i]; + + F_FLOAT rhotmp = 0.0; + F_FLOAT mutmp[3] = {0.0,0.0,0.0}; + F_FLOAT lambdatmp[6] = {0.0,0.0,0.0,0.0,0.0,0.0}; + int d_type_ji; + 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); + const int jtype = type(j); + const F_FLOAT rsq = delx*delx + dely*dely + delz*delz; + + if (rsq < cutforcesq) { + F_FLOAT p = sqrt(rsq)*rdr + 1.0; + int m = static_cast (p); + m = MIN(m,nr-1); + p -= m; + p = MIN(p,1.0); + + d_type_ji = d_type2rhor(jtype,itype); + rhotmp += ((d_rhor_spline(d_type_ji,m,3)*p + d_rhor_spline(d_type_ji,m,4))*p + + d_rhor_spline(d_type_ji,m,5))*p + d_rhor_spline(d_type_ji,m,6); + + d_type_ji = d_type2u2r(jtype,itype); + F_FLOAT u2 = ((d_u2r_spline(d_type_ji,m,3)*p + d_u2r_spline(d_type_ji,m,4))*p + + d_u2r_spline(d_type_ji,m,5))*p + d_u2r_spline(d_type_ji,m,6); + mutmp[0] += u2*delx; + mutmp[1] += u2*dely; + mutmp[2] += u2*delz; + + d_type_ji = d_type2w2r(jtype,itype); + F_FLOAT w2 = ((d_w2r_spline(d_type_ji,m,3)*p + d_w2r_spline(d_type_ji,m,4))*p + + d_w2r_spline(d_type_ji,m,5))*p + d_w2r_spline(d_type_ji,m,6); + lambdatmp[0] += w2*delx*delx; + lambdatmp[1] += w2*dely*dely; + lambdatmp[2] += w2*delz*delz; + lambdatmp[3] += w2*dely*delz; + lambdatmp[4] += w2*delx*delz; + lambdatmp[5] += w2*delx*dely; + + if (NEWTON_PAIR || j < nlocal) { + const int d_type2rhor_ij = d_type2rhor(itype,jtype); + a_rho[j] += ((d_rhor_spline(d_type2rhor_ij,m,3)*p + d_rhor_spline(d_type2rhor_ij,m,4))*p + + d_rhor_spline(d_type2rhor_ij,m,5))*p + d_rhor_spline(d_type2rhor_ij,m,6); + + const int d_type2u2r_ij = d_type2u2r(itype, jtype); + u2 = ((d_u2r_spline(d_type2u2r_ij,m,3)*p + d_u2r_spline(d_type2u2r_ij,m,4))*p + + d_u2r_spline(d_type2u2r_ij,m,5))*p + d_u2r_spline(d_type2u2r_ij,m,6); + a_mu(j,0) -= u2*delx; + a_mu(j,1) -= u2*dely; + a_mu(j,2) -= u2*delz; + + const int d_type2w2r_ij = d_type2w2r(itype, jtype); + w2 = ((d_w2r_spline(d_type2w2r_ij,m,3)*p + d_w2r_spline(d_type2w2r_ij,m,4))*p + + d_w2r_spline(d_type2w2r_ij,m,5))*p + d_w2r_spline(d_type2w2r_ij,m,6); + a_lambda(j, 0) += w2*delx*delx; + a_lambda(j, 1) += w2*dely*dely; + a_lambda(j, 2) += w2*delz*delz; + a_lambda(j, 3) += w2*dely*delz; + a_lambda(j, 4) += w2*delx*delz; + a_lambda(j, 5) += w2*delx*dely; + + } + } + + } + a_rho[i] += rhotmp; + a_mu(i, 0) += mutmp[0]; + a_mu(i, 1) += mutmp[1]; + a_mu(i, 2) += mutmp[2]; + a_lambda(i, 0) += lambdatmp[0]; + a_lambda(i, 1) += lambdatmp[1]; + a_lambda(i, 2) += lambdatmp[2]; + a_lambda(i, 3) += lambdatmp[3]; + a_lambda(i, 4) += lambdatmp[4]; + a_lambda(i, 5) += lambdatmp[5]; +} + +/* ---------------------------------------------------------------------- */ + +////Specialisation for Neighborlist types Half, HalfThread, Full +template +template +KOKKOS_INLINE_FUNCTION +void PairADPKokkos::operator()(TagPairADPKernelB, const int &ii, EV_FLOAT& ev) const { + // fp = derivative of embedding energy at each atom + // phi = embedding energy at each atom + + const int i = d_ilist[ii]; + const int itype = type(i); + + F_FLOAT p = d_rho[i]*rdrho + 1.0; + int m = static_cast (p); + m = MAX(1,MIN(m,nrho-1)); + p -= m; + p = MIN(p,1.0); + const int d_type2frho_i = d_type2frho[itype]; + d_fp[i] = (d_frho_spline(d_type2frho_i,m,0)*p + d_frho_spline(d_type2frho_i,m,1))*p + d_frho_spline(d_type2frho_i,m,2); + if (EFLAG) { + F_FLOAT phi = ((d_frho_spline(d_type2frho_i,m,3)*p + d_frho_spline(d_type2frho_i,m,4))*p + + d_frho_spline(d_type2frho_i,m,5))*p + d_frho_spline(d_type2frho_i,m,6); + phi += 0.5*(d_mu(i,0)*d_mu(i,0)+d_mu(i,1)*d_mu(i,1)+d_mu(i,2)*d_mu(i,2)); + phi += 0.5*(d_lambda(i,0)*d_lambda(i,0)+d_lambda(i,1)* + d_lambda(i,1)+d_lambda(i,2)*d_lambda(i,2)); + phi += 1.0*(d_lambda(i,3)*d_lambda(i,3)+d_lambda(i,4)* + d_lambda(i,4)+d_lambda(i,5)*d_lambda(i,5)); + phi -= 1.0/6.0*(d_lambda(i,0)+d_lambda(i,1)+d_lambda(i,2))* + (d_lambda(i,0)+d_lambda(i,1)+d_lambda(i,2)); + if (eflag_global) ev.evdwl += phi; + if (eflag_atom) d_eatom[i] += phi; + } +} + +template +template +KOKKOS_INLINE_FUNCTION +void PairADPKokkos::operator()(TagPairADPKernelB, const int &ii) const { + EV_FLOAT ev; + this->template operator()(TagPairADPKernelB(), ii, ev); +} + +/* ---------------------------------------------------------------------- */ + +////Specialisation for Neighborlist types Half, HalfThread, Full +template +template +KOKKOS_INLINE_FUNCTION +void PairADPKokkos::operator()(TagPairADPKernelAB, const int &ii, EV_FLOAT& ev) const { + + // rho = density at each atom + // loop over neighbors of my atoms + + const int i = d_ilist[ii]; + const X_FLOAT xtmp = x(i,0); + const X_FLOAT ytmp = x(i,1); + const X_FLOAT ztmp = x(i,2); + const int itype = type(i); + + const int jnum = d_numneigh[i]; + + F_FLOAT rhotmp = 0.0; + F_FLOAT mutmp[3] = {0.0,0.0,0.0}; + F_FLOAT lambdatmp[6] = {0.0,0.0,0.0,0.0,0.0,0.0}; + int d_type_ji; + + 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); + const int jtype = type(j); + const F_FLOAT rsq = delx*delx + dely*dely + delz*delz; + + if (rsq < cutforcesq) { + F_FLOAT p = sqrt(rsq)*rdr + 1.0; + int m = static_cast (p); + m = MIN(m,nr-1); + p -= m; + p = MIN(p,1.0); + d_type_ji = d_type2rhor(jtype,itype); + rhotmp += ((d_rhor_spline(d_type_ji,m,3)*p + d_rhor_spline(d_type_ji,m,4))*p + + d_rhor_spline(d_type_ji,m,5))*p + d_rhor_spline(d_type_ji,m,6); + + d_type_ji = d_type2u2r(jtype,itype); + F_FLOAT u2 = ((d_u2r_spline(d_type_ji,m,3)*p + d_u2r_spline(d_type_ji,m,4))*p + + d_u2r_spline(d_type_ji,m,5))*p + d_u2r_spline(d_type_ji,m,6); + mutmp[0] += u2*delx; + mutmp[1] += u2*dely; + mutmp[2] += u2*delz; + + d_type_ji = d_type2w2r(jtype,itype); + F_FLOAT w2 = ((d_w2r_spline(d_type_ji,m,3)*p + d_w2r_spline(d_type_ji,m,4))*p + + d_w2r_spline(d_type_ji,m,5))*p + d_w2r_spline(d_type_ji,m,6); + lambdatmp[0] += w2*delx*delx; + lambdatmp[1] += w2*dely*dely; + lambdatmp[2] += w2*delz*delz; + lambdatmp[3] += w2*dely*delz; + lambdatmp[4] += w2*delx*delz; + lambdatmp[5] += w2*delx*dely; + } + + } + d_rho[i] += rhotmp; + + d_mu(i, 0) += mutmp[0]; + d_mu(i, 1) += mutmp[1]; + d_mu(i, 2) += mutmp[2]; + + d_lambda(i, 0) += lambdatmp[0]; + d_lambda(i, 1) += lambdatmp[1]; + d_lambda(i, 2) += lambdatmp[2]; + d_lambda(i, 3) += lambdatmp[3]; + d_lambda(i, 4) += lambdatmp[4]; + d_lambda(i, 5) += lambdatmp[5]; + + // fp = derivative of embedding energy at each atom + // phi = embedding energy at each atom + + F_FLOAT p = d_rho[i]*rdrho + 1.0; + int m = static_cast (p); + m = MAX(1,MIN(m,nrho-1)); + p -= m; + p = MIN(p,1.0); + const int d_type2frho_i = d_type2frho[itype]; + d_fp[i] = (d_frho_spline(d_type2frho_i,m,0)*p + d_frho_spline(d_type2frho_i,m,1))*p + d_frho_spline(d_type2frho_i,m,2); + if (EFLAG) { + F_FLOAT phi = ((d_frho_spline(d_type2frho_i,m,3)*p + d_frho_spline(d_type2frho_i,m,4))*p + + d_frho_spline(d_type2frho_i,m,5))*p + d_frho_spline(d_type2frho_i,m,6); + + phi += 0.5*(d_mu(i,0)*d_mu(i,0)+d_mu(i,1)*d_mu(i,1)+d_mu(i,2)*d_mu(i,2)); + phi += 0.5*(d_lambda(i,0)*d_lambda(i,0)+d_lambda(i,1)* + d_lambda(i,1)+d_lambda(i,2)*d_lambda(i,2)); + phi += 1.0*(d_lambda(i,3)*d_lambda(i,3)+d_lambda(i,4)* + d_lambda(i,4)+d_lambda(i,5)*d_lambda(i,5)); + phi -= 1.0/6.0*(d_lambda(i,0)+d_lambda(i,1)+d_lambda(i,2))* + (d_lambda(i,0)+d_lambda(i,1)+d_lambda(i,2)); + if (eflag_global) ev.evdwl += phi; + if (eflag_atom) d_eatom[i] += phi; + } + +} + +template +template +KOKKOS_INLINE_FUNCTION +void PairADPKokkos::operator()(TagPairADPKernelAB, const int &ii) const { + EV_FLOAT ev; + this->template operator()(TagPairADPKernelAB(), ii, ev); +} + +/* ---------------------------------------------------------------------- */ + +////Specialisation for Neighborlist types Half, HalfThread, Full +template +template +KOKKOS_INLINE_FUNCTION +void PairADPKokkos::operator()(TagPairADPKernelC, const int &ii, EV_FLOAT& ev) const { + + // The f array is duplicated for OpenMP, atomic for CUDA, and neither for Serial + + auto v_f = ScatterViewHelper::value,decltype(dup_f),decltype(ndup_f)>::get(dup_f,ndup_f); + auto a_f = v_f.template access::value>(); + + const int i = d_ilist[ii]; + const X_FLOAT xtmp = x(i,0); + const X_FLOAT ytmp = x(i,1); + const X_FLOAT ztmp = x(i,2); + const int itype = type(i); + + const int jnum = d_numneigh[i]; + + F_FLOAT fxtmp = 0.0; + F_FLOAT fytmp = 0.0; + F_FLOAT fztmp = 0.0; + + 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); + const int jtype = type(j); + const F_FLOAT rsq = delx*delx + dely*dely + delz*delz; + + if (rsq < cutforcesq) { + const F_FLOAT r = sqrt(rsq); + F_FLOAT p = r*rdr + 1.0; + int m = static_cast (p); + m = MIN(m,nr-1); + p -= m; + p = MIN(p,1.0); + + // rhoip = derivative of (density at atom j due to atom i) + // rhojp = derivative of (density at atom i due to atom j) + // phi = pair potential energy + // phip = phi' + // z2 = phi * r + // z2p = (phi * r)' = (phi' r) + phi + // u2 = u + // u2p = u' + // w2 = w + // w2p = w' + // psip needs both fp[i] and fp[j] terms since r_ij appears in two + // terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji) + // hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip + + const int d_type2rhor_ij = d_type2rhor(itype,jtype); + const F_FLOAT rhoip = (d_rhor_spline(d_type2rhor_ij,m,0)*p + d_rhor_spline(d_type2rhor_ij,m,1))*p + + d_rhor_spline(d_type2rhor_ij,m,2); + const int d_type2rhor_ji = d_type2rhor(jtype,itype); + const F_FLOAT rhojp = (d_rhor_spline(d_type2rhor_ji,m,0)*p + d_rhor_spline(d_type2rhor_ji,m,1))*p + + d_rhor_spline(d_type2rhor_ji,m,2); + const int d_type2z2r_ij = d_type2z2r(itype,jtype); + const F_FLOAT z2p = (d_z2r_spline(d_type2z2r_ij,m,0)*p+d_z2r_spline(d_type2z2r_ij,m,1))*p + d_z2r_spline(d_type2z2r_ij,m,2); + const F_FLOAT z2 = ((d_z2r_spline(d_type2z2r_ij,m,3)*p + d_z2r_spline(d_type2z2r_ij,m,4))*p + + d_z2r_spline(d_type2z2r_ij,m,5))*p+d_z2r_spline(d_type2z2r_ij,m,6); + + const int d_type2u2r_ij = d_type2u2r(itype,jtype); + const F_FLOAT u2p = (d_u2r_spline(d_type2u2r_ij,m,0)*p + d_u2r_spline(d_type2u2r_ij,m,1))*p + + d_u2r_spline(d_type2u2r_ij,m,2); + + const F_FLOAT u2 = ((d_u2r_spline(d_type2u2r_ij,m,3)*p + d_u2r_spline(d_type2u2r_ij,m,4))*p + + d_u2r_spline(d_type2u2r_ij,m,5))*p + d_u2r_spline(d_type2u2r_ij,m,6); + + + const int d_type2w2r_ij = d_type2w2r(itype,jtype); + const F_FLOAT w2p = (d_w2r_spline(d_type2w2r_ij,m,0)*p + d_w2r_spline(d_type2w2r_ij,m,1))*p + + d_w2r_spline(d_type2w2r_ij,m,2); + + const F_FLOAT w2 = ((d_w2r_spline(d_type2w2r_ij,m,3)*p + d_w2r_spline(d_type2w2r_ij,m,4))*p + + d_w2r_spline(d_type2w2r_ij,m,5))*p + d_w2r_spline(d_type2w2r_ij,m,6); + + + + const F_FLOAT recip = 1.0/r; + const F_FLOAT phi = z2*recip; + const F_FLOAT phip = z2p*recip - phi*recip; + const F_FLOAT psip = d_fp[i]*rhojp + d_fp[j]*rhoip + phip; + const F_FLOAT fpair = -psip*recip; + + const F_FLOAT delmux = d_mu(i, 0)-d_mu(j,0); + const F_FLOAT delmuy = d_mu(i, 1)-d_mu(j,1); + const F_FLOAT delmuz = d_mu(i, 2)-d_mu(j,2); + const F_FLOAT trdelmu = delmux*delx+delmuy*dely+delmuz*delz; + const F_FLOAT sumlamxx = d_lambda(i,0)+d_lambda(j,0); + const F_FLOAT sumlamyy = d_lambda(i,1)+d_lambda(j,1); + const F_FLOAT sumlamzz = d_lambda(i,2)+d_lambda(j,2); + const F_FLOAT sumlamyz = d_lambda(i,3)+d_lambda(j,3); + const F_FLOAT sumlamxz = d_lambda(i,4)+d_lambda(j,4); + const F_FLOAT sumlamxy = d_lambda(i,5)+d_lambda(j,5); + + const F_FLOAT tradellam = sumlamxx*delx*delx+sumlamyy*dely*dely+ + sumlamzz*delz*delz+2.0*sumlamxy*delx*dely+ + 2.0*sumlamxz*delx*delz+2.0*sumlamyz*dely*delz; + const F_FLOAT nu = sumlamxx+sumlamyy+sumlamzz; + + const F_FLOAT adpx = -1.0*(delmux*u2 + trdelmu*u2p*delx*recip + + 2.0*w2*(sumlamxx*delx+sumlamxy*dely+sumlamxz*delz) + + w2p*delx*recip*tradellam - 1.0/3.0*nu*(w2p*r+2.0*w2)*delx); + const F_FLOAT adpy = -1.0*(delmuy*u2 + trdelmu*u2p*dely*recip + + 2.0*w2*(sumlamxy*delx+sumlamyy*dely+sumlamyz*delz) + + w2p*dely*recip*tradellam - 1.0/3.0*nu*(w2p*r+2.0*w2)*dely); + const F_FLOAT adpz = -1.0*(delmuz*u2 + trdelmu*u2p*delz*recip + + 2.0*w2*(sumlamxz*delx+sumlamyz*dely+sumlamzz*delz) + + w2p*delz*recip*tradellam - 1.0/3.0*nu*(w2p*r+2.0*w2)*delz); + + F_FLOAT fx = delx*fpair + adpx; + F_FLOAT fy = dely*fpair + adpy; + F_FLOAT fz = delz*fpair + adpz; + + fxtmp += fx; + fytmp += fy; + fztmp += fz; + + if ((NEIGHFLAG==HALF || NEIGHFLAG==HALFTHREAD) && (NEWTON_PAIR || j < nlocal)) { + a_f(j,0) -= fx; + a_f(j,1) -= fy; + a_f(j,2) -= fz; + } + + if (EVFLAG) { + if (eflag) { + ev.evdwl += (((NEIGHFLAG==HALF || NEIGHFLAG==HALFTHREAD)&&(NEWTON_PAIR||(jtemplate ev_tally_xyz(ev,i,j,phi,fx,fy,fz,delx,dely,delz); + } + + } + } + + a_f(i,0) += fxtmp; + a_f(i,1) += fytmp; + a_f(i,2) += fztmp; +} + +template +template +KOKKOS_INLINE_FUNCTION +void PairADPKokkos::operator()(TagPairADPKernelC, const int &ii) const { + EV_FLOAT ev; + this->template operator()(TagPairADPKernelC(), ii, ev); +} + +/* ---------------------------------------------------------------------- */ + +template +template +KOKKOS_INLINE_FUNCTION +void PairADPKokkos::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 +{ + const int EFLAG = eflag; + const int VFLAG = vflag_either; + + // The eatom and vatom arrays are duplicated for OpenMP, atomic for CUDA, and neither for Serial + + auto v_eatom = ScatterViewHelper,decltype(dup_eatom),decltype(ndup_eatom)>::get(dup_eatom,ndup_eatom); + auto a_eatom = v_eatom.template access>(); + + auto v_vatom = ScatterViewHelper,decltype(dup_vatom),decltype(ndup_vatom)>::get(dup_vatom,ndup_vatom); + auto a_vatom = v_vatom.template access>(); + + if (EFLAG) { + if (eflag_atom) { + const E_FLOAT epairhalf = 0.5 * epair; + if (NEIGHFLAG!=FULL) { + if (NEWTON_PAIR || i < nlocal) a_eatom[i] += epairhalf; + if (NEWTON_PAIR || j < nlocal) a_eatom[j] += epairhalf; + } else { + a_eatom[i] += epairhalf; + } + } + } + + if (VFLAG) { + const E_FLOAT v0 = delx*fx; + const E_FLOAT v1 = dely*fy; + const E_FLOAT v2 = delz*fz; + const E_FLOAT v3 = delx*fy; + const E_FLOAT v4 = delx*fz; + const E_FLOAT v5 = dely*fz; + + if (vflag_global) { + if (NEIGHFLAG!=FULL) { + if (NEWTON_PAIR || i < nlocal) { + ev.v[0] += 0.5*v0; + ev.v[1] += 0.5*v1; + ev.v[2] += 0.5*v2; + ev.v[3] += 0.5*v3; + ev.v[4] += 0.5*v4; + ev.v[5] += 0.5*v5; + } + if (NEWTON_PAIR || j < nlocal) { + ev.v[0] += 0.5*v0; + ev.v[1] += 0.5*v1; + ev.v[2] += 0.5*v2; + ev.v[3] += 0.5*v3; + ev.v[4] += 0.5*v4; + ev.v[5] += 0.5*v5; + } + } else { + ev.v[0] += 0.5*v0; + ev.v[1] += 0.5*v1; + ev.v[2] += 0.5*v2; + ev.v[3] += 0.5*v3; + ev.v[4] += 0.5*v4; + ev.v[5] += 0.5*v5; + } + } + + if (vflag_atom) { + if (NEIGHFLAG!=FULL) { + if (NEWTON_PAIR || i < nlocal) { + a_vatom(i,0) += 0.5*v0; + a_vatom(i,1) += 0.5*v1; + a_vatom(i,2) += 0.5*v2; + a_vatom(i,3) += 0.5*v3; + a_vatom(i,4) += 0.5*v4; + a_vatom(i,5) += 0.5*v5; + } + if (NEWTON_PAIR || j < nlocal) { + a_vatom(j,0) += 0.5*v0; + a_vatom(j,1) += 0.5*v1; + a_vatom(j,2) += 0.5*v2; + a_vatom(j,3) += 0.5*v3; + a_vatom(j,4) += 0.5*v4; + a_vatom(j,5) += 0.5*v5; + } + } else { + a_vatom(i,0) += 0.5*v0; + a_vatom(i,1) += 0.5*v1; + a_vatom(i,2) += 0.5*v2; + a_vatom(i,3) += 0.5*v3; + a_vatom(i,4) += 0.5*v4; + a_vatom(i,5) += 0.5*v5; + } + } + } +} + +namespace LAMMPS_NS { +template class PairADPKokkos; +#ifdef LMP_KOKKOS_GPU +template class PairADPKokkos; +#endif +} diff --git a/src/KOKKOS/pair_adp_kokkos.h b/src/KOKKOS/pair_adp_kokkos.h new file mode 100644 index 0000000000..ce163f6120 --- /dev/null +++ b/src/KOKKOS/pair_adp_kokkos.h @@ -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); +PairStyle(adp/kk/device,PairADPKokkos); +PairStyle(adp/kk/host,PairADPKokkos); +// 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 +struct TagPairADPKernelA{}; + +template +struct TagPairADPKernelB{}; + +template +struct TagPairADPKernelAB{}; + +template +struct TagPairADPKernelC{}; + +template +class PairADPKokkos : public PairADP, public KokkosBase +{ + public: + enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF}; + enum {COUL_FLAG=0}; + typedef DeviceType device_type; + typedef ArrayTypes 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 + KOKKOS_INLINE_FUNCTION + void operator()(TagPairADPKernelA, const int&) const; + + template + KOKKOS_INLINE_FUNCTION + void operator()(TagPairADPKernelB, const int&, EV_FLOAT&) const; + + template + KOKKOS_INLINE_FUNCTION + void operator()(TagPairADPKernelB, const int&) const; + + template + KOKKOS_INLINE_FUNCTION + void operator()(TagPairADPKernelAB, const int&, EV_FLOAT&) const; + + template + KOKKOS_INLINE_FUNCTION + void operator()(TagPairADPKernelAB, const int&) const; + + template + KOKKOS_INLINE_FUNCTION + void operator()(TagPairADPKernelC, const int&, EV_FLOAT&) const; + + template + KOKKOS_INLINE_FUNCTION + void operator()(TagPairADPKernelC, const int&) const; + + template + 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::value; + + template + using DupScatterView = KKScatterView; + + template + using NonDupScatterView = KKScatterView; + + DupScatterView dup_rho; + DupScatterView dup_mu; + DupScatterView dup_lambda; + DupScatterView dup_f; + DupScatterView dup_eatom; + DupScatterView dup_vatom; + NonDupScatterView ndup_rho; + NonDupScatterView ndup_mu; + NonDupScatterView ndup_lambda; + NonDupScatterView ndup_f; + NonDupScatterView ndup_eatom; + NonDupScatterView 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 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*); +}; + +} +#endif +#endif + diff --git a/src/KOKKOS/pair_eam_kokkos.cpp b/src/KOKKOS/pair_eam_kokkos.cpp index ac3b0d3f07..26b3ec6880 100644 --- a/src/KOKKOS/pair_eam_kokkos.cpp +++ b/src/KOKKOS/pair_eam_kokkos.cpp @@ -54,10 +54,10 @@ PairEAMKokkos::PairEAMKokkos(LAMMPS *lmp) : PairEAM(lmp) template PairEAMKokkos::~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::operator()(TagPairEAMKernelAB, 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); diff --git a/src/KOKKOS/pair_eam_kokkos.h b/src/KOKKOS/pair_eam_kokkos.h index b2ab3a6ac2..4bc7a55196 100644 --- a/src/KOKKOS/pair_eam_kokkos.h +++ b/src/KOKKOS/pair_eam_kokkos.h @@ -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 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 diff --git a/src/MANYBODY/pair_adp.cpp b/src/MANYBODY/pair_adp.cpp index 819a57efa2..b07f3e602d 100644 --- a/src/MANYBODY/pair_adp.cpp +++ b/src/MANYBODY/pair_adp.cpp @@ -76,6 +76,8 @@ PairADP::PairADP(LAMMPS *lmp) : Pair(lmp) PairADP::~PairADP() { + if (copymode) return; + memory->destroy(rho); memory->destroy(fp); memory->destroy(mu); diff --git a/src/MANYBODY/pair_adp.h b/src/MANYBODY/pair_adp.h index 416ad26fc3..a82028687c 100644 --- a/src/MANYBODY/pair_adp.h +++ b/src/MANYBODY/pair_adp.h @@ -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 diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp index 018f658c61..c8426673ab 100644 --- a/src/compute_reduce.cpp +++ b/src/compute_reduce.cpp @@ -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; diff --git a/src/compute_reduce.h b/src/compute_reduce.h index a7590ccd52..dc4ee1ef2c 100644 --- a/src/compute_reduce.h +++ b/src/compute_reduce.h @@ -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 **); diff --git a/src/delete_atoms.cpp b/src/delete_atoms.cpp index 8524d44b07..3646da36d4 100644 --- a/src/delete_atoms.cpp +++ b/src/delete_atoms.cpp @@ -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 #include @@ -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( atom->style_match("ellipsoid")); - auto avec_line = dynamic_cast( atom->style_match("line")); - auto avec_tri = dynamic_cast( atom->style_match("tri")); - auto avec_body = dynamic_cast( atom->style_match("body")); + auto avec_ellipsoid = dynamic_cast(atom->style_match("ellipsoid")); + auto avec_line = dynamic_cast(atom->style_match("line")); + auto avec_tri = dynamic_cast(atom->style_match("tri")); + auto avec_body = dynamic_cast(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(); + hash = new std::map(); // 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(); + hash = new std::map(); 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::iterator pos; + std::map::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 *hash = daptr->hash; + std::map *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 *hash = daptr->hash; + std::map *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]); } } diff --git a/src/delete_atoms.h b/src/delete_atoms.h index 33333a1ff4..633cb5a9a0 100644 --- a/src/delete_atoms.h +++ b/src/delete_atoms.h @@ -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(); diff --git a/unittest/commands/CMakeLists.txt b/unittest/commands/CMakeLists.txt index f6ef263f5f..deeee00f71 100644 --- a/unittest/commands/CMakeLists.txt +++ b/unittest/commands/CMakeLists.txt @@ -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) diff --git a/unittest/commands/test_delete_atoms.cpp b/unittest/commands/test_delete_atoms.cpp new file mode 100644 index 0000000000..d7c37d7311 --- /dev/null +++ b/unittest/commands/test_delete_atoms.cpp @@ -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 +#include + +// 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 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; +}