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

Merge pull request #3198 from bathmatt/pair-dpd

Browse Source 
DPD-BASIC kokkosification
This commit is contained in:
Axel Kohlmeyer
2022-04-08 17:48:28 -04:00
committed by GitHub
parent fc24cf15ce caa977e869
commit fbc55937fe
38 changed files with 3825 additions and 3182 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/DPD-BASIC/pair_dpd.cpp
View File

@ -46,6 +46,8 @@ PairDPD::PairDPD(LAMMPS *lmp) : Pair(lmp)
PairDPD::~PairDPD()
{
if (copymode) return;
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);

2
src/DPD-BASIC/pair_dpd.h
View File

@ -49,7 +49,7 @@ class PairDPD : public Pair {
double **sigma;
class RanMars *random;
void allocate();
virtual void allocate();
};
} // namespace LAMMPS_NS

2
src/DPD-BASIC/pair_dpd_ext.cpp
View File

@ -47,6 +47,8 @@ PairDPDExt::PairDPDExt(LAMMPS *lmp) : Pair(lmp)
PairDPDExt::~PairDPDExt()
{
if (copymode) return;
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);

4
src/DPD-BASIC/pair_dpd_ext.h
View File

@ -45,12 +45,12 @@ class PairDPDExt : public Pair {
double cut_global, temperature;
int seed;
double **cut;
double **a0, **gamma, **gammaII, **gammaT;
double **a0, **gamma, **gammaT;
double **sigma, **sigmaT;
double **ws, **wsT;
class RanMars *random;
void allocate();
virtual void allocate();
};
} // namespace LAMMPS_NS

1
src/Depend.sh
View File

@ -73,6 +73,7 @@ fi
if (test $1 = "DPD-BASIC") then
depend GPU
depend KOKKOS
depend OPENMP
depend INTEL
fi

8
src/KOKKOS/Install.sh
View File

@ -230,6 +230,14 @@ action pair_coul_long_kokkos.cpp pair_coul_long.cpp
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_kokkos.h pair_dpd.h
action pair_dpd_kokkos.cpp pair_dpd.cpp
action pair_dpd_ext_kokkos.cpp pair_dpd_ext.cpp
action pair_dpd_ext_kokkos.h pair_dpd_ext.h
action pair_dpd_ext_tstat_kokkos.h pair_dpd_ext_tstat.h
action pair_dpd_ext_tstat_kokkos.cpp pair_dpd_ext_tstat.cpp
action pair_dpd_tstat_kokkos.h pair_dpd_tstat.h
action pair_dpd_tstat_kokkos.cpp pair_dpd_tstat.cpp
action pair_dpd_fdt_energy_kokkos.cpp pair_dpd_fdt_energy.cpp
action pair_dpd_fdt_energy_kokkos.h pair_dpd_fdt_energy.h
action pair_eam_kokkos.cpp pair_eam.cpp

463
src/KOKKOS/pair_dpd_ext_kokkos.cpp Normal file
View File

