Homec4science
Diffusion lammps (master) Browse

lammps/src/KOKKOS/pair_multi_lucy_rx_kokkos.cpp
6f51c3b75c49master

Last Change
Raw File

pair_multi_lucy_rx_kokkos.cpp
View Options

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, 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:
Stan Moore (Sandia)
Please cite the related publications:
J.D. Moore, B.C. Barnes, S. Izvekov, M. Lisal, M.S. Sellers, D.E. Taylor & J.K. Brennan
"A coarse-grain force field for RDX: Density dependent and energy conserving"
The Journal of Chemical Physics, 2016, 144, 104501.
------------------------------------------------------------------------------------------- */
#include <mpi.h>
#include <math.h>
#include "math_const.h"
#include <stdlib.h>
#include <string.h>
#include "pair_multi_lucy_rx_kokkos.h"
#include "atom_kokkos.h"
#include "force.h"
#include "comm.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "citeme.h"
#include "modify.h"
#include "fix.h"
#include "atom_masks.h"
#include "neigh_request.h"
using namespace LAMMPS_NS;
enum{NONE,RLINEAR,RSQ};
#define MAXLINE 1024
#define oneFluidParameter (-1)
#define isOneFluid(_site) ( (_site) == oneFluidParameter )
static const char cite_pair_multi_lucy_rx[] =
"pair_style multi/lucy/rx command:\n\n"
"@Article{Moore16,\n"
" author = {J.D. Moore, B.C. Barnes, S. Izvekov, M. Lisal, M.S. Sellers, D.E. Taylor and J. K. Brennan},\n"
" title = {A coarse-grain force field for RDX: Density dependent and energy conserving},\n"
" journal = {J. Chem. Phys.},\n"
" year = 2016,\n"
" volume = 144\n"
" pages = {104501}\n"
"}\n\n";
/* ---------------------------------------------------------------------- */
template<class DeviceType>
PairMultiLucyRXKokkos<DeviceType>::PairMultiLucyRXKokkos(LAMMPS *lmp) : PairMultiLucyRX(lmp)
{
respa_enable = 0;
atomKK = (AtomKokkos *) atom;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = EMPTY_MASK;
datamask_modify = EMPTY_MASK;
update_table = 0;
h_table = new TableHost();
d_table = new TableDevice();
k_error_flag = DAT::tdual_int_scalar("pair:error_flag");
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
PairMultiLucyRXKokkos<DeviceType>::~PairMultiLucyRXKokkos()
{
if (copymode) return;
delete h_table;
delete d_table;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairMultiLucyRXKokkos<DeviceType>::init_style()
{
PairMultiLucyRX::init_style();
// irequest = neigh request made by parent class
neighflag = lmp->kokkos->neighflag;
int irequest = neighbor->nrequest - 1;
neighbor->requests[irequest]->
kokkos_host = Kokkos::Impl::is_same<DeviceType,LMPHostType>::value &&
!Kokkos::Impl::is_same<DeviceType,LMPDeviceType>::value;
neighbor->requests[irequest]->
kokkos_device = Kokkos::Impl::is_same<DeviceType,LMPDeviceType>::value;
if (neighflag == FULL) {
neighbor->requests[irequest]->full = 1;
neighbor->requests[irequest]->half = 0;
} else if (neighflag == HALF || neighflag == HALFTHREAD) {
neighbor->requests[irequest]->full = 0;
neighbor->requests[irequest]->half = 1;
} else {
error->all(FLERR,"Cannot use chosen neighbor list style with multi/lucy/rx/kk");
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairMultiLucyRXKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
{
if (update_table)
create_kokkos_tables();
if (tabstyle == LOOKUP)
compute_style<LOOKUP>(eflag_in,vflag_in);
else if(tabstyle == LINEAR)
compute_style<LINEAR>(eflag_in,vflag_in);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
template<int TABSTYLE>
void PairMultiLucyRXKokkos<DeviceType>::compute_style(int eflag_in, int vflag_in)
{
copymode = 1;
eflag = eflag_in;
vflag = vflag_in;
double evdwl,evdwlOld;
evdwlOld = 0.0;
evdwl = 0.0;
if (neighflag == FULL) no_virial_fdotr_compute = 1;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
// reallocate per-atom arrays if necessary
if (eflag_atom) {
memory->destroy_kokkos(k_eatom,eatom);
memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
d_eatom = k_eatom.d_view;
}
if (vflag_atom) {
memory->destroy_kokkos(k_vatom,vatom);
memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
d_vatom = k_vatom.d_view;
}
x = atomKK->k_x.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
rho = atomKK->k_rho.view<DeviceType>();
uCG = atomKK->k_uCG.view<DeviceType>();
uCGnew = atomKK->k_uCGnew.view<DeviceType>();
dvector = atomKK->k_dvector.view<DeviceType>();
atomKK->sync(execution_space,X_MASK | F_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK | DPDRHO_MASK | UCG_MASK | UCGNEW_MASK | DVECTOR_MASK);
if (evflag) atomKK->modified(execution_space,F_MASK | ENERGY_MASK | VIRIAL_MASK | UCG_MASK | UCGNEW_MASK);
else atomKK->modified(execution_space,F_MASK | UCG_MASK | UCGNEW_MASK);
nlocal = atom->nlocal;
int nghost = atom->nghost;
int newton_pair = force->newton_pair;
{
const int ntotal = nlocal + nghost;
d_fractionOld1 = typename AT::t_float_1d("PairMultiLucyRX::fractionOld1",ntotal);
d_fractionOld2 = typename AT::t_float_1d("PairMultiLucyRX::fractionOld2",ntotal);
d_fraction1 = typename AT::t_float_1d("PairMultiLucyRX::fraction1",ntotal);
d_fraction2 = typename AT::t_float_1d("PairMultiLucyRX::fraction2",ntotal);
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairMultiLucyRXgetParams>(0,ntotal),*this);
}
const int inum = list->inum;
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;
computeLocalDensity();
// loop over neighbors of my atoms
EV_FLOAT ev;
if (neighflag == HALF) {
if (newton_pair) {
if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagPairMultiLucyRXCompute<HALF,1,1,TABSTYLE> >(0,inum),*this,ev);
else Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairMultiLucyRXCompute<HALF,1,0,TABSTYLE> >(0,inum),*this);
} else {
if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagPairMultiLucyRXCompute<HALF,0,1,TABSTYLE> >(0,inum),*this,ev);
else Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairMultiLucyRXCompute<HALF,0,0,TABSTYLE> >(0,inum),*this);
}
} else if (neighflag == HALFTHREAD) {
if (newton_pair) {
if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagPairMultiLucyRXCompute<HALFTHREAD,1,1,TABSTYLE> >(0,inum),*this,ev);
else Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairMultiLucyRXCompute<HALFTHREAD,1,0,TABSTYLE> >(0,inum),*this);
} else {
if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagPairMultiLucyRXCompute<HALFTHREAD,0,1,TABSTYLE> >(0,inum),*this,ev);
else Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairMultiLucyRXCompute<HALFTHREAD,0,0,TABSTYLE> >(0,inum),*this);
}
}
k_error_flag.template modify<DeviceType>();
k_error_flag.template sync<LMPHostType>();
if (k_error_flag.h_view() == 1)
error->one(FLERR,"Density < table inner cutoff");
else if (k_error_flag.h_view() == 2)
error->one(FLERR,"Density > table outer cutoff");
else if (k_error_flag.h_view() == 3)
error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy/rx");
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<DeviceType>();
k_eatom.template sync<LMPHostType>();
}
if (vflag_atom) {
k_vatom.template modify<DeviceType>();
k_vatom.template sync<LMPHostType>();
}
copymode = 0;
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairMultiLucyRXKokkos<DeviceType>::operator()(TagPairMultiLucyRXgetParams, const int &i) const {
getParams(i, d_fractionOld1[i], d_fractionOld2[i], d_fraction1[i], d_fraction2[i]);
}
template<class DeviceType>
template<int NEIGHFLAG, int NEWTON_PAIR, int EVFLAG, int TABSTYLE>
KOKKOS_INLINE_FUNCTION
void PairMultiLucyRXKokkos<DeviceType>::operator()(TagPairMultiLucyRXCompute<NEIGHFLAG,NEWTON_PAIR,EVFLAG,TABSTYLE>, const int &ii, EV_FLOAT& ev) const {
// The f array is atomic for Half/Thread neighbor style
Kokkos::View<F_FLOAT*[3], typename DAT::t_f_array::array_layout,DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > a_f = f;
int i,j,jj,inum,jnum,itype,jtype,itable;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,evdwlOld,fpair;
double rsq;
double fractionOld1_i,fractionOld1_j;
double fractionOld2_i,fractionOld2_j;
double fraction1_i;
double pi = MathConst::MY_PI;
double A_i, A_j;
double fraction_i,fraction_j;
int jtable;
int tlm1 = tablength - 1;
i = d_ilist[ii];
xtmp = x(i,0);
ytmp = x(i,1);
ztmp = x(i,2);
itype = type[i];
jnum = d_numneigh[i];
double fx_i = 0.0;
double fy_i = 0.0;
double fz_i = 0.0;
fractionOld1_i = d_fractionOld1[i];
fractionOld2_i = d_fractionOld2[i];
fraction1_i = d_fraction1[i];
for (jj = 0; jj < jnum; jj++) {
int j = d_neighbors(i,jj);
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)) { // optimize
fpair = 0.0;
fractionOld1_j = d_fractionOld1[j];
fractionOld2_j = d_fractionOld2[j];
//tb = &tables[tabindex[itype][jtype]];
const int tidx = d_table_const.tabindex(itype,jtype);
//if (rho[i]*rho[i] < tb->innersq || rho[j]*rho[j] < tb->innersq){
if (rho[i]*rho[i] < d_table_const.innersq(tidx) || rho[j]*rho[j] < d_table_const.innersq(tidx)){
k_error_flag.d_view() = 1;
}
if (TABSTYLE == LOOKUP) {
//itable = static_cast<int> (((rho[i]*rho[i]) - tb->innersq) * tb->invdelta);
itable = static_cast<int> (((rho[i]*rho[i]) - d_table_const.innersq(tidx)) * d_table_const.invdelta(tidx));
//jtable = static_cast<int> (((rho[j]*rho[j]) - tb->innersq) * tb->invdelta);
jtable = static_cast<int> (((rho[j]*rho[j]) - d_table_const.innersq(tidx)) * d_table_const.invdelta(tidx));
if (itable >= tlm1 || jtable >= tlm1){
k_error_flag.d_view() = 2;
}
//A_i = tb->f[itable];
A_i = d_table_const.f(tidx,itable);
//A_j = tb->f[jtable];
A_j = d_table_const.f(tidx,jtable);
const double rfactor = 1.0-sqrt(rsq/d_cutsq(itype,jtype));
fpair = 0.5*(A_i + A_j)*(4.0-3.0*rfactor)*rfactor*rfactor*rfactor;
fpair /= sqrt(rsq);
} else if (TABSTYLE == LINEAR) {
//itable = static_cast<int> ((rho[i]*rho[i] - tb->innersq) * tb->invdelta);
itable = static_cast<int> ((rho[i]*rho[i] - d_table_const.innersq(tidx)) * d_table_const.invdelta(tidx));
//jtable = static_cast<int> (((rho[j]*rho[j]) - tb->innersq) * tb->invdelta);
jtable = static_cast<int> ((rho[j]*rho[j] - d_table_const.innersq(tidx)) * d_table_const.invdelta(tidx));
if (itable >= tlm1 || jtable >= tlm1){
k_error_flag.d_view() = 2;
}
if(itable<0) itable=0;
if(itable>=tlm1) itable=tlm1;
if(jtable<0) jtable=0;
if(jtable>=tlm1)jtable=tlm1;
//fraction_i = (((rho[i]*rho[i]) - tb->rsq[itable]) * tb->invdelta);
fraction_i = (((rho[i]*rho[i]) - d_table_const.rsq(tidx,itable)) * d_table_const.invdelta(tidx));
//fraction_j = (((rho[j]*rho[j]) - tb->rsq[jtable]) * tb->invdelta);
fraction_j = (((rho[j]*rho[j]) - d_table_const.rsq(tidx,jtable)) * d_table_const.invdelta(tidx));
if(itable==0) fraction_i=0.0;
if(itable==tlm1) fraction_i=0.0;
if(jtable==0) fraction_j=0.0;
if(jtable==tlm1) fraction_j=0.0;
//A_i = tb->f[itable] + fraction_i*tb->df[itable];
A_i = d_table_const.f(tidx,itable) + fraction_i*d_table_const.df(tidx,itable);
//A_j = tb->f[jtable] + fraction_j*tb->df[jtable];
A_j = d_table_const.f(tidx,jtable) + fraction_j*d_table_const.df(tidx,jtable);
const double rfactor = 1.0-sqrt(rsq/d_cutsq(itype,jtype));
fpair = 0.5*(A_i + A_j)*(4.0-3.0*rfactor)*rfactor*rfactor*rfactor;
fpair /= sqrt(rsq);
} else k_error_flag.d_view() = 3;
if (isite1 == isite2) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair;
else fpair = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*fpair;
fx_i += delx*fpair;
fy_i += dely*fpair;
fz_i += delz*fpair;
if (NEWTON_PAIR || j < nlocal) {
a_f(j,0) -= delx*fpair;
a_f(j,1) -= dely*fpair;
a_f(j,2) -= delz*fpair;
}
//if (evflag) ev_tally(i,j,nlocal,newton_pair,0.0,0.0,fpair,delx,dely,delz);
if (EVFLAG) this->template ev_tally<NEIGHFLAG,NEWTON_PAIR>(ev,i,j,0.0,fpair,delx,dely,delz);
}
}
a_f(i,0) += fx_i;
a_f(i,1) += fy_i;
a_f(i,2) += fz_i;
//tb = &tables[tabindex[itype][itype]];
const int tidx = d_table_const.tabindex(itype,itype);
//itable = static_cast<int> (((rho[i]*rho[i]) - tb->innersq) * tb->invdelta);
itable = static_cast<int> (((rho[i]*rho[i]) - d_table_const.innersq(tidx)) * d_table_const.invdelta(tidx));
//if (TABSTYLE == LOOKUP) evdwl = tb->e[itable];
if (TABSTYLE == LOOKUP) evdwl = d_table_const.e(tidx,itable);
else if (TABSTYLE == LINEAR){
if (itable >= tlm1){
k_error_flag.d_view() = 2;
}
if(itable==0) fraction_i=0.0;
//else fraction_i = (((rho[i]*rho[i]) - tb->rsq[itable]) * tb->invdelta);
else fraction_i = (((rho[i]*rho[i]) - d_table_const.rsq(tidx,itable)) * d_table_const.invdelta(tidx));
//evdwl = tb->e[itable] + fraction_i*tb->de[itable];
evdwl = d_table_const.e(tidx,itable); + fraction_i*d_table_const.de(tidx,itable);
} else k_error_flag.d_view() = 3;
evdwl *=(pi*d_cutsq(itype,itype)*d_cutsq(itype,itype))/84.0;
evdwlOld = fractionOld1_i*evdwl;
evdwl = fraction1_i*evdwl;
uCG[i] += evdwlOld;
uCGnew[i] += evdwl;
evdwl = evdwlOld;
//if (evflag) ev_tally(0,0,nlocal,newton_pair,evdwl,0.0,0.0,0.0,0.0,0.0);
if (EVFLAG) ev.evdwl += evdwl;
}
template<class DeviceType>
template<int NEIGHFLAG, int NEWTON_PAIR, int EVFLAG, int TABSTYLE>
KOKKOS_INLINE_FUNCTION
void PairMultiLucyRXKokkos<DeviceType>::operator()(TagPairMultiLucyRXCompute<NEIGHFLAG,NEWTON_PAIR,EVFLAG,TABSTYLE>, const int &ii) const {
EV_FLOAT ev;
this->template operator()<NEIGHFLAG,NEWTON_PAIR,EVFLAG,TABSTYLE>(TagPairMultiLucyRXCompute<NEIGHFLAG,NEWTON_PAIR,EVFLAG,TABSTYLE>(), ii, ev);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairMultiLucyRXKokkos<DeviceType>::computeLocalDensity()
{
copymode = 1;
x = atomKK->k_x.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
rho = atomKK->k_rho.view<DeviceType>();
h_rho = atomKK->k_rho.h_view;
nlocal = atom->nlocal;
atomKK->sync(execution_space,X_MASK | TYPE_MASK | DPDRHO_MASK);
atomKK->modified(execution_space,DPDRHO_MASK);
const int inum = list->inum;
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;
const double pi = MathConst::MY_PI;
const bool newton_pair = force->newton_pair;
one_type = (atom->ntypes == 1);
// Special cut-off values for when there's only one type.
cutsq_type11 = cutsq[1][1];
rcut_type11 = sqrt(cutsq_type11);
factor_type11 = 84.0/(5.0*pi*rcut_type11*rcut_type11*rcut_type11);
// zero out density
int m = nlocal;
if (newton_pair) m += atom->nghost;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairMultiLucyRXZero>(0,m),*this);
// rho = density at each atom
// loop over neighbors of my atoms
if (neighflag == HALF) {
if (newton_pair)
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairMultiLucyRXComputeLocalDensity<HALF,1> >(0,inum),*this);
else
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairMultiLucyRXComputeLocalDensity<HALF,0> >(0,inum),*this);
} else if (neighflag == HALFTHREAD) {
if (newton_pair)
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairMultiLucyRXComputeLocalDensity<HALFTHREAD,1> >(0,inum),*this);
else
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairMultiLucyRXComputeLocalDensity<HALFTHREAD,0> >(0,inum),*this);
}
// communicate and sum densities (on the host)
if (newton_pair) {
atomKK->modified(execution_space,DPDRHO_MASK);
atomKK->sync(Host,DPDRHO_MASK);
comm->reverse_comm_pair(this);
atomKK->modified(Host,DPDRHO_MASK);
atomKK->sync(execution_space,DPDRHO_MASK);
}
comm->forward_comm_pair(this);
copymode = 0;
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairMultiLucyRXKokkos<DeviceType>::operator()(TagPairMultiLucyRXZero, const int &i) const {
rho[i] = 0.0;
}
template<class DeviceType>
template<int NEIGHFLAG, int NEWTON_PAIR>
KOKKOS_INLINE_FUNCTION
void PairMultiLucyRXKokkos<DeviceType>::operator()(TagPairMultiLucyRXComputeLocalDensity<NEIGHFLAG,NEWTON_PAIR>, const int &ii) const {
// The rho array is atomic for Half/Thread neighbor style
Kokkos::View<E_FLOAT*, typename DAT::t_efloat_1d::array_layout,DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > a_rho = rho;
const int i = d_ilist[ii];
const double xtmp = x(i,0);
const double ytmp = x(i,1);
const double ztmp = x(i,2);
double rho_i = rho[i];
const int itype = type[i];
const int jnum = d_numneigh[i];
const double pi = MathConst::MY_PI;
for (int jj = 0; jj < jnum; jj++){
const int j = (d_neighbors(i,jj) & NEIGHMASK);
const int jtype = type[j];
const double delx = xtmp - x(j,0);
const double dely = ytmp - x(j,1);
const double delz = ztmp - x(j,2);
const double rsq = delx*delx + dely*dely + delz*delz;
if (one_type) {
if (rsq < cutsq_type11) {
const double rcut = rcut_type11;
const double r_over_rcut = sqrt(rsq) / rcut;
const double tmpFactor = 1.0 - r_over_rcut;
const double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor;
const double factor = factor_type11*(1.0 + 1.5*r_over_rcut)*tmpFactor4;
rho_i += factor;
if (NEWTON_PAIR || j < nlocal)
a_rho[j] += factor;
} else if (rsq < d_cutsq(itype,jtype)) {
const double rcut = sqrt(d_cutsq(itype,jtype));
const double tmpFactor = 1.0-sqrt(rsq)/rcut;
const double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor;
const double factor = (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4;
rho_i += factor;
if (NEWTON_PAIR || j < nlocal)
a_rho[j] += factor;
}
}
}
a_rho[i] = rho_i;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairMultiLucyRXKokkos<DeviceType>::getParams(int id, double &fractionOld1, double &fractionOld2, double &fraction1, double &fraction2) const
{
double fractionOld, fraction;
double nTotal, nTotalOld;
nTotal = 0.0;
nTotalOld = 0.0;
for (int ispecies = 0; ispecies < nspecies; ispecies++){
nTotal += dvector(ispecies,id);
nTotalOld += dvector(ispecies+nspecies,id);
}
if (isOneFluid(isite1) == false){
fractionOld1 = dvector(isite1+nspecies,id)/nTotalOld;
fraction1 = dvector(isite1,id)/nTotal;
}
if (isOneFluid(isite2) == false){
fractionOld2 = dvector(isite2+nspecies,id)/nTotalOld;
fraction2 = dvector(isite2,id)/nTotal;
}
if (isOneFluid(isite1) || isOneFluid(isite2)){
fractionOld = 0.0;
fraction = 0.0;
for (int ispecies = 0; ispecies < nspecies; ispecies++){
if (isite1 == ispecies || isite2 == ispecies) continue;
fractionOld += dvector(ispecies+nspecies,id) / nTotalOld;
fraction += dvector(ispecies,id) / nTotal;
}
if (isOneFluid(isite1)){
fractionOld1 = fractionOld;
fraction1 = fraction;
}
if (isOneFluid(isite2)){
fractionOld2 = fractionOld;
fraction2 = fraction;
}
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
int PairMultiLucyRXKokkos<DeviceType>::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<DeviceType>();
iswap = iswap_in;
v_buf = buf.view<DeviceType>();
Kokkos::parallel_for(Kokkos::RangePolicy<LMPDeviceType, TagPairMultiLucyRXPackForwardComm>(0,n),*this);
DeviceType::fence();
return n;
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairMultiLucyRXKokkos<DeviceType>::operator()(TagPairMultiLucyRXPackForwardComm, const int &i) const {
int j = d_sendlist(iswap, i);
v_buf[i] = rho[j];
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairMultiLucyRXKokkos<DeviceType>::unpack_forward_comm_kokkos(int n, int first_in, DAT::tdual_xfloat_1d &buf)
{
first = first_in;
v_buf = buf.view<DeviceType>();
Kokkos::parallel_for(Kokkos::RangePolicy<LMPDeviceType, TagPairMultiLucyRXUnpackForwardComm>(0,n),*this);
DeviceType::fence();
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairMultiLucyRXKokkos<DeviceType>::operator()(TagPairMultiLucyRXUnpackForwardComm, const int &i) const {
rho[i + first] = v_buf[i];
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
int PairMultiLucyRXKokkos<DeviceType>::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = h_rho[j];
}
return m;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairMultiLucyRXKokkos<DeviceType>::unpack_forward_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) h_rho[i] = buf[m++];
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
int PairMultiLucyRXKokkos<DeviceType>::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) buf[m++] = h_rho[i];
return m;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairMultiLucyRXKokkos<DeviceType>::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
h_rho[j] += buf[m++];
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG, int NEWTON_PAIR>
KOKKOS_INLINE_FUNCTION
void PairMultiLucyRXKokkos<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
{
const int EFLAG = eflag;
const int VFLAG = vflag_either;
// The eatom and vatom arrays are atomic for Half/Thread neighbor style
Kokkos::View<E_FLOAT*, typename DAT::t_efloat_1d::array_layout,DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > v_eatom = k_eatom.view<DeviceType>();
Kokkos::View<F_FLOAT*[6], typename DAT::t_virial_array::array_layout,DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > v_vatom = k_vatom.view<DeviceType>();
if (EFLAG) {
if (eflag_atom) {
const E_FLOAT epairhalf = 0.5 * epair;
if (NEIGHFLAG!=FULL) {
if (NEWTON_PAIR || i < nlocal) v_eatom[i] += epairhalf;
if (NEWTON_PAIR || j < nlocal) v_eatom[j] += epairhalf;
} else {
v_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) {
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;
}
} else {
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;
}
}
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairMultiLucyRXKokkos<DeviceType>::create_kokkos_tables()
{
const int tlm1 = tablength-1;
memory->create_kokkos(d_table->innersq,h_table->innersq,ntables,"Table::innersq");
memory->create_kokkos(d_table->invdelta,h_table->invdelta,ntables,"Table::invdelta");
memory->create_kokkos(d_table->deltasq6,h_table->deltasq6,ntables,"Table::deltasq6");
if(tabstyle == LOOKUP) {
memory->create_kokkos(d_table->e,h_table->e,ntables,tlm1,"Table::e");
memory->create_kokkos(d_table->f,h_table->f,ntables,tlm1,"Table::f");
}
if(tabstyle == LINEAR) {
memory->create_kokkos(d_table->rsq,h_table->rsq,ntables,tablength,"Table::rsq");
memory->create_kokkos(d_table->e,h_table->e,ntables,tablength,"Table::e");
memory->create_kokkos(d_table->f,h_table->f,ntables,tablength,"Table::f");
memory->create_kokkos(d_table->de,h_table->de,ntables,tlm1,"Table::de");
memory->create_kokkos(d_table->df,h_table->df,ntables,tlm1,"Table::df");
}
for(int i=0; i < ntables; i++) {
Table* tb = &tables[i];
h_table->innersq[i] = tb->innersq;
h_table->invdelta[i] = tb->invdelta;
h_table->deltasq6[i] = tb->deltasq6;
for(int j = 0; j<h_table->rsq.dimension_1(); j++)
h_table->rsq(i,j) = tb->rsq[j];
for(int j = 0; j<h_table->drsq.dimension_1(); j++)
h_table->drsq(i,j) = tb->drsq[j];
for(int j = 0; j<h_table->e.dimension_1(); j++)
h_table->e(i,j) = tb->e[j];
for(int j = 0; j<h_table->de.dimension_1(); j++)
h_table->de(i,j) = tb->de[j];
for(int j = 0; j<h_table->f.dimension_1(); j++)
h_table->f(i,j) = tb->f[j];
for(int j = 0; j<h_table->df.dimension_1(); j++)
h_table->df(i,j) = tb->df[j];
for(int j = 0; j<h_table->e2.dimension_1(); j++)
h_table->e2(i,j) = tb->e2[j];
for(int j = 0; j<h_table->f2.dimension_1(); j++)
h_table->f2(i,j) = tb->f2[j];
}
Kokkos::deep_copy(d_table->innersq,h_table->innersq);
Kokkos::deep_copy(d_table->invdelta,h_table->invdelta);
Kokkos::deep_copy(d_table->deltasq6,h_table->deltasq6);
Kokkos::deep_copy(d_table->rsq,h_table->rsq);
Kokkos::deep_copy(d_table->drsq,h_table->drsq);
Kokkos::deep_copy(d_table->e,h_table->e);
Kokkos::deep_copy(d_table->de,h_table->de);
Kokkos::deep_copy(d_table->f,h_table->f);
Kokkos::deep_copy(d_table->df,h_table->df);
Kokkos::deep_copy(d_table->e2,h_table->e2);
Kokkos::deep_copy(d_table->f2,h_table->f2);
Kokkos::deep_copy(d_table->tabindex,h_table->tabindex);
d_table_const.innersq = d_table->innersq;
d_table_const.invdelta = d_table->invdelta;
d_table_const.deltasq6 = d_table->deltasq6;
d_table_const.rsq = d_table->rsq;
d_table_const.drsq = d_table->drsq;
d_table_const.e = d_table->e;
d_table_const.de = d_table->de;
d_table_const.f = d_table->f;
d_table_const.df = d_table->df;
d_table_const.e2 = d_table->e2;
d_table_const.f2 = d_table->f2;
Kokkos::deep_copy(d_table->cutsq,h_table->cutsq);
update_table = 0;
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
template<class DeviceType>
void PairMultiLucyRXKokkos<DeviceType>::allocate()
{
allocated = 1;
const int nt = atom->ntypes + 1;
memory->create(setflag,nt,nt,"pair:setflag");
memory->create_kokkos(d_table->cutsq,h_table->cutsq,cutsq,nt,nt,"pair:cutsq");
memory->create_kokkos(d_table->tabindex,h_table->tabindex,tabindex,nt,nt,"pair:tabindex");
d_table_const.cutsq = d_table->cutsq;
d_table_const.tabindex = d_table->tabindex;
memset(&setflag[0][0],0,nt*nt*sizeof(int));
memset(&cutsq[0][0],0,nt*nt*sizeof(double));
memset(&tabindex[0][0],0,nt*nt*sizeof(int));
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
template<class DeviceType>
void PairMultiLucyRXKokkos<DeviceType>::settings(int narg, char **arg)
{
if (narg < 2) error->all(FLERR,"Illegal pair_style command");
// new settings
if (strcmp(arg[0],"lookup") == 0) tabstyle = LOOKUP;
else if (strcmp(arg[0],"linear") == 0) tabstyle = LINEAR;
else error->all(FLERR,"Unknown table style in pair_style command");
tablength = force->inumeric(FLERR,arg[1]);
if (tablength < 2) error->all(FLERR,"Illegal number of pair table entries");
// delete old tables, since cannot just change settings
for (int m = 0; m < ntables; m++) free_table(&tables[m]);
memory->sfree(tables);
if (allocated) {
memory->destroy(setflag);
d_table_const.tabindex = d_table->tabindex = typename ArrayTypes<DeviceType>::t_int_2d();
h_table->tabindex = typename ArrayTypes<LMPHostType>::t_int_2d();
d_table_const.cutsq = d_table->cutsq = typename ArrayTypes<DeviceType>::t_ffloat_2d();
h_table->cutsq = typename ArrayTypes<LMPHostType>::t_ffloat_2d();
}
allocated = 0;
ntables = 0;
tables = NULL;
}
/* ---------------------------------------------------------------------- */
namespace LAMMPS_NS {
template class PairMultiLucyRXKokkos<LMPDeviceType>;
#ifdef KOKKOS_HAVE_CUDA
template class PairMultiLucyRXKokkos<LMPHostType>;
#endif
}
c4science · Help