From 6929603eef93cf633fd8e731406782ccf68dcdd3 Mon Sep 17 00:00:00 2001 From: Anders Hafreager Date: Tue, 25 Oct 2016 14:15:34 +0200 Subject: [PATCH 1/8] Added KOKKOS vashishta (cherry picked from commit 5edc474bf0be574ddba96d00bb63894edf400ddb) --- src/KOKKOS/Install.sh | 2 + src/KOKKOS/pair_vashishta_kokkos.cpp | 915 +++++++++++++++++++++++++++ src/KOKKOS/pair_vashishta_kokkos.h | 154 +++++ 3 files changed, 1071 insertions(+) create mode 100644 src/KOKKOS/pair_vashishta_kokkos.cpp create mode 100644 src/KOKKOS/pair_vashishta_kokkos.h diff --git a/src/KOKKOS/Install.sh b/src/KOKKOS/Install.sh index c45be8c973..fdc87b8302 100644 --- a/src/KOKKOS/Install.sh +++ b/src/KOKKOS/Install.sh @@ -163,6 +163,8 @@ action pair_reax_c_kokkos.cpp pair_reax_c.cpp action pair_reax_c_kokkos.h pair_reax_c.h action pair_sw_kokkos.cpp pair_sw.cpp action pair_sw_kokkos.h pair_sw.h +action pair_vashishta_kokkos.cpp pair_vashishta.cpp +action pair_vashishta_kokkos.h pair_vashishta.h action pair_table_kokkos.cpp action pair_table_kokkos.h action pair_tersoff_kokkos.cpp pair_tersoff.cpp diff --git a/src/KOKKOS/pair_vashishta_kokkos.cpp b/src/KOKKOS/pair_vashishta_kokkos.cpp new file mode 100644 index 0000000000..b5fe84faf6 --- /dev/null +++ b/src/KOKKOS/pair_vashishta_kokkos.cpp @@ -0,0 +1,915 @@ +/* ---------------------------------------------------------------------- + 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: Stan Moore (SNL) + Anders Hafreager (UiO), andershaf@gmail.com +------------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include "pair_vashishta_kokkos.h" +#include "kokkos.h" +#include "pair_kokkos.h" +#include "atom_kokkos.h" +#include "neighbor.h" +#include "neigh_request.h" +#include "force.h" +#include "comm.h" +#include "memory.h" +#include "neighbor.h" +#include "neigh_list_kokkos.h" +#include "memory.h" +#include "error.h" +#include "atom_masks.h" +#include "math_const.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define MAXLINE 1024 +#define DELTA 4 + +/* ---------------------------------------------------------------------- */ + +template +PairVashishtaKokkos::PairVashishtaKokkos(LAMMPS *lmp) : PairVashishta(lmp) +{ + respa_enable = 0; + + atomKK = (AtomKokkos *) atom; + execution_space = ExecutionSpaceFromDevice::space; + datamask_read = X_MASK | F_MASK | TAG_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK; + datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK; +} + +/* ---------------------------------------------------------------------- + check if allocated, since class can be destructed when incomplete +------------------------------------------------------------------------- */ + +template +PairVashishtaKokkos::~PairVashishtaKokkos() +{ + if (!copymode) { + memory->destroy_kokkos(k_eatom,eatom); + memory->destroy_kokkos(k_vatom,vatom); + eatom = NULL; + vatom = NULL; + } +} + +/* ---------------------------------------------------------------------- */ + +template +void PairVashishtaKokkos::compute(int eflag_in, int vflag_in) +{ + eflag = eflag_in; + vflag = vflag_in; + + 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; + } + + 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(); + f = atomKK->k_f.view(); + tag = atomKK->k_tag.view(); + type = atomKK->k_type.view(); + nlocal = atom->nlocal; + newton_pair = force->newton_pair; + nall = atom->nlocal + atom->nghost; + + const int inum = list->inum; + const int ignum = inum + list->gnum; + NeighListKokkos* k_list = static_cast*>(list); + d_ilist = k_list->d_ilist; + d_numneigh = k_list->d_numneigh; + d_neighbors = k_list->d_neighbors; + + k_list->clean_copy(); + copymode = 1; + + EV_FLOAT ev; + EV_FLOAT ev_all; + + // loop over neighbor list of my atoms + + if (neighflag == HALF) { + if (evflag) + Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + else + Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + ev_all += ev; + } else if (neighflag == HALFTHREAD) { + if (evflag) + Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + else + Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + ev_all += ev; + } else if (neighflag == FULL) { + if (evflag) + Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,inum),*this,ev); + else + Kokkos::parallel_for(Kokkos::RangePolicy >(0,inum),*this); + ev_all += ev; + + if (evflag) + Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,ignum),*this,ev); + else + Kokkos::parallel_for(Kokkos::RangePolicy >(0,ignum),*this); + ev_all += ev; + } + + if (eflag_global) eng_vdwl += ev_all.evdwl; + if (vflag_global) { + virial[0] += ev_all.v[0]; + virial[1] += ev_all.v[1]; + virial[2] += ev_all.v[2]; + virial[3] += ev_all.v[3]; + virial[4] += ev_all.v[4]; + virial[5] += ev_all.v[5]; + } + + if (vflag_fdotr) pair_virial_fdotr_compute(this); + + if (eflag_atom) { + k_eatom.template modify(); + k_eatom.template sync(); + } + + if (vflag_atom) { + k_vatom.template modify(); + k_vatom.template sync(); + } + + copymode = 0; +} + +template +template +KOKKOS_INLINE_FUNCTION +void PairVashishtaKokkos::operator()(TagPairVashishtaComputeHalf, const int &ii, EV_FLOAT& ev) const { + + // The f array is atomic + + Kokkos::View::value> > a_f = f; + + F_FLOAT delr1[3],delr2[3],fj[3],fk[3]; + F_FLOAT evdwl = 0.0; + F_FLOAT fpair = 0.0; + + const int i = d_ilist[ii]; + const tagint itag = tag[i]; + const int itype = d_map[type[i]]; + const X_FLOAT xtmp = x(i,0); + const X_FLOAT ytmp = x(i,1); + const X_FLOAT ztmp = x(i,2); + + // two-body interactions, skip half of them + + const int jnum = d_numneigh[i]; + + F_FLOAT fxtmpi = 0.0; + F_FLOAT fytmpi = 0.0; + F_FLOAT fztmpi = 0.0; + + for (int jj = 0; jj < jnum; jj++) { + int j = d_neighbors(i,jj); + j &= NEIGHMASK; + const tagint jtag = tag[j]; + + if (itag > jtag) { + if ((itag+jtag) % 2 == 0) continue; + } else if (itag < jtag) { + if ((itag+jtag) % 2 == 1) continue; + } else { + if (x(j,2) < ztmp) continue; + if (x(j,2) == ztmp && x(j,1) < ytmp) continue; + if (x(j,2) == ztmp && x(j,1) == ytmp && x(j,0) < xtmp) continue; + } + + const int jtype = d_map[type[j]]; + + const X_FLOAT delx = xtmp - x(j,0); + const X_FLOAT dely = ytmp - x(j,1); + const X_FLOAT delz = ztmp - x(j,2); + const F_FLOAT rsq = delx*delx + dely*dely + delz*delz; + + const int ijparam = d_elem2param(itype,jtype,jtype); + if (rsq >= d_params[ijparam].cutsq) continue; + + twobody(d_params[ijparam],rsq,fpair,eflag,evdwl); + + fxtmpi += delx*fpair; + fytmpi += dely*fpair; + fztmpi += delz*fpair; + a_f(j,0) -= delx*fpair; + a_f(j,1) -= dely*fpair; + a_f(j,2) -= delz*fpair; + + if (EVFLAG) { + if (eflag) ev.evdwl += evdwl; + if (vflag_either || eflag_atom) this->template ev_tally(ev,i,j,evdwl,fpair,delx,dely,delz); + } + } + + const int jnumm1 = jnum - 1; + + for (int jj = 0; jj < jnumm1; jj++) { + int j = d_neighbors(i,jj); + j &= NEIGHMASK; + const int jtype = d_map[type[j]]; + const int ijparam = d_elem2param(itype,jtype,jtype); + delr1[0] = x(j,0) - xtmp; + delr1[1] = x(j,1) - ytmp; + delr1[2] = x(j,2) - ztmp; + const F_FLOAT rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; + if (rsq1 >= d_params[ijparam].cutsq2) continue; + + F_FLOAT fxtmpj = 0.0; + F_FLOAT fytmpj = 0.0; + F_FLOAT fztmpj = 0.0; + + for (int kk = jj+1; kk < jnum; kk++) { + int k = d_neighbors(i,kk); + k &= NEIGHMASK; + const int ktype = d_map[type[k]]; + const int ikparam = d_elem2param(itype,ktype,ktype); + const int ijkparam = d_elem2param(itype,jtype,ktype); + + delr2[0] = x(k,0) - xtmp; + delr2[1] = x(k,1) - ytmp; + delr2[2] = x(k,2) - ztmp; + const F_FLOAT rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; + + if (rsq2 >= d_params[ikparam].cutsq2) continue; + + threebody(d_params[ijparam],d_params[ikparam],d_params[ijkparam], + rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl); + + fxtmpi -= fj[0] + fk[0]; + fytmpi -= fj[1] + fk[1]; + fztmpi -= fj[2] + fk[2]; + fxtmpj += fj[0]; + fytmpj += fj[1]; + fztmpj += fj[2]; + a_f(k,0) += fk[0]; + a_f(k,1) += fk[1]; + a_f(k,2) += fk[2]; + + if (EVFLAG) { + if (eflag) ev.evdwl += evdwl; + if (vflag_either || eflag_atom) this->template ev_tally3(ev,i,j,k,evdwl,0.0,fj,fk,delr1,delr2); + } + } + + a_f(j,0) += fxtmpj; + a_f(j,1) += fytmpj; + a_f(j,2) += fztmpj; + } + + a_f(i,0) += fxtmpi; + a_f(i,1) += fytmpi; + a_f(i,2) += fztmpi; +} + +template +template +KOKKOS_INLINE_FUNCTION +void PairVashishtaKokkos::operator()(TagPairVashishtaComputeHalf, const int &ii) const { + EV_FLOAT ev; + this->template operator()(TagPairVashishtaComputeHalf(), ii, ev); +} + +/* ---------------------------------------------------------------------- */ + +template +template +KOKKOS_INLINE_FUNCTION +void PairVashishtaKokkos::operator()(TagPairVashishtaComputeFullA, const int &ii, EV_FLOAT& ev) const { + + F_FLOAT delr1[3],delr2[3],fj[3],fk[3]; + F_FLOAT evdwl = 0.0; + F_FLOAT fpair = 0.0; + + const int i = d_ilist[ii]; + + const tagint itag = tag[i]; + const int itype = d_map[type[i]]; + const X_FLOAT xtmp = x(i,0); + const X_FLOAT ytmp = x(i,1); + const X_FLOAT ztmp = x(i,2); + + // two-body interactions + + const int jnum = d_numneigh[i]; + + F_FLOAT fxtmpi = 0.0; + F_FLOAT fytmpi = 0.0; + F_FLOAT fztmpi = 0.0; + + for (int jj = 0; jj < jnum; jj++) { + int j = d_neighbors(i,jj); + j &= NEIGHMASK; + const tagint jtag = tag[j]; + + const int jtype = d_map[type[j]]; + + const X_FLOAT delx = xtmp - x(j,0); + const X_FLOAT dely = ytmp - x(j,1); + const X_FLOAT delz = ztmp - x(j,2); + const F_FLOAT rsq = delx*delx + dely*dely + delz*delz; + + const int ijparam = d_elem2param(itype,jtype,jtype); + + if (rsq >= d_params[ijparam].cutsq) continue; + + twobody(d_params[ijparam],rsq,fpair,eflag,evdwl); + + fxtmpi += delx*fpair; + fytmpi += dely*fpair; + fztmpi += delz*fpair; + + if (EVFLAG) { + if (eflag) ev.evdwl += 0.5*evdwl; + if (vflag_either || eflag_atom) this->template ev_tally(ev,i,j,evdwl,fpair,delx,dely,delz); + } + } + + const int jnumm1 = jnum - 1; + + for (int jj = 0; jj < jnumm1; jj++) { + int j = d_neighbors(i,jj); + j &= NEIGHMASK; + const int jtype = d_map[type[j]]; + const int ijparam = d_elem2param(itype,jtype,jtype); + delr1[0] = x(j,0) - xtmp; + delr1[1] = x(j,1) - ytmp; + delr1[2] = x(j,2) - ztmp; + const F_FLOAT rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; + + if (rsq1 >= d_params[ijparam].cutsq2) continue; + + for (int kk = jj+1; kk < jnum; kk++) { + int k = d_neighbors(i,kk); + k &= NEIGHMASK; + const int ktype = d_map[type[k]]; + const int ikparam = d_elem2param(itype,ktype,ktype); + const int ijkparam = d_elem2param(itype,jtype,ktype); + + delr2[0] = x(k,0) - xtmp; + delr2[1] = x(k,1) - ytmp; + delr2[2] = x(k,2) - ztmp; + const F_FLOAT rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; + + if (rsq2 >= d_params[ikparam].cutsq2) continue; + + threebody(d_params[ijparam],d_params[ikparam],d_params[ijkparam], + rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl); + + fxtmpi -= fj[0] + fk[0]; + fytmpi -= fj[1] + fk[1]; + fztmpi -= fj[2] + fk[2]; + + if (EVFLAG) { + if (eflag) ev.evdwl += evdwl; + if (vflag_either || eflag_atom) this->template ev_tally3(ev,i,j,k,evdwl,0.0,fj,fk,delr1,delr2); + } + } + } + + f(i,0) += fxtmpi; + f(i,1) += fytmpi; + f(i,2) += fztmpi; +} + +template +template +KOKKOS_INLINE_FUNCTION +void PairVashishtaKokkos::operator()(TagPairVashishtaComputeFullA, const int &ii) const { + EV_FLOAT ev; + this->template operator()(TagPairVashishtaComputeFullA(), ii, ev); +} + +/* ---------------------------------------------------------------------- */ + +template +template +KOKKOS_INLINE_FUNCTION +void PairVashishtaKokkos::operator()(TagPairVashishtaComputeFullB, const int &ii, EV_FLOAT& ev) const { + + F_FLOAT delr1[3],delr2[3],fj[3],fk[3]; + F_FLOAT evdwl = 0.0; + F_FLOAT fpair = 0.0; + + const int i = d_ilist[ii]; + + const int itype = d_map[type[i]]; + const X_FLOAT xtmpi = x(i,0); + const X_FLOAT ytmpi = x(i,1); + const X_FLOAT ztmpi = x(i,2); + + const int jnum = d_numneigh[i]; + + F_FLOAT fxtmpi = 0.0; + F_FLOAT fytmpi = 0.0; + F_FLOAT fztmpi = 0.0; + + for (int jj = 0; jj < jnum; jj++) { + int j = d_neighbors(i,jj); + j &= NEIGHMASK; + if (j >= nlocal) continue; + const int jtype = d_map[type[j]]; + const int jiparam = d_elem2param(jtype,itype,itype); + const X_FLOAT xtmpj = x(j,0); + const X_FLOAT ytmpj = x(j,1); + const X_FLOAT ztmpj = x(j,2); + + delr1[0] = xtmpi - xtmpj; + delr1[1] = ytmpi - ytmpj; + delr1[2] = ztmpi - ztmpj; + const F_FLOAT rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; + + if (rsq1 >= d_params[jiparam].cutsq2) continue; + + const int j_jnum = d_numneigh[j]; + + for (int kk = 0; kk < j_jnum; kk++) { + int k = d_neighbors(j,kk); + k &= NEIGHMASK; + if (k == i) continue; + const int ktype = d_map[type[k]]; + const int jkparam = d_elem2param(jtype,ktype,ktype); + const int jikparam = d_elem2param(jtype,itype,ktype); + + delr2[0] = x(k,0) - xtmpj; + delr2[1] = x(k,1) - ytmpj; + delr2[2] = x(k,2) - ztmpj; + const F_FLOAT rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; + + if (rsq2 >= d_params[jkparam].cutsq2) continue; + + if (vflag_atom) + threebody(d_params[jiparam],d_params[jkparam],d_params[jikparam], + rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl); + else + threebodyj(d_params[jiparam],d_params[jkparam],d_params[jikparam], + rsq1,rsq2,delr1,delr2,fj); + + fxtmpi += fj[0]; + fytmpi += fj[1]; + fztmpi += fj[2]; + + if (EVFLAG) + if (vflag_atom || eflag_atom) ev_tally3_atom(ev,i,evdwl,0.0,fj,fk,delr1,delr2); + } + } + + f(i,0) += fxtmpi; + f(i,1) += fytmpi; + f(i,2) += fztmpi; +} + +template +template +KOKKOS_INLINE_FUNCTION +void PairVashishtaKokkos::operator()(TagPairVashishtaComputeFullB, const int &ii) const { + EV_FLOAT ev; + this->template operator()(TagPairVashishtaComputeFullB(), ii, ev); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +template +void PairVashishtaKokkos::coeff(int narg, char **arg) +{ + PairVashishta::coeff(narg,arg); + + // sync map + + int n = atom->ntypes; + + DAT::tdual_int_1d k_map = DAT::tdual_int_1d("pair:map",n+1); + HAT::t_int_1d h_map = k_map.h_view; + + for (int i = 1; i <= n; i++) + h_map[i] = map[i]; + + k_map.template modify(); + k_map.template sync(); + + d_map = k_map.d_view; +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +template +void PairVashishtaKokkos::init_style() +{ + PairVashishta::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::value && + !Kokkos::Impl::is_same::value; + neighbor->requests[irequest]-> + kokkos_device = Kokkos::Impl::is_same::value; + + // always request a full neighbor list + + if (neighflag == FULL || neighflag == HALF || neighflag == HALFTHREAD) { + neighbor->requests[irequest]->full = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full_cluster = 0; + if (neighflag == FULL) + neighbor->requests[irequest]->ghost = 1; + else + neighbor->requests[irequest]->ghost = 0; + } else { + error->all(FLERR,"Cannot use chosen neighbor list style with pair vashishta/kk"); + } +} + +/* ---------------------------------------------------------------------- */ + +template +void PairVashishtaKokkos::setup_params() +{ + PairVashishta::setup_params(); + + // sync elem2param and params + + tdual_int_3d k_elem2param = tdual_int_3d("pair:elem2param",nelements,nelements,nelements); + t_host_int_3d h_elem2param = k_elem2param.h_view; + + tdual_param_1d k_params = tdual_param_1d("pair:params",nparams); + t_host_param_1d h_params = k_params.h_view; + + for (int i = 0; i < nelements; i++) + for (int j = 0; j < nelements; j++) + for (int k = 0; k < nelements; k++) + h_elem2param(i,j,k) = elem2param[i][j][k]; + + for (int m = 0; m < nparams; m++) + h_params[m] = params[m]; + + k_elem2param.template modify(); + k_elem2param.template sync(); + k_params.template modify(); + k_params.template sync(); + + d_elem2param = k_elem2param.d_view; + d_params = k_params.d_view; +} + +/* ---------------------------------------------------------------------- */ + +template +KOKKOS_INLINE_FUNCTION +void PairVashishtaKokkos::twobody(const Param& param, const F_FLOAT& rsq, F_FLOAT& fforce, + const int& eflag, F_FLOAT& eng) const +{ + F_FLOAT r,rinvsq,r4inv,r6inv,reta,lam1r,lam4r,vc2,vc3; + r = sqrt(rsq); + rinvsq = 1.0/rsq; + r4inv = rinvsq*rinvsq; + r6inv = rinvsq*r4inv; + reta = pow(r,-param.eta); + lam1r = r*param.lam1inv; + lam4r = r*param.lam4inv; + vc2 = param.zizj * exp(-lam1r)/r; + vc3 = param.mbigd * r4inv*exp(-lam4r); + + fforce = (param.dvrc*r + - (4.0*vc3 + lam4r*vc3+param.big6w*r6inv + - param.heta*reta - vc2 - lam1r*vc2) + ) * rinvsq; + + if (eflag) eng = param.bigh*reta + vc2 - vc3 - param.bigw*r6inv - r*param.dvrc + param.c0; +} + +/* ---------------------------------------------------------------------- */ + +template +KOKKOS_INLINE_FUNCTION +void PairVashishtaKokkos::threebody(const Param& paramij, const Param& paramik, const Param& paramijk, + const F_FLOAT& rsq1, const F_FLOAT& rsq2, + F_FLOAT *delr1, F_FLOAT *delr2, + F_FLOAT *fj, F_FLOAT *fk, const int& eflag, F_FLOAT& eng) const +{ + F_FLOAT r1,rinvsq1,rainv1,gsrainv1,gsrainvsq1,expgsrainv1; + F_FLOAT r2,rinvsq2,rainv2,gsrainv2,gsrainvsq2,expgsrainv2; + F_FLOAT rinv12,cs,delcs,delcssq,facexp,facrad,frad1,frad2,pcsinv,pcsinvsq,pcs; + F_FLOAT facang,facang12,csfacang,csfac1,csfac2; + + r1 = sqrt(rsq1); + rinvsq1 = 1.0/rsq1; + rainv1 = 1.0/(r1 - paramij.r0); + gsrainv1 = paramij.gamma * rainv1; + gsrainvsq1 = gsrainv1*rainv1/r1; + expgsrainv1 = exp(gsrainv1); + + r2 = sqrt(rsq2); + rinvsq2 = 1.0/rsq2; + rainv2 = 1.0/(r2 - paramik.r0); + gsrainv2 = paramik.gamma * rainv2; + gsrainvsq2 = gsrainv2*rainv2/r2; + expgsrainv2 = exp(gsrainv2); + + rinv12 = 1.0/(r1*r2); + cs = (delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]) * rinv12; + delcs = cs - paramijk.costheta; + delcssq = delcs*delcs; + pcsinv = paramijk.bigc*delcssq + 1.0; + pcsinvsq = pcsinv*pcsinv; + pcs = delcssq/pcsinv; + + facexp = expgsrainv1*expgsrainv2; + + facrad = paramijk.bigb * facexp * pcs; + frad1 = facrad*gsrainvsq1; + frad2 = facrad*gsrainvsq2; + facang = paramijk.big2b * facexp * delcs/pcsinvsq; + facang12 = rinv12*facang; + csfacang = cs*facang; + csfac1 = rinvsq1*csfacang; + + fj[0] = delr1[0]*(frad1+csfac1)-delr2[0]*facang12; + fj[1] = delr1[1]*(frad1+csfac1)-delr2[1]*facang12; + fj[2] = delr1[2]*(frad1+csfac1)-delr2[2]*facang12; + + csfac2 = rinvsq2*csfacang; + + fk[0] = delr2[0]*(frad2+csfac2)-delr1[0]*facang12; + fk[1] = delr2[1]*(frad2+csfac2)-delr1[1]*facang12; + fk[2] = delr2[2]*(frad2+csfac2)-delr1[2]*facang12; + + if (eflag) eng = facrad; +} + +/* ---------------------------------------------------------------------- */ + +template +KOKKOS_INLINE_FUNCTION +void PairVashishtaKokkos::threebodyj(const Param& paramij, const Param& paramik, const Param& paramijk, + const F_FLOAT& rsq1, const F_FLOAT& rsq2, F_FLOAT *delr1, F_FLOAT *delr2, F_FLOAT *fj) const +{ + F_FLOAT r1,rinvsq1,rainv1,gsrainv1,gsrainvsq1,expgsrainv1; + F_FLOAT r2,rinvsq2,rainv2,gsrainv2,gsrainvsq2,expgsrainv2; + F_FLOAT rinv12,cs,delcs,delcssq,facexp,facrad,frad1,frad2,pcsinv,pcsinvsq,pcs; + F_FLOAT facang,facang12,csfacang,csfac1,csfac2; + + r1 = sqrt(rsq1); + rinvsq1 = 1.0/rsq1; + rainv1 = 1.0/(r1 - paramij.r0); + gsrainv1 = paramij.gamma * rainv1; + gsrainvsq1 = gsrainv1*rainv1/r1; + expgsrainv1 = exp(gsrainv1); + + r2 = sqrt(rsq2); + rinvsq2 = 1.0/rsq2; + rainv2 = 1.0/(r2 - paramik.r0); + gsrainv2 = paramik.gamma * rainv2; + gsrainvsq2 = gsrainv2*rainv2/r2; + expgsrainv2 = exp(gsrainv2); + + rinv12 = 1.0/(r1*r2); + cs = (delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]) * rinv12; + delcs = cs - paramijk.costheta; + delcssq = delcs*delcs; + pcsinv = paramijk.bigc*delcssq + 1.0; + pcsinvsq = pcsinv*pcsinv; + pcs = delcssq/pcsinv; + + facexp = expgsrainv1*expgsrainv2; + + facrad = paramijk.bigb * facexp * pcs; + frad1 = facrad*gsrainvsq1; + frad2 = facrad*gsrainvsq2; + facang = paramijk.big2b * facexp * delcs/pcsinvsq; + facang12 = rinv12*facang; + csfacang = cs*facang; + csfac1 = rinvsq1*csfacang; + + fj[0] = delr1[0]*(frad1+csfac1)-delr2[0]*facang12; + fj[1] = delr1[1]*(frad1+csfac1)-delr2[1]*facang12; + fj[2] = delr1[2]*(frad1+csfac1)-delr2[2]*facang12; +} + +/* ---------------------------------------------------------------------- */ + +template +template +KOKKOS_INLINE_FUNCTION +void PairVashishtaKokkos::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 VFLAG = vflag_either; + + // The eatom and vatom arrays are atomic for half/thread neighbor list + + Kokkos::View::value> > v_eatom = k_eatom.view(); + Kokkos::View::value> > v_vatom = k_vatom.view(); + + + if (eflag_atom) { + const E_FLOAT epairhalf = 0.5 * epair; + v_eatom[i] += epairhalf; + if (NEIGHFLAG != FULL) + v_eatom[j] += 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) { + 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 { + 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) { + 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 (NEIGHFLAG != FULL) { + 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; + } + } + } +} + +/* ---------------------------------------------------------------------- + tally eng_vdwl and virial into global and per-atom accumulators + called by SW and hbond potentials, newton_pair is always on + virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drji*fj + drki*fk + ------------------------------------------------------------------------- */ + +template +template +KOKKOS_INLINE_FUNCTION +void PairVashishtaKokkos::ev_tally3(EV_FLOAT &ev, const int &i, const int &j, int &k, + const F_FLOAT &evdwl, const F_FLOAT &ecoul, + F_FLOAT *fj, F_FLOAT *fk, F_FLOAT *drji, F_FLOAT *drki) const +{ + F_FLOAT epairthird,v[6]; + + const int VFLAG = vflag_either; + +// The eatom and vatom arrays are atomic for half/thread neighbor list + + Kokkos::View::value> > v_eatom = k_eatom.view(); + Kokkos::View::value> > v_vatom = k_vatom.view(); + + if (eflag_atom) { + epairthird = THIRD * (evdwl + ecoul); + v_eatom[i] += epairthird; + if (NEIGHFLAG != FULL) { + v_eatom[j] += epairthird; + v_eatom[k] += epairthird; + } + } + + if (VFLAG) { + v[0] = drji[0]*fj[0] + drki[0]*fk[0]; + v[1] = drji[1]*fj[1] + drki[1]*fk[1]; + v[2] = drji[2]*fj[2] + drki[2]*fk[2]; + v[3] = drji[0]*fj[1] + drki[0]*fk[1]; + v[4] = drji[0]*fj[2] + drki[0]*fk[2]; + v[5] = drji[1]*fj[2] + drki[1]*fk[2]; + + if (vflag_global) { + ev.v[0] += v[0]; + ev.v[1] += v[1]; + ev.v[2] += v[2]; + ev.v[3] += v[3]; + ev.v[4] += v[4]; + ev.v[5] += v[5]; + } + + if (vflag_atom) { + v_vatom(i,0) += THIRD*v[0]; v_vatom(i,1) += THIRD*v[1]; + v_vatom(i,2) += THIRD*v[2]; v_vatom(i,3) += THIRD*v[3]; + v_vatom(i,4) += THIRD*v[4]; v_vatom(i,5) += THIRD*v[5]; + + if (NEIGHFLAG != FULL) { + v_vatom(j,0) += THIRD*v[0]; v_vatom(j,1) += THIRD*v[1]; + v_vatom(j,2) += THIRD*v[2]; v_vatom(j,3) += THIRD*v[3]; + v_vatom(j,4) += THIRD*v[4]; v_vatom(j,5) += THIRD*v[5]; + + v_vatom(k,0) += THIRD*v[0]; v_vatom(k,1) += THIRD*v[1]; + v_vatom(k,2) += THIRD*v[2]; v_vatom(k,3) += THIRD*v[3]; + v_vatom(k,4) += THIRD*v[4]; v_vatom(k,5) += THIRD*v[5]; + } + } + } +} + +/* ---------------------------------------------------------------------- + tally eng_vdwl and virial into global and per-atom accumulators + called by SW and hbond potentials, newton_pair is always on + virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drji*fj + drki*fk + ------------------------------------------------------------------------- */ + +template +KOKKOS_INLINE_FUNCTION +void PairVashishtaKokkos::ev_tally3_atom(EV_FLOAT &ev, const int &i, + const F_FLOAT &evdwl, const F_FLOAT &ecoul, + F_FLOAT *fj, F_FLOAT *fk, F_FLOAT *drji, F_FLOAT *drki) const +{ + F_FLOAT epairthird,v[6]; + + const int VFLAG = vflag_either; + + if (eflag_atom) { + epairthird = THIRD * (evdwl + ecoul); + d_eatom[i] += epairthird; + } + + if (VFLAG) { + v[0] = drji[0]*fj[0] + drki[0]*fk[0]; + v[1] = drji[1]*fj[1] + drki[1]*fk[1]; + v[2] = drji[2]*fj[2] + drki[2]*fk[2]; + v[3] = drji[0]*fj[1] + drki[0]*fk[1]; + v[4] = drji[0]*fj[2] + drki[0]*fk[2]; + v[5] = drji[1]*fj[2] + drki[1]*fk[2]; + + if (vflag_atom) { + d_vatom(i,0) += THIRD*v[0]; d_vatom(i,1) += THIRD*v[1]; + d_vatom(i,2) += THIRD*v[2]; d_vatom(i,3) += THIRD*v[3]; + d_vatom(i,4) += THIRD*v[4]; d_vatom(i,5) += THIRD*v[5]; + } + } +} + +namespace LAMMPS_NS { +template class PairVashishtaKokkos; +#ifdef KOKKOS_HAVE_CUDA +template class PairVashishtaKokkos; +#endif +} + diff --git a/src/KOKKOS/pair_vashishta_kokkos.h b/src/KOKKOS/pair_vashishta_kokkos.h new file mode 100644 index 0000000000..87452a952f --- /dev/null +++ b/src/KOKKOS/pair_vashishta_kokkos.h @@ -0,0 +1,154 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(vashishta/kk,PairVashishtaKokkos) +PairStyle(vashishta/kk/device,PairVashishtaKokkos) +PairStyle(vashishta/kk/host,PairVashishtaKokkos) + +#else + +#ifndef LMP_PAIR_VASHISHTA_KOKKOS_H +#define LMP_PAIR_VASHISHTA_KOKKOS_H + +#include "pair_vashishta.h" +#include "pair_kokkos.h" + +template +struct TagPairVashishtaComputeHalf{}; + +template +struct TagPairVashishtaComputeFullA{}; + +template +struct TagPairVashishtaComputeFullB{}; + +namespace LAMMPS_NS { + +template +class PairVashishtaKokkos : public PairVashishta { + public: + enum {EnabledNeighFlags=FULL}; + enum {COUL_FLAG=0}; + typedef DeviceType device_type; + typedef ArrayTypes AT; + typedef EV_FLOAT value_type; + + PairVashishtaKokkos(class LAMMPS *); + virtual ~PairVashishtaKokkos(); + virtual void compute(int, int); + virtual void coeff(int, char **); + virtual void init_style(); + + template + KOKKOS_INLINE_FUNCTION + void operator()(TagPairVashishtaComputeHalf, const int&, EV_FLOAT&) const; + + template + KOKKOS_INLINE_FUNCTION + void operator()(TagPairVashishtaComputeHalf, const int&) const; + + template + KOKKOS_INLINE_FUNCTION + void operator()(TagPairVashishtaComputeFullA, const int&, EV_FLOAT&) const; + + template + KOKKOS_INLINE_FUNCTION + void operator()(TagPairVashishtaComputeFullA, const int&) const; + + template + KOKKOS_INLINE_FUNCTION + void operator()(TagPairVashishtaComputeFullB, const int&, EV_FLOAT&) const; + + template + KOKKOS_INLINE_FUNCTION + void operator()(TagPairVashishtaComputeFullB, const int&) const; + + template + KOKKOS_INLINE_FUNCTION + void ev_tally(EV_FLOAT &ev, const int &i, const int &j, + const F_FLOAT &epair, const F_FLOAT &fpair, const F_FLOAT &delx, + const F_FLOAT &dely, const F_FLOAT &delz) const; + + template + KOKKOS_INLINE_FUNCTION + void ev_tally3(EV_FLOAT &ev, const int &i, const int &j, int &k, + const F_FLOAT &evdwl, const F_FLOAT &ecoul, + F_FLOAT *fj, F_FLOAT *fk, F_FLOAT *drji, F_FLOAT *drki) const; + + KOKKOS_INLINE_FUNCTION + void ev_tally3_atom(EV_FLOAT &ev, const int &i, + const F_FLOAT &evdwl, const F_FLOAT &ecoul, + F_FLOAT *fj, F_FLOAT *fk, F_FLOAT *drji, F_FLOAT *drki) const; + + protected: + typedef Kokkos::DualView tdual_int_3d; + typedef typename tdual_int_3d::t_dev_const_randomread t_int_3d_randomread; + typedef typename tdual_int_3d::t_host t_host_int_3d; + + t_int_3d_randomread d_elem2param; + DAT::t_int_1d_randomread d_map; + + typedef Kokkos::DualView tdual_param_1d; + typedef typename tdual_param_1d::t_dev t_param_1d; + typedef typename tdual_param_1d::t_host t_host_param_1d; + + t_param_1d d_params; + + virtual void setup_params(); + void twobody(const Param&, const F_FLOAT&, F_FLOAT&, const int&, F_FLOAT&) const; + void threebody(const Param&, const Param&, const Param&, const F_FLOAT&, const F_FLOAT&, F_FLOAT *, F_FLOAT *, + F_FLOAT *, F_FLOAT *, const int&, F_FLOAT&) const; + void threebodyj(const Param&, const Param&, const Param&, const F_FLOAT&, const F_FLOAT&, F_FLOAT *, F_FLOAT *, + F_FLOAT *) const; + + typename ArrayTypes::t_x_array_randomread x; + typename ArrayTypes::t_f_array f; + typename ArrayTypes::t_tagint_1d tag; + typename ArrayTypes::t_int_1d_randomread type; + + DAT::tdual_efloat_1d k_eatom; + DAT::tdual_virial_array k_vatom; + DAT::t_efloat_1d d_eatom; + DAT::t_virial_array d_vatom; + + DAT::t_int_1d_randomread d_type2frho; + DAT::t_int_2d_randomread d_type2rhor; + DAT::t_int_2d_randomread d_type2z2r; + + typename ArrayTypes::t_neighbors_2d d_neighbors; + typename ArrayTypes::t_int_1d_randomread d_ilist; + typename ArrayTypes::t_int_1d_randomread d_numneigh; + //NeighListKokkos k_list; + + int neighflag,newton_pair; + int nlocal,nall,eflag,vflag; + + int inum; + + friend void pair_virial_fdotr_compute(PairVashishtaKokkos*); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Cannot use chosen neighbor list style with pair vashishta/kk + +Self-explanatory. + +*/ From 039bda9b61db398171f990c1cd0bd6c4d951dfd3 Mon Sep 17 00:00:00 2001 From: Anders Hafreager Date: Tue, 25 Oct 2016 15:02:00 +0200 Subject: [PATCH 2/8] Added updated vashishta for KOKKOS support (cherry picked from commit 96089a42547f625e70aa2ac3933d248d2731b731) --- src/MANYBODY/pair_vashishta.cpp | 2 ++ src/MANYBODY/pair_vashishta.h | 3 +-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/MANYBODY/pair_vashishta.cpp b/src/MANYBODY/pair_vashishta.cpp index 19f6017907..86391dc39d 100644 --- a/src/MANYBODY/pair_vashishta.cpp +++ b/src/MANYBODY/pair_vashishta.cpp @@ -60,6 +60,8 @@ PairVashishta::PairVashishta(LAMMPS *lmp) : Pair(lmp) PairVashishta::~PairVashishta() { + if (copymode) return; + if (elements) for (int i = 0; i < nelements; i++) delete [] elements[i]; delete [] elements; diff --git a/src/MANYBODY/pair_vashishta.h b/src/MANYBODY/pair_vashishta.h index 87077011e6..5e91fb245d 100644 --- a/src/MANYBODY/pair_vashishta.h +++ b/src/MANYBODY/pair_vashishta.h @@ -34,7 +34,6 @@ class PairVashishta : public Pair { double init_one(int, int); void init_style(); - protected: struct Param { double bigb,gamma,r0,bigc,costheta; double bigh,eta,zi,zj; @@ -45,7 +44,7 @@ class PairVashishta : public Pair { double lam1rc,lam4rc,vrcc2,vrcc3,vrc,dvrc,c0; int ielement,jelement,kelement; }; - + protected: double cutmax; // max cutoff for all elements int nelements; // # of unique elements char **elements; // names of unique elements From 21619f6a2f808847ce4edc4d54abb38d4e593eec Mon Sep 17 00:00:00 2001 From: stamoor Date: Wed, 19 Oct 2016 19:02:08 +0000 Subject: [PATCH 3/8] Recommitting reverted change git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15794 f3b2605a-c512-4ea7-a41b-209d697bcdaa (cherry picked from commit c0b98f5299e0b13802ff7ec32f162f477f407682) --- src/region.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/region.h b/src/region.h index c813772a82..0d406a3e65 100644 --- a/src/region.h +++ b/src/region.h @@ -106,11 +106,12 @@ class Region : protected Pointers { void options(int, char **); void point_on_line_segment(double *, double *, double *, double *); void forward_transform(double &, double &, double &); + double point[3],runit[3]; private: char *xstr,*ystr,*zstr,*tstr; int xvar,yvar,zvar,tvar; - double axis[3],point[3],runit[3]; + double axis[3]; void inverse_transform(double &, double &, double &); void rotate(double &, double &, double &, double); From e37d2b5c9431098d9f932bf6df7db7a074ea2944 Mon Sep 17 00:00:00 2001 From: Anders Hafreager Date: Sat, 29 Oct 2016 22:20:37 +0200 Subject: [PATCH 4/8] Calculating short neighbor lists --- src/KOKKOS/pair_vashishta_kokkos.cpp | 59 +++++++++++++++++++++++++++- src/KOKKOS/pair_vashishta_kokkos.h | 10 ++++- src/MANYBODY/pair_vashishta.cpp | 3 +- src/MANYBODY/pair_vashishta.h | 1 + 4 files changed, 70 insertions(+), 3 deletions(-) diff --git a/src/KOKKOS/pair_vashishta_kokkos.cpp b/src/KOKKOS/pair_vashishta_kokkos.cpp index b5fe84faf6..c710320c67 100644 --- a/src/KOKKOS/pair_vashishta_kokkos.cpp +++ b/src/KOKKOS/pair_vashishta_kokkos.cpp @@ -108,7 +108,7 @@ void PairVashishtaKokkos::compute(int eflag_in, int vflag_in) newton_pair = force->newton_pair; nall = atom->nlocal + atom->nghost; - const int inum = list->inum; + inum = list->inum; const int ignum = inum + list->gnum; NeighListKokkos* k_list = static_cast*>(list); d_ilist = k_list->d_ilist; @@ -121,6 +121,29 @@ void PairVashishtaKokkos::compute(int eflag_in, int vflag_in) EV_FLOAT ev; EV_FLOAT ev_all; + int max_neighs = d_neighbors.dimension_1(); + + if ((d_neighbors_short_2body.dimension_1() != max_neighs) || + (d_neighbors_short_2body.dimension_0() != ignum)) { + d_neighbors_short_2body = Kokkos::View("Vashishta::neighbors_short_2body",ignum,max_neighs); + } + if (d_numneigh_short_2body.dimension_0()!=ignum) { + d_numneigh_short_2body = Kokkos::View("Vashishta::numneighs_short_2body",ignum); + } + + if ((d_neighbors_short_3body.dimension_1() != max_neighs) || + (d_neighbors_short_3body.dimension_0() != ignum)) { + d_neighbors_short_3body = Kokkos::View("Vashishta::neighbors_short_3body",ignum,max_neighs); + } + + if (d_numneigh_short_3body.dimension_0()!=ignum) { + d_numneigh_short_3body = Kokkos::View("Vashishta::numneighs_short_3body",ignum); + } + + Kokkos::parallel_for(Kokkos::RangePolicy(0,neighflag==FULL?ignum:inum), *this); + + + // loop over neighbor list of my atoms if (neighflag == HALF) { @@ -174,6 +197,40 @@ void PairVashishtaKokkos::compute(int eflag_in, int vflag_in) copymode = 0; } +template +KOKKOS_INLINE_FUNCTION +void PairVashishtaKokkos::operator()(TagPairVashishtaComputeShortNeigh, const int& ii) const { + const int i = d_ilist[ii]; + const X_FLOAT xtmp = x(i,0); + const X_FLOAT ytmp = x(i,1); + const X_FLOAT ztmp = x(i,2); + + const int jnum = d_numneigh[i]; + int inside_2body = 0; + int inside_3body = 0; + for (int jj = 0; jj < jnum; jj++) { + int j = d_neighbors(i,jj); + j &= NEIGHMASK; + + const X_FLOAT delx = xtmp - x(j,0); + const X_FLOAT dely = ytmp - x(j,1); + const X_FLOAT delz = ztmp - x(j,2); + const F_FLOAT rsq = delx*delx + dely*dely + delz*delz; + + if (rsq < cutmax*cutmax) { + d_neighbors_short_2body(i,inside_2body) = j; + inside_2body++; + } + + if (rsq < cutmax_3body*cutmax_3body) { + d_neighbors_short_3body(i,inside_3body) = j; + inside_3body++; + } + } + d_numneigh_short_2body(i) = inside_2body; + d_numneigh_short_3body(i) = inside_3body; +} + template template KOKKOS_INLINE_FUNCTION diff --git a/src/KOKKOS/pair_vashishta_kokkos.h b/src/KOKKOS/pair_vashishta_kokkos.h index 87452a952f..8bc4d621c2 100644 --- a/src/KOKKOS/pair_vashishta_kokkos.h +++ b/src/KOKKOS/pair_vashishta_kokkos.h @@ -34,6 +34,8 @@ struct TagPairVashishtaComputeFullA{}; template struct TagPairVashishtaComputeFullB{}; +struct TagPairVashishtaComputeShortNeigh{}; + namespace LAMMPS_NS { template @@ -75,6 +77,9 @@ class PairVashishtaKokkos : public PairVashishta { KOKKOS_INLINE_FUNCTION void operator()(TagPairVashishtaComputeFullB, const int&) const; + KOKKOS_INLINE_FUNCTION + void operator()(TagPairVashishtaComputeShortNeigh, const int&) const; + template KOKKOS_INLINE_FUNCTION void ev_tally(EV_FLOAT &ev, const int &i, const int &j, @@ -136,7 +141,10 @@ class PairVashishtaKokkos : public PairVashishta { int nlocal,nall,eflag,vflag; int inum; - + Kokkos::View d_neighbors_short_2body; + Kokkos::View d_numneigh_short_2body; + Kokkos::View d_neighbors_short_3body; + Kokkos::View d_numneigh_short_3body; friend void pair_virial_fdotr_compute(PairVashishtaKokkos*); }; diff --git a/src/MANYBODY/pair_vashishta.cpp b/src/MANYBODY/pair_vashishta.cpp index 86391dc39d..cc28632f60 100644 --- a/src/MANYBODY/pair_vashishta.cpp +++ b/src/MANYBODY/pair_vashishta.cpp @@ -532,11 +532,12 @@ void PairVashishta::setup_params() } // set cutmax to max of all params - + cutmax_3body = 0.0; cutmax = 0.0; for (m = 0; m < nparams; m++) { if (params[m].cut > cutmax) cutmax = params[m].cut; if (params[m].r0 > cutmax) cutmax = params[m].r0; + if (params[m].r0 > cutmax_3body) cutmax_3body = params[m].r0; } } diff --git a/src/MANYBODY/pair_vashishta.h b/src/MANYBODY/pair_vashishta.h index 5e91fb245d..3fa5a3e185 100644 --- a/src/MANYBODY/pair_vashishta.h +++ b/src/MANYBODY/pair_vashishta.h @@ -46,6 +46,7 @@ class PairVashishta : public Pair { }; protected: double cutmax; // max cutoff for all elements + double cutmax_3body; // max cutoff for all elements int nelements; // # of unique elements char **elements; // names of unique elements int ***elem2param; // mapping from element triplets to parameters From 51e2313fac2cc424e230efe9f7245e35d67bc553 Mon Sep 17 00:00:00 2001 From: Anders Hafreager Date: Sat, 29 Oct 2016 22:35:29 +0200 Subject: [PATCH 5/8] Using short neighbor lists --- src/KOKKOS/pair_vashishta_kokkos.cpp | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/KOKKOS/pair_vashishta_kokkos.cpp b/src/KOKKOS/pair_vashishta_kokkos.cpp index c710320c67..ea1bd12da1 100644 --- a/src/KOKKOS/pair_vashishta_kokkos.cpp +++ b/src/KOKKOS/pair_vashishta_kokkos.cpp @@ -253,14 +253,16 @@ void PairVashishtaKokkos::operator()(TagPairVashishtaComputeHalf::operator()(TagPairVashishtaComputeHalf::operator()(TagPairVashishtaComputeHalf Date: Sat, 29 Oct 2016 22:54:43 +0200 Subject: [PATCH 6/8] Using short neighborlists in neigh full --- src/KOKKOS/pair_vashishta_kokkos.cpp | 24 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/src/KOKKOS/pair_vashishta_kokkos.cpp b/src/KOKKOS/pair_vashishta_kokkos.cpp index ea1bd12da1..3e4142bba4 100644 --- a/src/KOKKOS/pair_vashishta_kokkos.cpp +++ b/src/KOKKOS/pair_vashishta_kokkos.cpp @@ -253,7 +253,6 @@ void PairVashishtaKokkos::operator()(TagPairVashishtaComputeHalf::operator()(TagPairVashishtaComputeHalf::operator()(TagPairVashishtaComputeFullA::operator()(TagPairVashishtaComputeFullA::operator()(TagPairVashishtaComputeFullA= d_params[ijparam].cutsq2) continue; - for (int kk = jj+1; kk < jnum; kk++) { - int k = d_neighbors(i,kk); + for (int kk = jj+1; kk < jnumm1; kk++) { + int k = d_neighbors_short_3body(i,kk); k &= NEIGHMASK; const int ktype = d_map[type[k]]; const int ikparam = d_elem2param(itype,ktype,ktype); @@ -497,14 +495,14 @@ void PairVashishtaKokkos::operator()(TagPairVashishtaComputeFullB= nlocal) continue; const int jtype = d_map[type[j]]; @@ -520,10 +518,10 @@ void PairVashishtaKokkos::operator()(TagPairVashishtaComputeFullB= d_params[jiparam].cutsq2) continue; - const int j_jnum = d_numneigh[j]; + const int j_jnum = d_numneigh_short_3body[j]; for (int kk = 0; kk < j_jnum; kk++) { - int k = d_neighbors(j,kk); + int k = d_neighbors_short_3body(j,kk); k &= NEIGHMASK; if (k == i) continue; const int ktype = d_map[type[k]]; From d20b32092ec129525c82488b6907ff4b4957f415 Mon Sep 17 00:00:00 2001 From: Anders Hafreager Date: Tue, 8 Nov 2016 18:57:27 +0100 Subject: [PATCH 7/8] Building correct shortlists and removed rsq test in force loops --- src/KOKKOS/pair_vashishta_kokkos.cpp | 27 ++++++++------------------- 1 file changed, 8 insertions(+), 19 deletions(-) diff --git a/src/KOKKOS/pair_vashishta_kokkos.cpp b/src/KOKKOS/pair_vashishta_kokkos.cpp index 3e4142bba4..73e4e04f98 100644 --- a/src/KOKKOS/pair_vashishta_kokkos.cpp +++ b/src/KOKKOS/pair_vashishta_kokkos.cpp @@ -201,6 +201,7 @@ template KOKKOS_INLINE_FUNCTION void PairVashishtaKokkos::operator()(TagPairVashishtaComputeShortNeigh, const int& ii) const { const int i = d_ilist[ii]; + const int itype = d_map[type[i]]; const X_FLOAT xtmp = x(i,0); const X_FLOAT ytmp = x(i,1); const X_FLOAT ztmp = x(i,2); @@ -211,18 +212,20 @@ void PairVashishtaKokkos::operator()(TagPairVashishtaComputeShortNei for (int jj = 0; jj < jnum; jj++) { int j = d_neighbors(i,jj); j &= NEIGHMASK; - + const int jtype = d_map[type[j]]; + const int ijparam = d_elem2param(itype,jtype,jtype); + const X_FLOAT delx = xtmp - x(j,0); const X_FLOAT dely = ytmp - x(j,1); const X_FLOAT delz = ztmp - x(j,2); const F_FLOAT rsq = delx*delx + dely*dely + delz*delz; - if (rsq < cutmax*cutmax) { + if (rsq < d_params[ijparam].cutsq) { d_neighbors_short_2body(i,inside_2body) = j; inside_2body++; } - if (rsq < cutmax_3body*cutmax_3body) { + if (rsq < d_params[ijparam].cutsq2) { d_neighbors_short_3body(i,inside_3body) = j; inside_3body++; } @@ -282,8 +285,7 @@ void PairVashishtaKokkos::operator()(TagPairVashishtaComputeHalf= d_params[ijparam].cutsq) continue; - + twobody(d_params[ijparam],rsq,fpair,eflag,evdwl); fxtmpi += delx*fpair; @@ -310,8 +312,7 @@ void PairVashishtaKokkos::operator()(TagPairVashishtaComputeHalf= d_params[ijparam].cutsq2) continue; - + F_FLOAT fxtmpj = 0.0; F_FLOAT fytmpj = 0.0; F_FLOAT fztmpj = 0.0; @@ -328,8 +329,6 @@ void PairVashishtaKokkos::operator()(TagPairVashishtaComputeHalf= d_params[ikparam].cutsq2) continue; - threebody(d_params[ijparam],d_params[ikparam],d_params[ijkparam], rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl); @@ -408,8 +407,6 @@ void PairVashishtaKokkos::operator()(TagPairVashishtaComputeFullA= d_params[ijparam].cutsq) continue; - twobody(d_params[ijparam],rsq,fpair,eflag,evdwl); fxtmpi += delx*fpair; @@ -434,8 +431,6 @@ void PairVashishtaKokkos::operator()(TagPairVashishtaComputeFullA= d_params[ijparam].cutsq2) continue; - for (int kk = jj+1; kk < jnumm1; kk++) { int k = d_neighbors_short_3body(i,kk); k &= NEIGHMASK; @@ -448,8 +443,6 @@ void PairVashishtaKokkos::operator()(TagPairVashishtaComputeFullA= d_params[ikparam].cutsq2) continue; - threebody(d_params[ijparam],d_params[ikparam],d_params[ijkparam], rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl); @@ -516,8 +509,6 @@ void PairVashishtaKokkos::operator()(TagPairVashishtaComputeFullB= d_params[jiparam].cutsq2) continue; - const int j_jnum = d_numneigh_short_3body[j]; for (int kk = 0; kk < j_jnum; kk++) { @@ -533,8 +524,6 @@ void PairVashishtaKokkos::operator()(TagPairVashishtaComputeFullB= d_params[jkparam].cutsq2) continue; - if (vflag_atom) threebody(d_params[jiparam],d_params[jkparam],d_params[jikparam], rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl); From 3e36ec37542c653b25d99d136b672f65a8a32af9 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 9 Nov 2016 16:17:46 -0500 Subject: [PATCH 8/8] remove unused class member --- src/MANYBODY/pair_vashishta.h | 1 - 1 file changed, 1 deletion(-) diff --git a/src/MANYBODY/pair_vashishta.h b/src/MANYBODY/pair_vashishta.h index e365ba2ecc..5ec47cf92d 100644 --- a/src/MANYBODY/pair_vashishta.h +++ b/src/MANYBODY/pair_vashishta.h @@ -46,7 +46,6 @@ class PairVashishta : public Pair { }; protected: double cutmax; // max cutoff for all elements - double cutmax_3body; // max cutoff for all elements int nelements; // # of unique elements char **elements; // names of unique elements int ***elem2param; // mapping from element triplets to parameters