From 5003c359637c94f81c1ea36caba9777c2aba7d5f Mon Sep 17 00:00:00 2001 From: Matt Bettencourt Date: Wed, 30 Mar 2022 16:22:26 +0200 Subject: [PATCH] have first cut of just pair-dpd --- src/DPD-BASIC/pair_dpd.cpp | 2 + src/DPD-BASIC/pair_dpd.h | 8 +- src/KOKKOS/Install.sh | 2 + src/KOKKOS/pair_dpd_kokkos.cpp | 451 +++++++++++++++++++++++++++++++++ src/KOKKOS/pair_dpd_kokkos.h | 126 +++++++++ 5 files changed, 585 insertions(+), 4 deletions(-) create mode 100644 src/KOKKOS/pair_dpd_kokkos.cpp create mode 100644 src/KOKKOS/pair_dpd_kokkos.h diff --git a/src/DPD-BASIC/pair_dpd.cpp b/src/DPD-BASIC/pair_dpd.cpp index caa8161573..e988f1521e 100644 --- a/src/DPD-BASIC/pair_dpd.cpp +++ b/src/DPD-BASIC/pair_dpd.cpp @@ -46,6 +46,8 @@ PairDPD::PairDPD(LAMMPS *lmp) : Pair(lmp) PairDPD::~PairDPD() { + if (copymode) return; + if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); diff --git a/src/DPD-BASIC/pair_dpd.h b/src/DPD-BASIC/pair_dpd.h index 796228878c..7a5f6708d7 100644 --- a/src/DPD-BASIC/pair_dpd.h +++ b/src/DPD-BASIC/pair_dpd.h @@ -28,11 +28,11 @@ class PairDPD : public Pair { public: PairDPD(class LAMMPS *); ~PairDPD() override; - void compute(int, int) override; + virtual void compute(int, int) override; void settings(int, char **) override; void coeff(int, char **) override; - void init_style() override; - double init_one(int, int) override; + virtual void init_style() override; + virtual double init_one(int, int) override; void write_restart(FILE *) override; void read_restart(FILE *) override; void write_restart_settings(FILE *) override; @@ -49,7 +49,7 @@ class PairDPD : public Pair { double **sigma; class RanMars *random; - void allocate(); + virtual void allocate(); }; } // namespace LAMMPS_NS diff --git a/src/KOKKOS/Install.sh b/src/KOKKOS/Install.sh index da6bd8f1e3..8c4a3ea9ae 100755 --- a/src/KOKKOS/Install.sh +++ b/src/KOKKOS/Install.sh @@ -231,6 +231,8 @@ action pair_coul_long_kokkos.h pair_coul_long.h action pair_coul_wolf_kokkos.cpp action pair_coul_wolf_kokkos.h action pair_dpd_ext_kokkos.cpp pair_dpd_ext.cpp +action pair_dpd_kokkos.h pair_dpd.h +action pair_dpd_kokkos.cpp pair_dpd.cpp action pair_dpd_ext_kokkos.h pair_dpd_ext.h action pair_dpd_fdt_energy_kokkos.cpp pair_dpd_fdt_energy.cpp action pair_dpd_fdt_energy_kokkos.h pair_dpd_fdt_energy.h diff --git a/src/KOKKOS/pair_dpd_kokkos.cpp b/src/KOKKOS/pair_dpd_kokkos.cpp new file mode 100644 index 0000000000..190acac98e --- /dev/null +++ b/src/KOKKOS/pair_dpd_kokkos.cpp @@ -0,0 +1,451 @@ +// 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 author: Matt Bettencourt (NVIDIA) +------------------------------------------------------------------------- */ + +#include "pair_dpd_kokkos.h" + +#include "atom.h" +#include "atom_kokkos.h" +#include "memory_kokkos.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "neighbor.h" +#include "random_mars.h" +#include "update.h" +#include "atom_masks.h" +#include "kokkos.h" + +#include + +using namespace LAMMPS_NS; + +#define EPSILON 1.0e-10 + + +template +PairDPDKokkos::PairDPDKokkos(class LAMMPS *lmp) : + PairDPD(lmp) , +#ifdef DPD_USE_RAN_MARS + rand_pool(0 /* unused */, lmp) +#else + rand_pool() +#endif +{ + kokkosable = 1; + atomKK = (AtomKokkos *) atom; + execution_space = ExecutionSpaceFromDevice::space; + + datamask_read = X_MASK | V_MASK | F_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK; + datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK; +} + +/* ---------------------------------------------------------------------- */ + +template +PairDPDKokkos::~PairDPDKokkos() { + if (copymode) return; + +#ifdef DPD_USE_RAN_MARS + rand_pool.destroy(); +#endif + + memoryKK->destroy_kokkos(k_eatom,eatom); + memoryKK->destroy_kokkos(k_vatom,vatom); + + memoryKK->destroy_kokkos(k_cutsq,cutsq); +} + +/* ---------------------------------------------------------------------- */ + +template +void PairDPDKokkos::init_style() +{ + PairDPD::init_style(); + +#ifdef DPD_USE_RAN_MARS + rand_pool.init(random,seed); +#else + typedef Kokkos::Experimental::UniqueToken< + DeviceType, Kokkos::Experimental::UniqueTokenScope::Global> unique_token_type; + unique_token_type unique_token; + rand_pool.init(seed + comm->me,unique_token.size()); +#endif + + neighflag = lmp->kokkos->neighflag; + + if (force->newton_pair == 0 || neighflag == FULL ) + error->all(FLERR,"Must use half neighbor list style and newton on with pair dpd/kk"); + + 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(); +} + +/* ---------------------------------------------------------------------- */ + +template +void PairDPDKokkos::compute(int eflagin, int vflagin) +{ + eflag = eflagin; vflag = vflagin; + if (neighflag == FULL) no_virial_fdotr_compute = 1; + + ev_init(eflag,vflag); + + if (eflag_atom) { + memoryKK->destroy_kokkos(k_eatom,eatom); + memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom"); + d_eatom = k_eatom.template view(); + } + if (vflag_atom) { + memoryKK->destroy_kokkos(k_vatom,vatom); + memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom"); + d_vatom = k_vatom.template view(); + } + + atomKK->sync(execution_space,X_MASK | V_MASK | F_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK); + + x = atomKK->k_x.view(); + v = atomKK->k_v.view(); + f = atomKK->k_f.view(); + type = atomKK->k_type.view(); + + k_cutsq.template sync(); + k_params.template sync(); + + special_lj[0] = force->special_lj[0]; + special_lj[1] = force->special_lj[1]; + special_lj[2] = force->special_lj[2]; + special_lj[3] = force->special_lj[3]; + + nlocal = atom->nlocal; + newton_pair = force->newton_pair; + dtinvsqrt = 1.0/sqrt(update->dt); + + NeighListKokkos* k_list = static_cast*>(list); + d_numneigh = k_list->d_numneigh; + d_neighbors = k_list->d_neighbors; + d_ilist = k_list->d_ilist; + + // loop over neighbors of my atoms + + int inum = list->inum; + EV_FLOAT ev; + copymode = 1; + if (neighflag == HALF) { + if (newton_pair) { + if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + else Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } else { + if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + else Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } + } else if (neighflag == HALFTHREAD) { + if (newton_pair) { + if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + else Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } else { + if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + else Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } + } else if (neighflag == FULL) { + if (newton_pair) { + if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + else Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } else { + if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + else Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + } + } + + 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) { + k_eatom.template modify(); + k_eatom.template sync(); + } + + if (vflag_atom) { + k_vatom.template modify(); + k_vatom.template sync(); + } + + copymode = 0; + + if (evflag) atomKK->modified(execution_space,F_MASK | ENERGY_MASK | VIRIAL_MASK); + else atomKK->modified(execution_space,F_MASK); +} + +/* ---------------------------------------------------------------------- */ + +template +template +KOKKOS_INLINE_FUNCTION +void PairDPDKokkos::operator() (TagDPDKokkos, const int &ii) const { + EV_FLOAT ev; + this->template operator()(TagDPDKokkos(), ii, ev); +} +template +template +KOKKOS_INLINE_FUNCTION +void PairDPDKokkos::operator() (TagDPDKokkos, const int &ii, EV_FLOAT &ev) const { + + Kokkos::View::value,Kokkos::MemoryTraits::value> > a_f = f; + + int i,j,jj,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,fpair; + double vxtmp,vytmp,vztmp,delvx,delvy,delvz; + double rsq,r,rinv,dot,wd,randnum,factor_dpd; + double fx = 0,fy = 0,fz = 0; + double evdwl = 0; + i = d_ilist[ii]; + xtmp = x(i,0); + ytmp = x(i,1); + ztmp = x(i,2); + vxtmp = v(i,0); + vytmp = v(i,1); + vztmp = v(i,2); + itype = type(i); + jnum = d_numneigh[i]; + rand_type rand_gen = rand_pool.get_state(); + for (jj = 0; jj < jnum; jj++) { + j = d_neighbors(i,jj); + factor_dpd = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x(j,0); + dely = ytmp - x(j,1); + delz = ztmp - x(j,2); + rsq = delx*delx + dely*dely + delz*delz; + jtype = type(j); + if (rsq < d_cutsq(itype,jtype)) { + r = sqrt(rsq); + if (r < EPSILON) continue; // r can be 0.0 in DPD systems + rinv = 1.0/r; + delvx = vxtmp - v(j,0); + delvy = vytmp - v(j,1); + delvz = vztmp - v(j,2); + dot = delx*delvx + dely*delvy + delz*delvz; + + wd = 1.0 - r/params(itype,jtype).cut; + + randnum = rand_gen.normal(); + + // conservative force + fpair = params(itype,jtype).a0*wd; + + // drag force - parallel + fpair -= params(itype,jtype).gamma*wd*wd*dot*rinv; + + // random force - parallel + fpair += params(itype,jtype).sigma*wd*randnum*dtinvsqrt; + fpair *= factor_dpd*rinv; + + fx += fpair*delx; + fy += fpair*dely; + fz += fpair*delz; + + if ((neighflag==HALF || neighflag==HALFTHREAD) && (NEWTON_PAIR || j < nlocal) ) { + a_f(j,0) -= fpair*delx; + a_f(j,1) -= fpair*dely; + a_f(j,2) -= fpair*delz; + } + + if (eflag) { + // unshifted eng of conservative term: + // evdwl = -a0[itype][jtype]*r * (1.0-0.5*r/cut[itype][jtype]); + // eng shifted to 0.0 at cutoff + evdwl = 0.5*params(itype,jtype).a0*params(itype,jtype).cut* wd*wd; + evdwl *= factor_dpd; + if (EVFLAG) + ev.evdwl += (((NEIGHFLAG==HALF || NEIGHFLAG==HALFTHREAD) && (NEWTON_PAIR||(jtemplate ev_tally(ev,i,j,evdwl,fpair,delx,dely,delz); + } + } + a_f(i,0) += fx; + a_f(i,1) += fy; + a_f(i,2) += fz; + rand_pool.free_state(rand_gen); +} + +/* ---------------------------------------------------------------------- */ + +template +template +KOKKOS_INLINE_FUNCTION +void PairDPDKokkos::ev_tally(EV_FLOAT &ev, const int &i, const int &j, + const F_FLOAT &epair, const F_FLOAT &fpair, 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 atomic for Half/Thread neighbor style + Kokkos::View::value,Kokkos::MemoryTraits::value> > a_eatom = k_eatom.view(); + Kokkos::View::value,Kokkos::MemoryTraits::value> > a_vatom = k_vatom.view(); + + 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*delx*fpair; + const E_FLOAT v1 = dely*dely*fpair; + const E_FLOAT v2 = delz*delz*fpair; + const E_FLOAT v3 = delx*dely*fpair; + const E_FLOAT v4 = delx*delz*fpair; + const E_FLOAT v5 = dely*delz*fpair; + + 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; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +template +void PairDPDKokkos::allocate() +{ + PairDPD::allocate(); + int n = atom->ntypes; + + memory->destroy(cutsq); + memoryKK->create_kokkos(k_cutsq,cutsq,n+1,n+1,"pair:cutsq"); + d_cutsq = k_cutsq.template view(); + + k_params = Kokkos::DualView("PairDPD::params",n+1,n+1); + params = k_params.template view(); +} + +/* ---------------------------------------------------------------------- */ + +template +KOKKOS_INLINE_FUNCTION +int PairDPDKokkos::sbmask(const int& j) const { + return j >> SBBITS & 3; +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +template +double PairDPDKokkos::init_one(int i, int j) +{ + double cutone = PairDPD::init_one(i,j); + + k_params.h_view(i,j).cut = cut[i][j]; + k_params.h_view(i,j).a0 = a0[i][j]; + k_params.h_view(i,j).gamma = gamma[i][j]; + k_params.h_view(i,j).sigma = sigma[i][j]; + k_params.h_view(j,i) = k_params.h_view(i,j); + + k_params.template modify(); + + k_cutsq.h_view(i,j) = cutone*cutone; + k_cutsq.h_view(j,i) = k_cutsq.h_view(i,j); + k_cutsq.template modify(); + + return cutone; +} + +namespace LAMMPS_NS { +template class PairDPDKokkos; +#ifdef LMP_KOKKOS_GPU +template class PairDPDKokkos; +#endif +} diff --git a/src/KOKKOS/pair_dpd_kokkos.h b/src/KOKKOS/pair_dpd_kokkos.h new file mode 100644 index 0000000000..da2558a6fb --- /dev/null +++ b/src/KOKKOS/pair_dpd_kokkos.h @@ -0,0 +1,126 @@ +/* -*- 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. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(dpd/kk,PairDPDKokkos); +PairStyle(dpd/kk/device,PairDPDKokkos); +PairStyle(dpd/kk/host,PairDPDKokkos); +// clang-format on +#else + +#ifndef LMP_PAIR_DPD_KOKKOS_H +#define LMP_PAIR_DPD_KOKKOS_H + +#include "pair_dpd.h" +#include "pair_kokkos.h" +#include "kokkos_type.h" + +#if !defined(DPD_USE_RAN_MARS) && !defined(DPD_USE_Random_XorShift64) && !defined(Random_XorShift1024) +#define DPD_USE_Random_XorShift64 +#endif + +#ifdef DPD_USE_RAN_MARS +#include "rand_pool_wrap_kokkos.h" +#else +#include "Kokkos_Random.hpp" +#endif + +namespace LAMMPS_NS { + +template +class PairDPDKokkos : public PairDPD { + public: + typedef DeviceType device_type; + typedef ArrayTypes AT; + typedef EV_FLOAT value_type; + + PairDPDKokkos(class LAMMPS*); + ~PairDPDKokkos() override; + + void allocate() override; + + void init_style() override; + double init_one(int i, int j) override; + void compute(int, int) override; + + struct params_dpd { + KOKKOS_INLINE_FUNCTION + params_dpd() {cut=a0=gamma=sigma=0;} + KOKKOS_INLINE_FUNCTION + params_dpd(int /*i*/) {cut=a0=gamma=sigma=0;} + F_FLOAT cutsq,cut,a0,gamma,sigma; + }; + + template + struct TagDPDKokkos{}; + + template + KOKKOS_INLINE_FUNCTION + void operator () (TagDPDKokkos, const int &i) const; + + template + KOKKOS_INLINE_FUNCTION + void operator () (TagDPDKokkos, const int &i, EV_FLOAT&) const; + + template + KOKKOS_INLINE_FUNCTION + void ev_tally(EV_FLOAT &ev, const int &i, const int &j, + const F_FLOAT &epair, const F_FLOAT &fpair, const F_FLOAT &delx, + const F_FLOAT &dely, const F_FLOAT &delz) const; + private: + + double special_lj[4]; + int eflag,vflag; + int neighflag, nlocal,newton_pair; + double dtinvsqrt; + +#ifdef DPD_USE_RAN_MARS + RandPoolWrap rand_pool; + typedef RandWrap rand_type; +#elif defined(DPD_USE_Random_XorShift64) + Kokkos::Random_XorShift64_Pool rand_pool; + typedef typename Kokkos::Random_XorShift64_Pool::generator_type rand_type; +#elif defined(DPD_USE_Random_XorShift1024) + Kokkos::Random_XorShift1024_Pool rand_pool; + typedef typename Kokkos::Random_XorShift1024_Pool::generator_type rand_type; +#endif + typename AT::t_x_array_randomread x; + typename AT::t_x_array_randomread v; + typename AT::t_f_array f; + typename AT::t_int_1d_randomread type; + + typename AT::t_neighbors_2d d_neighbors; + typename AT::t_int_1d_randomread d_ilist; + typename AT::t_int_1d_randomread d_numneigh; + + typename AT::tdual_ffloat_2d k_cutsq; + typename AT::t_ffloat_2d d_cutsq; + + Kokkos::DualView k_params; + typename Kokkos::DualView::t_dev_const_um params; + + 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; + + KOKKOS_INLINE_FUNCTION + int sbmask(const int& j) const; + friend void pair_virial_fdotr_compute(PairDPDKokkos*); + +}; +} +#endif +#endif