diff --git a/lib/gpu/lal_lj_cubic.cpp b/lib/gpu/lal_lj_cubic.cpp new file mode 100644 index 0000000000..25f83166e1 --- /dev/null +++ b/lib/gpu/lal_lj_cubic.cpp @@ -0,0 +1,159 @@ +/*************************************************************************** + lj_cubic.cpp + ------------------- + Trung Dac Nguyen + + Class for acceleration of the lj/cubic pair style. + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : + email : ndactrung@gmail.com + ***************************************************************************/ + +#if defined(USE_OPENCL) +#include "lj_cubic_cl.h" +#elif defined(USE_CUDART) +const char *lj_cubic=0; +#else +#include "lj_cubic_cubin.h" +#endif + +#include "lal_lj_cubic.h" +#include +using namespace LAMMPS_AL; +#define LJCubicT LJCubic + +extern Device device; + +template +LJCubicT::LJCubic() : BaseAtomic(), _allocated(false) { +} + +template +LJCubicT::~LJCubic() { + clear(); +} + +template +int LJCubicT::bytes_per_atom(const int max_nbors) const { + return this->bytes_per_atom_atomic(max_nbors); +} + +template +int LJCubicT::init(const int ntypes, + double **host_cutsq, double **host_cut_inner_sq, + double **host_cut_inner, double **host_sigma, + double **host_epsilon, double **host_lj1, + double **host_lj2, double **host_lj3, double **host_lj4, + double *host_special_lj, const int nlocal, + const int nall, const int max_nbors, + const int maxspecial, const double cell_size, + const double gpu_split, FILE *_screen) { + int success; + success=this->init_atomic(nlocal,nall,max_nbors,maxspecial,cell_size,gpu_split, + _screen,lj_cubic,"k_lj_cubic"); + if (success!=0) + return success; + + // If atom type constants fit in shared memory use fast kernel + int lj_types=ntypes; + shared_types=false; + int max_shared_types=this->device->max_shared_types(); + if (lj_types<=max_shared_types && this->_block_size>=max_shared_types) { + lj_types=max_shared_types; + shared_types=true; + } + _lj_types=lj_types; + + // Allocate a host write buffer for data initialization + UCL_H_Vec host_write(lj_types*lj_types*32,*(this->ucl_device), + UCL_WRITE_ONLY); + + for (int i=0; iucl_device),UCL_READ_ONLY); + this->atom->type_pack4(ntypes,lj_types,lj1,host_write,host_lj1,host_lj2, + host_cutsq); + + lj2.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY); + this->atom->type_pack4(ntypes,lj_types,lj2,host_write,host_cut_inner_sq, + host_cut_inner,host_sigma,host_epsilon); + + lj3.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY); + this->atom->type_pack2(ntypes,lj_types,lj3,host_write,host_lj3,host_lj4); + + UCL_H_Vec dview; + sp_lj.alloc(4,*(this->ucl_device),UCL_READ_ONLY); + dview.view(host_special_lj,4,*(this->ucl_device)); + ucl_copy(sp_lj,dview,false); + + _allocated=true; + this->_max_bytes=lj1.row_bytes()+lj2.row_bytes()+lj3.row_bytes() + +sp_lj.row_bytes(); + return 0; +} + +template +void LJCubicT::clear() { + if (!_allocated) + return; + _allocated=false; + + lj1.clear(); + lj2.clear(); + lj3.clear(); + sp_lj.clear(); + this->clear_atomic(); +} + +template +double LJCubicT::host_memory_usage() const { + return this->host_memory_usage_atomic()+sizeof(LJCubic); +} + +// --------------------------------------------------------------------------- +// Calculate energies, forces, and torques +// --------------------------------------------------------------------------- +template +void LJCubicT::loop(const bool _eflag, const bool _vflag) { + // Compute the block size and grid size to keep all cores busy + const int BX=this->block_size(); + int eflag, vflag; + if (_eflag) + eflag=1; + else + eflag=0; + + if (_vflag) + vflag=1; + else + vflag=0; + + int GX=static_cast(ceil(static_cast(this->ans->inum())/ + (BX/this->_threads_per_atom))); + + int ainum=this->ans->inum(); + int nbor_pitch=this->nbor->nbor_pitch(); + this->time_pair.start(); + if (shared_types) { + this->k_pair_fast.set_size(GX,BX); + this->k_pair_fast.run(&this->atom->x, &lj1, &lj2, &lj3, &sp_lj, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &this->ans->force, &this->ans->engv, &eflag, + &vflag, &ainum, &nbor_pitch, + &this->_threads_per_atom); + } else { + this->k_pair.set_size(GX,BX); + this->k_pair.run(&this->atom->x, &lj1, &lj2, &lj3, &_lj_types, &sp_lj, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &this->ans->force, &this->ans->engv, &eflag, &vflag, + &ainum, &nbor_pitch, &this->_threads_per_atom); + } + this->time_pair.stop(); +} + +template class LJCubic; diff --git a/lib/gpu/lal_lj_cubic.h b/lib/gpu/lal_lj_cubic.h new file mode 100644 index 0000000000..0fefc727eb --- /dev/null +++ b/lib/gpu/lal_lj_cubic.h @@ -0,0 +1,81 @@ +/*************************************************************************** + lj_cubic.h + ------------------- + Trung Dac Nguyen + + Class for acceleration of the lj/cubic pair style. + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : + email : ndactrung@gmail.com + ***************************************************************************/ + +#ifndef LAL_LJ_H +#define LAL_LJ_H + +#include "lal_base_atomic.h" + +namespace LAMMPS_AL { + +template +class LJCubic : public BaseAtomic { + public: + LJCubic(); + ~LJCubic(); + + /// Clear any previous data and set up for a new LAMMPS run + /** \param max_nbors initial number of rows in the neighbor matrix + * \param cell_size cutoff + skin + * \param gpu_split fraction of particles handled by device + * + * Returns: + * - 0 if successfull + * - -1 if fix gpu not found + * - -3 if there is an out of memory error + * - -4 if the GPU library was not compiled for GPU + * - -5 Double precision is not supported on card **/ + int init(const int ntypes, double **host_cutsq, double **host_cut_inner_sq, + double **host_cut_inner, double **host_sigma, double **host_epsilon, + double **host_lj1, double **host_lj2, double **host_lj3, + double **host_lj4, double *host_special_lj, const int nlocal, + const int nall, const int max_nbors, const int maxspecial, + const double cell_size, const double gpu_split, FILE *screen); + + /// Clear all host and device data + /** \note This is called at the beginning of the init() routine **/ + void clear(); + + /// Returns memory usage on device per atom + int bytes_per_atom(const int max_nbors) const; + + /// Total host memory used by library for pair style + double host_memory_usage() const; + + // --------------------------- TYPE DATA -------------------------- + + /// lj1.x = lj1, lj1.y = lj2, lj1.z = cutsq + UCL_D_Vec lj1; + /// lj2.x = cut_inner_sq, lj2.y = cut_inner, lj2.z = sigma, lj2.w = epsilon + UCL_D_Vec lj2; + /// lj3.x = lj3, lj3.y = lj4 + UCL_D_Vec lj3; + /// Special LJ values + UCL_D_Vec sp_lj; + + /// If atom type constants fit in shared memory, use fast kernels + bool shared_types; + + /// Number of atom types + int _lj_types; + + private: + bool _allocated; + void loop(const bool _eflag, const bool _vflag); +}; + +} + +#endif diff --git a/lib/gpu/lal_lj_cubic_ext.cpp b/lib/gpu/lal_lj_cubic_ext.cpp new file mode 100644 index 0000000000..518f706781 --- /dev/null +++ b/lib/gpu/lal_lj_cubic_ext.cpp @@ -0,0 +1,124 @@ +/*************************************************************************** + lj_cubic_ext.cpp + ------------------- + Trung Dac Nguyen + + Functions for LAMMPS access to lj/cubic acceleration routines. + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : + email : ndactrung@gmail.com + ***************************************************************************/ + +#include +#include +#include + +#include "lal_lj_cubic.h" + +using namespace std; +using namespace LAMMPS_AL; + +static LJCubic LJCubicLMF; + +// --------------------------------------------------------------------------- +// Allocate memory on host and device and copy constants to device +// --------------------------------------------------------------------------- +int ljcb_gpu_init(const int ntypes, double **cutsq, double **cut_inner_sq, + double **cut_inner, double **sigma, double **epsilon, + double **host_lj1, double **host_lj2, double **host_lj3, + double **host_lj4, double *special_lj, + const int inum, const int nall, const int max_nbors, + const int maxspecial, const double cell_size, + int &gpu_mode, FILE *screen) { + LJCubicLMF.clear(); + gpu_mode=LJCubicLMF.device->gpu_mode(); + double gpu_split=LJCubicLMF.device->particle_split(); + int first_gpu=LJCubicLMF.device->first_device(); + int last_gpu=LJCubicLMF.device->last_device(); + int world_me=LJCubicLMF.device->world_me(); + int gpu_rank=LJCubicLMF.device->gpu_rank(); + int procs_per_gpu=LJCubicLMF.device->procs_per_gpu(); + + LJCubicLMF.device->init_message(screen,"lj/cubic",first_gpu,last_gpu); + + bool message=false; + if (LJCubicLMF.device->replica_me()==0 && screen) + message=true; + + if (message) { + fprintf(screen,"Initializing Device and compiling on process 0..."); + fflush(screen); + } + + int init_ok=0; + if (world_me==0) + init_ok=LJCubicLMF.init(ntypes, cutsq, cut_inner_sq, cut_inner, sigma, + epsilon, host_lj1, host_lj2, host_lj3, host_lj4, + special_lj, inum, nall, 300, maxspecial, + cell_size, gpu_split, screen); + + LJCubicLMF.device->world_barrier(); + if (message) + fprintf(screen,"Done.\n"); + + for (int i=0; igpu_barrier(); + if (message) + fprintf(screen,"Done.\n"); + } + if (message) + fprintf(screen,"\n"); + + if (init_ok==0) + LJCubicLMF.estimate_gpu_overhead(); + return init_ok; +} + +void ljcb_gpu_clear() { + LJCubicLMF.clear(); +} + +int ** ljcb_gpu_compute_n(const int ago, const int inum_full, + const int nall, double **host_x, int *host_type, + double *sublo, double *subhi, tagint *tag, int **nspecial, + tagint **special, const bool eflag, const bool vflag, + const bool eatom, const bool vatom, int &host_start, + int **ilist, int **jnum, const double cpu_time, + bool &success) { + return LJCubicLMF.compute(ago, inum_full, nall, host_x, host_type, sublo, + subhi, tag, nspecial, special, eflag, vflag, eatom, + vatom, host_start, ilist, jnum, cpu_time, success); +} + +void ljcb_gpu_compute(const int ago, const int inum_full, const int nall, + double **host_x, int *host_type, int *ilist, int *numj, + int **firstneigh, const bool eflag, const bool vflag, + const bool eatom, const bool vatom, int &host_start, + const double cpu_time, bool &success) { + LJCubicLMF.compute(ago,inum_full,nall,host_x,host_type,ilist,numj, + firstneigh,eflag,vflag,eatom,vatom,host_start,cpu_time,success); +} + +double ljcb_gpu_bytes() { + return LJCubicLMF.host_memory_usage(); +} + + diff --git a/lib/gpu/lal_tersoff.cpp b/lib/gpu/lal_tersoff.cpp new file mode 100644 index 0000000000..ebf3b2fcb1 --- /dev/null +++ b/lib/gpu/lal_tersoff.cpp @@ -0,0 +1,451 @@ +/*************************************************************************** + tersoff.cpp + ------------------- + Trung Dac Nguyen + + Class for acceleration of the tersoff pair style. + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : Thu April 17, 2014 + email : ndactrung@gmail.com + ***************************************************************************/ + +#if defined(USE_OPENCL) +#include "tersoff_cl.h" +#elif defined(USE_CUDART) +const char *tersoff=0; +#else +#include "tersoff_cubin.h" +#endif + +#include "lal_tersoff.h" +#include +using namespace LAMMPS_AL; +#define TersoffT Tersoff + +extern Device device; + +template +TersoffT::Tersoff() : BaseThree(), _allocated(false) { +} + +template +TersoffT::~Tersoff() { + clear(); +} + +template +int TersoffT::bytes_per_atom(const int max_nbors) const { + return this->bytes_per_atom_atomic(max_nbors); +} + +template +int TersoffT::init(const int ntypes, const int nlocal, const int nall, const int max_nbors, + const double cell_size, const double gpu_split, FILE *_screen, + int* host_map, const int nelements, int*** host_elem2param, const int nparams, + const double* lam1, const double* lam2, const double* lam3,const double* powermint, + const double* biga, const double* bigb, const double* bigr, const double* bigd, + const double* c1, const double* c2, const double* c3, const double* c4, + const double* c, const double* d, const double* h, const double* gamma, + const double* beta, const double* powern, const double* host_cutsq) +{ + int success; + success=this->init_three(nlocal,nall,max_nbors,0,cell_size,gpu_split, + _screen,tersoff,"k_tersoff_repulsive", + "k_tersoff_three_center", "k_tersoff_three_end"); + if (success!=0) + return success; + + int ef_nall=nall; + if (ef_nall==0) + ef_nall=2000; + + _max_zij_size=static_cast(static_cast(ef_nall)*1.10); + _zetaij.alloc(_max_zij_size*max_nbors,*(this->ucl_device),UCL_READ_WRITE); + zeta_tex.get_texture(*(this->pair_program),"zeta_tex"); + zeta_tex.bind_float(_zetaij,1); + + k_zeta.set_function(*(this->pair_program),"k_tersoff_zeta"); + + // If atom type constants fit in shared memory use fast kernel + int lj_types=ntypes; + shared_types=false; + int max_shared_types=this->device->max_shared_types(); + if (lj_types<=max_shared_types && this->_block_size>=max_shared_types) { + lj_types=max_shared_types; + shared_types=true; + } + _lj_types=lj_types; + + _nparams = nparams; + _nelements = nelements; + + UCL_H_Vec dview(nparams,*(this->ucl_device), + UCL_WRITE_ONLY); + + for (int i=0; iucl_device),UCL_READ_ONLY); + + for (int i=0; i(lam1[i]); + dview[i].y=static_cast(lam2[i]); + dview[i].z=static_cast(lam3[i]); + dview[i].w=static_cast(powermint[i]); + } + + ucl_copy(ts1,dview,false); + ts1_tex.get_texture(*(this->pair_program),"ts1_tex"); + ts1_tex.bind_float(ts1,4); + + ts2.alloc(nparams,*(this->ucl_device),UCL_READ_ONLY); + + for (int i=0; i(biga[i]); + dview[i].y=static_cast(bigb[i]); + dview[i].z=static_cast(bigr[i]); + dview[i].w=static_cast(bigd[i]); + } + + ucl_copy(ts2,dview,false); + ts2_tex.get_texture(*(this->pair_program),"ts2_tex"); + ts2_tex.bind_float(ts2,4); + + ts3.alloc(nparams,*(this->ucl_device),UCL_READ_ONLY); + + for (int i=0; i(c1[i]); + dview[i].y=static_cast(c2[i]); + dview[i].z=static_cast(c3[i]); + dview[i].w=static_cast(c4[i]); + } + + ucl_copy(ts3,dview,false); + ts3_tex.get_texture(*(this->pair_program),"ts3_tex"); + ts3_tex.bind_float(ts3,4); + + ts4.alloc(nparams,*(this->ucl_device),UCL_READ_ONLY); + + for (int i=0; i(c[i]); + dview[i].y=static_cast(d[i]); + dview[i].z=static_cast(h[i]); + dview[i].w=static_cast(gamma[i]); + } + + ucl_copy(ts4,dview,false); + ts4_tex.get_texture(*(this->pair_program),"ts4_tex"); + ts4_tex.bind_float(ts4,4); + + ts5.alloc(nparams,*(this->ucl_device),UCL_READ_ONLY); + + for (int i=0; i(beta[i]); + dview[i].y=static_cast(powern[i]); + dview[i].z=(numtyp)0; + dview[i].w=(numtyp)0; + } + + ucl_copy(ts5,dview,false); + ts5_tex.get_texture(*(this->pair_program),"ts5_tex"); + ts5_tex.bind_float(ts5,4); + + UCL_H_Vec cutsq_view(nparams,*(this->ucl_device), + UCL_WRITE_ONLY); + for (int i=0; i(host_cutsq[i]); + cutsq.alloc(nparams,*(this->ucl_device),UCL_READ_ONLY); + ucl_copy(cutsq,cutsq_view,false); + + UCL_H_Vec dview_elem2param(nelements*nelements*nelements, + *(this->ucl_device), UCL_WRITE_ONLY); + + elem2param.alloc(nelements*nelements*nelements,*(this->ucl_device), + UCL_READ_ONLY); + + for (int i = 0; i < nelements; i++) + for (int j = 0; j < nelements; j++) + for (int k = 0; k < nelements; k++) { + int idx = i*nelements*nelements+j*nelements+k; + dview_elem2param[idx] = host_elem2param[i][j][k]; + } + + ucl_copy(elem2param,dview_elem2param,false); + + UCL_H_Vec dview_map(lj_types, *(this->ucl_device), UCL_WRITE_ONLY); + for (int i = 0; i < lj_types; i++) + dview_map[i] = host_map[i]; + + map.alloc(lj_types,*(this->ucl_device), UCL_READ_ONLY); + ucl_copy(map,dview_map,false); + + _allocated=true; + this->_max_bytes=ts1.row_bytes()+ts2.row_bytes()+ts3.row_bytes()+ + ts4.row_bytes()+ts5.row_bytes()+cutsq.row_bytes()+ + map.row_bytes()+elem2param.row_bytes()+_zetaij.row_bytes(); + return 0; +} + +template +void TersoffT::clear() { + if (!_allocated) + return; + _allocated=false; + + ts1.clear(); + ts2.clear(); + ts3.clear(); + ts4.clear(); + ts5.clear(); + cutsq.clear(); + map.clear(); + elem2param.clear(); + _zetaij.clear(); + + k_zeta.clear(); + + this->clear_atomic(); +} + +template +double TersoffT::host_memory_usage() const { + return this->host_memory_usage_atomic()+sizeof(Tersoff); +} + +#define KTHREADS this->_threads_per_atom +#define JTHREADS this->_threads_per_atom +// --------------------------------------------------------------------------- +// Copy nbor list from host if necessary and then calculate forces, virials,.. +// --------------------------------------------------------------------------- +template +void TersoffT::compute(const int f_ago, const int nlocal, const int nall, + const int nlist, double **host_x, int *host_type, + int *ilist, int *numj, int **firstneigh, + const bool eflag, const bool vflag, const bool eatom, + const bool vatom, int &host_start, + const double cpu_time, bool &success) { + this->acc_timers(); + if (nlist==0) { + host_start=0; + // Make sure textures are correct if realloc by a different hybrid style + this->resize_atom(0,nall,success); + this->zero_timers(); + return; + } + + int ago=this->hd_balancer.ago_first(f_ago); + int inum=this->hd_balancer.balance(ago,nlocal,cpu_time); + this->ans->inum(inum); + #ifdef THREE_CONCURRENT + this->ans2->inum(inum); + #endif + host_start=inum; + + if (ago==0) { + this->reset_nbors(nall, inum, nlist, ilist, numj, firstneigh, success); + if (!success) + return; + } + + this->atom->cast_x_data(host_x,host_type); + this->hd_balancer.start_timer(); + this->atom->add_x_data(host_x,host_type); + + // re-allocate zetaij if necessary + if (nall>_max_zij_size) { + _max_nbors = this->nbor->max_nbor_loop(nlist,numj,ilist); + _max_zij_size=static_cast(static_cast(nall)*1.10); + _zetaij.resize(_max_nbors*_max_zij_size); + zeta_tex.bind_float(_zetaij,1); + } + + int ainum=nall; + int nbor_pitch=this->nbor->nbor_pitch(); + int BX=this->block_pair(); + int GX=static_cast(ceil(static_cast(ainum)/ + (BX/this->_threads_per_atom))); + + this->k_zeta.set_size(GX,BX); + this->k_zeta.run(&this->atom->x, &ts1, &ts2, &ts4, &cutsq, + &map, &elem2param, &_nelements, &_nparams, &_zetaij, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &nall, &ainum, &nbor_pitch, &this->_threads_per_atom); + + int evatom=0; + if (eatom || vatom) + evatom=1; + #ifdef THREE_CONCURRENT + this->ucl_device->sync(); + #endif + loop(eflag,vflag,evatom); + this->ans->copy_answers(eflag,vflag,eatom,vatom,ilist); + this->device->add_ans_object(this->ans); + #ifdef THREE_CONCURRENT + this->ans2->copy_answers(eflag,vflag,eatom,vatom,ilist); + this->device->add_ans_object(this->ans2); + #endif + this->hd_balancer.stop_timer(); +} + +// --------------------------------------------------------------------------- +// Reneighbor on GPU if necessary and then compute forces, virials, energies +// --------------------------------------------------------------------------- +template +int ** TersoffT::compute(const int ago, const int inum_full, + const int nall, double **host_x, int *host_type, + double *sublo, double *subhi, tagint *tag, + int **nspecial, tagint **special, const bool eflag, + const bool vflag, const bool eatom, + const bool vatom, int &host_start, + int **ilist, int **jnum, + const double cpu_time, bool &success) { + this->acc_timers(); + + if (inum_full==0) { + host_start=0; + // Make sure textures are correct if realloc by a different hybrid style + this->resize_atom(0,nall,success); + this->zero_timers(); + return NULL; + } + + this->hd_balancer.balance(cpu_time); + int inum=this->hd_balancer.get_gpu_count(ago,inum_full); + this->ans->inum(inum); + #ifdef THREE_CONCURRENT + this->ans2->inum(inum); + #endif + host_start=inum; + + // Build neighbor list on GPU if necessary + if (ago==0) { + _max_nbors = this->build_nbor_list(inum, inum_full-inum, nall, host_x, host_type, + sublo, subhi, tag, nspecial, special, success); + if (!success) + return NULL; + this->hd_balancer.start_timer(); + } else { + this->atom->cast_x_data(host_x,host_type); + this->hd_balancer.start_timer(); + this->atom->add_x_data(host_x,host_type); + } + *ilist=this->nbor->host_ilist.begin(); + *jnum=this->nbor->host_acc.begin(); + + // re-allocate zetaij if necessary + if (nall>_max_zij_size) { + _max_zij_size=static_cast(static_cast(nall)*1.10); + _zetaij.resize(_max_nbors*_max_zij_size); + zeta_tex.bind_float(_zetaij,1); + } + + int ainum=nall; + int nbor_pitch=this->nbor->nbor_pitch(); + int BX=this->block_pair(); + int GX=static_cast(ceil(static_cast(ainum)/ + (BX/this->_threads_per_atom))); + + this->k_zeta.set_size(GX,BX); + this->k_zeta.run(&this->atom->x, &ts1, &ts2, &ts4, &cutsq, &map, + &elem2param, &_nelements, &_nparams, &_zetaij, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &nall, &ainum, &nbor_pitch, &this->_threads_per_atom); + + int evatom=0; + if (eatom || vatom) + evatom=1; + #ifdef THREE_CONCURRENT + this->ucl_device->sync(); + #endif + loop(eflag,vflag,evatom); + this->ans->copy_answers(eflag,vflag,eatom,vatom); + this->device->add_ans_object(this->ans); + #ifdef THREE_CONCURRENT + this->ans2->copy_answers(eflag,vflag,eatom,vatom); + this->device->add_ans_object(this->ans2); + #endif + this->hd_balancer.stop_timer(); + + return this->nbor->host_jlist.begin()-host_start; +} + +// --------------------------------------------------------------------------- +// Calculate energies, forces, and torques +// --------------------------------------------------------------------------- +template +void TersoffT::loop(const bool _eflag, const bool _vflag, const int evatom) { + // Compute the block size and grid size to keep all cores busy + int BX=this->block_pair(); + int eflag, vflag; + if (_eflag) + eflag=1; + else + eflag=0; + + if (_vflag) + vflag=1; + else + vflag=0; + + int ainum=this->ans->inum(); + int nbor_pitch=this->nbor->nbor_pitch(); + int GX=static_cast(ceil(static_cast(this->ans->inum())/ + (BX/this->_threads_per_atom))); + + this->time_pair.start(); + this->k_pair.set_size(GX,BX); + this->k_pair.run(&this->atom->x, &ts1, &ts2, &cutsq, + &map, &elem2param, &_nelements, &_nparams, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &this->ans->force, &this->ans->engv, + &eflag, &vflag, &ainum, &nbor_pitch, + &this->_threads_per_atom); + + BX=this->block_size(); + GX=static_cast(ceil(static_cast(this->ans->inum())/ + (BX/(KTHREADS*JTHREADS)))); + this->k_three_center.set_size(GX,BX); + this->k_three_center.run(&this->atom->x, &ts1, &ts2, &ts3, &ts4, &ts5, &cutsq, + &map, &elem2param, &_nelements, &_nparams, &_zetaij, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &this->ans->force, &this->ans->engv, &eflag, &vflag, &ainum, + &nbor_pitch, &this->_threads_per_atom, &evatom); + + Answer *end_ans; + #ifdef THREE_CONCURRENT + end_ans=this->ans2; + #else + end_ans=this->ans; + #endif + if (evatom!=0) { + this->k_three_end_vatom.set_size(GX,BX); + this->k_three_end_vatom.run(&this->atom->x, &ts1, &ts2, &ts3, &ts4, &ts5, &cutsq, + &map, &elem2param, &_nelements, &_nparams, &_zetaij, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum, + &nbor_pitch, &this->_threads_per_atom); + + } else { + this->k_three_end.set_size(GX,BX); + this->k_three_end.run(&this->atom->x, &ts1, &ts2, &ts3, &ts4, &ts5, &cutsq, + &map, &elem2param, &_nelements, &_nparams, &_zetaij, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum, + &nbor_pitch, &this->_threads_per_atom); + + } + + this->time_pair.stop(); +} + +template class Tersoff; + diff --git a/lib/gpu/lal_tersoff.h b/lib/gpu/lal_tersoff.h new file mode 100644 index 0000000000..5e760297ba --- /dev/null +++ b/lib/gpu/lal_tersoff.h @@ -0,0 +1,116 @@ +/*************************************************************************** + tersoff.h + ------------------- + Trung Dac Nguyen + + Class for acceleration of the tersoff pair style. + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : Thu April 17, 2014 + email : ndactrung@gmail.com + ***************************************************************************/ + +#ifndef LAL_TERSOFF_H +#define LAL_TERSOFF_H + +#include "lal_base_three.h" + +namespace LAMMPS_AL { + +template +class Tersoff : public BaseThree { + public: + Tersoff(); + ~Tersoff(); + + /// Clear any previous data and set up for a new LAMMPS run for generic systems + /** \param max_nbors initial number of rows in the neighbor matrix + * \param cell_size cutoff + skin + * \param gpu_split fraction of particles handled by device + * + * Returns: + * - 0 if successfull + * - -1 if fix gpu not found + * - -3 if there is an out of memory error + * - -4 if the GPU library was not compiled for GPU + * - -5 Double precision is not supported on card **/ + int init(const int ntypes, const int nlocal, const int nall, const int max_nbors, + const double cell_size, const double gpu_split, FILE *screen, + int* host_map, const int nelements, int*** host_elem2param, const int nparams, + const double* lam1, const double* lam2, const double* lam3, + const double* powermint, const double* biga, const double* bigb, + const double* bigr, const double* bigd, const double* c1, const double* c2, + const double* c3, const double* c4, const double* c, const double* d, + const double* h, const double* gamma, const double* beta, + const double* powern, const double* cutsq); + + /// Pair loop with host neighboring + void compute(const int f_ago, const int inum_full, const int nall, + const int nlist, double **host_x, int *host_type, + int *ilist, int *numj, int **firstneigh, const bool eflag, + const bool vflag, const bool eatom, const bool vatom, + int &host_start, const double cpu_time, bool &success); + + /// Pair loop with device neighboring + int ** compute(const int ago, const int inum_full, + const int nall, double **host_x, int *host_type, double *sublo, + double *subhi, tagint *tag, int **nspecial, + tagint **special, const bool eflag, const bool vflag, + const bool eatom, const bool vatom, int &host_start, + int **ilist, int **numj, const double cpu_time, bool &success); + + /// Clear all host and device data + /** \note This is called at the beginning of the init() routine **/ + void clear(); + + /// Returns memory usage on device per atom + int bytes_per_atom(const int max_nbors) const; + + /// Total host memory used by library for pair style + double host_memory_usage() const; + + // --------------------------- TYPE DATA -------------------------- + + /// If atom type constants fit in shared memory, use fast kernels + bool shared_types; + + /// Number of atom types + int _lj_types; + + /// ts1.x = lam1, ts1.y = lam2, ts1.z = lam3, ts1.w = powermint + UCL_D_Vec ts1; + /// ts2.x = biga, ts2.y = bigb, ts2.z = bigr, ts2.w = bigd + UCL_D_Vec ts2; + /// ts3.x = c1, ts3.y = c2, ts3.z = c3, ts3.w = c4 + UCL_D_Vec ts3; + /// ts4.x = c, ts4.y = d, ts4.z = h, ts4.w = gamma + UCL_D_Vec ts4; + /// ts5.x = beta, ts5.y = powern + UCL_D_Vec ts5; + + UCL_D_Vec cutsq; + + UCL_D_Vec elem2param; + UCL_D_Vec map; + int _nparams,_nelements; + + /// Per-atom arrays + UCL_D_Vec _zetaij; + + UCL_Kernel k_zeta; + UCL_Texture ts1_tex, ts2_tex, ts3_tex, ts4_tex, ts5_tex, zeta_tex; + + int _max_zij_size, _max_nbors; + + private: + bool _allocated; + void loop(const bool _eflag, const bool _vflag, const int evatom); +}; + +} + +#endif + diff --git a/lib/gpu/lal_tersoff_ext.cpp b/lib/gpu/lal_tersoff_ext.cpp new file mode 100644 index 0000000000..c33ea6615e --- /dev/null +++ b/lib/gpu/lal_tersoff_ext.cpp @@ -0,0 +1,139 @@ +/*************************************************************************** + tersoff_ext.cpp + ------------------- + Trung Dac Nguyen + + Functions for LAMMPS access to tersoff acceleration routines. + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : Thu April 17, 2014 + email : ndactrung@gmail.com + ***************************************************************************/ + +#include +#include +#include + +#include "lal_tersoff.h" + +using namespace std; +using namespace LAMMPS_AL; + +static Tersoff TSMF; + +// --------------------------------------------------------------------------- +// Allocate memory on host and device and copy constants to device +// --------------------------------------------------------------------------- +int tersoff_gpu_init(const int ntypes, const int inum, const int nall, const int max_nbors, + const double cell_size, int &gpu_mode, FILE *screen, + int* host_map, const int nelements, int*** host_elem2param, const int nparams, + const double* ts_lam1, const double* ts_lam2, const double* ts_lam3, + const double* ts_powermint, const double* ts_biga, const double* ts_bigb, + const double* ts_bigr, const double* ts_bigd, + const double* ts_c1, const double* ts_c2, const double* ts_c3, const double* ts_c4, + const double* ts_c, const double* ts_d, const double* ts_h, + const double* ts_gamma, const double* ts_beta, + const double* ts_powern, const double* ts_cutsq) { + TSMF.clear(); + gpu_mode=TSMF.device->gpu_mode(); + double gpu_split=TSMF.device->particle_split(); + int first_gpu=TSMF.device->first_device(); + int last_gpu=TSMF.device->last_device(); + int world_me=TSMF.device->world_me(); + int gpu_rank=TSMF.device->gpu_rank(); + int procs_per_gpu=TSMF.device->procs_per_gpu(); + + // disable host/device split for now + if (gpu_split != 1.0) + return -8; + + // disable multiple threads per atom for now + if (TSMF.device->threads_per_atom() != 1) + return -10; + + TSMF.device->init_message(screen,"tersoff/gpu",first_gpu,last_gpu); + + bool message=false; + if (TSMF.device->replica_me()==0 && screen) + message=true; + + if (message) { + fprintf(screen,"Initializing Device and compiling on process 0..."); + fflush(screen); + } + + int init_ok=0; + if (world_me==0) + init_ok=TSMF.init(ntypes, inum, nall, 300, cell_size, gpu_split, screen, + host_map, nelements, host_elem2param, nparams, + ts_lam1, ts_lam2, ts_lam3, ts_powermint, + ts_biga, ts_bigb, ts_bigr, ts_bigd, + ts_c1, ts_c2, ts_c3, ts_c4, ts_c, ts_d, ts_h, + ts_gamma, ts_beta, ts_powern, ts_cutsq); + + TSMF.device->world_barrier(); + if (message) + fprintf(screen,"Done.\n"); + + for (int i=0; igpu_barrier(); + if (message) + fprintf(screen,"Done.\n"); + } + if (message) + fprintf(screen,"\n"); + + if (init_ok==0) + TSMF.estimate_gpu_overhead(); + return init_ok; +} + +void tersoff_gpu_clear() { + TSMF.clear(); +} + +int ** tersoff_gpu_compute_n(const int ago, const int inum_full, + const int nall, double **host_x, int *host_type, + double *sublo, double *subhi, tagint *tag, int **nspecial, + tagint **special, const bool eflag, const bool vflag, + const bool eatom, const bool vatom, int &host_start, + int **ilist, int **jnum, const double cpu_time, + bool &success) { + return TSMF.compute(ago, inum_full, nall, host_x, host_type, sublo, + subhi, tag, nspecial, special, eflag, vflag, eatom, + vatom, host_start, ilist, jnum, cpu_time, success); +} + +void tersoff_gpu_compute(const int ago, const int nlocal, const int nall, + const int nlist, double **host_x, int *host_type, + int *ilist, int *numj, int **firstneigh, const bool eflag, + const bool vflag, const bool eatom, const bool vatom, + int &host_start, const double cpu_time, bool &success) { + TSMF.compute(ago,nlocal,nall,nlist,host_x,host_type,ilist,numj, + firstneigh,eflag,vflag,eatom,vatom,host_start,cpu_time,success); +} + +double tersoff_gpu_bytes() { + return TSMF.host_memory_usage(); +} + + diff --git a/lib/gpu/lal_tersoff_extra.h b/lib/gpu/lal_tersoff_extra.h new file mode 100644 index 0000000000..7a8fc5872a --- /dev/null +++ b/lib/gpu/lal_tersoff_extra.h @@ -0,0 +1,612 @@ +/// ************************************************************************** +// tersoff_extra.h +// ------------------- +// Trung Dac Nguyen +// +// Device code for Tersoff math routines +// +// __________________________________________________________________________ +// This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) +// __________________________________________________________________________ +// +// begin : +// email : ndactrung@gmail.com +// ***************************************************************************/* + +#ifndef LAL_TERSOFF_EXTRA_H +#define LAL_TERSOFF_EXTRA_H + +#ifdef NV_KERNEL +#include "lal_aux_fun1.h" +#else +#endif + +#define MY_PI2 (numtyp)1.57079632679489661923 +#define MY_PI4 (numtyp)0.78539816339744830962 + +/* ---------------------------------------------------------------------- */ + +ucl_inline numtyp vec3_dot(const numtyp x[3], const numtyp y[3]) +{ + return (x[0]*y[0] + x[1]*y[1] + x[2]*y[2]); +} + +ucl_inline void vec3_add(const numtyp x[3], const numtyp y[3], numtyp z[3]) +{ + z[0] = x[0]+y[0]; z[1] = x[1]+y[1]; z[2] = x[2]+y[2]; +} + +ucl_inline void vec3_scale(const numtyp k, const numtyp x[3], numtyp y[3]) +{ + y[0] = k*x[0]; y[1] = k*x[1]; y[2] = k*x[2]; +} + +ucl_inline void vec3_scaleadd(const numtyp k, const numtyp x[3], + const numtyp y[3], numtyp z[3]) +{ + z[0] = k*x[0]+y[0]; z[1] = k*x[1]+y[1]; z[2] = k*x[2]+y[2]; +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline numtyp ters_gijk(const numtyp costheta, + const numtyp param_c, + const numtyp param_d, + const numtyp param_h, + const numtyp param_gamma) +{ + const numtyp ters_c = param_c * param_c; + const numtyp ters_d = param_d * param_d; + const numtyp hcth = param_h - costheta; + return param_gamma*((numtyp)1.0 + ters_c*ucl_recip(ters_d) - + ters_c *ucl_recip(ters_d + hcth*hcth)); +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline numtyp ters_gijk_d(const numtyp costheta, + const numtyp param_c, + const numtyp param_d, + const numtyp param_h, + const numtyp param_gamma) +{ + const numtyp ters_c = param_c * param_c; + const numtyp ters_d = param_d * param_d; + const numtyp hcth = param_h - costheta; + const numtyp numerator = (numtyp)-2.0 * ters_c * hcth; + const numtyp denominator = ucl_recip(ters_d + hcth*hcth); + return param_gamma*numerator*denominator*denominator; +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline void costheta_d(const numtyp rij_hat[3], + const numtyp rij, + const numtyp rik_hat[3], + const numtyp rik, + numtyp *dri, + numtyp *drj, + numtyp *drk) +{ + // first element is derivative wrt Ri, second wrt Rj, third wrt Rk + + numtyp cos_theta = vec3_dot(rij_hat,rik_hat); + + vec3_scaleadd(-cos_theta,rij_hat,rik_hat,drj); + vec3_scale(ucl_recip(rij),drj,drj); + vec3_scaleadd(-cos_theta,rik_hat,rij_hat,drk); + vec3_scale(ucl_recip(rik),drk,drk); + vec3_add(drj,drk,dri); + vec3_scale((numtyp)-1.0,dri,dri); +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline numtyp ters_fc(const numtyp r, + const numtyp param_bigr, + const numtyp param_bigd) +{ + if (r < param_bigr-param_bigd) return (numtyp)1.0; + if (r > param_bigr+param_bigd) return (numtyp)0.0; + return (numtyp)0.5*((numtyp)1.0 - sin(MY_PI2*(r - param_bigr)/param_bigd)); +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline numtyp ters_fc_d(const numtyp r, + const numtyp param_bigr, + const numtyp param_bigd) +{ + if (r < param_bigr-param_bigd) return (numtyp)0.0; + if (r > param_bigr+param_bigd) return (numtyp)0.0; + return -(MY_PI4/param_bigd) * cos(MY_PI2*(r - param_bigr)/param_bigd); +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline numtyp ters_fa(const numtyp r, + const numtyp param_bigb, + const numtyp param_bigr, + const numtyp param_bigd, + const numtyp param_lam2) +{ + if (r > param_bigr + param_bigd) return (numtyp)0.0; + return -param_bigb * ucl_exp(-param_lam2 * r) * + ters_fc(r,param_bigr,param_bigd); +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline numtyp ters_fa_d(const numtyp r, + const numtyp param_bigb, + const numtyp param_bigr, + const numtyp param_bigd, + const numtyp param_lam2) +{ + if (r > param_bigr + param_bigd) return (numtyp)0.0; + return param_bigb * ucl_exp(-param_lam2 * r) * (param_lam2 * + ters_fc(r,param_bigr,param_bigd) - ters_fc_d(r,param_bigr,param_bigd)); +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline numtyp ters_bij(const numtyp zeta, + const numtyp param_beta, + const numtyp param_powern, + const numtyp param_c1, + const numtyp param_c2, + const numtyp param_c3, + const numtyp param_c4) +{ + numtyp tmp = param_beta * zeta; + if (tmp > param_c1) return ucl_rsqrt(tmp); + if (tmp > param_c2) + return ((numtyp)1.0 - ucl_powr(tmp,-param_powern) / + ((numtyp)2.0*param_powern))*ucl_rsqrt(tmp); + if (tmp < param_c4) return (numtyp)1.0; + if (tmp < param_c3) + return (numtyp)1.0 - ucl_powr(tmp,param_powern)/((numtyp)2.0*param_powern); + return ucl_powr((numtyp)1.0 + ucl_powr(tmp,param_powern), + (numtyp)-1.0/((numtyp)2.0*param_powern)); +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline numtyp ters_bij_d(const numtyp zeta, + const numtyp param_beta, + const numtyp param_powern, + const numtyp param_c1, + const numtyp param_c2, + const numtyp param_c3, + const numtyp param_c4) +{ + numtyp tmp = param_beta * zeta; + if (tmp > param_c1) + return param_beta * (numtyp)-0.5*ucl_powr(tmp,(numtyp)-1.5); + if (tmp > param_c2) + return param_beta * ((numtyp)-0.5*ucl_powr(tmp,(numtyp)-1.5) * + ((numtyp)1.0 - (numtyp)0.5 * ((numtyp)1.0 + (numtyp)1.0 / + ((numtyp)2.0 * param_powern)) * ucl_powr(tmp,-param_powern))); + if (tmp < param_c4) return (numtyp)0.0; + if (tmp < param_c3) + return (numtyp)-0.5*param_beta * ucl_powr(tmp,param_powern-(numtyp)1.0); + + numtyp tmp_n = ucl_powr(tmp,param_powern); + return (numtyp)-0.5 * ucl_powr((numtyp)1.0+tmp_n, (numtyp) - + (numtyp)1.0-((numtyp)1.0 / ((numtyp)2.0 * param_powern)))*tmp_n / zeta; +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline void ters_zetaterm_d(const numtyp prefactor, + const numtyp rij_hat[3], + const numtyp rij, + const numtyp rik_hat[3], + const numtyp rik, + const numtyp param_bigr, + const numtyp param_bigd, + const numtyp param_powermint, + const numtyp param_lam3, + const numtyp param_c, + const numtyp param_d, + const numtyp param_h, + const numtyp param_gamma, + numtyp dri[3], + numtyp drj[3], + numtyp drk[3]) +{ + numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp; + numtyp dcosdri[3],dcosdrj[3],dcosdrk[3]; + + fc = ters_fc(rik,param_bigr,param_bigd); + dfc = ters_fc_d(rik,param_bigr,param_bigd); + + numtyp t = param_lam3*(rij-rik); + if ((int)param_powermint == 3) tmp = t*t*t; + else tmp = t; + + if (tmp > (numtyp)69.0776) ex_delr = (acctyp)1.e30; + else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0; + else ex_delr = ucl_exp(tmp); + + if ((int)param_powermint == 3) + ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr; + else ex_delr_d = param_lam3 * ex_delr; + + cos_theta = vec3_dot(rij_hat,rik_hat); + gijk = ters_gijk(cos_theta,param_c,param_d,param_h,param_gamma); + gijk_d = ters_gijk_d(cos_theta,param_c,param_d,param_h,param_gamma); + costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk); + + // compute the derivative wrt Ri + // dri = -dfc*gijk*ex_delr*rik_hat; + // dri += fc*gijk_d*ex_delr*dcosdri; + // dri += fc*gijk*ex_delr_d*(rik_hat - rij_hat); + + vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri); + vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri); + vec3_scaleadd(fc*gijk*ex_delr_d,rik_hat,dri,dri); + vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri); + vec3_scale(prefactor,dri,dri); + + // compute the derivative wrt Rj + // drj = fc*gijk_d*ex_delr*dcosdrj; + // drj += fc*gijk*ex_delr_d*rij_hat; + + vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj); + vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj); + vec3_scale(prefactor,drj,drj); + + // compute the derivative wrt Rk + // drk = dfc*gijk*ex_delr*rik_hat; + // drk += fc*gijk_d*ex_delr*dcosdrk; + // drk += -fc*gijk*ex_delr_d*rik_hat; + + vec3_scale(dfc*gijk*ex_delr,rik_hat,drk); + vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk); + vec3_scaleadd(-fc*gijk*ex_delr_d,rik_hat,drk,drk); + vec3_scale(prefactor,drk,drk); +} + +ucl_inline void ters_zetaterm_d_fi(const numtyp prefactor, + const numtyp rij_hat[3], + const numtyp rij, + const numtyp rik_hat[3], + const numtyp rik, + const numtyp param_bigr, + const numtyp param_bigd, + const numtyp param_powermint, + const numtyp param_lam3, + const numtyp param_c, + const numtyp param_d, + const numtyp param_h, + const numtyp param_gamma, + numtyp dri[3]) +{ + numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp; + numtyp dcosdri[3],dcosdrj[3],dcosdrk[3]; + + fc = ters_fc(rik,param_bigr,param_bigd); + dfc = ters_fc_d(rik,param_bigr,param_bigd); + + numtyp t = param_lam3*(rij-rik); + if ((int)param_powermint == 3) tmp = t*t*t; + else tmp = t; + + if (tmp > (numtyp)69.0776) ex_delr = (acctyp)1.e30; + else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0; + else ex_delr = ucl_exp(tmp); + + if ((int)param_powermint == 3) + ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr; + else ex_delr_d = param_lam3 * ex_delr; + + cos_theta = vec3_dot(rij_hat,rik_hat); + gijk = ters_gijk(cos_theta,param_c,param_d,param_h,param_gamma); + gijk_d = ters_gijk_d(cos_theta,param_c,param_d,param_h,param_gamma); + costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk); + + // compute the derivative wrt Ri + // dri = -dfc*gijk*ex_delr*rik_hat; + // dri += fc*gijk_d*ex_delr*dcosdri; + // dri += fc*gijk*ex_delr_d*(rik_hat - rij_hat); + + vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri); + vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri); + vec3_scaleadd(fc*gijk*ex_delr_d,rik_hat,dri,dri); + vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri); + vec3_scale(prefactor,dri,dri); +} + +ucl_inline void ters_zetaterm_d_fj(const numtyp prefactor, + const numtyp rij_hat[3], + const numtyp rij, + const numtyp rik_hat[3], + const numtyp rik, + const numtyp param_bigr, + const numtyp param_bigd, + const numtyp param_powermint, + const numtyp param_lam3, + const numtyp param_c, + const numtyp param_d, + const numtyp param_h, + const numtyp param_gamma, + numtyp drj[3]) +{ + numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,cos_theta,tmp; + numtyp dcosdri[3],dcosdrj[3],dcosdrk[3]; + + fc = ters_fc(rik,param_bigr,param_bigd); + + numtyp t = param_lam3*(rij-rik); + if ((int)param_powermint == 3) tmp = t*t*t; + else tmp = t; + + if (tmp > (numtyp)69.0776) ex_delr = (acctyp)1.e30; + else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0; + else ex_delr = ucl_exp(tmp); + + if ((int)param_powermint == 3) + ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr; + else ex_delr_d = param_lam3 * ex_delr; + + cos_theta = vec3_dot(rij_hat,rik_hat); + gijk = ters_gijk(cos_theta,param_c,param_d,param_h,param_gamma); + gijk_d = ters_gijk_d(cos_theta,param_c,param_d,param_h,param_gamma); + costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk); + + // compute the derivative wrt Rj + // drj = fc*gijk_d*ex_delr*dcosdrj; + // drj += fc*gijk*ex_delr_d*rij_hat; + + vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj); + vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj); + vec3_scale(prefactor,drj,drj); +} + +ucl_inline void ters_zetaterm_d_fk(const numtyp prefactor, + const numtyp rij_hat[3], + const numtyp rij, + const numtyp rik_hat[3], + const numtyp rik, + const numtyp param_bigr, + const numtyp param_bigd, + const numtyp param_powermint, + const numtyp param_lam3, + const numtyp param_c, + const numtyp param_d, + const numtyp param_h, + const numtyp param_gamma, + numtyp drk[3]) +{ + numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp; + numtyp dcosdri[3],dcosdrj[3],dcosdrk[3]; + + fc = ters_fc(rik,param_bigr,param_bigd); + dfc = ters_fc_d(rik,param_bigr,param_bigd); + + numtyp t = param_lam3*(rij-rik); + if ((int)param_powermint == 3) tmp = t*t*t; + else tmp = t; + + if (tmp > (numtyp)69.0776) ex_delr = (acctyp)1.e30; + else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0; + else ex_delr = ucl_exp(tmp); + + if ((int)param_powermint == 3) + ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr; + else ex_delr_d = param_lam3 * ex_delr; + + cos_theta = vec3_dot(rij_hat,rik_hat); + gijk = ters_gijk(cos_theta,param_c,param_d,param_h,param_gamma); + gijk_d = ters_gijk_d(cos_theta,param_c,param_d,param_h,param_gamma); + costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk); + + // compute the derivative wrt Rk + // drk = dfc*gijk*ex_delr*rik_hat; + // drk += fc*gijk_d*ex_delr*dcosdrk; + // drk += -fc*gijk*ex_delr_d*rik_hat; + + vec3_scale(dfc*gijk*ex_delr,rik_hat,drk); + vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk); + vec3_scaleadd(-fc*gijk*ex_delr_d,rik_hat,drk,drk); + vec3_scale(prefactor,drk,drk); +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline void repulsive(const numtyp param_bigr, + const numtyp param_bigd, + const numtyp param_lam1, + const numtyp param_biga, + const numtyp rsq, + const int eflag, + numtyp *ans) +{ + numtyp r,tmp_fc,tmp_fc_d,tmp_exp; + r = ucl_sqrt(rsq); + tmp_fc = ters_fc(r,param_bigr,param_bigd); + tmp_fc_d = ters_fc_d(r,param_bigr,param_bigd); + tmp_exp = ucl_exp(-param_lam1 * r); + // fforce + ans[0] = -param_biga*tmp_exp*(tmp_fc_d - tmp_fc*param_lam1)*ucl_recip(r); + // eng + if (eflag) ans[1] = tmp_fc * param_biga * tmp_exp; +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline numtyp zeta(const numtyp param_powermint, + const numtyp param_lam3, + const numtyp param_bigr, + const numtyp param_bigd, + const numtyp param_c, + const numtyp param_d, + const numtyp param_h, + const numtyp param_gamma, + const numtyp rsqij, + const numtyp rsqik, + const numtyp4 delrij, + const numtyp4 delrik) +{ + numtyp rij,rik,costheta,arg,ex_delr; + + rij = ucl_sqrt(rsqij); + rik = ucl_sqrt(rsqik); + costheta = (delrij.x*delrik.x + delrij.y*delrik.y + + delrij.z*delrik.z) / (rij*rik); + + numtyp t = param_lam3*(rij-rik); + if ((int)param_powermint == 3) arg = t*t*t; + else arg = t; + + if (arg > (numtyp)69.0776) ex_delr = (numtyp)1.e30; + else if (arg < (numtyp)-69.0776) ex_delr = (numtyp)0.0; + else ex_delr = ucl_exp(arg); + + return ters_fc(rik,param_bigr,param_bigd) * + ters_gijk(costheta,param_c, param_d, param_h, param_gamma) * ex_delr; +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline void force_zeta(const numtyp param_bigb, + const numtyp param_bigr, + const numtyp param_bigd, + const numtyp param_lam2, + const numtyp param_beta, + const numtyp param_powern, + const numtyp param_c1, + const numtyp param_c2, + const numtyp param_c3, + const numtyp param_c4, + const numtyp rsq, + const numtyp zeta_ij, + const int eflag, + numtyp fpfeng[4]) +{ + numtyp r,fa,fa_d,bij; + + r = ucl_sqrt(rsq); + fa = ters_fa(r,param_bigb,param_bigr,param_bigd,param_lam2); + fa_d = ters_fa_d(r,param_bigb,param_bigr,param_bigd,param_lam2); + bij = ters_bij(zeta_ij,param_beta,param_powern, + param_c1,param_c2, param_c3, param_c4); + fpfeng[0] = (numtyp)0.5*bij*fa_d * ucl_recip(r); // fforce + fpfeng[1] = (numtyp)-0.5*fa * ters_bij_d(zeta_ij,param_beta, param_powern, + param_c1,param_c2, param_c3, param_c4); // prefactor + if (eflag) fpfeng[2] = (numtyp)0.5*bij*fa; // eng +} + +/* ---------------------------------------------------------------------- + attractive term + use param_ij cutoff for rij test + use param_ijk cutoff for rik test +------------------------------------------------------------------------- */ + +ucl_inline void attractive(const numtyp param_bigr, + const numtyp param_bigd, + const numtyp param_powermint, + const numtyp param_lam3, + const numtyp param_c, + const numtyp param_d, + const numtyp param_h, + const numtyp param_gamma, + const numtyp prefactor, + const numtyp rij, + const numtyp rijinv, + const numtyp rik, + const numtyp rikinv, + const numtyp delrij[3], + const numtyp delrik[3], + numtyp fi[3], + numtyp fj[3], + numtyp fk[3]) +{ + numtyp rij_hat[3],rik_hat[3]; + vec3_scale(rijinv,delrij,rij_hat); + vec3_scale(rikinv,delrik,rik_hat); + ters_zetaterm_d(prefactor,rij_hat,rij,rik_hat,rik, + param_bigr, param_bigd, param_powermint, param_lam3, + param_c, param_d, param_h, param_gamma, fi, fj, fk); +} + +ucl_inline void attractive_fi(const numtyp param_bigr, + const numtyp param_bigd, + const numtyp param_powermint, + const numtyp param_lam3, + const numtyp param_c, + const numtyp param_d, + const numtyp param_h, + const numtyp param_gamma, + const numtyp prefactor, + const numtyp rij, + const numtyp rijinv, + const numtyp rik, + const numtyp rikinv, + const numtyp delrij[3], + const numtyp delrik[3], + numtyp fi[3]) +{ + numtyp rij_hat[3],rik_hat[3]; + vec3_scale(rijinv,delrij,rij_hat); + vec3_scale(rikinv,delrik,rik_hat); + ters_zetaterm_d_fi(prefactor,rij_hat,rij,rik_hat,rik, + param_bigr, param_bigd, param_powermint, param_lam3, + param_c, param_d, param_h, param_gamma, fi); +} + +ucl_inline void attractive_fj(const numtyp param_bigr, + const numtyp param_bigd, + const numtyp param_powermint, + const numtyp param_lam3, + const numtyp param_c, + const numtyp param_d, + const numtyp param_h, + const numtyp param_gamma, + const numtyp prefactor, + const numtyp rij, + const numtyp rijinv, + const numtyp rik, + const numtyp rikinv, + const numtyp delrij[3], + const numtyp delrik[3], + numtyp fj[3]) +{ + numtyp rij_hat[3],rik_hat[3]; + vec3_scale(rijinv,delrij,rij_hat); + vec3_scale(rikinv,delrik,rik_hat); + ters_zetaterm_d_fj(prefactor,rij_hat,rij,rik_hat,rik, + param_bigr, param_bigd, param_powermint, param_lam3, + param_c, param_d, param_h, param_gamma, fj); +} + +ucl_inline void attractive_fk(const numtyp param_bigr, + const numtyp param_bigd, + const numtyp param_powermint, + const numtyp param_lam3, + const numtyp param_c, + const numtyp param_d, + const numtyp param_h, + const numtyp param_gamma, + const numtyp prefactor, + const numtyp rij, + const numtyp rijinv, + const numtyp rik, + const numtyp rikinv, + const numtyp delrij[3], + const numtyp delrik[3], + numtyp fk[3]) +{ + numtyp rij_hat[3],rik_hat[3]; + vec3_scale(rijinv,delrij,rij_hat); + vec3_scale(rikinv,delrik,rik_hat); + ters_zetaterm_d_fk(prefactor,rij_hat,rij,rik_hat,rik, + param_bigr, param_bigd, param_powermint, param_lam3, + param_c, param_d, param_h, param_gamma, fk); +} + + +#endif + + diff --git a/lib/gpu/lal_zbl.cpp b/lib/gpu/lal_zbl.cpp new file mode 100644 index 0000000000..203def347a --- /dev/null +++ b/lib/gpu/lal_zbl.cpp @@ -0,0 +1,159 @@ +/*************************************************************************** + zbl.cpp + ------------------- + Trung Dac Nguyen + + Class for acceleration of the zbl pair style. + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : + email : ndactrung@gmail.com + ***************************************************************************/ + +#ifdef USE_OPENCL +#include "zbl_cl.h" +#elif defined(USE_CUDART) +const char *zbl=0; +#else +#include "zbl_cubin.h" +#endif + +#include "lal_zbl.h" +#include +using namespace LAMMPS_AL; +#define ZBLT ZBL + +extern Device device; + +template +ZBLT::ZBL() : BaseAtomic(), _allocated(false) { +} + +template +ZBLT::~ZBL() { + clear(); +} + +template +int ZBLT::bytes_per_atom(const int max_nbors) const { + return this->bytes_per_atom_atomic(max_nbors); +} + +template +int ZBLT::init(const int ntypes, double **host_cutsq, + double **host_sw1, double **host_sw2, + double **host_sw3, double **host_sw4, + double **host_sw5, + double **host_d1a, double **host_d2a, + double **host_d3a, double **host_d4a, + double **host_zze, double cut_globalsq, + double cut_innersq, double cut_inner, + const int nlocal, const int nall, const int max_nbors, + const int maxspecial, const double cell_size, + const double gpu_split, FILE *_screen) { + int success; + success=this->init_atomic(nlocal,nall,max_nbors,maxspecial,cell_size,gpu_split, + _screen,zbl,"k_zbl"); + if (success!=0) + return success; + + // If atom type constants fit in shared memory use fast kernel + int lj_types=ntypes; + shared_types=false; + int max_shared_types=this->device->max_shared_types(); + if (lj_types<=max_shared_types && this->_block_size>=max_shared_types) { + lj_types=max_shared_types; + shared_types=true; + } + _lj_types=lj_types; + + // Allocate a host write buffer for data initialization + UCL_H_Vec host_write(lj_types*lj_types*32,*(this->ucl_device), + UCL_WRITE_ONLY); + + for (int i=0; iucl_device),UCL_READ_ONLY); + this->atom->type_pack4(ntypes,lj_types,coeff1,host_write,host_sw1,host_sw2, + host_zze, host_cutsq); + + coeff2.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY); + this->atom->type_pack4(ntypes,lj_types,coeff2,host_write,host_d1a,host_d2a, + host_d3a,host_d4a); + + coeff3.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY); + this->atom->type_pack4(ntypes,lj_types,coeff3,host_write,host_sw3,host_sw4,host_sw5); + + _cut_globalsq = cut_globalsq; + _cut_innersq = cut_innersq; + _cut_inner = cut_inner; + printf("params: %f %f %f\n", _cut_globalsq, _cut_innersq, _cut_inner); + _allocated=true; + this->_max_bytes=coeff1.row_bytes()+coeff2.row_bytes()+coeff3.row_bytes(); + return 0; +} + +template +void ZBLT::clear() { + if (!_allocated) + return; + _allocated=false; + + coeff1.clear(); + coeff2.clear(); + coeff3.clear(); + this->clear_atomic(); +} + +template +double ZBLT::host_memory_usage() const { + return this->host_memory_usage_atomic()+sizeof(ZBL); +} + +// --------------------------------------------------------------------------- +// Calculate energies, forces, and torques +// --------------------------------------------------------------------------- +template +void ZBLT::loop(const bool _eflag, const bool _vflag) { + // Compute the block size and grid size to keep all cores busy + const int BX=this->block_size(); + int eflag, vflag; + if (_eflag) + eflag=1; + else + eflag=0; + + if (_vflag) + vflag=1; + else + vflag=0; + + int GX=static_cast(ceil(static_cast(this->ans->inum())/ + (BX/this->_threads_per_atom))); + + int ainum=this->ans->inum(); + int nbor_pitch=this->nbor->nbor_pitch(); + this->time_pair.start(); + if (shared_types) { + this->k_pair_fast.set_size(GX,BX); + this->k_pair_fast.run(&this->atom->x, &coeff1, &coeff2, &coeff3, + &_cut_globalsq, &_cut_innersq, &_cut_inner, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &this->ans->force, &this->ans->engv, &eflag, &vflag, + &ainum, &nbor_pitch, &this->_threads_per_atom); + } else { + this->k_pair.set_size(GX,BX); + this->k_pair.run(&this->atom->x, &coeff1, &coeff2, &coeff3, + &_cut_globalsq, &_cut_innersq, &_cut_inner, &_lj_types, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &this->ans->force, &this->ans->engv, &eflag, &vflag, + &ainum, &nbor_pitch, &this->_threads_per_atom); + } + this->time_pair.stop(); +} + +template class ZBL; diff --git a/lib/gpu/lal_zbl.h b/lib/gpu/lal_zbl.h new file mode 100644 index 0000000000..2996d90a5c --- /dev/null +++ b/lib/gpu/lal_zbl.h @@ -0,0 +1,84 @@ +/*************************************************************************** + zbl.h + ------------------- + Trung Dac Nguyen (ORNL) + + Class for acceleration of the zbl pair style. + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : + email : ndactrung@gmail.com + ***************************************************************************/ + +#ifndef LAL_ZBL_H +#define LAL_ZBL_H + +#include "lal_base_atomic.h" + +namespace LAMMPS_AL { + +template +class ZBL : public BaseAtomic { + public: + ZBL(); + ~ZBL(); + + /// Clear any previous data and set up for a new LAMMPS run + /** \param max_nbors initial number of rows in the neighbor matrix + * \param cell_size cutoff + skin + * \param gpu_split fraction of particles handled by device + * + * Returns: + * - 0 if successfull + * - -1 if fix gpu not found + * - -3 if there is an out of memory error + * - -4 if the GPU library was not compiled for GPU + * - -5 Double precision is not supported on card **/ + int init(const int ntypes, double **host_cutsq, double **host_sw1, + double **host_sw2, double **host_sw3, double **host_sw4, double **host_sw5, + double **host_d1a, double **host_d2a, double **host_d3a, double **host_d4a, + double **host_zze, double cut_globalsq, double cut_innersq, double cut_inner, + const int nlocal, const int nall, const int max_nbors, + const int maxspecial, const double cell_size, + const double gpu_split, FILE *screen); + + /// Clear all host and device data + /** \note This is called at the beginning of the init() routine **/ + void clear(); + + /// Returns memory usage on device per atom + int bytes_per_atom(const int max_nbors) const; + + /// Total host memory used by library for pair style + double host_memory_usage() const; + + // --------------------------- TYPE DATA -------------------------- + + /// coeff1.x = sw1, coeff1.y = sw2, coeff1.z = zze, coeff1.w = cutsq + UCL_D_Vec coeff1; + /// coeff2.x = d1a, coeff2.y = d2a, coeff2.z = d3a, coeff2.w = d4a + UCL_D_Vec coeff2; + /// coeff3.x = sw3, coeff3.y = sw4, coeff3.z = sw5; + UCL_D_Vec coeff3; + + /// If atom type constants fit in shared memory, use fast kernels + bool shared_types; + + double _cut_globalsq; + double _cut_innersq; + double _cut_inner; + + /// Number of atom types + int _lj_types; + + private: + bool _allocated; + void loop(const bool _eflag, const bool _vflag); +}; + +} + +#endif diff --git a/lib/gpu/lal_zbl_ext.cpp b/lib/gpu/lal_zbl_ext.cpp new file mode 100644 index 0000000000..ddce858076 --- /dev/null +++ b/lib/gpu/lal_zbl_ext.cpp @@ -0,0 +1,123 @@ +/*************************************************************************** + zbl_ext.cpp + ------------------- + Trung Dac Nguyen + + Functions for LAMMPS access to zbl acceleration routines. + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : + email : ndactrung@gmail.com + ***************************************************************************/ + +#include +#include +#include + +#include "lal_zbl.h" + +using namespace std; +using namespace LAMMPS_AL; + +static ZBL ZBLMF; + +// --------------------------------------------------------------------------- +// Allocate memory on host and device and copy constants to device +// --------------------------------------------------------------------------- +int zbl_gpu_init(const int ntypes, double **cutsq, double **host_sw1, + double **host_sw2, double **host_sw3, double **host_sw4, double **host_sw5, + double **host_d1a, double **host_d2a, double **host_d3a, double **host_d4a, + double **host_zze, double cut_globalsq, double cut_innersq, double cut_inner, + const int inum, const int nall, const int max_nbors, + const int maxspecial, const double cell_size, int &gpu_mode, FILE *screen) { + ZBLMF.clear(); + gpu_mode=ZBLMF.device->gpu_mode(); + double gpu_split=ZBLMF.device->particle_split(); + int first_gpu=ZBLMF.device->first_device(); + int last_gpu=ZBLMF.device->last_device(); + int world_me=ZBLMF.device->world_me(); + int gpu_rank=ZBLMF.device->gpu_rank(); + int procs_per_gpu=ZBLMF.device->procs_per_gpu(); + + ZBLMF.device->init_message(screen,"zbl",first_gpu,last_gpu); + + bool message=false; + if (ZBLMF.device->replica_me()==0 && screen) + message=true; + + if (message) { + fprintf(screen,"Initializing Device and compiling on process 0..."); + fflush(screen); + } + + int init_ok=0; + if (world_me==0) + init_ok=ZBLMF.init(ntypes, cutsq, host_sw1, host_sw2, host_sw3, host_sw4, + host_sw5, host_d1a, host_d2a, host_d3a, host_d4a, host_zze, + cut_globalsq, cut_innersq, cut_inner, + inum, nall, 300, maxspecial, cell_size, gpu_split, screen); + + ZBLMF.device->world_barrier(); + if (message) + fprintf(screen,"Done.\n"); + + for (int i=0; igpu_barrier(); + if (message) + fprintf(screen,"Done.\n"); + } + if (message) + fprintf(screen,"\n"); + + if (init_ok==0) + ZBLMF.estimate_gpu_overhead(); + return init_ok; +} + +void zbl_gpu_clear() { + ZBLMF.clear(); +} + +int ** zbl_gpu_compute_n(const int ago, const int inum_full, + const int nall, double **host_x, int *host_type, + double *sublo, double *subhi, tagint *tag, int **nspecial, + tagint **special, const bool eflag, const bool vflag, + const bool eatom, const bool vatom, int &host_start, + int **ilist, int **jnum, const double cpu_time, + bool &success) { + return ZBLMF.compute(ago, inum_full, nall, host_x, host_type, sublo, + subhi, tag, nspecial, special, eflag, vflag, eatom, + vatom, host_start, ilist, jnum, cpu_time, success); +} + +void zbl_gpu_compute(const int ago, const int inum_full, const int nall, + double **host_x, int *host_type, int *ilist, int *numj, + int **firstneigh, const bool eflag, const bool vflag, + const bool eatom, const bool vatom, int &host_start, + const double cpu_time, bool &success) { + ZBLMF.compute(ago,inum_full,nall,host_x,host_type,ilist,numj, + firstneigh,eflag,vflag,eatom,vatom,host_start,cpu_time,success); +} + +double zbl_gpu_bytes() { + return ZBLMF.host_memory_usage(); +} + +