/* ---------------------------------------------------------------------- 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: Christian Trott (SNL) ------------------------------------------------------------------------- */ #include #include #include #include #include "pair_table_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_kokkos.h" #include "error.h" #include "atom_masks.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ template PairTableKokkos::PairTableKokkos(LAMMPS *lmp) : PairTable(lmp) { update_table = 0; atomKK = (AtomKokkos *) atom; execution_space = ExecutionSpaceFromDevice::space; datamask_read = X_MASK | F_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK; datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK; h_table = new TableHost(); d_table = new TableDevice(); } /* ---------------------------------------------------------------------- */ template PairTableKokkos::~PairTableKokkos() { if (copymode) return; delete h_table; h_table = nullptr; delete d_table; d_table = nullptr; copymode = true; //prevents base class destructor from running } /* ---------------------------------------------------------------------- */ template void PairTableKokkos::compute(int eflag_in, int vflag_in) { if(update_table) create_kokkos_tables(); if(tabstyle == LOOKUP) compute_style(eflag_in,vflag_in); if(tabstyle == LINEAR) compute_style(eflag_in,vflag_in); if(tabstyle == SPLINE) compute_style(eflag_in,vflag_in); if(tabstyle == BITMAP) compute_style(eflag_in,vflag_in); } template template void PairTableKokkos::compute_style(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,0); else evflag = vflag_fdotr = 0; // reallocate per-atom arrays if necessary if (eflag_atom) { memoryKK->destroy_kokkos(k_eatom,eatom); memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom"); d_eatom = k_eatom.view(); } if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom,vatom); memoryKK->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom"); d_vatom = k_vatom.view(); } atomKK->sync(execution_space,datamask_read); //k_cutsq.template sync(); //k_params.template sync(); if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); else atomKK->modified(execution_space,F_MASK); x = c_x = atomKK->k_x.view(); f = atomKK->k_f.view(); type = atomKK->k_type.view(); nlocal = atom->nlocal; nall = atom->nlocal + atom->nghost; special_lj[0] = force->special_lj[0]; special_lj[1] = force->special_lj[1]; special_lj[2] = force->special_lj[2]; special_lj[3] = force->special_lj[3]; newton_pair = force->newton_pair; d_cutsq = d_table->cutsq; // loop over neighbors of my atoms EV_FLOAT ev; if(atom->ntypes > MAX_TYPES_STACKPARAMS) { if (neighflag == FULL) { PairComputeFunctor,FULL,false,S_TableCompute > ff(this,(NeighListKokkos*) list); if (eflag || vflag) Kokkos::parallel_reduce(list->inum,ff,ev); else Kokkos::parallel_for(list->inum,ff); ff.contribute(); } else if (neighflag == HALFTHREAD) { PairComputeFunctor,HALFTHREAD,false,S_TableCompute > ff(this,(NeighListKokkos*) list); if (eflag || vflag) Kokkos::parallel_reduce(list->inum,ff,ev); else Kokkos::parallel_for(list->inum,ff); ff.contribute(); } else if (neighflag == HALF) { PairComputeFunctor,HALF,false,S_TableCompute > f(this,(NeighListKokkos*) list); if (eflag || vflag) Kokkos::parallel_reduce(list->inum,f,ev); else Kokkos::parallel_for(list->inum,f); f.contribute(); } else if (neighflag == N2) { PairComputeFunctor,N2,false,S_TableCompute > f(this,(NeighListKokkos*) list); if (eflag || vflag) Kokkos::parallel_reduce(list->inum,f,ev); else Kokkos::parallel_for(list->inum,f); f.contribute(); } } else { if (neighflag == FULL) { PairComputeFunctor,FULL,true,S_TableCompute > f(this,(NeighListKokkos*) list); if (eflag || vflag) Kokkos::parallel_reduce(list->inum,f,ev); else Kokkos::parallel_for(list->inum,f); f.contribute(); } else if (neighflag == HALFTHREAD) { PairComputeFunctor,HALFTHREAD,true,S_TableCompute > f(this,(NeighListKokkos*) list); if (eflag || vflag) Kokkos::parallel_reduce(list->inum,f,ev); else Kokkos::parallel_for(list->inum,f); f.contribute(); } else if (neighflag == HALF) { PairComputeFunctor,HALF,true,S_TableCompute > f(this,(NeighListKokkos*) list); if (eflag || vflag) Kokkos::parallel_reduce(list->inum,f,ev); else Kokkos::parallel_for(list->inum,f); f.contribute(); } else if (neighflag == N2) { PairComputeFunctor,N2,true,S_TableCompute > f(this,(NeighListKokkos*) list); if (eflag || vflag) Kokkos::parallel_reduce(list->inum,f,ev); else Kokkos::parallel_for(list->inum,f); f.contribute(); } } 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 (eflag_atom) { k_eatom.template modify(); k_eatom.template sync(); } if (vflag_atom) { k_vatom.template modify(); k_vatom.template sync(); } if (vflag_fdotr) pair_virial_fdotr_compute(this); } template template KOKKOS_INLINE_FUNCTION F_FLOAT PairTableKokkos:: compute_fpair(const F_FLOAT& rsq, const int& i, const int&j, const int& itype, const int& jtype) const { (void) i; (void) j; union_int_float_t rsq_lookup; double fpair; const int tidx = d_table_const.tabindex(itype,jtype); if (Specialisation::TabStyle == LOOKUP) { const int itable = static_cast ((rsq - d_table_const.innersq(tidx)) * d_table_const.invdelta(tidx)); fpair = d_table_const.f(tidx,itable); } else if (Specialisation::TabStyle == LINEAR) { const int itable = static_cast ((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 (Specialisation::TabStyle == SPLINE) { const int itable = static_cast ((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 template KOKKOS_INLINE_FUNCTION F_FLOAT PairTableKokkos:: compute_evdwl(const F_FLOAT& rsq, const int& i, const int&j, const int& itype, const int& jtype) const { (void) i; (void) j; double evdwl; union_int_float_t rsq_lookup; const int tidx = d_table_const.tabindex(itype,jtype); if (Specialisation::TabStyle == LOOKUP) { const int itable = static_cast ((rsq - d_table_const.innersq(tidx)) * d_table_const.invdelta(tidx)); evdwl = d_table_const.e(tidx,itable); } else if (Specialisation::TabStyle == LINEAR) { const int itable = static_cast ((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 (Specialisation::TabStyle == SPLINE) { const int itable = static_cast ((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 void PairTableKokkos::create_kokkos_tables() { const int tlm1 = tablength-1; memoryKK->create_kokkos(d_table->nshiftbits,h_table->nshiftbits,ntables,"Table::nshiftbits"); memoryKK->create_kokkos(d_table->nmask,h_table->nmask,ntables,"Table::nmask"); memoryKK->create_kokkos(d_table->innersq,h_table->innersq,ntables,"Table::innersq"); memoryKK->create_kokkos(d_table->invdelta,h_table->invdelta,ntables,"Table::invdelta"); memoryKK->create_kokkos(d_table->deltasq6,h_table->deltasq6,ntables,"Table::deltasq6"); if(tabstyle == LOOKUP) { memoryKK->create_kokkos(d_table->e,h_table->e,ntables,tlm1,"Table::e"); memoryKK->create_kokkos(d_table->f,h_table->f,ntables,tlm1,"Table::f"); } if(tabstyle == LINEAR) { memoryKK->create_kokkos(d_table->rsq,h_table->rsq,ntables,tablength,"Table::rsq"); memoryKK->create_kokkos(d_table->e,h_table->e,ntables,tablength,"Table::e"); memoryKK->create_kokkos(d_table->f,h_table->f,ntables,tablength,"Table::f"); memoryKK->create_kokkos(d_table->de,h_table->de,ntables,tlm1,"Table::de"); memoryKK->create_kokkos(d_table->df,h_table->df,ntables,tlm1,"Table::df"); } if(tabstyle == SPLINE) { memoryKK->create_kokkos(d_table->rsq,h_table->rsq,ntables,tablength,"Table::rsq"); memoryKK->create_kokkos(d_table->e,h_table->e,ntables,tablength,"Table::e"); memoryKK->create_kokkos(d_table->f,h_table->f,ntables,tablength,"Table::f"); memoryKK->create_kokkos(d_table->e2,h_table->e2,ntables,tablength,"Table::e2"); memoryKK->create_kokkos(d_table->f2,h_table->f2,ntables,tablength,"Table::f2"); } if(tabstyle == BITMAP) { int ntable = 1 << tablength; memoryKK->create_kokkos(d_table->rsq,h_table->rsq,ntables,ntable,"Table::rsq"); memoryKK->create_kokkos(d_table->e,h_table->e,ntables,ntable,"Table::e"); memoryKK->create_kokkos(d_table->f,h_table->f,ntables,ntable,"Table::f"); memoryKK->create_kokkos(d_table->de,h_table->de,ntables,ntable,"Table::de"); memoryKK->create_kokkos(d_table->df,h_table->df,ntables,ntable,"Table::df"); memoryKK->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; jrsq.extent(1); j++) h_table->rsq(i,j) = tb->rsq[j]; for(int j = 0; jdrsq.extent(1); j++) h_table->drsq(i,j) = tb->drsq[j]; for(int j = 0; je.extent(1); j++) h_table->e(i,j) = tb->e[j]; for(int j = 0; jde.extent(1); j++) h_table->de(i,j) = tb->de[j]; for(int j = 0; jf.extent(1); j++) h_table->f(i,j) = tb->f[j]; for(int j = 0; jdf.extent(1); j++) h_table->df(i,j) = tb->df[j]; for(int j = 0; je2.extent(1); j++) h_table->e2(i,j) = tb->e2[j]; for(int j = 0; jf2.extent(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 void PairTableKokkos::allocate() { allocated = 1; const int nt = atom->ntypes + 1; memory->create(setflag,nt,nt,"pair:setflag"); memoryKK->create_kokkos(d_table->cutsq,h_table->cutsq,cutsq,nt,nt,"pair:cutsq"); memoryKK->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 void PairTableKokkos::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 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); if (allocated) { memory->destroy(setflag); d_table_const.tabindex = d_table->tabindex = typename ArrayTypes::t_int_2d(); h_table->tabindex = typename ArrayTypes::t_int_2d(); d_table_const.cutsq = d_table->cutsq = typename ArrayTypes::t_ffloat_2d(); h_table->cutsq = typename ArrayTypes::t_ffloat_2d(); } allocated = 0; ntables = 0; tables = NULL; } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ template double PairTableKokkos::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 void PairTableKokkos::compute_table(Table *tb) { update_table = 1; PairTable::compute_table(tb); } template void PairTableKokkos::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::value && !Kokkos::Impl::is_same::value; neighbor->requests[irequest]-> kokkos_device = Kokkos::Impl::is_same::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 void PairTableKokkos::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 PairTableKokkos; #ifdef KOKKOS_HAVE_CUDA template class PairTableKokkos; #endif }