Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78001165
pair_table_rx_kokkos.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Aug 17, 18:47
Size
45 KB
Mime Type
text/x-c++
Expires
Mon, Aug 19, 18:47 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
19964835
Attached To
rLAMMPS lammps
pair_table_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 author: Dan Ibanez (SNL)
------------------------------------------------------------------------- */
#include <mpi.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "pair_table_rx_kokkos.h"
#include "kokkos.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "memory.h"
#include "error.h"
#include "atom_masks.h"
#include "fix.h"
#include "kokkos_few.h"
#include <cassert>
using namespace LAMMPS_NS;
enum{NONE,RLINEAR,RSQ,BMP};
#ifdef DBL_EPSILON
#define MY_EPSILON (10.0*DBL_EPSILON)
#else
#define MY_EPSILON (10.0*2.220446049250313e-16)
#endif
#define OneFluidValue (-1)
#define isOneFluid(_site_) ( (_site_) == OneFluidValue )
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void getMixingWeights(
typename ArrayTypes<DeviceType>::t_float_2d_randomread dvector,
int nspecies,
int isite1, int isite2,
bool fractionalWeighting,
int id,
double &mixWtSite1old, double &mixWtSite2old,
double &mixWtSite1, double &mixWtSite2) {
double fractionOFAold, fractionOFA;
double fractionOld1, fraction1;
double fractionOld2, fraction2;
double nMoleculesOFAold, nMoleculesOFA;
double nMoleculesOld1, nMolecules1;
double nMoleculesOld2, nMolecules2;
double nTotal, nTotalOld;
nTotal = 0.0;
nTotalOld = 0.0;
assert(id >= 0);
assert(id < dvector.dimension_1());
for (int ispecies = 0; ispecies < nspecies; ++ispecies){
assert(ispecies < dvector.dimension_0());
nTotal += dvector(ispecies,id);
assert(ispecies+nspecies < dvector.dimension_0());
nTotalOld += dvector(ispecies+nspecies,id);
}
assert(isite1 >= 0);
assert(isite1 < nspecies);
assert(isite2 >= 0);
assert(isite2 < nspecies);
if (isOneFluid(isite1) == false){
nMoleculesOld1 = dvector(isite1+nspecies,id);
nMolecules1 = dvector(isite1,id);
fractionOld1 = nMoleculesOld1/nTotalOld;
fraction1 = nMolecules1/nTotal;
}
if (isOneFluid(isite2) == false){
nMoleculesOld2 = dvector(isite2+nspecies,id);
nMolecules2 = dvector(isite2,id);
fractionOld2 = nMoleculesOld2/nTotalOld;
fraction2 = nMolecules2/nTotal;
}
if (isOneFluid(isite1) || isOneFluid(isite2)){
nMoleculesOFAold = 0.0;
nMoleculesOFA = 0.0;
fractionOFAold = 0.0;
fractionOFA = 0.0;
for (int ispecies = 0; ispecies < nspecies; ispecies++){
if (isite1 == ispecies || isite2 == ispecies) continue;
nMoleculesOFAold += dvector(ispecies+nspecies,id);
nMoleculesOFA += dvector(ispecies,id);
fractionOFAold += dvector(ispecies+nspecies,id)/nTotalOld;
fractionOFA += dvector(ispecies,id)/nTotal;
}
if(isOneFluid(isite1)){
nMoleculesOld1 = 1.0-(nTotalOld-nMoleculesOFAold);
nMolecules1 = 1.0-(nTotal-nMoleculesOFA);
fractionOld1 = fractionOFAold;
fraction1 = fractionOFA;
}
if(isOneFluid(isite2)){
nMoleculesOld2 = 1.0-(nTotalOld-nMoleculesOFAold);
nMolecules2 = 1.0-(nTotal-nMoleculesOFA);
fractionOld2 = fractionOFAold;
fraction2 = fractionOFA;
}
}
if(fractionalWeighting){
mixWtSite1old = fractionOld1;
mixWtSite1 = fraction1;
mixWtSite2old = fractionOld2;
mixWtSite2 = fraction2;
} else {
mixWtSite1old = nMoleculesOld1;
mixWtSite1 = nMolecules1;
mixWtSite2old = nMoleculesOld2;
mixWtSite2 = nMolecules2;
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
PairTableRXKokkos<DeviceType>::PairTableRXKokkos(LAMMPS *lmp) : PairTable(lmp)
{
update_table = 0;
atomKK = (AtomKokkos *) atom;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = X_MASK | F_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK |
DVECTOR_MASK | UCG_MASK | UCGNEW_MASK;
datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK | UCG_MASK | UCGNEW_MASK;
h_table = new TableHost();
d_table = new TableDevice();
fractionalWeighting = true;
site1 = nullptr;
site2 = nullptr;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
PairTableRXKokkos<DeviceType>::~PairTableRXKokkos()
{
if (copymode) return;
delete [] site1;
delete [] site2;
memory->destroy_kokkos(k_eatom,eatom);
memory->destroy_kokkos(k_vatom,vatom);
if (allocated) {
memory->destroy_kokkos(d_table->cutsq, cutsq);
memory->destroy_kokkos(d_table->tabindex, tabindex);
}
delete h_table;
h_table = nullptr;
delete d_table;
d_table = nullptr;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairTableRXKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
{
if(update_table)
create_kokkos_tables();
if(tabstyle == LOOKUP)
compute_style<LOOKUP>(eflag_in,vflag_in);
if(tabstyle == LINEAR)
compute_style<LINEAR>(eflag_in,vflag_in);
if(tabstyle == SPLINE)
compute_style<SPLINE>(eflag_in,vflag_in);
if(tabstyle == BITMAP)
compute_style<BITMAP>(eflag_in,vflag_in);
}
KOKKOS_INLINE_FUNCTION static int sbmask(const int& j)
{
return j >> SBBITS & 3;
}
template <class DeviceType, int TABSTYLE>
KOKKOS_INLINE_FUNCTION
static F_FLOAT
compute_fpair(F_FLOAT rsq,
int itype, int jtype,
typename PairTableRXKokkos<DeviceType>::TableDeviceConst const& d_table_const
) {
Pair::union_int_float_t rsq_lookup;
double fpair;
const int tidx = d_table_const.tabindex(itype,jtype);
if (TABSTYLE == PairTable::LOOKUP) {
const int itable = static_cast<int> ((rsq - d_table_const.innersq(tidx)) * d_table_const.invdelta(tidx));
fpair = d_table_const.f(tidx,itable);
} else if (TABSTYLE == PairTable::LINEAR) {
const int itable = static_cast<int> ((rsq - d_table_const.innersq(tidx)) * d_table_const.invdelta(tidx));
const double fraction = (rsq - d_table_const.rsq(tidx,itable)) * d_table_const.invdelta(tidx);
fpair = d_table_const.f(tidx,itable) + fraction*d_table_const.df(tidx,itable);
} else if (TABSTYLE == PairTable::SPLINE) {
const int itable = static_cast<int> ((rsq - d_table_const.innersq(tidx)) * d_table_const.invdelta(tidx));
const double b = (rsq - d_table_const.rsq(tidx,itable)) * d_table_const.invdelta(tidx);
const double a = 1.0 - b;
fpair = a * d_table_const.f(tidx,itable) + b * d_table_const.f(tidx,itable+1) +
((a*a*a-a)*d_table_const.f2(tidx,itable) + (b*b*b-b)*d_table_const.f2(tidx,itable+1)) *
d_table_const.deltasq6(tidx);
} else {
rsq_lookup.f = rsq;
int itable = rsq_lookup.i & d_table_const.nmask(tidx);
itable >>= d_table_const.nshiftbits(tidx);
const double fraction = (rsq_lookup.f - d_table_const.rsq(tidx,itable)) * d_table_const.drsq(tidx,itable);
fpair = d_table_const.f(tidx,itable) + fraction*d_table_const.df(tidx,itable);
}
return fpair;
}
template<class DeviceType, int TABSTYLE>
KOKKOS_INLINE_FUNCTION
static F_FLOAT
compute_evdwl(
F_FLOAT rsq,
int itype, int jtype,
typename PairTableRXKokkos<DeviceType>::TableDeviceConst const& d_table_const
) {
double evdwl;
Pair::union_int_float_t rsq_lookup;
const int tidx = d_table_const.tabindex(itype,jtype);
if (TABSTYLE == PairTable::LOOKUP) {
const int itable = static_cast<int> ((rsq - d_table_const.innersq(tidx)) * d_table_const.invdelta(tidx));
evdwl = d_table_const.e(tidx,itable);
} else if (TABSTYLE == PairTable::LINEAR) {
const int itable = static_cast<int> ((rsq - d_table_const.innersq(tidx)) * d_table_const.invdelta(tidx));
const double fraction = (rsq - d_table_const.rsq(tidx,itable)) * d_table_const.invdelta(tidx);
evdwl = d_table_const.e(tidx,itable) + fraction*d_table_const.de(tidx,itable);
} else if (TABSTYLE == PairTable::SPLINE) {
const int itable = static_cast<int> ((rsq - d_table_const.innersq(tidx)) * d_table_const.invdelta(tidx));
const double b = (rsq - d_table_const.rsq(tidx,itable)) * d_table_const.invdelta(tidx);
const double a = 1.0 - b;
evdwl = a * d_table_const.e(tidx,itable) + b * d_table_const.e(tidx,itable+1) +
((a*a*a-a)*d_table_const.e2(tidx,itable) + (b*b*b-b)*d_table_const.e2(tidx,itable+1)) *
d_table_const.deltasq6(tidx);
} else {
rsq_lookup.f = rsq;
int itable = rsq_lookup.i & d_table_const.nmask(tidx);
itable >>= d_table_const.nshiftbits(tidx);
const double fraction = (rsq_lookup.f - d_table_const.rsq(tidx,itable)) * d_table_const.drsq(tidx,itable);
evdwl = d_table_const.e(tidx,itable) + fraction*d_table_const.de(tidx,itable);
}
return evdwl;
}
template<class DeviceType, int NEIGHFLAG, int TABSTYLE, int NEWTON_PAIR>
KOKKOS_INLINE_FUNCTION
void
ev_tally(
int eflag,
int eflag_atom,
int vflag,
int vflag_global,
int vflag_atom,
int nlocal,
int i, int j,
EV_FLOAT& ev,
F_FLOAT epair, F_FLOAT fpair,
F_FLOAT delx, F_FLOAT dely, F_FLOAT delz,
Kokkos::View<F_FLOAT*[6],
typename ArrayTypes<DeviceType>::t_virial_array::array_layout,
DeviceType,
Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > const& v_vatom,
Kokkos::View<E_FLOAT*,
typename ArrayTypes<DeviceType>::t_efloat_1d::array_layout,
DeviceType,
Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > const& v_eatom)
{
if (eflag) {
if (eflag_atom) {
auto 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) {
auto v0 = delx*delx*fpair;
auto v1 = dely*dely*fpair;
auto v2 = delz*delz*fpair;
auto v3 = delx*dely*fpair;
auto v4 = delx*delz*fpair;
auto v5 = dely*delz*fpair;
if (vflag_global) {
if (NEIGHFLAG!=FULL) {
if (NEWTON_PAIR) {
ev.v[0] += v0;
ev.v[1] += v1;
ev.v[2] += v2;
ev.v[3] += v3;
ev.v[4] += v4;
ev.v[5] += v5;
} else {
if (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 (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, int NEIGHFLAG, bool STACKPARAMS, int TABSTYLE,
int EVFLAG, int NEWTON_PAIR>
KOKKOS_INLINE_FUNCTION
static EV_FLOAT
compute_item(
int ii,
int nlocal,
typename ArrayTypes<DeviceType>::t_int_1d_const const& d_ilist,
typename ArrayTypes<DeviceType>::t_neighbors_2d_const const& d_neighbors,
typename ArrayTypes<DeviceType>::t_int_1d_const const& d_numneigh,
typename ArrayTypes<DeviceType>::t_x_array_randomread const& x,
typename ArrayTypes<DeviceType>::t_int_1d_randomread const& type,
Kokkos::View<double*, DeviceType> const& mixWtSite1old,
Kokkos::View<double*, DeviceType> const& mixWtSite2old,
Kokkos::View<double*, DeviceType> const& mixWtSite1,
Kokkos::View<double*, DeviceType> const& mixWtSite2,
Few<int, 4> const& special_lj,
Few<Few<F_FLOAT, MAX_TYPES_STACKPARAMS+1>, MAX_TYPES_STACKPARAMS+1> const& m_cutsq,
typename ArrayTypes<DeviceType>::t_ffloat_2d const& d_cutsq,
Kokkos::View<F_FLOAT*[3],
typename ArrayTypes<DeviceType>::t_f_array::array_layout,
DeviceType,
Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > const& f,
Kokkos::View<E_FLOAT*,
typename ArrayTypes<DeviceType>::t_efloat_1d::array_layout,
DeviceType,
Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > const& uCG,
Kokkos::View<E_FLOAT*,
typename ArrayTypes<DeviceType>::t_efloat_1d::array_layout,
DeviceType,
Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > const& uCGnew,
int isite1, int isite2,
typename PairTableRXKokkos<DeviceType>::TableDeviceConst const& d_table_const,
int eflag,
int eflag_atom,
int vflag,
int vflag_global,
int vflag_atom,
Kokkos::View<F_FLOAT*[6],
typename ArrayTypes<DeviceType>::t_virial_array::array_layout,
DeviceType,
Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > const& v_vatom,
Kokkos::View<E_FLOAT*,
typename ArrayTypes<DeviceType>::t_efloat_1d::array_layout,
DeviceType,
Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > const& v_eatom) {
EV_FLOAT ev;
auto i = d_ilist(ii);
auto xtmp = x(i,0);
auto ytmp = x(i,1);
auto ztmp = x(i,2);
auto itype = type(i);
auto jlist = NeighListKokkos<DeviceType>::static_neighbors_const(i,
d_neighbors, d_numneigh);
auto jnum = d_numneigh(i);
double uCG_i = 0.0;
double uCGnew_i = 0.0;
double fx_i = 0.0, fy_i = 0.0, fz_i = 0.0;
auto mixWtSite1old_i = mixWtSite1old(i);
auto mixWtSite2old_i = mixWtSite2old(i);
auto mixWtSite1_i = mixWtSite1(i);
auto mixWtSite2_i = mixWtSite2(i);
for (int jj = 0; jj < jnum; jj++) {
auto j = jlist(jj);
const F_FLOAT factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;
auto delx = xtmp - x(j,0);
auto dely = ytmp - x(j,1);
auto delz = ztmp - x(j,2);
auto rsq = delx*delx + dely*dely + delz*delz;
auto jtype = type(j);
if(rsq < (STACKPARAMS ? m_cutsq[itype][jtype] : d_cutsq(itype,jtype))) {
auto mixWtSite1old_j = mixWtSite1old(j);
auto mixWtSite2old_j = mixWtSite2old(j);
auto mixWtSite1_j = mixWtSite1(j);
auto mixWtSite2_j = mixWtSite2(j);
auto fpair = factor_lj * compute_fpair<DeviceType,TABSTYLE>(
rsq,itype,jtype,d_table_const);
if (isite1 == isite2) fpair *= sqrt(mixWtSite1old_i * mixWtSite2old_j);
else fpair *= (sqrt(mixWtSite1old_i * mixWtSite2old_j) +
sqrt(mixWtSite2old_i * mixWtSite1old_j));
fx_i += delx*fpair;
fy_i += dely*fpair;
fz_i += delz*fpair;
auto do_half = (NEIGHFLAG==HALF || NEIGHFLAG==HALFTHREAD) &&
(NEWTON_PAIR || j < nlocal);
if (do_half) {
f(j,0) -= delx*fpair;
f(j,1) -= dely*fpair;
f(j,2) -= delz*fpair;
}
auto evdwl = compute_evdwl<DeviceType,TABSTYLE>(
rsq,itype,jtype,d_table_const);
double evdwlOld;
if (isite1 == isite2) {
evdwlOld = sqrt(mixWtSite1old_i*mixWtSite2old_j)*evdwl;
evdwl = sqrt(mixWtSite1_i*mixWtSite2_j)*evdwl;
} else {
evdwlOld = (sqrt(mixWtSite1old_i*mixWtSite2old_j) +
sqrt(mixWtSite2old_i*mixWtSite1old_j))*evdwl;
evdwl = (sqrt(mixWtSite1_i*mixWtSite2_j) +
sqrt(mixWtSite2_i*mixWtSite1_j))*evdwl;
}
evdwlOld *= factor_lj;
evdwl *= factor_lj;
uCG_i += 0.5*evdwlOld;
if (do_half) uCG(j) += 0.5*evdwlOld;
uCGnew_i += 0.5*evdwl;
if (do_half) uCGnew(j) += 0.5*evdwl;
evdwl = evdwlOld;
ev.evdwl += (do_half ? 1.0 : 0.5)*evdwl;
if (EVFLAG) {
ev_tally<DeviceType,NEIGHFLAG,TABSTYLE,NEWTON_PAIR>(
eflag,eflag_atom,
vflag,vflag_global,vflag_atom,
nlocal,i,j,ev,evdwl,fpair,delx,dely,delz,
v_vatom, v_eatom);
}
}
}
uCG(i) += uCG_i;
uCGnew(i) += uCGnew_i;
f(i,0) += fx_i;
f(i,1) += fy_i;
f(i,2) += fz_i;
return ev;
}
template<class DeviceType, int NEIGHFLAG, bool STACKPARAMS, int TABSTYLE, int NEWTON_PAIR>
static void compute_all_items(
EV_FLOAT& ev,
int nlocal,
int inum,
typename ArrayTypes<DeviceType>::t_int_1d_const d_ilist,
typename ArrayTypes<DeviceType>::t_neighbors_2d_const d_neighbors,
typename ArrayTypes<DeviceType>::t_int_1d_const d_numneigh,
typename ArrayTypes<DeviceType>::t_x_array_randomread x,
typename ArrayTypes<DeviceType>::t_int_1d_randomread type,
Kokkos::View<double*, DeviceType> mixWtSite1old,
Kokkos::View<double*, DeviceType> mixWtSite2old,
Kokkos::View<double*, DeviceType> mixWtSite1,
Kokkos::View<double*, DeviceType> mixWtSite2,
Few<int, 4> special_lj,
Few<Few<F_FLOAT, MAX_TYPES_STACKPARAMS+1>, MAX_TYPES_STACKPARAMS+1> m_cutsq,
typename ArrayTypes<DeviceType>::t_ffloat_2d d_cutsq,
Kokkos::View<F_FLOAT*[3],
typename ArrayTypes<DeviceType>::t_f_array::array_layout,
DeviceType,
Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > f,
Kokkos::View<E_FLOAT*,
typename ArrayTypes<DeviceType>::t_efloat_1d::array_layout,
DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > uCG,
Kokkos::View<E_FLOAT*,
typename ArrayTypes<DeviceType>::t_efloat_1d::array_layout,
DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > uCGnew,
int isite1, int isite2,
typename PairTableRXKokkos<DeviceType>::TableDeviceConst d_table_const,
int eflag,
int eflag_atom,
int vflag,
int vflag_global,
int vflag_atom,
Kokkos::View<F_FLOAT*[6],
typename ArrayTypes<DeviceType>::t_virial_array::array_layout,
DeviceType,
Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > v_vatom,
Kokkos::View<E_FLOAT*,
typename ArrayTypes<DeviceType>::t_efloat_1d::array_layout,
DeviceType,
Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > v_eatom) {
if (eflag || vflag) {
Kokkos::parallel_reduce(inum,
LAMMPS_LAMBDA(int i, EV_FLOAT& energy_virial) {
energy_virial +=
compute_item<DeviceType,NEIGHFLAG,STACKPARAMS,TABSTYLE,1,NEWTON_PAIR>(
i, nlocal, d_ilist, d_neighbors, d_numneigh, x, type,
mixWtSite1old, mixWtSite2old, mixWtSite1, mixWtSite2,
special_lj, m_cutsq, d_cutsq, f, uCG, uCGnew, isite1, isite2,
d_table_const, eflag, eflag_atom,
vflag, vflag_global, vflag_atom, v_vatom, v_eatom);
}, ev);
} else {
Kokkos::parallel_for(inum,
LAMMPS_LAMBDA(int i) {
compute_item<DeviceType,NEIGHFLAG,STACKPARAMS,TABSTYLE,0,NEWTON_PAIR>(
i, nlocal, d_ilist, d_neighbors, d_numneigh, x, type,
mixWtSite1old, mixWtSite2old, mixWtSite1, mixWtSite2,
special_lj, m_cutsq, d_cutsq, f, uCG, uCGnew, isite1, isite2,
d_table_const, eflag, eflag_atom,
vflag, vflag_global, vflag_atom, v_vatom, v_eatom);
});
}
}
template<class DeviceType>
static void getAllMixingWeights(
int ntotal,
typename ArrayTypes<DeviceType>::t_float_2d_randomread dvector,
int nspecies,
int isite1, int isite2,
bool fractionalWeighting,
Kokkos::View<double*, DeviceType> mixWtSite1old,
Kokkos::View<double*, DeviceType> mixWtSite2old,
Kokkos::View<double*, DeviceType> mixWtSite1,
Kokkos::View<double*, DeviceType> mixWtSite2) {
Kokkos::parallel_for(ntotal,
LAMMPS_LAMBDA(int i) {
getMixingWeights<DeviceType>(dvector,nspecies,isite1,isite2,fractionalWeighting,
i, mixWtSite1old(i), mixWtSite2old(i), mixWtSite1(i), mixWtSite2(i));
});
}
template<class DeviceType>
template<int TABSTYLE>
void PairTableRXKokkos<DeviceType>::compute_style(int eflag_in, int vflag_in)
{
auto eflag = eflag_in;
auto vflag = vflag_in;
if (neighflag == FULL) no_virial_fdotr_compute = 1;
if (eflag || vflag) ev_setup(eflag,vflag,0);
else evflag = vflag_fdotr = 0;
if (eflag_atom) {
memory->destroy_kokkos(k_eatom,eatom);
memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
d_eatom = k_eatom.template view<DeviceType>();
}
if (vflag_atom) {
memory->destroy_kokkos(k_vatom,vatom);
memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
d_vatom = k_vatom.template view<DeviceType>();
}
atomKK->sync(execution_space,datamask_read);
if (eflag || vflag) atomKK->modified(execution_space,datamask_modify);
else atomKK->modified(execution_space,F_MASK);
x = atomKK->k_x.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
auto type = atomKK->k_type.view<DeviceType>();
auto uCG = atomKK->k_uCG.view<DeviceType>();
auto uCGnew = atomKK->k_uCGnew.view<DeviceType>();
auto nlocal = atom->nlocal;
Few<int, 4> special_lj_local;
special_lj_local[0] = force->special_lj[0];
special_lj_local[1] = force->special_lj[1];
special_lj_local[2] = force->special_lj[2];
special_lj_local[3] = force->special_lj[3];
auto newton_pair = force->newton_pair;
d_cutsq = d_table->cutsq;
// loop over neighbors of my atoms
const int ntotal = atom->nlocal + atom->nghost;
auto mixWtSite1old = Kokkos::View<double*, DeviceType>("PairTableRXKokkos::mixWtSite1old", ntotal);
auto mixWtSite2old = Kokkos::View<double*, DeviceType>("PairTableRXKokkos::mixWtSite2old", ntotal);
auto mixWtSite1 = Kokkos::View<double*, DeviceType>("PairTableRXKokkos::mixWtSite1", ntotal);
auto mixWtSite2 = Kokkos::View<double*, DeviceType>("PairTableRXKokkos::mixWtSite2", ntotal);
getAllMixingWeights(ntotal, atomKK->k_dvector.template view<DeviceType>(),
nspecies, isite1, isite2, fractionalWeighting,
mixWtSite1old, mixWtSite2old, mixWtSite1, mixWtSite2);
if (neighflag == N2) error->all(FLERR,"pair table/rx/kk can't handle N2 yet\n");
NeighListKokkos<DeviceType>* l =
dynamic_cast<NeighListKokkos<DeviceType>*>(list);
EV_FLOAT ev;
if(atom->ntypes > MAX_TYPES_STACKPARAMS) {
if (neighflag == HALFTHREAD) {
if (newton_pair) {
compute_all_items<DeviceType,HALFTHREAD,false,TABSTYLE,1>(
ev, nlocal, l->inum, l->d_ilist, l->d_neighbors, l->d_numneigh,
x, type, mixWtSite1old, mixWtSite2old, mixWtSite1, mixWtSite2,
special_lj_local, m_cutsq, d_cutsq, f, uCG, uCGnew, isite1, isite2,
d_table_const, eflag, eflag_atom,
vflag, vflag_global, vflag_atom, d_vatom, d_eatom);
} else {
compute_all_items<DeviceType,HALFTHREAD,false,TABSTYLE,0>(
ev, nlocal, l->inum, l->d_ilist, l->d_neighbors, l->d_numneigh,
x, type, mixWtSite1old, mixWtSite2old, mixWtSite1, mixWtSite2,
special_lj_local, m_cutsq, d_cutsq, f, uCG, uCGnew, isite1, isite2,
d_table_const, eflag, eflag_atom,
vflag, vflag_global, vflag_atom, d_vatom, d_eatom);
}
} else if (neighflag == HALF) {
if (newton_pair) {
compute_all_items<DeviceType,HALF,false,TABSTYLE,1>(
ev, nlocal, l->inum, l->d_ilist, l->d_neighbors, l->d_numneigh,
x, type, mixWtSite1old, mixWtSite2old, mixWtSite1, mixWtSite2,
special_lj_local, m_cutsq, d_cutsq, f, uCG, uCGnew, isite1, isite2,
d_table_const, eflag, eflag_atom,
vflag, vflag_global, vflag_atom, d_vatom, d_eatom);
} else {
compute_all_items<DeviceType,HALF,false,TABSTYLE,0>(
ev, nlocal, l->inum, l->d_ilist, l->d_neighbors, l->d_numneigh,
x, type, mixWtSite1old, mixWtSite2old, mixWtSite1, mixWtSite2,
special_lj_local, m_cutsq, d_cutsq, f, uCG, uCGnew, isite1, isite2,
d_table_const, eflag, eflag_atom,
vflag, vflag_global, vflag_atom, d_vatom, d_eatom);
}
} else if (neighflag == FULL) {
if (newton_pair) {
compute_all_items<DeviceType,FULL,false,TABSTYLE,1>(
ev, nlocal, l->inum, l->d_ilist, l->d_neighbors, l->d_numneigh,
x, type, mixWtSite1old, mixWtSite2old, mixWtSite1, mixWtSite2,
special_lj_local, m_cutsq, d_cutsq, f, uCG, uCGnew, isite1, isite2,
d_table_const, eflag, eflag_atom,
vflag, vflag_global, vflag_atom, d_vatom, d_eatom);
} else {
compute_all_items<DeviceType,FULL,false,TABSTYLE,0>(
ev, nlocal, l->inum, l->d_ilist, l->d_neighbors, l->d_numneigh,
x, type, mixWtSite1old, mixWtSite2old, mixWtSite1, mixWtSite2,
special_lj_local, m_cutsq, d_cutsq, f, uCG, uCGnew, isite1, isite2,
d_table_const, eflag, eflag_atom,
vflag, vflag_global, vflag_atom, d_vatom, d_eatom);
}
}
} else {
if (neighflag == HALFTHREAD) {
if (newton_pair) {
compute_all_items<DeviceType,HALFTHREAD,true,TABSTYLE,1>(
ev, nlocal, l->inum, l->d_ilist, l->d_neighbors, l->d_numneigh,
x, type, mixWtSite1old, mixWtSite2old, mixWtSite1, mixWtSite2,
special_lj_local, m_cutsq, d_cutsq, f, uCG, uCGnew, isite1, isite2,
d_table_const, eflag, eflag_atom,
vflag, vflag_global, vflag_atom, d_vatom, d_eatom);
} else {
compute_all_items<DeviceType,HALFTHREAD,true,TABSTYLE,0>(
ev, nlocal, l->inum, l->d_ilist, l->d_neighbors, l->d_numneigh,
x, type, mixWtSite1old, mixWtSite2old, mixWtSite1, mixWtSite2,
special_lj_local, m_cutsq, d_cutsq, f, uCG, uCGnew, isite1, isite2,
d_table_const, eflag, eflag_atom,
vflag, vflag_global, vflag_atom, d_vatom, d_eatom);
}
} else if (neighflag == HALF) {
if (newton_pair) {
compute_all_items<DeviceType,HALF,true,TABSTYLE,1>(
ev, nlocal, l->inum, l->d_ilist, l->d_neighbors, l->d_numneigh,
x, type, mixWtSite1old, mixWtSite2old, mixWtSite1, mixWtSite2,
special_lj_local, m_cutsq, d_cutsq, f, uCG, uCGnew, isite1, isite2,
d_table_const, eflag, eflag_atom,
vflag, vflag_global, vflag_atom, d_vatom, d_eatom);
} else {
compute_all_items<DeviceType,HALF,true,TABSTYLE,0>(
ev, nlocal, l->inum, l->d_ilist, l->d_neighbors, l->d_numneigh,
x, type, mixWtSite1old, mixWtSite2old, mixWtSite1, mixWtSite2,
special_lj_local, m_cutsq, d_cutsq, f, uCG, uCGnew, isite1, isite2,
d_table_const, eflag, eflag_atom,
vflag, vflag_global, vflag_atom, d_vatom, d_eatom);
}
} else if (neighflag == FULL) {
if (newton_pair) {
compute_all_items<DeviceType,FULL,true,TABSTYLE,1>(
ev, nlocal, l->inum, l->d_ilist, l->d_neighbors, l->d_numneigh,
x, type, mixWtSite1old, mixWtSite2old, mixWtSite1, mixWtSite2,
special_lj_local, m_cutsq, d_cutsq, f, uCG, uCGnew, isite1, isite2,
d_table_const, eflag, eflag_atom,
vflag, vflag_global, vflag_atom, d_vatom, d_eatom);
} else {
compute_all_items<DeviceType,FULL,true,TABSTYLE,0>(
ev, nlocal, l->inum, l->d_ilist, l->d_neighbors, l->d_numneigh,
x, type, mixWtSite1old, mixWtSite2old, mixWtSite1, mixWtSite2,
special_lj_local, m_cutsq, d_cutsq, f, uCG, uCGnew, isite1, isite2,
d_table_const, eflag, eflag_atom,
vflag, vflag_global, vflag_atom, d_vatom, d_eatom);
}
}
}
if (eflag) 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>();
}
}
template<class DeviceType>
void PairTableRXKokkos<DeviceType>::create_kokkos_tables()
{
const int tlm1 = tablength-1;
memory->create_kokkos(d_table->nshiftbits,h_table->nshiftbits,ntables,"Table::nshiftbits");
memory->create_kokkos(d_table->nmask,h_table->nmask,ntables,"Table::nmask");
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");
}
if(tabstyle == SPLINE) {
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->e2,h_table->e2,ntables,tablength,"Table::e2");
memory->create_kokkos(d_table->f2,h_table->f2,ntables,tablength,"Table::f2");
}
if(tabstyle == BITMAP) {
int ntable = 1 << tablength;
memory->create_kokkos(d_table->rsq,h_table->rsq,ntables,ntable,"Table::rsq");
memory->create_kokkos(d_table->e,h_table->e,ntables,ntable,"Table::e");
memory->create_kokkos(d_table->f,h_table->f,ntables,ntable,"Table::f");
memory->create_kokkos(d_table->de,h_table->de,ntables,ntable,"Table::de");
memory->create_kokkos(d_table->df,h_table->df,ntables,ntable,"Table::df");
memory->create_kokkos(d_table->drsq,h_table->drsq,ntables,ntable,"Table::drsq");
}
for(int i=0; i < ntables; i++) {
Table* tb = &tables[i];
h_table->nshiftbits[i] = tb->nshiftbits;
h_table->nmask[i] = tb->nmask;
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->nshiftbits,h_table->nshiftbits);
d_table_const.nshiftbits = d_table->nshiftbits;
Kokkos::deep_copy(d_table->nmask,h_table->nmask);
d_table_const.nmask = d_table->nmask;
Kokkos::deep_copy(d_table->innersq,h_table->innersq);
d_table_const.innersq = d_table->innersq;
Kokkos::deep_copy(d_table->invdelta,h_table->invdelta);
d_table_const.invdelta = d_table->invdelta;
Kokkos::deep_copy(d_table->deltasq6,h_table->deltasq6);
d_table_const.deltasq6 = d_table->deltasq6;
if(tabstyle == LOOKUP) {
Kokkos::deep_copy(d_table->e,h_table->e);
d_table_const.e = d_table->e;
Kokkos::deep_copy(d_table->f,h_table->f);
d_table_const.f = d_table->f;
}
if(tabstyle == LINEAR) {
Kokkos::deep_copy(d_table->rsq,h_table->rsq);
d_table_const.rsq = d_table->rsq;
Kokkos::deep_copy(d_table->e,h_table->e);
d_table_const.e = d_table->e;
Kokkos::deep_copy(d_table->f,h_table->f);
d_table_const.f = d_table->f;
Kokkos::deep_copy(d_table->de,h_table->de);
d_table_const.de = d_table->de;
Kokkos::deep_copy(d_table->df,h_table->df);
d_table_const.df = d_table->df;
}
if(tabstyle == SPLINE) {
Kokkos::deep_copy(d_table->rsq,h_table->rsq);
d_table_const.rsq = d_table->rsq;
Kokkos::deep_copy(d_table->e,h_table->e);
d_table_const.e = d_table->e;
Kokkos::deep_copy(d_table->f,h_table->f);
d_table_const.f = d_table->f;
Kokkos::deep_copy(d_table->e2,h_table->e2);
d_table_const.e2 = d_table->e2;
Kokkos::deep_copy(d_table->f2,h_table->f2);
d_table_const.f2 = d_table->f2;
}
if(tabstyle == BITMAP) {
Kokkos::deep_copy(d_table->rsq,h_table->rsq);
d_table_const.rsq = d_table->rsq;
Kokkos::deep_copy(d_table->e,h_table->e);
d_table_const.e = d_table->e;
Kokkos::deep_copy(d_table->f,h_table->f);
d_table_const.f = d_table->f;
Kokkos::deep_copy(d_table->de,h_table->de);
d_table_const.de = d_table->de;
Kokkos::deep_copy(d_table->df,h_table->df);
d_table_const.df = d_table->df;
Kokkos::deep_copy(d_table->drsq,h_table->drsq);
d_table_const.drsq = d_table->drsq;
}
Kokkos::deep_copy(d_table->cutsq,h_table->cutsq);
d_table_const.cutsq = d_table->cutsq;
Kokkos::deep_copy(d_table->tabindex,h_table->tabindex);
d_table_const.tabindex = d_table->tabindex;
update_table = 0;
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
template<class DeviceType>
void PairTableRXKokkos<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 PairTableRXKokkos<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 if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE;
else if (strcmp(arg[0],"bitmap") == 0) tabstyle = BITMAP;
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");
// optional keywords
// assert the tabulation is compatible with a specific long-range solver
int iarg = 2;
while (iarg < narg) {
if (strcmp(arg[iarg],"ewald") == 0) ewaldflag = 1;
else if (strcmp(arg[iarg],"pppm") == 0) pppmflag = 1;
else if (strcmp(arg[iarg],"msm") == 0) msmflag = 1;
else if (strcmp(arg[iarg],"dispersion") == 0) dispersionflag = 1;
else if (strcmp(arg[iarg],"tip4p") == 0) tip4pflag = 1;
else if (strcmp(arg[iarg],"fractional") == 0) fractionalWeighting = true;
else if (strcmp(arg[iarg],"molecular") == 0) fractionalWeighting = false;
else error->all(FLERR,"Illegal pair_style command");
iarg++;
}
// delete old tables, since cannot just change settings
for (int m = 0; m < ntables; m++) free_table(&tables[m]);
memory->sfree(tables);
ntables = 0;
tables = NULL;
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;
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
template<class DeviceType>
void PairTableRXKokkos<DeviceType>::coeff(int narg, char **arg)
{
if (narg != 6 && narg != 7) error->all(FLERR,"Illegal pair_coeff command");
if (!allocated) allocate();
bool rx_flag = false;
for (int i = 0; i < modify->nfix; i++)
if (strncmp(modify->fix[i]->style,"rx",2) == 0) rx_flag = true;
if (!rx_flag) error->all(FLERR,"PairTableRX requires a fix rx command.");
int ilo,ihi,jlo,jhi;
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
int me;
MPI_Comm_rank(world,&me);
tables = (Table *)
memory->srealloc(tables,(ntables+1)*sizeof(Table),"pair:tables");
Table *tb = &tables[ntables];
null_table(tb);
if (me == 0) read_table(tb,arg[2],arg[3]);
bcast_table(tb);
nspecies = atom->nspecies_dpd;
if(nspecies==0) error->all(FLERR,"There are no rx species specified.");
int n;
n = strlen(arg[4]) + 1;
site1 = new char[n];
strcpy(site1,arg[4]);
int ispecies;
for (ispecies = 0; ispecies < nspecies; ispecies++){
if (strcmp(site1,&atom->dname[ispecies][0]) == 0) break;
}
if (ispecies == nspecies && strcmp(site1,"1fluid") != 0)
error->all(FLERR,"Site1 name not recognized in pair coefficients");
n = strlen(arg[5]) + 1;
site2 = new char[n];
strcpy(site2,arg[5]);
for (ispecies = 0; ispecies < nspecies; ispecies++){
if (strcmp(site2,&atom->dname[ispecies][0]) == 0) break;
}
if (ispecies == nspecies && strcmp(site2,"1fluid") != 0)
error->all(FLERR,"Site2 name not recognized in pair coefficients");
// set table cutoff
if (narg == 7) tb->cut = force->numeric(FLERR,arg[6]);
else if (tb->rflag) tb->cut = tb->rhi;
else tb->cut = tb->rfile[tb->ninput-1];
// error check on table parameters
// insure cutoff is within table
// for BITMAP tables, file values can be in non-ascending order
if (tb->ninput <= 1) error->one(FLERR,"Invalid pair table length");
double rlo,rhi;
if (tb->rflag == 0) {
rlo = tb->rfile[0];
rhi = tb->rfile[tb->ninput-1];
} else {
rlo = tb->rlo;
rhi = tb->rhi;
}
if (tb->cut <= rlo || tb->cut > rhi)
error->all(FLERR,"Invalid pair table cutoff");
if (rlo <= 0.0) error->all(FLERR,"Invalid pair table cutoff");
// match = 1 if don't need to spline read-in tables
// this is only the case if r values needed by final tables
// exactly match r values read from file
// for tabstyle SPLINE, always need to build spline tables
tb->match = 0;
if (tabstyle == LINEAR && tb->ninput == tablength &&
tb->rflag == RSQ && tb->rhi == tb->cut) tb->match = 1;
if (tabstyle == BITMAP && tb->ninput == 1 << tablength &&
tb->rflag == BMP && tb->rhi == tb->cut) tb->match = 1;
if (tb->rflag == BMP && tb->match == 0)
error->all(FLERR,"Bitmapped table in file does not match requested table");
// spline read-in values and compute r,e,f vectors within table
if (tb->match == 0) spline_table(tb);
compute_table(tb);
// store ptr to table in tabindex
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
tabindex[i][j] = ntables;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Illegal pair_coeff command");
ntables++;
{
if ( strcmp(site1,"1fluid") == 0 )
isite1 = OneFluidValue;
else {
isite1 = nspecies;
for (int k = 0; k < nspecies; k++){
if (strcmp(site1, atom->dname[k]) == 0){
isite1 = k;
break;
}
}
if (isite1 == nspecies) error->all(FLERR,"isite1 == nspecies");
}
if ( strcmp(site2,"1fluid") == 0 )
isite2 = OneFluidValue;
else {
isite2 = nspecies;
for (int k = 0; k < nspecies; k++){
if (strcmp(site2, atom->dname[k]) == 0){
isite2 = ispecies;
break;
}
}
if (isite2 == nspecies)
error->all(FLERR,"isite2 == nspecies");
}
}
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
template<class DeviceType>
double PairTableRXKokkos<DeviceType>::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
tabindex[j][i] = tabindex[i][j];
if(i<MAX_TYPES_STACKPARAMS+1 && j<MAX_TYPES_STACKPARAMS+1) {
m_cutsq[j][i] = m_cutsq[i][j] = tables[tabindex[i][j]].cut*tables[tabindex[i][j]].cut;
}
return tables[tabindex[i][j]].cut;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
double PairTableRXKokkos<DeviceType>::single(int i, int j, int itype, int jtype, double rsq,
double factor_coul, double factor_lj,
double &fforce)
{
int itable;
double fraction,value,a,b,phi;
int tlm1 = tablength - 1;
Table *tb = &tables[tabindex[itype][jtype]];
double mixWtSite1_i, mixWtSite1_j;
double mixWtSite2_i, mixWtSite2_j;
double mixWtSite1old_i, mixWtSite1old_j;
double mixWtSite2old_i, mixWtSite2old_j;
fraction = 0.0;
a = 0.0;
b = 0.0;
atomKK->k_dvector.template sync<LMPHostType>();
typename ArrayTypes<LMPHostType>::t_float_2d_randomread h_dvector =
atomKK->k_dvector.view<LMPHostType>();
getMixingWeights<LMPHostType>(h_dvector,
nspecies, isite1, isite2, fractionalWeighting,
i,mixWtSite1old_i,mixWtSite2old_i,
mixWtSite1_i,mixWtSite2_i);
getMixingWeights<LMPHostType>(h_dvector,
nspecies, isite1, isite2, fractionalWeighting,
j,mixWtSite1old_j,mixWtSite2old_j,
mixWtSite1_j,mixWtSite2_j);
if (rsq < tb->innersq) error->one(FLERR,"Pair distance < table inner cutoff");
if (tabstyle == LOOKUP) {
itable = static_cast<int> ((rsq-tb->innersq) * tb->invdelta);
if (itable >= tlm1) error->one(FLERR,"Pair distance > table outer cutoff");
fforce = factor_lj * tb->f[itable];
} else if (tabstyle == LINEAR) {
itable = static_cast<int> ((rsq-tb->innersq) * tb->invdelta);
if (itable >= tlm1) error->one(FLERR,"Pair distance > table outer cutoff");
fraction = (rsq - tb->rsq[itable]) * tb->invdelta;
value = tb->f[itable] + fraction*tb->df[itable];
fforce = factor_lj * value;
} else if (tabstyle == SPLINE) {
itable = static_cast<int> ((rsq-tb->innersq) * tb->invdelta);
if (itable >= tlm1) error->one(FLERR,"Pair distance > table outer cutoff");
b = (rsq - tb->rsq[itable]) * tb->invdelta;
a = 1.0 - b;
value = a * tb->f[itable] + b * tb->f[itable+1] +
((a*a*a-a)*tb->f2[itable] + (b*b*b-b)*tb->f2[itable+1]) *
tb->deltasq6;
fforce = factor_lj * value;
} else {
Pair::union_int_float_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & tb->nmask;
itable >>= tb->nshiftbits;
fraction = (rsq_lookup.f - tb->rsq[itable]) * tb->drsq[itable];
value = tb->f[itable] + fraction*tb->df[itable];
fforce = factor_lj * value;
}
if (isite1 == isite2) fforce = sqrt(mixWtSite1_i*mixWtSite2_j)*fforce;
else fforce = (sqrt(mixWtSite1_i*mixWtSite2_j) + sqrt(mixWtSite2_i*mixWtSite1_j))*fforce;
if (tabstyle == LOOKUP)
phi = tb->e[itable];
else if (tabstyle == LINEAR || tabstyle == BITMAP)
phi = tb->e[itable] + fraction*tb->de[itable];
else
phi = a * tb->e[itable] + b * tb->e[itable+1] +
((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) * tb->deltasq6;
if (isite1 == isite2) phi = sqrt(mixWtSite1_i*mixWtSite2_j)*phi;
else phi = (sqrt(mixWtSite1_i*mixWtSite2_j) + sqrt(mixWtSite2_i*mixWtSite1_j))*phi;
return factor_lj*phi;
}
/* ----------------------------------------------------------------------
compute r,e,f vectors from splined values
------------------------------------------------------------------------- */
template<class DeviceType>
void PairTableRXKokkos<DeviceType>::compute_table(Table *tb)
{
update_table = 1;
PairTable::compute_table(tb);
}
template<class DeviceType>
void PairTableRXKokkos<DeviceType>::init_style()
{
neighbor->request(this,instance_me);
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 if (neighflag == N2) {
neighbor->requests[irequest]->full = 0;
neighbor->requests[irequest]->half = 0;
} else {
error->all(FLERR,"Cannot use chosen neighbor list style with lj/cut/kk");
}
}
template<class DeviceType>
void PairTableRXKokkos<DeviceType>::cleanup_copy() {
// WHY needed: this prevents parent copy from deallocating any arrays
allocated = 0;
cutsq = NULL;
eatom = NULL;
vatom = NULL;
h_table=NULL; d_table=NULL;
}
namespace LAMMPS_NS {
template class PairTableRXKokkos<LMPDeviceType>;
#ifdef KOKKOS_HAVE_CUDA
template class PairTableRXKokkos<LMPHostType>;
#endif
}
Event Timeline
Log In to Comment