@ -0,0 +1,463 @@
// 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_ext_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 <cmath>
using namespace LAMMPS_NS;
#define EPSILON 1.0e-10
template<class DeviceType>
PairDPDExtKokkos<DeviceType>::PairDPDExtKokkos(class LAMMPS *lmp) :
PairDPDExt(lmp) ,
#ifdef DPD_USE_RAN_MARS
rand_pool(0 /* unused */, lmp)
#else
rand_pool()
#endif
{
kokkosable = 1;
atomKK = (AtomKokkos *) atom;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = EMPTY_MASK;
datamask_modify = EMPTY_MASK;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
PairDPDExtKokkos<DeviceType>::~PairDPDExtKokkos() {
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<class DeviceType>
void PairDPDExtKokkos<DeviceType>::init_style()
{
PairDPDExt::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/ext/kk");
auto request = neighbor->find_request(this);
request->set_kokkos_host(std::is_same<DeviceType,LMPHostType>::value &&
!std::is_same<DeviceType,LMPDeviceType>::value);
request->set_kokkos_device(std::is_same<DeviceType,LMPDeviceType>::value);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairDPDExtKokkos<DeviceType>::compute(int eflagin, int vflagin)
{
eflag = eflagin; vflag = vflagin;
ev_init(eflag,vflag,0);
if (eflag_atom) {
memoryKK->destroy_kokkos(k_eatom,eatom);
memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
d_eatom = k_eatom.template view<DeviceType>();
}
if (vflag_atom) {
memoryKK->destroy_kokkos(k_vatom,vatom);
memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom");
d_vatom = k_vatom.template view<DeviceType>();
}
atomKK->sync(execution_space,X_MASK | V_MASK | F_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK);
x = atomKK->k_x.view<DeviceType>();
v = atomKK->k_v.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
k_cutsq.template sync<DeviceType>();
k_params.template sync<DeviceType>();
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;
dtinvsqrt = 1.0/sqrt(update->dt);
NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(list);
d_numneigh = k_list->d_numneigh;
d_neighbors = k_list->d_neighbors;
d_ilist = k_list->d_ilist;
need_dup = lmp->kokkos->need_dup<DeviceType>();
if (need_dup) {
dup_f = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterDuplicated>(f);
dup_eatom = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterDuplicated>(d_eatom);
dup_vatom = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterDuplicated>(d_vatom);
} else {
ndup_f = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterNonDuplicated>(f);
ndup_eatom = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterNonDuplicated>(d_eatom);
ndup_vatom = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterNonDuplicated>(d_vatom);
}
// loop over neighbors of my atoms
int inum = list->inum;
EV_FLOAT ev;
copymode = 1;
if (neighflag == HALF) {
if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagDPDExtKokkos<HALF,1> >(0,inum),*this,ev);
else Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagDPDExtKokkos<HALF,0> >(0,inum),*this);
} else if (neighflag == HALFTHREAD) {
if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagDPDExtKokkos<HALFTHREAD,1> >(0,inum),*this,ev);
else Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagDPDExtKokkos<HALFTHREAD,0> >(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<DeviceType>();
k_eatom.template sync<LMPHostType>();
}
if (vflag_atom) {
if (need_dup)
Kokkos::Experimental::contribute(d_vatom, dup_vatom);
k_vatom.template modify<DeviceType>();
k_vatom.template sync<LMPHostType>();
}
copymode = 0;
if (evflag) atomKK->modified(execution_space,F_MASK | ENERGY_MASK | VIRIAL_MASK);
else atomKK->modified(execution_space,F_MASK);
// free duplicated memory
if (need_dup) {
dup_f = decltype(dup_f)();
dup_eatom = decltype(dup_eatom)();
dup_vatom = decltype(dup_vatom)();
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void PairDPDExtKokkos<DeviceType>::operator() (TagDPDExtKokkos<NEIGHFLAG,EVFLAG>, const int &ii) const {
EV_FLOAT ev;
this->template operator()<NEIGHFLAG,EVFLAG>(TagDPDExtKokkos<NEIGHFLAG,EVFLAG>(), ii, ev);
}
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void PairDPDExtKokkos<DeviceType>::operator() (TagDPDExtKokkos<NEIGHFLAG,EVFLAG>, 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<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_f),decltype(ndup_f)>::get(dup_f,ndup_f);
auto a_f = v_f.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
int i,j,jj,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,fpairx,fpairy,fpairz,fpair;
double vxtmp,vytmp,vztmp,delvx,delvy,delvz;
double rsq,r,rinv,dot,wd,wdPar,wdPerp,randnum,randnumx,randnumy,randnumz,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++) {
double P[3][3];
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;
P[0][0] = 1.0 - delx*delx*rinv*rinv;
P[0][1] = - delx*dely*rinv*rinv;
P[0][2] = - delx*delz*rinv*rinv;
P[1][0] = P[0][1];
P[1][1] = 1.0 - dely*dely*rinv*rinv;
P[1][2] = - dely*delz*rinv*rinv;
P[2][0] = P[0][2];
P[2][1] = P[1][2];
P[2][2] = 1.0 - delz*delz*rinv*rinv;
wd = 1.0 - r/params(itype,jtype).cut;
wdPar = pow(wd,params(itype,jtype).ws);
wdPerp = pow(wd,params(itype,jtype).wsT);
randnum = rand_gen.normal();
randnumx = rand_gen.normal();
randnumy = rand_gen.normal();
randnumz = rand_gen.normal();
// conservative force
fpair = params(itype,jtype).a0*wd;
// drag force - parallel
fpair -= params(itype,jtype).gamma*wdPar*wdPar*dot*rinv;
// random force - parallel
fpair += params(itype,jtype).sigma*wdPar*randnum*dtinvsqrt;
fpairx = fpair*rinv*delx;
fpairy = fpair*rinv*dely;
fpairz = fpair*rinv*delz;
// drag force - perpendicular
fpairx -= params(itype,jtype).gammaT*wdPerp*wdPerp*
(P[0][0]*delvx + P[0][1]*delvy + P[0][2]*delvz);
fpairy -= params(itype,jtype).gammaT*wdPerp*wdPerp*
(P[1][0]*delvx + P[1][1]*delvy + P[1][2]*delvz);
fpairz -= params(itype,jtype).gammaT*wdPerp*wdPerp*
(P[2][0]*delvx + P[2][1]*delvy + P[2][2]*delvz);
// random force - perpendicular
fpairx += params(itype,jtype).sigmaT*wdPerp*
(P[0][0]*randnumx + P[0][1]*randnumy + P[0][2]*randnumz)*dtinvsqrt;
fpairy += params(itype,jtype).sigmaT*wdPerp*
(P[1][0]*randnumx + P[1][1]*randnumy + P[1][2]*randnumz)*dtinvsqrt;
fpairz += params(itype,jtype).sigmaT*wdPerp*
(P[2][0]*randnumx + P[2][1]*randnumy + P[2][2]*randnumz)*dtinvsqrt;
fpairx *= factor_dpd;
fpairy *= factor_dpd;
fpairz *= factor_dpd;
fx += fpairx;
fy += fpairy;
fz += fpairz;
a_f(j,0) -= fpairx;
a_f(j,1) -= fpairy;
a_f(j,2) -= fpairz;
if (EVFLAG && 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 && eflag_global)
ev.evdwl += evdwl;
}
if (EVFLAG && (eflag_atom || vflag_either))
this->template ev_tally_xyz<NEIGHFLAG>(ev,i,j,evdwl,fpairx,fpairy,fpairz,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<class DeviceType>
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void PairDPDExtKokkos<DeviceType>::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
{
// The eatom and vatom arrays are duplicated for OpenMP, atomic for CUDA, and neither for Serial
auto v_eatom = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_eatom),decltype(ndup_eatom)>::get(dup_eatom,ndup_eatom);
auto a_eatom = v_eatom.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
auto v_vatom = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_vatom),decltype(ndup_vatom)>::get(dup_vatom,ndup_vatom);
auto a_vatom = v_vatom.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
if (eflag_atom) {
const E_FLOAT epairhalf = 0.5 * epair;
a_eatom[i] += epairhalf;
a_eatom[j] += epairhalf;
}
if (vflag_either) {
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) {
ev.v[0] += v0;
ev.v[1] += v1;
ev.v[2] += v2;
ev.v[3] += v3;
ev.v[4] += v4;
ev.v[5] += v5;
}
if (vflag_atom) {
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;
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;
}
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairDPDExtKokkos<DeviceType>::allocate()
{
PairDPDExt::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<DeviceType>();
k_params = Kokkos::DualView<params_dpd**,Kokkos::LayoutRight,DeviceType>("PairDPDExt::params",n+1,n+1);
params = k_params.template view<DeviceType>();
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
int PairDPDExtKokkos<DeviceType>::sbmask(const int& j) const {
return j >> SBBITS & 3;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
template<class DeviceType>
double PairDPDExtKokkos<DeviceType>::init_one(int i, int j)
{
double cutone = PairDPDExt::init_one(i,j);
k_params.h_view(i,j).cut = cut[i][j];
k_params.h_view(i,j).ws = ws[i][j];
k_params.h_view(i,j).wsT = wsT[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(i,j).gammaT = gammaT[i][j];
k_params.h_view(i,j).sigmaT = sigmaT[i][j];
k_params.h_view(j,i) = k_params.h_view(i,j);
k_params.template modify<LMPHostType>();
k_cutsq.h_view(i,j) = cutone*cutone;
k_cutsq.h_view(j,i) = k_cutsq.h_view(i,j);
k_cutsq.template modify<LMPHostType>();
return cutone;
}
namespace LAMMPS_NS {
template class PairDPDExtKokkos<LMPDeviceType>;
#ifdef LMP_KOKKOS_GPU
template class PairDPDExtKokkos<LMPHostType>;
#endif
}

143
src/KOKKOS/pair_dpd_ext_kokkos.h Normal file
View File

@ -0,0 +1,143 @@
/* -*- 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/ext/kk,PairDPDExtKokkos<LMPDeviceType>);
PairStyle(dpd/ext/kk/device,PairDPDExtKokkos<LMPDeviceType>);
PairStyle(dpd/ext/kk/host,PairDPDExtKokkos<LMPHostType>);
// clang-format on
#else
#ifndef LMP_PAIR_DPD_EXT_KOKKOS_H
#define LMP_PAIR_DPD_EXT_KOKKOS_H
#include "pair_dpd_ext.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 DeviceType>
class PairDPDExtKokkos : public PairDPDExt {
public:
typedef DeviceType device_type;
typedef ArrayTypes<DeviceType> AT;
typedef EV_FLOAT value_type;
PairDPDExtKokkos(class LAMMPS*);
~PairDPDExtKokkos() 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=ws=wsT=a0=gamma=sigma=gammaT=sigmaT=0;}
KOKKOS_INLINE_FUNCTION
params_dpd(int /*i*/) {cut=ws=wsT=a0=gamma=sigma=gammaT=sigmaT=0;}
F_FLOAT cut,ws,wsT,a0,gamma,sigma,gammaT,sigmaT;
};
template<int NEIGHFLAG, int EVFLAG>
struct TagDPDExtKokkos{};
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator () (TagDPDExtKokkos<NEIGHFLAG,EVFLAG>, const int &i) const;
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator () (TagDPDExtKokkos<NEIGHFLAG,EVFLAG>, const int &i, EV_FLOAT&) const;
template<int NEIGHFLAG>
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;
private:
double special_lj[4];
int eflag,vflag;
int neighflag,nlocal;
double dtinvsqrt;
int need_dup;
using KKDeviceType = typename KKDevice<DeviceType>::value;
template<typename DataType, typename Layout>
using DupScatterView = KKScatterView<DataType, Layout, KKDeviceType, KKScatterSum, KKScatterDuplicated>;
template<typename DataType, typename Layout>
using NonDupScatterView = KKScatterView<DataType, Layout, KKDeviceType, KKScatterSum, KKScatterNonDuplicated>;
DupScatterView<F_FLOAT*[3], typename DAT::t_f_array::array_layout> dup_f;
DupScatterView<E_FLOAT*, typename DAT::t_efloat_1d::array_layout> dup_eatom;
DupScatterView<F_FLOAT*[6], typename DAT::t_virial_array::array_layout> dup_vatom;
NonDupScatterView<F_FLOAT*[3], typename DAT::t_f_array::array_layout> ndup_f;
NonDupScatterView<E_FLOAT*, typename DAT::t_efloat_1d::array_layout> ndup_eatom;
NonDupScatterView<F_FLOAT*[6], typename DAT::t_virial_array::array_layout> ndup_vatom;
#ifdef DPD_USE_RAN_MARS
RandPoolWrap rand_pool;
typedef RandWrap rand_type;
#elif defined(DPD_USE_Random_XorShift64)
Kokkos::Random_XorShift64_Pool<DeviceType> rand_pool;
typedef typename Kokkos::Random_XorShift64_Pool<DeviceType>::generator_type rand_type;
#elif defined(DPD_USE_Random_XorShift1024)
Kokkos::Random_XorShift1024_Pool<DeviceType> rand_pool;
typedef typename Kokkos::Random_XorShift1024_Pool<DeviceType>::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<params_dpd**,Kokkos::LayoutRight,DeviceType> k_params;
typename Kokkos::DualView<params_dpd**,
Kokkos::LayoutRight,DeviceType>::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<PairDPDExtKokkos>(PairDPDExtKokkos*);
};
}
#endif
#endif

444
src/KOKKOS/pair_dpd_ext_tstat_kokkos.cpp Normal file
View File

@ -0,0 +1,444 @@
// 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_ext_tstat_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 <cmath>
using namespace LAMMPS_NS;
#define EPSILON 1.0e-10
template<class DeviceType>
PairDPDExtTstatKokkos<DeviceType>::PairDPDExtTstatKokkos(class LAMMPS *lmp) :
PairDPDExtTstat(lmp) ,
#ifdef DPD_USE_RAN_MARS
rand_pool(0 /* unused */, lmp)
#else
rand_pool()
#endif
{
kokkosable = 1;
atomKK = (AtomKokkos *) atom;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = EMPTY_MASK;
datamask_modify = EMPTY_MASK;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
PairDPDExtTstatKokkos<DeviceType>::~PairDPDExtTstatKokkos() {
if (copymode) return;
#ifdef DPD_USE_RAN_MARS
rand_pool.destroy();
#endif
memoryKK->destroy_kokkos(k_vatom,vatom);
memoryKK->destroy_kokkos(k_cutsq,cutsq);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairDPDExtTstatKokkos<DeviceType>::init_style()
{
PairDPDExt::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/ext/kk");
auto request = neighbor->find_request(this);
request->set_kokkos_host(std::is_same<DeviceType,LMPHostType>::value &&
!std::is_same<DeviceType,LMPDeviceType>::value);
request->set_kokkos_device(std::is_same<DeviceType,LMPDeviceType>::value);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairDPDExtTstatKokkos<DeviceType>::compute(int eflagin, int vflagin)
{
eflag = eflagin; vflag = vflagin;
ev_init(eflag,vflag,0);
// adjust sigma if target T is changing
if (t_start != t_stop) {
double delta = update->ntimestep - update->beginstep;
if (delta != 0.0) delta /= update->endstep - update->beginstep;
temperature = t_start + delta * (t_stop-t_start);
double boltz = force->boltz;
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++) {
k_params.h_view(i,j).sigma = k_params.h_view(j,i).sigma =
sqrt(2.0*boltz*temperature*gamma[i][j]);
}
}
k_params.template modify<LMPHostType>();
if (eflag_atom) {
maxeatom = atom->nmax;
memory->destroy(eatom);
memory->create(eatom,maxeatom,"pair:eatom");
memset(&eatom[0], 0, maxeatom * sizeof(double));
}
if (vflag_atom) {
memoryKK->destroy_kokkos(k_vatom,vatom);
memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom");
d_vatom = k_vatom.template view<DeviceType>();
}
atomKK->sync(execution_space,X_MASK | V_MASK | F_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK);
x = atomKK->k_x.view<DeviceType>();
v = atomKK->k_v.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
k_cutsq.template sync<DeviceType>();
k_params.template sync<DeviceType>();
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;
dtinvsqrt = 1.0/sqrt(update->dt);
NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(list);
d_numneigh = k_list->d_numneigh;
d_neighbors = k_list->d_neighbors;
d_ilist = k_list->d_ilist;
need_dup = lmp->kokkos->need_dup<DeviceType>();
if (need_dup) {
dup_f = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterDuplicated>(f);
dup_vatom = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterDuplicated>(d_vatom);
} else {
ndup_f = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterNonDuplicated>(f);
ndup_vatom = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterNonDuplicated>(d_vatom);
}
// loop over neighbors of my atoms
int inum = list->inum;
EV_FLOAT ev;
copymode = 1;
if (neighflag == HALF) {
if (vflag_either) Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagDPDExtTstatKokkos<HALF,1> >(0,inum),*this,ev);
else Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagDPDExtTstatKokkos<HALF,0> >(0,inum),*this);
} else if (neighflag == HALFTHREAD) {
if (vflag_either) Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagDPDExtTstatKokkos<HALFTHREAD,1> >(0,inum),*this,ev);
else Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagDPDExtTstatKokkos<HALFTHREAD,0> >(0,inum),*this);
}
if (need_dup)
Kokkos::Experimental::contribute(f, dup_f);
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 (vflag_atom) {
if (need_dup)
Kokkos::Experimental::contribute(d_vatom, dup_vatom);
k_vatom.template modify<DeviceType>();
k_vatom.template sync<LMPHostType>();
}
copymode = 0;
if (evflag) atomKK->modified(execution_space,F_MASK | ENERGY_MASK | VIRIAL_MASK);
else atomKK->modified(execution_space,F_MASK);
// free duplicated memory
if (need_dup) {
dup_f = decltype(dup_f)();
dup_vatom = decltype(dup_vatom)();
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG, int VFLAG>
KOKKOS_INLINE_FUNCTION
void PairDPDExtTstatKokkos<DeviceType>::operator() (TagDPDExtTstatKokkos<NEIGHFLAG,VFLAG>, const int &ii) const {
EV_FLOAT ev;
this->template operator()<NEIGHFLAG,VFLAG>(TagDPDExtTstatKokkos<NEIGHFLAG,VFLAG>(), ii, ev);
}
template<class DeviceType>
template<int NEIGHFLAG, int VFLAG>
KOKKOS_INLINE_FUNCTION
void PairDPDExtTstatKokkos<DeviceType>::operator() (TagDPDExtTstatKokkos<NEIGHFLAG,VFLAG>, 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<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_f),decltype(ndup_f)>::get(dup_f,ndup_f);
auto a_f = v_f.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
int i,j,jj,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,fpairx,fpairy,fpairz,fpair;
double vxtmp,vytmp,vztmp,delvx,delvy,delvz;
double rsq,r,rinv,dot,wd,wdPar,wdPerp,randnum,randnumx,randnumy,randnumz,factor_dpd;
double fx = 0,fy = 0,fz = 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++) {
double P[3][3];
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;
P[0][0] = 1.0 - delx*delx*rinv*rinv;
P[0][1] = - delx*dely*rinv*rinv;
P[0][2] = - delx*delz*rinv*rinv;
P[1][0] = P[0][1];
P[1][1] = 1.0 - dely*dely*rinv*rinv;
P[1][2] = - dely*delz*rinv*rinv;
P[2][0] = P[0][2];
P[2][1] = P[1][2];
P[2][2] = 1.0 - delz*delz*rinv*rinv;
wd = 1.0 - r/params(itype,jtype).cut;
wdPar = pow(wd,params(itype,jtype).ws);
wdPerp = pow(wd,params(itype,jtype).wsT);
randnum = rand_gen.normal();
randnumx = rand_gen.normal();
randnumy = rand_gen.normal();
randnumz = rand_gen.normal();
// drag force - parallel
fpair = -params(itype,jtype).gamma*wdPar*wdPar*dot*rinv;
// random force - parallel
fpair += params(itype,jtype).sigma*wdPar*randnum*dtinvsqrt;
fpairx = fpair*rinv*delx;
fpairy = fpair*rinv*dely;
fpairz = fpair*rinv*delz;
// drag force - perpendicular
fpairx -= params(itype,jtype).gammaT*wdPerp*wdPerp*
(P[0][0]*delvx + P[0][1]*delvy + P[0][2]*delvz);
fpairy -= params(itype,jtype).gammaT*wdPerp*wdPerp*
(P[1][0]*delvx + P[1][1]*delvy + P[1][2]*delvz);
fpairz -= params(itype,jtype).gammaT*wdPerp*wdPerp*
(P[2][0]*delvx + P[2][1]*delvy + P[2][2]*delvz);
// random force - perpendicular
fpairx += params(itype,jtype).sigmaT*wdPerp*
(P[0][0]*randnumx + P[0][1]*randnumy + P[0][2]*randnumz)*dtinvsqrt;
fpairy += params(itype,jtype).sigmaT*wdPerp*
(P[1][0]*randnumx + P[1][1]*randnumy + P[1][2]*randnumz)*dtinvsqrt;
fpairz += params(itype,jtype).sigmaT*wdPerp*
(P[2][0]*randnumx + P[2][1]*randnumy + P[2][2]*randnumz)*dtinvsqrt;
fpairx *= factor_dpd;
fpairy *= factor_dpd;
fpairz *= factor_dpd;
fx += fpairx;
fy += fpairy;
fz += fpairz;
a_f(j,0) -= fpairx;
a_f(j,1) -= fpairy;
a_f(j,2) -= fpairz;
if (VFLAG)
this->template v_tally_xyz<NEIGHFLAG>(ev,i,j,fpairx,fpairy,fpairz,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<class DeviceType>
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void PairDPDExtTstatKokkos<DeviceType>::v_tally_xyz(EV_FLOAT &ev, const int &i, const int &j,
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
{
// The vatom array is duplicated for OpenMP, atomic for CUDA, and neither for Serial
auto v_vatom = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_vatom),decltype(ndup_vatom)>::get(dup_vatom,ndup_vatom);
auto a_vatom = v_vatom.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
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) {
ev.v[0] += v0;
ev.v[1] += v1;
ev.v[2] += v2;
ev.v[3] += v3;
ev.v[4] += v4;
ev.v[5] += v5;
}
if (vflag_atom) {
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;
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;
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairDPDExtTstatKokkos<DeviceType>::allocate()
{
PairDPDExt::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<DeviceType>();
k_params = Kokkos::DualView<params_dpd**,Kokkos::LayoutRight,DeviceType>("PairDPDExt::params",n+1,n+1);
params = k_params.template view<DeviceType>();
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
int PairDPDExtTstatKokkos<DeviceType>::sbmask(const int& j) const {
return j >> SBBITS & 3;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
template<class DeviceType>
double PairDPDExtTstatKokkos<DeviceType>::init_one(int i, int j)
{
double cutone = PairDPDExt::init_one(i,j);
k_params.h_view(i,j).cut = cut[i][j];
k_params.h_view(i,j).ws = ws[i][j];
k_params.h_view(i,j).wsT = wsT[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(i,j).gammaT = gammaT[i][j];
k_params.h_view(i,j).sigmaT = sigmaT[i][j];
k_params.h_view(j,i) = k_params.h_view(i,j);
k_params.template modify<LMPHostType>();
k_cutsq.h_view(i,j) = cutone*cutone;
k_cutsq.h_view(j,i) = k_cutsq.h_view(i,j);
k_cutsq.template modify<LMPHostType>();
return cutone;
}
namespace LAMMPS_NS {
template class PairDPDExtTstatKokkos<LMPDeviceType>;
#ifdef LMP_KOKKOS_GPU
template class PairDPDExtTstatKokkos<LMPHostType>;
#endif
}

138
src/KOKKOS/pair_dpd_ext_tstat_kokkos.h Normal file
View File

@ -0,0 +1,138 @@
/* -*- 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/ext/tstat/kk,PairDPDExtTstatKokkos<LMPDeviceType>);
PairStyle(dpd/ext/tstat/kk/device,PairDPDExtTstatKokkos<LMPDeviceType>);
PairStyle(dpd/ext/tstat/kk/host,PairDPDExtTstatKokkos<LMPHostType>);
// clang-format on
#else
#ifndef LMP_PAIR_DPD_EXT_TSTAT_KOKKOS_H
#define LMP_PAIR_DPD_EXT_TSTAT_KOKKOS_H
#include "pair_dpd_ext_tstat.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 DeviceType>
class PairDPDExtTstatKokkos : public PairDPDExtTstat {
public:
typedef DeviceType device_type;
typedef ArrayTypes<DeviceType> AT;
typedef EV_FLOAT value_type;
PairDPDExtTstatKokkos(class LAMMPS*);
~PairDPDExtTstatKokkos() 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=ws=wsT=gamma=sigma=gammaT=sigmaT=0;}
KOKKOS_INLINE_FUNCTION
params_dpd(int /*i*/) {cut=ws=wsT=gamma=sigma=gammaT=sigmaT=0;}
F_FLOAT cut,ws,wsT,gamma,sigma,gammaT,sigmaT;
};
template<int NEIGHFLAG, int VFLAG>
struct TagDPDExtTstatKokkos{};
template<int NEIGHFLAG, int VFLAG>
KOKKOS_INLINE_FUNCTION
void operator () (TagDPDExtTstatKokkos<NEIGHFLAG,VFLAG>, const int &i) const;
template<int NEIGHFLAG, int VFLAG>
KOKKOS_INLINE_FUNCTION
void operator () (TagDPDExtTstatKokkos<NEIGHFLAG,VFLAG>, const int &i, EV_FLOAT&) const;
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void v_tally_xyz(EV_FLOAT &ev, const int &i, const int &j,
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;
private:
double special_lj[4];
int eflag,vflag;
int neighflag,nlocal;
double dtinvsqrt;
int need_dup;
using KKDeviceType = typename KKDevice<DeviceType>::value;
template<typename DataType, typename Layout>
using DupScatterView = KKScatterView<DataType, Layout, KKDeviceType, KKScatterSum, KKScatterDuplicated>;
template<typename DataType, typename Layout>
using NonDupScatterView = KKScatterView<DataType, Layout, KKDeviceType, KKScatterSum, KKScatterNonDuplicated>;
DupScatterView<F_FLOAT*[3], typename DAT::t_f_array::array_layout> dup_f;
DupScatterView<F_FLOAT*[6], typename DAT::t_virial_array::array_layout> dup_vatom;
NonDupScatterView<F_FLOAT*[3], typename DAT::t_f_array::array_layout> ndup_f;
NonDupScatterView<F_FLOAT*[6], typename DAT::t_virial_array::array_layout> ndup_vatom;
#ifdef DPD_USE_RAN_MARS
RandPoolWrap rand_pool;
typedef RandWrap rand_type;
#elif defined(DPD_USE_Random_XorShift64)
Kokkos::Random_XorShift64_Pool<DeviceType> rand_pool;
typedef typename Kokkos::Random_XorShift64_Pool<DeviceType>::generator_type rand_type;
#elif defined(DPD_USE_Random_XorShift1024)
Kokkos::Random_XorShift1024_Pool<DeviceType> rand_pool;
typedef typename Kokkos::Random_XorShift1024_Pool<DeviceType>::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<params_dpd**,Kokkos::LayoutRight,DeviceType> k_params;
typename Kokkos::DualView<params_dpd**,
Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
DAT::tdual_virial_array k_vatom;
typename AT::t_virial_array d_vatom;
KOKKOS_INLINE_FUNCTION
int sbmask(const int& j) const;
friend void pair_virial_fdotr_compute<PairDPDExtTstatKokkos>(PairDPDExtTstatKokkos*);
};
}
#endif
#endif

26
src/KOKKOS/pair_dpd_fdt_energy_kokkos.cpp
View File

@ -729,21 +729,19 @@ void PairDPDfdtEnergyKokkos<DeviceType>::ev_tally(EV_FLOAT &ev, const int &i, co
if (vflag_atom) {
if (NEIGHFLAG!=FULL) {
if (NEWTON_PAIR || i < nlocal) {
v_vatom(i,0) += 0.5*v0;
v_vatom(i,1) += 0.5*v1;
v_vatom(i,2) += 0.5*v2;
v_vatom(i,3) += 0.5*v3;
v_vatom(i,4) += 0.5*v4;
v_vatom(i,5) += 0.5*v5;
}
v_vatom(i,0) += 0.5*v0;
v_vatom(i,1) += 0.5*v1;
v_vatom(i,2) += 0.5*v2;
v_vatom(i,3) += 0.5*v3;
v_vatom(i,4) += 0.5*v4;
v_vatom(i,5) += 0.5*v5;
if (NEWTON_PAIR || j < nlocal) {
v_vatom(j,0) += 0.5*v0;
v_vatom(j,1) += 0.5*v1;
v_vatom(j,2) += 0.5*v2;
v_vatom(j,3) += 0.5*v3;
v_vatom(j,4) += 0.5*v4;
v_vatom(j,5) += 0.5*v5;
v_vatom(j,0) += 0.5*v0;
v_vatom(j,1) += 0.5*v1;
v_vatom(j,2) += 0.5*v2;
v_vatom(j,3) += 0.5*v3;
v_vatom(j,4) += 0.5*v4;
v_vatom(j,5) += 0.5*v5;
}
} else {
v_vatom(i,0) += 0.5*v0;

417
src/KOKKOS/pair_dpd_kokkos.cpp Normal file
View File

@ -0,0 +1,417 @@
// 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 <cmath>
using namespace LAMMPS_NS;
#define EPSILON 1.0e-10
template<class DeviceType>
PairDPDKokkos<DeviceType>::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<DeviceType>::space;
datamask_read = EMPTY_MASK;
datamask_modify = EMPTY_MASK;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
PairDPDKokkos<DeviceType>::~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<class DeviceType>
void PairDPDKokkos<DeviceType>::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<DeviceType,LMPHostType>::value &&
!std::is_same<DeviceType,LMPDeviceType>::value);
request->set_kokkos_device(std::is_same<DeviceType,LMPDeviceType>::value);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairDPDKokkos<DeviceType>::compute(int eflagin, int vflagin)
{
eflag = eflagin; vflag = vflagin;
ev_init(eflag,vflag,0);
if (eflag_atom) {
memoryKK->destroy_kokkos(k_eatom,eatom);
memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
d_eatom = k_eatom.template view<DeviceType>();
}
if (vflag_atom) {
memoryKK->destroy_kokkos(k_vatom,vatom);
memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom");
d_vatom = k_vatom.template view<DeviceType>();
}
atomKK->sync(execution_space,X_MASK | V_MASK | F_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK);
x = atomKK->k_x.view<DeviceType>();
v = atomKK->k_v.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
k_cutsq.template sync<DeviceType>();
k_params.template sync<DeviceType>();
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;
dtinvsqrt = 1.0/sqrt(update->dt);
NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(list);
d_numneigh = k_list->d_numneigh;
d_neighbors = k_list->d_neighbors;
d_ilist = k_list->d_ilist;
need_dup = lmp->kokkos->need_dup<DeviceType>();
if (need_dup) {
dup_f = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterDuplicated>(f);
dup_eatom = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterDuplicated>(d_eatom);
dup_vatom = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterDuplicated>(d_vatom);
} else {
ndup_f = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterNonDuplicated>(f);
ndup_eatom = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterNonDuplicated>(d_eatom);
ndup_vatom = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterNonDuplicated>(d_vatom);
}
// loop over neighbors of my atoms
int inum = list->inum;
EV_FLOAT ev;
copymode = 1;
if (neighflag == HALF) {
if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagDPDKokkos<HALF,1> >(0,inum),*this,ev);
else Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagDPDKokkos<HALF,0> >(0,inum),*this);
} else if (neighflag == HALFTHREAD) {
if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagDPDKokkos<HALFTHREAD,1> >(0,inum),*this,ev);
else Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagDPDKokkos<HALFTHREAD,0> >(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<DeviceType>();
k_eatom.template sync<LMPHostType>();
}
if (vflag_atom) {
if (need_dup)
Kokkos::Experimental::contribute(d_vatom, dup_vatom);
k_vatom.template modify<DeviceType>();
k_vatom.template sync<LMPHostType>();
}
copymode = 0;
if (evflag) atomKK->modified(execution_space,F_MASK | ENERGY_MASK | VIRIAL_MASK);
else atomKK->modified(execution_space,F_MASK);
// free duplicated memory
if (need_dup) {
dup_f = decltype(dup_f)();
dup_eatom = decltype(dup_eatom)();
dup_vatom = decltype(dup_vatom)();
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void PairDPDKokkos<DeviceType>::operator() (TagDPDKokkos<NEIGHFLAG,EVFLAG>, const int &ii) const {
EV_FLOAT ev;
this->template operator()<NEIGHFLAG,EVFLAG>(TagDPDKokkos<NEIGHFLAG,EVFLAG>(), ii, ev);
}
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void PairDPDKokkos<DeviceType>::operator() (TagDPDKokkos<NEIGHFLAG,EVFLAG>, 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<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_f),decltype(ndup_f)>::get(dup_f,ndup_f);
auto a_f = v_f.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
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;
a_f(j,0) -= fpair*delx;
a_f(j,1) -= fpair*dely;
a_f(j,2) -= fpair*delz;
if (EVFLAG && eflag_global) {
// 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;
ev.evdwl += evdwl;
}
if (EVFLAG && (eflag_atom || vflag_either))
this->template ev_tally<NEIGHFLAG>(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<class DeviceType>
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void PairDPDKokkos<DeviceType>::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
{
// The eatom and vatom arrays are duplicated for OpenMP, atomic for CUDA, and neither for Serial
auto v_eatom = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_eatom),decltype(ndup_eatom)>::get(dup_eatom,ndup_eatom);
auto a_eatom = v_eatom.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
auto v_vatom = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_vatom),decltype(ndup_vatom)>::get(dup_vatom,ndup_vatom);
auto a_vatom = v_vatom.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
if (eflag_atom) {
const E_FLOAT epairhalf = 0.5 * epair;
a_eatom[i] += epairhalf;
a_eatom[j] += epairhalf;
}
if (vflag_either) {
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) {
ev.v[0] += v0;
ev.v[1] += v1;
ev.v[2] += v2;
ev.v[3] += v3;
ev.v[4] += v4;
ev.v[5] += v5;
}
if (vflag_atom) {
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;
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;
}
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairDPDKokkos<DeviceType>::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<DeviceType>();
k_params = Kokkos::DualView<params_dpd**,Kokkos::LayoutRight,DeviceType>("PairDPD::params",n+1,n+1);
params = k_params.template view<DeviceType>();
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
int PairDPDKokkos<DeviceType>::sbmask(const int& j) const {
return j >> SBBITS & 3;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
template<class DeviceType>
double PairDPDKokkos<DeviceType>::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<LMPHostType>();
k_cutsq.h_view(i,j) = cutone*cutone;
k_cutsq.h_view(j,i) = k_cutsq.h_view(i,j);
k_cutsq.template modify<LMPHostType>();
return cutone;
}
namespace LAMMPS_NS {
template class PairDPDKokkos<LMPDeviceType>;
#ifdef LMP_KOKKOS_GPU
template class PairDPDKokkos<LMPHostType>;
#endif
}

142
src/KOKKOS/pair_dpd_kokkos.h Normal file
View File

@ -0,0 +1,142 @@
/* -*- 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<LMPDeviceType>);
PairStyle(dpd/kk/device,PairDPDKokkos<LMPDeviceType>);
PairStyle(dpd/kk/host,PairDPDKokkos<LMPHostType>);
// 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 DeviceType>
class PairDPDKokkos : public PairDPD {
public:
typedef DeviceType device_type;
typedef ArrayTypes<DeviceType> 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 cut,a0,gamma,sigma;
};
template<int NEIGHFLAG, int EVFLAG>
struct TagDPDKokkos{};
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator () (TagDPDKokkos<NEIGHFLAG,EVFLAG>, const int &i) const;
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator () (TagDPDKokkos<NEIGHFLAG,EVFLAG>, const int &i, EV_FLOAT&) const;
template<int NEIGHFLAG>
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;
double dtinvsqrt;
int need_dup;
using KKDeviceType = typename KKDevice<DeviceType>::value;
template<typename DataType, typename Layout>
using DupScatterView = KKScatterView<DataType, Layout, KKDeviceType, KKScatterSum, KKScatterDuplicated>;
template<typename DataType, typename Layout>
using NonDupScatterView = KKScatterView<DataType, Layout, KKDeviceType, KKScatterSum, KKScatterNonDuplicated>;
DupScatterView<F_FLOAT*[3], typename DAT::t_f_array::array_layout> dup_f;
DupScatterView<E_FLOAT*, typename DAT::t_efloat_1d::array_layout> dup_eatom;
DupScatterView<F_FLOAT*[6], typename DAT::t_virial_array::array_layout> dup_vatom;
NonDupScatterView<F_FLOAT*[3], typename DAT::t_f_array::array_layout> ndup_f;
NonDupScatterView<E_FLOAT*, typename DAT::t_efloat_1d::array_layout> ndup_eatom;
NonDupScatterView<F_FLOAT*[6], typename DAT::t_virial_array::array_layout> ndup_vatom;
#ifdef DPD_USE_RAN_MARS
RandPoolWrap rand_pool;
typedef RandWrap rand_type;
#elif defined(DPD_USE_Random_XorShift64)
Kokkos::Random_XorShift64_Pool<DeviceType> rand_pool;
typedef typename Kokkos::Random_XorShift64_Pool<DeviceType>::generator_type rand_type;
#elif defined(DPD_USE_Random_XorShift1024)
Kokkos::Random_XorShift1024_Pool<DeviceType> rand_pool;
typedef typename Kokkos::Random_XorShift1024_Pool<DeviceType>::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<params_dpd**,Kokkos::LayoutRight,DeviceType> k_params;
typename Kokkos::DualView<params_dpd**,
Kokkos::LayoutRight,DeviceType>::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>(PairDPDKokkos*);
};
}
#endif
#endif

399
src/KOKKOS/pair_dpd_tstat_kokkos.cpp Normal file
View File

@ -0,0 +1,399 @@
// 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_tstat_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 <cmath>
using namespace LAMMPS_NS;
#define EPSILON 1.0e-10
template<class DeviceType>
PairDPDTstatKokkos<DeviceType>::PairDPDTstatKokkos(class LAMMPS *lmp) :
PairDPDTstat(lmp) ,
#ifdef DPD_USE_RAN_MARS
rand_pool(0 /* unused */, lmp)
#else
rand_pool()
#endif
{
kokkosable = 1;
atomKK = (AtomKokkos *) atom;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = EMPTY_MASK;
datamask_modify = EMPTY_MASK;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
PairDPDTstatKokkos<DeviceType>::~PairDPDTstatKokkos() {
if (copymode) return;
#ifdef DPD_USE_RAN_MARS
rand_pool.destroy();
#endif
memoryKK->destroy_kokkos(k_vatom,vatom);
memoryKK->destroy_kokkos(k_cutsq,cutsq);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairDPDTstatKokkos<DeviceType>::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<DeviceType,LMPHostType>::value &&
!std::is_same<DeviceType,LMPDeviceType>::value);
request->set_kokkos_device(std::is_same<DeviceType,LMPDeviceType>::value);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairDPDTstatKokkos<DeviceType>::compute(int eflagin, int vflagin)
{
eflag = eflagin; vflag = vflagin;
ev_init(eflag,vflag,0);
// adjust sigma if target T is changing
if (t_start != t_stop) {
double delta = update->ntimestep - update->beginstep;
if (delta != 0.0) delta /= update->endstep - update->beginstep;
temperature = t_start + delta * (t_stop-t_start);
double boltz = force->boltz;
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++) {
k_params.h_view(i,j).sigma = k_params.h_view(j,i).sigma =
sqrt(2.0*boltz*temperature*gamma[i][j]);
}
}
k_params.template modify<LMPHostType>();
if (eflag_atom) {
maxeatom = atom->nmax;
memory->destroy(eatom);
memory->create(eatom,maxeatom,"pair:eatom");
memset(&eatom[0], 0, maxeatom * sizeof(double));
}
if (vflag_atom) {
memoryKK->destroy_kokkos(k_vatom,vatom);
memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom");
d_vatom = k_vatom.template view<DeviceType>();
}
atomKK->sync(execution_space,X_MASK | V_MASK | F_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK);
x = atomKK->k_x.view<DeviceType>();
v = atomKK->k_v.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
k_cutsq.template sync<DeviceType>();
k_params.template sync<DeviceType>();
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;
dtinvsqrt = 1.0/sqrt(update->dt);
NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(list);
d_numneigh = k_list->d_numneigh;
d_neighbors = k_list->d_neighbors;
d_ilist = k_list->d_ilist;
need_dup = lmp->kokkos->need_dup<DeviceType>();
if (need_dup) {
dup_f = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterDuplicated>(f);
dup_vatom = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterDuplicated>(d_vatom);
} else {
ndup_f = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterNonDuplicated>(f);
ndup_vatom = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterNonDuplicated>(d_vatom);
}
// loop over neighbors of my atoms
int inum = list->inum;
EV_FLOAT ev;
copymode = 1;
if (neighflag == HALF) {
if (vflag_either) Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagDPDTstatKokkos<HALF,1> >(0,inum),*this,ev);
else Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagDPDTstatKokkos<HALF,0> >(0,inum),*this);
} else if (neighflag == HALFTHREAD) {
if (vflag_either) Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagDPDTstatKokkos<HALFTHREAD,1> >(0,inum),*this,ev);
else Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagDPDTstatKokkos<HALFTHREAD,0> >(0,inum),*this);
}
if (need_dup)
Kokkos::Experimental::contribute(f, dup_f);
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 (vflag_atom) {
if (need_dup)
Kokkos::Experimental::contribute(d_vatom, dup_vatom);
k_vatom.template modify<DeviceType>();
k_vatom.template sync<LMPHostType>();
}
copymode = 0;
if (evflag) atomKK->modified(execution_space,F_MASK | ENERGY_MASK | VIRIAL_MASK);
else atomKK->modified(execution_space,F_MASK);
// free duplicated memory
if (need_dup) {
dup_f = decltype(dup_f)();
dup_vatom = decltype(dup_vatom)();
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG, int VFLAG>
KOKKOS_INLINE_FUNCTION
void PairDPDTstatKokkos<DeviceType>::operator() (TagDPDTstatKokkos<NEIGHFLAG,VFLAG>, const int &ii) const {
EV_FLOAT ev;
this->template operator()<NEIGHFLAG,VFLAG>(TagDPDTstatKokkos<NEIGHFLAG,VFLAG>(), ii, ev);
}
template<class DeviceType>
template<int NEIGHFLAG, int VFLAG>
KOKKOS_INLINE_FUNCTION
void PairDPDTstatKokkos<DeviceType>::operator() (TagDPDTstatKokkos<NEIGHFLAG,VFLAG>, 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<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_f),decltype(ndup_f)>::get(dup_f,ndup_f);
auto a_f = v_f.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
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;
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();
// 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;
a_f(j,0) -= fpair*delx;
a_f(j,1) -= fpair*dely;
a_f(j,2) -= fpair*delz;
if (VFLAG)
this->template v_tally<NEIGHFLAG>(ev,i,j,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<class DeviceType>
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void PairDPDTstatKokkos<DeviceType>::v_tally(EV_FLOAT &ev, const int &i, const int &j,
const F_FLOAT &fpair, const F_FLOAT &delx,
const F_FLOAT &dely, const F_FLOAT &delz) const
{
// The vatom array is duplicated for OpenMP, atomic for CUDA, and neither for Serial
auto v_vatom = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_vatom),decltype(ndup_vatom)>::get(dup_vatom,ndup_vatom);
auto a_vatom = v_vatom.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
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) {
ev.v[0] += v0;
ev.v[1] += v1;
ev.v[2] += v2;
ev.v[3] += v3;
ev.v[4] += v4;
ev.v[5] += v5;
}
if (vflag_atom) {
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;
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;
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairDPDTstatKokkos<DeviceType>::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<DeviceType>();
k_params = Kokkos::DualView<params_dpd**,Kokkos::LayoutRight,DeviceType>("PairDPD::params",n+1,n+1);
params = k_params.template view<DeviceType>();
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
int PairDPDTstatKokkos<DeviceType>::sbmask(const int& j) const {
return j >> SBBITS & 3;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
template<class DeviceType>
double PairDPDTstatKokkos<DeviceType>::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).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<LMPHostType>();
k_cutsq.h_view(i,j) = cutone*cutone;
k_cutsq.h_view(j,i) = k_cutsq.h_view(i,j);
k_cutsq.template modify<LMPHostType>();
return cutone;
}
namespace LAMMPS_NS {
template class PairDPDTstatKokkos<LMPDeviceType>;
#ifdef LMP_KOKKOS_GPU
template class PairDPDTstatKokkos<LMPHostType>;
#endif
}

138
src/KOKKOS/pair_dpd_tstat_kokkos.h Normal file
View File

@ -0,0 +1,138 @@
/* -*- 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/tstat/kk,PairDPDTstatKokkos<LMPDeviceType>);
PairStyle(dpd/tstat/kk/device,PairDPDTstatKokkos<LMPDeviceType>);
PairStyle(dpd/tstat/kk/host,PairDPDTstatKokkos<LMPHostType>);
// clang-format on
#else
#ifndef LMP_PAIR_DPD_TSTAT_KOKKOS_H
#define LMP_PAIR_DPD_TSTAT_KOKKOS_H
#include "pair_dpd_tstat.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 DeviceType>
class PairDPDTstatKokkos : public PairDPDTstat {
public:
typedef DeviceType device_type;
typedef ArrayTypes<DeviceType> AT;
typedef EV_FLOAT value_type;
PairDPDTstatKokkos(class LAMMPS*);
~PairDPDTstatKokkos() 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=gamma=sigma=0;}
KOKKOS_INLINE_FUNCTION
params_dpd(int /*i*/) {cut=gamma=sigma=0;}
F_FLOAT cut,gamma,sigma;
};
template<int NEIGHFLAG, int VFLAG>
struct TagDPDTstatKokkos{};
template<int NEIGHFLAG, int VFLAG>
KOKKOS_INLINE_FUNCTION
void operator () (TagDPDTstatKokkos<NEIGHFLAG,VFLAG>, const int &i) const;
template<int NEIGHFLAG, int VFLAG>
KOKKOS_INLINE_FUNCTION
void operator () (TagDPDTstatKokkos<NEIGHFLAG,VFLAG>, const int &i, EV_FLOAT&) const;
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void v_tally(EV_FLOAT &ev, const int &i, const int &j,
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;
double dtinvsqrt;
int need_dup;
using KKDeviceType = typename KKDevice<DeviceType>::value;
template<typename DataType, typename Layout>
using DupScatterView = KKScatterView<DataType, Layout, KKDeviceType, KKScatterSum, KKScatterDuplicated>;
template<typename DataType, typename Layout>
using NonDupScatterView = KKScatterView<DataType, Layout, KKDeviceType, KKScatterSum, KKScatterNonDuplicated>;
DupScatterView<F_FLOAT*[3], typename DAT::t_f_array::array_layout> dup_f;
DupScatterView<F_FLOAT*[6], typename DAT::t_virial_array::array_layout> dup_vatom;
NonDupScatterView<F_FLOAT*[3], typename DAT::t_f_array::array_layout> ndup_f;
NonDupScatterView<F_FLOAT*[6], typename DAT::t_virial_array::array_layout> ndup_vatom;
#ifdef DPD_USE_RAN_MARS
RandPoolWrap rand_pool;
typedef RandWrap rand_type;
#elif defined(DPD_USE_Random_XorShift64)
Kokkos::Random_XorShift64_Pool<DeviceType> rand_pool;
typedef typename Kokkos::Random_XorShift64_Pool<DeviceType>::generator_type rand_type;
#elif defined(DPD_USE_Random_XorShift1024)
Kokkos::Random_XorShift1024_Pool<DeviceType> rand_pool;
typedef typename Kokkos::Random_XorShift1024_Pool<DeviceType>::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<params_dpd**,Kokkos::LayoutRight,DeviceType> k_params;
typename Kokkos::DualView<params_dpd**,
Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
DAT::tdual_virial_array k_vatom;
typename AT::t_virial_array d_vatom;
KOKKOS_INLINE_FUNCTION
int sbmask(const int& j) const;
friend void pair_virial_fdotr_compute<PairDPDTstatKokkos>(PairDPDTstatKokkos*);
};
}
#endif
#endif

2
src/KOKKOS/pair_reaxff_kokkos.cpp
View File

@ -694,7 +694,7 @@ void PairReaxFFKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
eflag = eflag_in;
vflag = vflag_in;
ev_init(eflag,vflag);
ev_init(eflag,vflag,0);
atomKK->sync(execution_space,datamask_read);
k_params_sing.template sync<DeviceType>();