From 008896a77d6cb6ea989ef81ef712c7ce0329c49f Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 7 Apr 2016 21:10:37 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14808 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- lib/gpu/lal_tersoff_mod.cpp | 456 +++++++++++++ lib/gpu/lal_tersoff_mod.cu | 1070 +++++++++++++++++++++++++++++++ lib/gpu/lal_tersoff_mod.h | 118 ++++ lib/gpu/lal_tersoff_mod_ext.cpp | 135 ++++ lib/gpu/lal_tersoff_mod_extra.h | 627 ++++++++++++++++++ lib/gpu/lal_tersoff_zbl.cpp | 482 ++++++++++++++ lib/gpu/lal_tersoff_zbl.cu | 1065 ++++++++++++++++++++++++++++++ lib/gpu/lal_tersoff_zbl.h | 123 ++++ lib/gpu/lal_tersoff_zbl_ext.cpp | 146 +++++ lib/gpu/lal_tersoff_zbl_extra.h | 690 ++++++++++++++++++++ 10 files changed, 4912 insertions(+) create mode 100644 lib/gpu/lal_tersoff_mod.cpp create mode 100644 lib/gpu/lal_tersoff_mod.cu create mode 100644 lib/gpu/lal_tersoff_mod.h create mode 100644 lib/gpu/lal_tersoff_mod_ext.cpp create mode 100644 lib/gpu/lal_tersoff_mod_extra.h create mode 100644 lib/gpu/lal_tersoff_zbl.cpp create mode 100644 lib/gpu/lal_tersoff_zbl.cu create mode 100644 lib/gpu/lal_tersoff_zbl.h create mode 100644 lib/gpu/lal_tersoff_zbl_ext.cpp create mode 100644 lib/gpu/lal_tersoff_zbl_extra.h diff --git a/lib/gpu/lal_tersoff_mod.cpp b/lib/gpu/lal_tersoff_mod.cpp new file mode 100644 index 0000000000..bfcc9c3bd3 --- /dev/null +++ b/lib/gpu/lal_tersoff_mod.cpp @@ -0,0 +1,456 @@ +/*************************************************************************** + tersoff_mod.cpp + ------------------- + Trung Dac Nguyen + + Class for acceleration of the tersoff pair style. + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : + email : ndactrung@gmail.com + ***************************************************************************/ + +#if defined(USE_OPENCL) +#include "tersoff_mod_cl.h" +#elif defined(USE_CUDART) +const char *tersoff_mod=0; +#else +#include "tersoff_mod_cubin.h" +#endif + +#include "lal_tersoff_mod.h" +#include +using namespace LAMMPS_AL; +#define TersoffMT TersoffMod + +extern Device device; + +template +TersoffMT::TersoffMod() : BaseThree(), _allocated(false) { +} + +template +TersoffMT::~TersoffMod() { + clear(); +} + +template +int TersoffMT::bytes_per_atom(const int max_nbors) const { + return this->bytes_per_atom_atomic(max_nbors); +} + +template +int TersoffMT::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* c5, const double* h, const double* beta, const double* powern, + const double* powern_del, const double* ca1, const double* host_cutsq) +{ + int success; + success=this->init_three(nlocal,nall,max_nbors,0,cell_size,gpu_split, + _screen,tersoff_mod,"k_tersoff_mod_repulsive", + "k_tersoff_mod_three_center", "k_tersoff_mod_three_end"); + if (success!=0) + return success; + + int ef_nall=nall; + if (ef_nall==0) + ef_nall=2000; + _zetaij.alloc(ef_nall*max_nbors,*(this->ucl_device),UCL_READ_WRITE); + + k_zeta.set_function(*(this->pair_program),"k_tersoff_mod_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(beta[i]); + dview[i].y=static_cast(powern[i]); + dview[i].z=static_cast(powern_del[i]); + dview[i].w=static_cast(ca1[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(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(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(c5[i]); + dview[i].y=static_cast(h[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 < ntypes; 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()+cutsq.row_bytes()+ + map.row_bytes()+elem2param.row_bytes()+_zetaij.row_bytes(); + return 0; +} + +template +void TersoffMT::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 TersoffMT::host_memory_usage() const { + return this->host_memory_usage_atomic()+sizeof(TersoffMod); +} + +#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 TersoffMT::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; + _max_nbors = this->nbor->max_nbor_loop(nlist,numj,ilist); + } + + 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_nbors > _zetaij.cols()) { + int _nmax=static_cast(static_cast(nall)*1.10); + _zetaij.resize(_max_nbors*_nmax); + } + + int _eflag; + if (eflag) + _eflag=1; + else + _eflag=0; + + int ainum=nall; + int nbor_pitch=this->nbor->nbor_pitch(); + int BX=this->block_pair(); + int GX=static_cast(ceil(static_cast(ainum)/ + (BX/(JTHREADS*KTHREADS)))); + + this->k_zeta.set_size(GX,BX); + this->k_zeta.run(&this->atom->x, &ts1, &ts2, &ts3, &ts4, &ts5, &cutsq, + &map, &elem2param, &_nelements, &_nparams, &_zetaij, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &_eflag, &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 ** TersoffMT::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_nbors > _zetaij.cols()) { + int _nmax=static_cast(static_cast(nall)*1.10); + _zetaij.resize(_max_nbors*_nmax); + } + + int _eflag; + if (eflag) + _eflag=1; + else + _eflag=0; + + int ainum=nall; + int nbor_pitch=this->nbor->nbor_pitch(); + int BX=this->block_pair(); + int GX=static_cast(ceil(static_cast(ainum)/ + (BX/(JTHREADS*KTHREADS)))); + + this->k_zeta.set_size(GX,BX); + this->k_zeta.run(&this->atom->x, &ts1, &ts2, &ts3, &ts4, &ts5, &cutsq, + &map, &elem2param, &_nelements, &_nparams, &_zetaij, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &_eflag, &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 TersoffMT::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, &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, &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, &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 TersoffMod; + diff --git a/lib/gpu/lal_tersoff_mod.cu b/lib/gpu/lal_tersoff_mod.cu new file mode 100644 index 0000000000..ba4ad32005 --- /dev/null +++ b/lib/gpu/lal_tersoff_mod.cu @@ -0,0 +1,1070 @@ +// ************************************************************************** +// tersoff_mod.cu +// ------------------- +// Trung Dac Nguyen +// +// Device code for acceleration of the tersoff pair style +// +// __________________________________________________________________________ +// This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) +// __________________________________________________________________________ +// +// begin : +// email : ndactrung@gmail.com +// ***************************************************************************/ + +#ifdef NV_KERNEL +#include "lal_tersoff_mod_extra.h" + +#ifndef _DOUBLE_DOUBLE +texture pos_tex; +texture ts1_tex; +texture ts2_tex; +texture ts3_tex; +texture ts4_tex; +texture ts5_tex; +#else +texture pos_tex; +texture ts1_tex; +texture ts2_tex; +texture ts3_tex; +texture ts4_tex; +texture ts5_tex; +#endif + +#else +#define pos_tex x_ +#define ts1_tex ts1 +#define ts2_tex ts2 +#define ts3_tex ts3 +#define ts4_tex ts4 +#define ts5_tex ts5 +#endif + +//#define THREE_CONCURRENT + +#define TWOTHIRD (numtyp)0.66666666666666666667 + +#define zeta_idx(nbor_mem, packed_mem, nbor_pitch, n_stride, t_per_atom, \ + i, nbor_j, offset_j, idx) \ + if (nbor_mem==packed_mem) { \ + int jj = (nbor_j-offset_j-2*nbor_pitch)/n_stride; \ + idx = jj*n_stride + i*t_per_atom + offset_j; \ + } else { \ + idx = nbor_j; \ + } + +#if (ARCH < 300) + +#define store_answers_p(f, energy, virial, ii, inum, tid, t_per_atom, \ + offset, eflag, vflag, ans, engv) \ + if (t_per_atom>1) { \ + __local acctyp red_acc[6][BLOCK_PAIR]; \ + red_acc[0][tid]=f.x; \ + red_acc[1][tid]=f.y; \ + red_acc[2][tid]=f.z; \ + red_acc[3][tid]=energy; \ + for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \ + if (offset < s) { \ + for (int r=0; r<4; r++) \ + red_acc[r][tid] += red_acc[r][tid+s]; \ + } \ + } \ + f.x=red_acc[0][tid]; \ + f.y=red_acc[1][tid]; \ + f.z=red_acc[2][tid]; \ + energy=red_acc[3][tid]; \ + if (vflag>0) { \ + for (int r=0; r<6; r++) \ + red_acc[r][tid]=virial[r]; \ + for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \ + if (offset < s) { \ + for (int r=0; r<6; r++) \ + red_acc[r][tid] += red_acc[r][tid+s]; \ + } \ + } \ + for (int r=0; r<6; r++) \ + virial[r]=red_acc[r][tid]; \ + } \ + } \ + if (offset==0) { \ + int ei=ii; \ + if (eflag>0) { \ + engv[ei]+=energy*(acctyp)0.5; \ + ei+=inum; \ + } \ + if (vflag>0) { \ + for (int i=0; i<6; i++) { \ + engv[ei]+=virial[i]*(acctyp)0.5; \ + ei+=inum; \ + } \ + } \ + acctyp4 old=ans[ii]; \ + old.x+=f.x; \ + old.y+=f.y; \ + old.z+=f.z; \ + ans[ii]=old; \ + } + +#define store_zeta(z, tid, t_per_atom, offset) \ + if (t_per_atom>1) { \ + __local acctyp red_acc[BLOCK_PAIR]; \ + red_acc[tid]=z; \ + for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \ + if (offset < s) { \ + red_acc[tid] += red_acc[tid+s]; \ + } \ + } \ + z=red_acc[tid]; \ + } + +#else + +#define store_answers_p(f, energy, virial, ii, inum, tid, t_per_atom, \ + offset, eflag, vflag, ans, engv) \ + if (t_per_atom>1) { \ + for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \ + f.x += shfl_xor(f.x, s, t_per_atom); \ + f.y += shfl_xor(f.y, s, t_per_atom); \ + f.z += shfl_xor(f.z, s, t_per_atom); \ + energy += shfl_xor(energy, s, t_per_atom); \ + } \ + if (vflag>0) { \ + for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \ + for (int r=0; r<6; r++) \ + virial[r] += shfl_xor(virial[r], s, t_per_atom); \ + } \ + } \ + } \ + if (offset==0) { \ + int ei=ii; \ + if (eflag>0) { \ + engv[ei]+=energy*(acctyp)0.5; \ + ei+=inum; \ + } \ + if (vflag>0) { \ + for (int i=0; i<6; i++) { \ + engv[ei]+=virial[i]*(acctyp)0.5; \ + ei+=inum; \ + } \ + } \ + acctyp4 old=ans[ii]; \ + old.x+=f.x; \ + old.y+=f.y; \ + old.z+=f.z; \ + ans[ii]=old; \ + } + +#define store_zeta(z, tid, t_per_atom, offset) \ + if (t_per_atom>1) { \ + for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \ + z += shfl_xor(z, s, t_per_atom); \ + } \ + } + +#endif + +// Tersoff is currently used for 3 elements at most: 3*3*3 = 27 entries +// while the block size should never be less than 32. +// SHARED_SIZE = 32 for now to reduce the pressure on the shared memory per block +// must be increased if there will be more than 3 elements in the future. + +#define SHARED_SIZE 32 + +__kernel void k_tersoff_mod_zeta(const __global numtyp4 *restrict x_, + const __global numtyp4 *restrict ts1_in, + const __global numtyp4 *restrict ts2_in, + const __global numtyp4 *restrict ts3_in, + const __global numtyp4 *restrict ts4_in, + const __global numtyp4 *restrict ts5_in, + const __global numtyp *restrict cutsq, + const __global int *restrict map, + const __global int *restrict elem2param, + const int nelements, const int nparams, + __global acctyp4 * zetaij, + const __global int * dev_nbor, + const __global int * dev_packed, + const int eflag, const int nall, const int inum, + const int nbor_pitch, const int t_per_atom) { + __local int tpa_sq,n_stride; + tpa_sq = fast_mul(t_per_atom,t_per_atom); + + int tid, ii, offset; + atom_info(tpa_sq,ii,tid,offset); + + // must be increased if there will be more than 3 elements in the future. + __local numtyp4 ts1[SHARED_SIZE]; + __local numtyp4 ts2[SHARED_SIZE]; + __local numtyp4 ts3[SHARED_SIZE]; + __local numtyp4 ts4[SHARED_SIZE]; + __local numtyp4 ts5[SHARED_SIZE]; + if (tid cutsq[ijparam]) continue; + + // compute zeta_ij + z = (numtyp)0; + + int nbor_k = nborj_start-offset_j+offset_k; + for ( ; nbor_k < nbor_end; nbor_k+=n_stride) { + int k=dev_packed[nbor_k]; + k &= NEIGHMASK; + + if (k == j) continue; + + numtyp4 kx; fetch4(kx,k,pos_tex); //x_[k]; + int ktype=kx.w; + ktype=map[ktype]; + int ijkparam=elem2param[itype*nelements*nelements+jtype*nelements+ktype]; + + // Compute rik + delr2.x = kx.x-ix.x; + delr2.y = kx.y-ix.y; + delr2.z = kx.z-ix.z; + numtyp rsq2 = delr2.x*delr2.x+delr2.y*delr2.y+delr2.z*delr2.z; + + if (rsq2 > cutsq[ijkparam]) continue; + + numtyp4 ts1_ijkparam = ts1[ijkparam]; //fetch4(ts1_ijkparam,ijkparam,ts1_tex); + numtyp ijkparam_lam3 = ts1_ijkparam.z; + numtyp ijkparam_powermint = ts1_ijkparam.w; + numtyp4 ts2_ijkparam = ts2[ijkparam]; //fetch4(ts2_ijkparam,ijkparam,ts2_tex); + numtyp ijkparam_bigr = ts2_ijkparam.z; + numtyp ijkparam_bigd = ts2_ijkparam.w; + numtyp4 ts4_ijkparam = ts4[ijkparam]; //fetch4(ts4_ijkparam,ijkparam,ts4_tex); + numtyp ijkparam_c1 = ts4_ijkparam.x; + numtyp ijkparam_c2 = ts4_ijkparam.y; + numtyp ijkparam_c3 = ts4_ijkparam.z; + numtyp ijkparam_c4 = ts4_ijkparam.w; + numtyp4 ts5_ijkparam = ts5[ijkparam]; //fetch4(ts4_ijkparam,ijkparam,ts4_tex); + numtyp ijkparam_c5 = ts5_ijkparam.x; + numtyp ijkparam_h = ts5_ijkparam.y; + z += zeta(ijkparam_powermint, ijkparam_lam3, ijkparam_bigr, ijkparam_bigd, + ijkparam_h, ijkparam_c1, ijkparam_c2, ijkparam_c3, ijkparam_c4, + ijkparam_c5, rsq1, rsq2, delr1, delr2); + } + + //int jj = (nbor_j-offset_j-2*nbor_pitch)/n_stride; + //int idx = jj*n_stride + i*t_per_atom + offset_j; + int idx; + zeta_idx(dev_nbor,dev_packed, nbor_pitch, n_stride, t_per_atom, + i, nbor_j, offset_j, idx); + store_zeta(z, tid, t_per_atom, offset_k); + + numtyp4 ts1_ijparam = ts1[ijparam]; //fetch4(ts1_ijparam,ijparam,ts1_tex); + numtyp ijparam_lam2 = ts1_ijparam.y; + numtyp4 ts2_ijparam = ts2[ijparam]; //fetch4(ts2_ijparam,ijparam,ts2_tex); + numtyp ijparam_bigb = ts2_ijparam.y; + numtyp ijparam_bigr = ts2_ijparam.z; + numtyp ijparam_bigd = ts2_ijparam.w; + numtyp4 ts3_ijparam = ts3[ijparam]; //fetch4(ts3_ijparam,ijparam,ts3_tex); + numtyp ijparam_beta = ts3_ijparam.x; + numtyp ijparam_powern = ts3_ijparam.y; + numtyp ijparam_powern_del = ts3_ijparam.z; + numtyp ijparam_ca1 = ts3_ijparam.w; + numtyp ijparam_ca4 = ucl_recip(ts3_ijparam.w); + + if (offset_k == 0) { + numtyp fpfeng[4]; + force_zeta(ijparam_bigb, ijparam_bigr, ijparam_bigd, ijparam_lam2, + ijparam_beta, ijparam_powern, ijparam_powern_del, ijparam_ca1, + ijparam_ca4, rsq1, z, eflag, fpfeng); + acctyp4 zij; + zij.x = fpfeng[0]; + zij.y = fpfeng[1]; + zij.z = fpfeng[2]; + zij.w = z; + zetaij[idx] = zij; + } + + } // for nbor + } // if ii +} + +__kernel void k_tersoff_mod_repulsive(const __global numtyp4 *restrict x_, + const __global numtyp4 *restrict ts1_in, + const __global numtyp4 *restrict ts2_in, + const __global numtyp *restrict cutsq, + const __global int *restrict map, + const __global int *restrict elem2param, + const int nelements, const int nparams, + const __global int * dev_nbor, + const __global int * dev_packed, + __global acctyp4 *restrict ans, + __global acctyp *restrict engv, + const int eflag, const int vflag, + const int inum, const int nbor_pitch, + const int t_per_atom) { + __local int n_stride; + int tid, ii, offset; + atom_info(t_per_atom,ii,tid,offset); + + __local numtyp4 ts1[SHARED_SIZE]; + __local numtyp4 ts2[SHARED_SIZE]; + if (tid0) + energy+=feng[1]; + if (vflag>0) { + virial[0] += delx*delx*force; + virial[1] += dely*dely*force; + virial[2] += delz*delz*force; + virial[3] += delx*dely*force; + virial[4] += delx*delz*force; + virial[5] += dely*delz*force; + } + } + } // for nbor + + store_answers(f,energy,virial,ii,inum,tid,t_per_atom,offset,eflag,vflag, + ans,engv); + } // if ii + +} + +__kernel void k_tersoff_mod_three_center(const __global numtyp4 *restrict x_, + const __global numtyp4 *restrict ts1_in, + const __global numtyp4 *restrict ts2_in, + const __global numtyp4 *restrict ts4_in, + const __global numtyp4 *restrict ts5_in, + const __global numtyp *restrict cutsq, + const __global int *restrict map, + const __global int *restrict elem2param, + const int nelements, const int nparams, + const __global acctyp4 *restrict zetaij, + const __global int * dev_nbor, + const __global int * dev_packed, + __global acctyp4 *restrict ans, + __global acctyp *restrict engv, + const int eflag, const int vflag, + const int inum, const int nbor_pitch, + const int t_per_atom, const int evatom) { + __local int tpa_sq, n_stride; + tpa_sq=fast_mul(t_per_atom,t_per_atom); + numtyp lam3, powermint, bigr, bigd, c1, c2, c3, c4, c5, h; + + int tid, ii, offset; + atom_info(tpa_sq,ii,tid,offset); // offset ranges from 0 to tpa_sq-1 + + __local numtyp4 ts1[SHARED_SIZE]; + __local numtyp4 ts2[SHARED_SIZE]; + __local numtyp4 ts4[SHARED_SIZE]; + __local numtyp4 ts5[SHARED_SIZE]; + if (tid cutsq[ijparam]) continue; + numtyp r1 = ucl_sqrt(rsq1); + numtyp r1inv = ucl_rsqrt(rsq1); + + // look up for zeta_ij + + //int jj = (nbor_j-offset_j-2*nbor_pitch) / n_stride; + //int idx = jj*n_stride + i*t_per_atom + offset_j; + int idx; + zeta_idx(dev_nbor,dev_packed, nbor_pitch, n_stride, t_per_atom, + i, nbor_j, offset_j, idx); + acctyp4 zeta_ij = zetaij[idx]; // fetch(zeta_ij,idx,zeta_tex); + numtyp force = zeta_ij.x*tpainv; + numtyp prefactor = zeta_ij.y; + f.x += delr1[0]*force; + f.y += delr1[1]*force; + f.z += delr1[2]*force; + + if (eflag>0) { + energy+=zeta_ij.z*tpainv; + } + if (vflag>0) { + numtyp mforce = -force; + virial[0] += delr1[0]*delr1[0]*mforce; + virial[1] += delr1[1]*delr1[1]*mforce; + virial[2] += delr1[2]*delr1[2]*mforce; + virial[3] += delr1[0]*delr1[1]*mforce; + virial[4] += delr1[0]*delr1[2]*mforce; + virial[5] += delr1[1]*delr1[2]*mforce; + } + + int nbor_k=nborj_start-offset_j+offset_k; + for ( ; nbor_k cutsq[ijkparam]) continue; + numtyp r2 = ucl_sqrt(rsq2); + numtyp r2inv = ucl_rsqrt(rsq2); + + numtyp fi[3], fj[3], fk[3]; + numtyp4 ts1_ijkparam = ts1[ijkparam]; //fetch4(ts1_ijkparam,ijkparam,ts1_tex); + lam3 = ts1_ijkparam.z; + powermint = ts1_ijkparam.w; + numtyp4 ts2_ijkparam = ts2[ijkparam]; //fetch4(ts2_ijkparam,ijkparam,ts2_tex); + bigr = ts2_ijkparam.z; + bigd = ts2_ijkparam.w; + numtyp4 ts4_ijkparam = ts4[ijkparam]; //fetch4(ts4_ijkparam,ijkparam,ts4_tex); + c1 = ts4_ijkparam.x; + c2 = ts4_ijkparam.y; + c3 = ts4_ijkparam.z; + c4 = ts4_ijkparam.w; + numtyp4 ts5_ijkparam = ts5[ijkparam]; //fetch4(ts5_ijkparam,ijkparam,ts5_tex); + c5 = ts5_ijkparam.x; + h = ts5_ijkparam.y; + if (vflag>0) + attractive(bigr, bigd, powermint, lam3, h, c1, c2, c3, c4, c5, + prefactor, r1, r1inv, r2, r2inv, delr1, delr2, fi, fj, fk); + else + attractive_fi(bigr, bigd, powermint, lam3, h, c1, c2, c3, c4, c5, + prefactor, r1, r1inv, r2, r2inv, delr1, delr2, fi); + f.x += fi[0]; + f.y += fi[1]; + f.z += fi[2]; + + if (vflag>0) { + acctyp v[6]; + numtyp pre = (numtyp)2.0; + if (evatom==1) pre = TWOTHIRD; + v[0] = pre*(delr1[0]*fj[0] + delr2[0]*fk[0]); + v[1] = pre*(delr1[1]*fj[1] + delr2[1]*fk[1]); + v[2] = pre*(delr1[2]*fj[2] + delr2[2]*fk[2]); + v[3] = pre*(delr1[0]*fj[1] + delr2[0]*fk[1]); + v[4] = pre*(delr1[0]*fj[2] + delr2[0]*fk[2]); + v[5] = pre*(delr1[1]*fj[2] + delr2[1]*fk[2]); + + virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2]; + virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5]; + } + } // nbor_k + } // for nbor_j + + store_answers_p(f,energy,virial,ii,inum,tid,tpa_sq, + offset,eflag,vflag,ans,engv); + } // if ii +} + +__kernel void k_tersoff_mod_three_end(const __global numtyp4 *restrict x_, + const __global numtyp4 *restrict ts1_in, + const __global numtyp4 *restrict ts2_in, + const __global numtyp4 *restrict ts4_in, + const __global numtyp4 *restrict ts5_in, + const __global numtyp *restrict cutsq, + const __global int *restrict map, + const __global int *restrict elem2param, + const int nelements, const int nparams, + const __global acctyp4 *restrict zetaij, + const __global int * dev_nbor, + const __global int * dev_packed, + __global acctyp4 *restrict ans, + __global acctyp *restrict engv, + const int eflag, const int vflag, + const int inum, const int nbor_pitch, + const int t_per_atom) { + __local int tpa_sq, n_stride; + tpa_sq=fast_mul(t_per_atom,t_per_atom); + numtyp lam3, powermint, bigr, bigd, c1, c2, c3, c4, c5, h; + + int tid, ii, offset; + atom_info(tpa_sq,ii,tid,offset); + + __local numtyp4 ts1[SHARED_SIZE]; + __local numtyp4 ts2[SHARED_SIZE]; + __local numtyp4 ts4[SHARED_SIZE]; + __local numtyp4 ts5[SHARED_SIZE]; + if (tid cutsq[ijparam]) continue; + + numtyp mdelr1[3]; + mdelr1[0] = -delr1[0]; + mdelr1[1] = -delr1[1]; + mdelr1[2] = -delr1[2]; + + int nbor_k=j+nbor_pitch; + int numk=dev_nbor[nbor_k]; + if (dev_nbor==dev_packed) { + nbor_k+=nbor_pitch+fast_mul(j,t_per_atom-1); + k_end=nbor_k+fast_mul(numk/t_per_atom,n_stride)+(numk & (t_per_atom-1)); + nbor_k+=offset_k; + } else { + nbor_k+=nbor_pitch; + nbor_k=dev_nbor[nbor_k]; + k_end=nbor_k+numk; + nbor_k+=offset_k; + } + int nbork_start = nbor_k; + + // look up for zeta_ji: find i in the j's neighbor list + int m = tid / t_per_atom; + int ijnum = -1; + for ( ; nbor_k= 0) { + offset_kf = offset_k; + } else { + ijnum = red_acc[2*m+0]; + offset_kf = red_acc[2*m+1]; + } + + //int iix = (ijnum - offset_kf - 2*nbor_pitch) / n_stride; + //int idx = iix*n_stride + j*t_per_atom + offset_kf; + int idx; + zeta_idx(dev_nbor,dev_packed, nbor_pitch, n_stride, t_per_atom, + j, ijnum, offset_kf, idx); + acctyp4 zeta_ji = zetaij[idx]; // fetch(zeta_ji,idx,zeta_tex); + numtyp force = zeta_ji.x*tpainv; + numtyp prefactor_ji = zeta_ji.y; + f.x += delr1[0]*force; + f.y += delr1[1]*force; + f.z += delr1[2]*force; + + if (eflag>0) { + energy+=zeta_ji.z*tpainv; + } + if (vflag>0) { + numtyp mforce = -force; + virial[0] += mdelr1[0]*mdelr1[0]*mforce; + virial[1] += mdelr1[1]*mdelr1[1]*mforce; + virial[2] += mdelr1[2]*mdelr1[2]*mforce; + virial[3] += mdelr1[0]*mdelr1[1]*mforce; + virial[4] += mdelr1[0]*mdelr1[2]*mforce; + virial[5] += mdelr1[1]*mdelr1[2]*mforce; + } + + // attractive forces + for (nbor_k = nbork_start ; nbor_k cutsq[jikparam]) continue; + numtyp r2 = ucl_sqrt(rsq2); + numtyp r2inv = ucl_rsqrt(rsq2); + numtyp4 ts1_param, ts2_param, ts4_param, ts5_param; + numtyp fi[3]; + + ts1_param = ts1[jikparam]; //fetch4(ts1_jikparam,jikparam,ts1_tex); + lam3 = ts1_param.z; + powermint = ts1_param.w; + ts2_param = ts2[jikparam]; //fetch4(ts2_jikparam,jikparam,ts2_tex); + bigr = ts2_param.z; + bigd = ts2_param.w; + ts4_param = ts4[jikparam]; //fetch4(ts4_jikparam,jikparam,ts4_tex); + c1 = ts4_param.x; + c2 = ts4_param.y; + c3 = ts4_param.z; + c4 = ts4_param.w; + ts5_param = ts5[jikparam]; //fetch4(ts5_jikparam,jikparam,ts5_tex); + c5 = ts5_param.x; + h = ts5_param.y; + attractive_fj(bigr, bigd, powermint, lam3, h, c1, c2, c3, c4, c5, + prefactor_ji, r1, r1inv, r2, r2inv, mdelr1, delr2, fi); + f.x += fi[0]; + f.y += fi[1]; + f.z += fi[2]; + + //int kk = (nbor_k - offset_k - 2*nbor_pitch) / n_stride; + //int idx = kk*n_stride + j*t_per_atom + offset_k; + int idx; + zeta_idx(dev_nbor,dev_packed, nbor_pitch, n_stride, t_per_atom, + j, nbor_k, offset_k, idx); + acctyp4 zeta_jk = zetaij[idx]; // fetch(zeta_jk,idx,zeta_tex); + numtyp prefactor_jk = zeta_jk.y; + int jkiparam=elem2param[jtype*nelements*nelements+ktype*nelements+itype]; + ts1_param = ts1[jkiparam]; //fetch4(ts1_jkiparam,jkiparam,ts1_tex); + lam3 = ts1_param.z; + powermint = ts1_param.w; + ts2_param = ts2[jkiparam]; //fetch4(ts2_jkiparam,jkiparam,ts2_tex); + bigr = ts2_param.z; + bigd = ts2_param.w; + ts4_param = ts4[jkiparam]; //fetch4(ts4_jkiparam,jkiparam,ts4_tex); + c1 = ts4_param.x; + c2 = ts4_param.y; + c3 = ts4_param.z; + c4 = ts4_param.w; + ts5_param = ts5[jkiparam]; //fetch4(ts5_ikiparam,jkiparam,ts5_tex); + c5 = ts5_param.x; + h = ts5_param.y; + attractive_fk(bigr, bigd, powermint, lam3, h, c1, c2, c3, c4, c5, + prefactor_jk, r2, r2inv, r1, r1inv, delr2, mdelr1, fi); + f.x += fi[0]; + f.y += fi[1]; + f.z += fi[2]; + } // for nbor_k + } // for nbor_j + + #ifdef THREE_CONCURRENT + store_answers(f,energy,virial,ii,inum,tid,tpa_sq,offset, + eflag,vflag,ans,engv); + #else + store_answers_p(f,energy,virial,ii,inum,tid,tpa_sq,offset, + eflag,vflag,ans,engv); + #endif + } // if ii +} + +__kernel void k_tersoff_mod_three_end_vatom(const __global numtyp4 *restrict x_, + const __global numtyp4 *restrict ts1_in, + const __global numtyp4 *restrict ts2_in, + const __global numtyp4 *restrict ts4_in, + const __global numtyp4 *restrict ts5_in, + const __global numtyp *restrict cutsq, + const __global int *restrict map, + const __global int *restrict elem2param, + const int nelements, const int nparams, + const __global acctyp4 *restrict zetaij, + const __global int * dev_nbor, + const __global int * dev_packed, + __global acctyp4 *restrict ans, + __global acctyp *restrict engv, + const int eflag, const int vflag, + const int inum, const int nbor_pitch, + const int t_per_atom) { + __local int tpa_sq, n_stride; + tpa_sq=fast_mul(t_per_atom,t_per_atom); + numtyp lam3, powermint, bigr, bigd, c1, c2, c3, c4, c5, h; + + int tid, ii, offset; + atom_info(tpa_sq,ii,tid,offset); + + __local numtyp4 ts1[SHARED_SIZE]; + __local numtyp4 ts2[SHARED_SIZE]; + __local numtyp4 ts4[SHARED_SIZE]; + __local numtyp4 ts5[SHARED_SIZE]; + if (tid cutsq[ijparam]) continue; + + numtyp mdelr1[3]; + mdelr1[0] = -delr1[0]; + mdelr1[1] = -delr1[1]; + mdelr1[2] = -delr1[2]; + + int nbor_k=j+nbor_pitch; + int numk=dev_nbor[nbor_k]; + if (dev_nbor==dev_packed) { + nbor_k+=nbor_pitch+fast_mul(j,t_per_atom-1); + k_end=nbor_k+fast_mul(numk/t_per_atom,n_stride)+(numk & (t_per_atom-1)); + nbor_k+=offset_k; + } else { + nbor_k+=nbor_pitch; + nbor_k=dev_nbor[nbor_k]; + k_end=nbor_k+numk; + nbor_k+=offset_k; + } + int nbork_start = nbor_k; + + // look up for zeta_ji + int m = tid / t_per_atom; + int ijnum = -1; + for ( ; nbor_k= 0) { + offset_kf = offset_k; + } else { + ijnum = red_acc[2*m+0]; + offset_kf = red_acc[2*m+1]; + } + + //int iix = (ijnum - offset_kf - 2*nbor_pitch) / n_stride; + //int idx = iix*n_stride + j*t_per_atom + offset_kf; + int idx; + zeta_idx(dev_nbor,dev_packed, nbor_pitch, n_stride, t_per_atom, + j, ijnum, offset_kf, idx); + acctyp4 zeta_ji = zetaij[idx]; // fetch(zeta_ji,idx,zeta_tex); + numtyp force = zeta_ji.x*tpainv; + numtyp prefactor_ji = zeta_ji.y; + f.x += delr1[0]*force; + f.y += delr1[1]*force; + f.z += delr1[2]*force; + + if (eflag>0) { + energy+=zeta_ji.z*tpainv; + } + if (vflag>0) { + numtyp mforce = -force; + virial[0] += mdelr1[0]*mdelr1[0]*mforce; + virial[1] += mdelr1[1]*mdelr1[1]*mforce; + virial[2] += mdelr1[2]*mdelr1[2]*mforce; + virial[3] += mdelr1[0]*mdelr1[1]*mforce; + virial[4] += mdelr1[0]*mdelr1[2]*mforce; + virial[5] += mdelr1[1]*mdelr1[2]*mforce; + } + + // attractive forces + for (nbor_k = nbork_start; nbor_k cutsq[jikparam]) continue; + numtyp r2 = ucl_sqrt(rsq2); + numtyp r2inv = ucl_rsqrt(rsq2); + + numtyp fi[3], fj[3], fk[3]; + numtyp4 ts1_param, ts2_param, ts4_param, ts5_param; + ts1_param = ts1[jikparam]; //fetch4(ts1_jikparam,jikparam,ts1_tex); + lam3 = ts1_param.z; + powermint = ts1_param.w; + ts2_param = ts2[jikparam]; //fetch4(ts2_jikparam,jikparam,ts2_tex); + bigr = ts2_param.z; + bigd = ts2_param.w; + ts4_param = ts4[jikparam]; //fetch4(ts4_jikparam,jikparam,ts4_tex); + c1 = ts4_param.x; + c2 = ts4_param.y; + c3 = ts4_param.z; + c4 = ts4_param.w; + ts5_param = ts5[jikparam]; //fetch4(ts5_jijparam,jikparam,ts5_tex); + c5 = ts5_param.x; + h = ts5_param.y; + attractive(bigr, bigd, powermint, lam3, h, c1, c2, c3, c4, c5, + prefactor_ji, r1, r1inv, r2, r2inv, mdelr1, delr2, fi, fj, fk); + f.x += fj[0]; + f.y += fj[1]; + f.z += fj[2]; + + virial[0] += TWOTHIRD*(mdelr1[0]*fj[0] + delr2[0]*fk[0]); + virial[1] += TWOTHIRD*(mdelr1[1]*fj[1] + delr2[1]*fk[1]); + virial[2] += TWOTHIRD*(mdelr1[2]*fj[2] + delr2[2]*fk[2]); + virial[3] += TWOTHIRD*(mdelr1[0]*fj[1] + delr2[0]*fk[1]); + virial[4] += TWOTHIRD*(mdelr1[0]*fj[2] + delr2[0]*fk[2]); + virial[5] += TWOTHIRD*(mdelr1[1]*fj[2] + delr2[1]*fk[2]); + + //int kk = (nbor_k - offset_k - 2*nbor_pitch) / n_stride; + //int idx = kk*n_stride + j*t_per_atom + offset_k; + int idx; + zeta_idx(dev_nbor,dev_packed, nbor_pitch, n_stride, t_per_atom, + j, nbor_k, offset_k, idx); + acctyp4 zeta_jk = zetaij[idx]; // fetch(zeta_jk,idx,zeta_tex); + numtyp prefactor_jk = zeta_jk.y; + + int jkiparam=elem2param[jtype*nelements*nelements+ktype*nelements+itype]; + ts1_param = ts1[jkiparam]; //fetch4(ts1_jkiparam,jkiparam,ts1_tex); + lam3 = ts1_param.z; + powermint = ts1_param.w; + ts2_param = ts2[jkiparam]; //fetch4(ts2_jkiparam,jkiparam,ts2_tex); + bigr = ts2_param.z; + bigd = ts2_param.w; + ts4_param = ts4[jkiparam]; //fetch4(ts4_jkiparam,jkiparam,ts4_tex); + c1 = ts4_param.x; + c2 = ts4_param.y; + c3 = ts4_param.z; + c4 = ts4_param.w; + ts5_param = ts5[jkiparam]; //fetch4(ts5_ikiparam,jkiparam,ts5_tex); + c5 = ts5_param.x; + h = ts5_param.y; + attractive(bigr, bigd, powermint, lam3, h, c1, c2, c3, c4, c5, + prefactor_jk, r2, r2inv, r1, r1inv, delr2, mdelr1, fi, fj, fk); + f.x += fk[0]; + f.y += fk[1]; + f.z += fk[2]; + + virial[0] += TWOTHIRD*(delr2[0]*fj[0] + mdelr1[0]*fk[0]); + virial[1] += TWOTHIRD*(delr2[1]*fj[1] + mdelr1[1]*fk[1]); + virial[2] += TWOTHIRD*(delr2[2]*fj[2] + mdelr1[2]*fk[2]); + virial[3] += TWOTHIRD*(delr2[0]*fj[1] + mdelr1[0]*fk[1]); + virial[4] += TWOTHIRD*(delr2[0]*fj[2] + mdelr1[0]*fk[2]); + virial[5] += TWOTHIRD*(delr2[1]*fj[2] + mdelr1[1]*fk[2]); + } + } // for nbor + + #ifdef THREE_CONCURRENT + store_answers(f,energy,virial,ii,inum,tid,tpa_sq,offset, + eflag,vflag,ans,engv); + #else + store_answers_p(f,energy,virial,ii,inum,tid,tpa_sq,offset, + eflag,vflag,ans,engv); + #endif + } // if ii +} + diff --git a/lib/gpu/lal_tersoff_mod.h b/lib/gpu/lal_tersoff_mod.h new file mode 100644 index 0000000000..9a05c66009 --- /dev/null +++ b/lib/gpu/lal_tersoff_mod.h @@ -0,0 +1,118 @@ +/*************************************************************************** + tersoff_mod.h + ------------------- + Trung Dac Nguyen + + Class for acceleration of the tersoff/mod pair style. + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : + email : ndactrung@gmail.com + ***************************************************************************/ + +#ifndef LAL_TERSOFF_MOD_H +#define LAL_TERSOFF_MOD_H + +#include "lal_base_three.h" + +namespace LAMMPS_AL { + +template +class TersoffMod : public BaseThree { + public: + TersoffMod(); + ~TersoffMod(); + + /// 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* c5, + const double* h, const double* beta, const double* powern, + const double* powern_del, const double* ca1, 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 = beta, ts3.y = powern, ts3.z = powern_del, ts3.w = ca1 + UCL_D_Vec ts3; + /// ts4.x = c1, ts4.y = c2, ts4.z = c3, ts4.w = c4 + UCL_D_Vec ts4; + /// ts5.x = c5, ts5.y = h + UCL_D_Vec ts5; + + UCL_D_Vec cutsq; + + UCL_D_Vec elem2param; + UCL_D_Vec map; + int _nparams,_nelements; + + /// Per-atom arrays: + /// zetaij.x = force, zetaij.y = prefactor, zetaij.z = evdwl, + /// zetaij.w = zetaij + UCL_D_Vec _zetaij; + + UCL_Kernel k_zeta; + UCL_Texture ts1_tex, ts2_tex, ts3_tex, ts4_tex, ts5_tex; + + int _max_nbors; + + private: + bool _allocated; + void loop(const bool _eflag, const bool _vflag, const int evatom); +}; + +} + +#endif + diff --git a/lib/gpu/lal_tersoff_mod_ext.cpp b/lib/gpu/lal_tersoff_mod_ext.cpp new file mode 100644 index 0000000000..7817e7d08d --- /dev/null +++ b/lib/gpu/lal_tersoff_mod_ext.cpp @@ -0,0 +1,135 @@ +/*************************************************************************** + tersoff_mod_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_mod.h" + +using namespace std; +using namespace LAMMPS_AL; + +static TersoffMod TSMMF; + +// --------------------------------------------------------------------------- +// Allocate memory on host and device and copy constants to device +// --------------------------------------------------------------------------- +int tersoff_mod_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_c5, const double* ts_h, const double* ts_beta, + const double* ts_powern, const double* ts_powern_del, + const double* ts_ca1, const double* ts_cutsq) { + TSMMF.clear(); + gpu_mode=TSMMF.device->gpu_mode(); + double gpu_split=TSMMF.device->particle_split(); + int first_gpu=TSMMF.device->first_device(); + int last_gpu=TSMMF.device->last_device(); + int world_me=TSMMF.device->world_me(); + int gpu_rank=TSMMF.device->gpu_rank(); + int procs_per_gpu=TSMMF.device->procs_per_gpu(); + + // disable host/device split for now + if (gpu_split != 1.0) + return -8; + + TSMMF.device->init_message(screen,"tersoff/mod/gpu",first_gpu,last_gpu); + + bool message=false; + if (TSMMF.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=TSMMF.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_c5, ts_h, ts_beta, ts_powern, + ts_powern_del, ts_ca1, ts_cutsq); + + TSMMF.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) + TSMMF.estimate_gpu_overhead(); + return init_ok; +} + +void tersoff_mod_gpu_clear() { + TSMMF.clear(); +} + +int ** tersoff_mod_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 TSMMF.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_mod_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) { + TSMMF.compute(ago,nlocal,nall,nlist,host_x,host_type,ilist,numj, + firstneigh,eflag,vflag,eatom,vatom,host_start,cpu_time,success); +} + +double tersoff_mod_gpu_bytes() { + return TSMMF.host_memory_usage(); +} + + diff --git a/lib/gpu/lal_tersoff_mod_extra.h b/lib/gpu/lal_tersoff_mod_extra.h new file mode 100644 index 0000000000..370aceb634 --- /dev/null +++ b/lib/gpu/lal_tersoff_mod_extra.h @@ -0,0 +1,627 @@ +/// ************************************************************************** +// tersoff_mod_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_MOD_EXTRA_H +#define LAL_TERSOFF_MOD_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_mod(const numtyp costheta, + const numtyp param_c1, + const numtyp param_c2, + const numtyp param_c3, + const numtyp param_c4, + const numtyp param_c5, + const numtyp param_h) +{ + const numtyp tmp_h = (param_h - costheta)*(param_h - costheta); + return param_c1 + (param_c2*tmp_h/(param_c3 + tmp_h)) * + ((numtyp)1.0 + param_c4*ucl_exp(-param_c5*tmp_h)); +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline numtyp ters_gijk_d_mod(const numtyp costheta, + const numtyp param_c2, + const numtyp param_c3, + const numtyp param_c4, + const numtyp param_c5, + const numtyp param_h) +{ + const numtyp tmp_h = (param_h - costheta)*(param_h - costheta); + const numtyp g1 = (param_h - costheta)/(param_c3 + tmp_h); + const numtyp g2 = ucl_exp(-param_c5*tmp_h); + return (numtyp)-2.0*param_c2*g1*((1 + param_c4*g2) * + (1 + g1*(costheta - param_h)) - tmp_h*param_c4*param_c5*g2); +} + +/* ---------------------------------------------------------------------- */ + +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 - + (numtyp)1.125*sin(MY_PI2*(r - param_bigr)/param_bigd) - + (numtyp)0.125*sin(3*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 -((numtyp)0.375*MY_PI4/param_bigd) * + ((numtyp)3*cos(MY_PI2*(r - param_bigr)/param_bigd) + + cos((numtyp)3*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_powern_del, + const numtyp param_ca1, + const numtyp param_ca4) +{ + numtyp tmp = param_beta * zeta; + if (tmp > param_ca1) + return ucl_powr(tmp, -param_powern/((numtyp)2.0*param_powern_del)); + if (tmp < param_ca4) return (numtyp)1.0; + return ucl_powr((numtyp)1.0 + ucl_powr(tmp,param_powern), + (numtyp)-1.0/((numtyp)2.0*param_powern_del)); +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline numtyp ters_bij_d(const numtyp zeta, + const numtyp param_beta, + const numtyp param_powern, + const numtyp param_powern_del, + const numtyp param_ca1, + const numtyp param_ca4) +{ + numtyp tmp = param_beta * zeta; + if (tmp > param_ca1) return (numtyp)-0.5*(param_powern/param_powern_del) * + ucl_powr(tmp,(numtyp)-0.5*(param_powern/param_powern_del)) / zeta; + if (tmp < param_ca4) return (numtyp)0.0; + + numtyp tmp_n = ucl_powr(tmp,param_powern); + return (numtyp)-0.5 *(param_powern/param_powern_del) * + ucl_powr((numtyp)1.0+tmp_n, (numtyp)-1.0-((numtyp)1.0 / + ((numtyp)2.0*param_powern_del)))*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_h, + const numtyp param_c1, + const numtyp param_c2, + const numtyp param_c3, + const numtyp param_c4, + const numtyp param_c5, + 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 = (numtyp)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_mod(cos_theta,param_c1,param_c2,param_c3,param_c4,param_c5,param_h); + gijk_d = ters_gijk_d_mod(cos_theta,param_c2,param_c3,param_c4,param_c5,param_h); + 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_h, + const numtyp param_c1, + const numtyp param_c2, + const numtyp param_c3, + const numtyp param_c4, + const numtyp param_c5, + 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 = (numtyp)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_mod(cos_theta,param_c1,param_c2,param_c3,param_c4,param_c5,param_h); + gijk_d = ters_gijk_d_mod(cos_theta,param_c2,param_c3,param_c4,param_c5,param_h); + 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_h, + const numtyp param_c1, + const numtyp param_c2, + const numtyp param_c3, + const numtyp param_c4, + const numtyp param_c5, + 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 = (numtyp)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_mod(cos_theta,param_c1,param_c2,param_c3,param_c4,param_c5,param_h); + gijk_d = ters_gijk_d_mod(cos_theta,param_c2,param_c3,param_c4,param_c5,param_h); + 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_h, + const numtyp param_c1, + const numtyp param_c2, + const numtyp param_c3, + const numtyp param_c4, + const numtyp param_c5, + 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 = (numtyp)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_mod(cos_theta,param_c1,param_c2,param_c3,param_c4,param_c5,param_h); + gijk_d = ters_gijk_d_mod(cos_theta,param_c2,param_c3,param_c4,param_c5,param_h); + 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_h, + const numtyp param_c1, + const numtyp param_c2, + const numtyp param_c3, + const numtyp param_c4, + const numtyp param_c5, + 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_mod(costheta,param_c1,param_c2,param_c3,param_c4,param_c5, + param_h) * 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_powern_del, + const numtyp param_ca1, + const numtyp param_ca4, + 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_powern_del,param_ca1,param_ca4); + 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_powern_del,param_ca1,param_ca4); // 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_h, + const numtyp param_c1, + const numtyp param_c2, + const numtyp param_c3, + const numtyp param_c4, + const numtyp param_c5, + 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_h, param_c1, param_c2, param_c3, param_c4, param_c5, + 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_h, + const numtyp param_c1, + const numtyp param_c2, + const numtyp param_c3, + const numtyp param_c4, + const numtyp param_c5, + 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_h, param_c1, param_c2, param_c3, param_c4, param_c5, + 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_h, + const numtyp param_c1, + const numtyp param_c2, + const numtyp param_c3, + const numtyp param_c4, + const numtyp param_c5, + 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_h, param_c1, param_c2, param_c3, param_c4, param_c5, + 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_h, + const numtyp param_c1, + const numtyp param_c2, + const numtyp param_c3, + const numtyp param_c4, + const numtyp param_c5, + 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_h, param_c1, param_c2, param_c3, param_c4, param_c5, + fk); +} + + +#endif + + diff --git a/lib/gpu/lal_tersoff_zbl.cpp b/lib/gpu/lal_tersoff_zbl.cpp new file mode 100644 index 0000000000..57688f53ab --- /dev/null +++ b/lib/gpu/lal_tersoff_zbl.cpp @@ -0,0 +1,482 @@ +/*************************************************************************** + tersoff_zbl.cpp + ------------------- + Trung Dac Nguyen + + Class for acceleration of the tersoff/zbl pair style. + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : + email : ndactrung@gmail.com + ***************************************************************************/ + +#if defined(USE_OPENCL) +#include "tersoff_zbl_cl.h" +#elif defined(USE_CUDART) +const char *tersoff_zbl=0; +#else +#include "tersoff_zbl_cubin.h" +#endif + +#include "lal_tersoff_zbl.h" +#include +using namespace LAMMPS_AL; +#define TersoffZT TersoffZBL + +extern Device device; + +template +TersoffZT::TersoffZBL() : BaseThree(), _allocated(false) { +} + +template +TersoffZT::~TersoffZBL() { + clear(); +} + +template +int TersoffZT::bytes_per_atom(const int max_nbors) const { + return this->bytes_per_atom_atomic(max_nbors); +} + +template +int TersoffZT::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* Z_i, + const double* Z_j, const double* ZBLcut, + const double* ZBLexpscale, const double global_e, + const double global_a_0, const double global_epsilon_0, + const double* host_cutsq) +{ + int success; + success=this->init_three(nlocal,nall,max_nbors,0,cell_size,gpu_split, + _screen,tersoff_zbl,"k_tersoff_zbl_repulsive", + "k_tersoff_zbl_three_center", "k_tersoff_zbl_three_end"); + if (success!=0) + return success; + + int ef_nall=nall; + if (ef_nall==0) + ef_nall=2000; + _zetaij.alloc(ef_nall*max_nbors,*(this->ucl_device),UCL_READ_WRITE); + + k_zeta.set_function(*(this->pair_program),"k_tersoff_zbl_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); + + ts6.alloc(nparams,*(this->ucl_device),UCL_READ_ONLY); + + for (int i=0; i(Z_i[i]); + dview[i].y=static_cast(Z_j[i]); + dview[i].z=static_cast(ZBLcut[i]); + dview[i].w=static_cast(ZBLexpscale[i]); + } + + ucl_copy(ts6,dview,false); + ts6_tex.get_texture(*(this->pair_program),"ts6_tex"); + ts6_tex.bind_float(ts6,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 < ntypes; i++) + dview_map[i] = host_map[i]; + + map.alloc(lj_types,*(this->ucl_device), UCL_READ_ONLY); + ucl_copy(map,dview_map,false); + + _global_e = global_e; + _global_a_0 = global_a_0; + _global_epsilon_0 = global_epsilon_0; + + _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 TersoffZT::clear() { + if (!_allocated) + return; + _allocated=false; + + ts1.clear(); + ts2.clear(); + ts3.clear(); + ts4.clear(); + ts5.clear(); + ts6.clear(); + cutsq.clear(); + map.clear(); + elem2param.clear(); + _zetaij.clear(); + + k_zeta.clear(); + + this->clear_atomic(); +} + +template +double TersoffZT::host_memory_usage() const { + return this->host_memory_usage_atomic()+sizeof(TersoffZBL); +} + +#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 TersoffZT::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; + _max_nbors = this->nbor->max_nbor_loop(nlist,numj,ilist); + } + + 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_nbors > _zetaij.cols()) { + int _nmax=static_cast(static_cast(nall)*1.10); + _zetaij.resize(_max_nbors*_nmax); + } + + int _eflag; + if (eflag) + _eflag=1; + else + _eflag=0; + + int ainum=nall; + int nbor_pitch=this->nbor->nbor_pitch(); + int BX=this->block_pair(); + int GX=static_cast(ceil(static_cast(ainum)/ + (BX/(JTHREADS*KTHREADS)))); + + this->k_zeta.set_size(GX,BX); + this->k_zeta.run(&this->atom->x, &ts1, &ts2, &ts3, &ts4, &ts5, &ts6, &cutsq, + &map, &elem2param, &_nelements, &_nparams, &_zetaij, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &_eflag, &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 ** TersoffZT::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_nbors > _zetaij.cols()) { + int _nmax=static_cast(static_cast(nall)*1.10); + _zetaij.resize(_max_nbors*_nmax); + } + + int _eflag; + if (eflag) + _eflag=1; + else + _eflag=0; + + int ainum=nall; + int nbor_pitch=this->nbor->nbor_pitch(); + int BX=this->block_pair(); + int GX=static_cast(ceil(static_cast(ainum)/ + (BX/(JTHREADS*KTHREADS)))); + + this->k_zeta.set_size(GX,BX); + this->k_zeta.run(&this->atom->x, &ts1, &ts2, &ts3, &ts4, &ts5, &ts6, &cutsq, + &map, &elem2param, &_nelements, &_nparams, &_zetaij, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &_eflag, &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 TersoffZT::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, &ts6, + &_global_e, &_global_a_0, &_global_epsilon_0, &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, &ts4, &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, &ts4, &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, &ts4, &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 TersoffZBL; + diff --git a/lib/gpu/lal_tersoff_zbl.cu b/lib/gpu/lal_tersoff_zbl.cu new file mode 100644 index 0000000000..0d6c5a38d6 --- /dev/null +++ b/lib/gpu/lal_tersoff_zbl.cu @@ -0,0 +1,1065 @@ +// ************************************************************************** +// tersoff_zbl.cu +// ------------------- +// Trung Dac Nguyen +// +// Device code for acceleration of the tersoff/zbl pair style +// +// __________________________________________________________________________ +// This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) +// __________________________________________________________________________ +// +// begin : +// email : ndactrung@gmail.com +// ***************************************************************************/ + +#ifdef NV_KERNEL +#include "lal_tersoff_zbl_extra.h" + +#ifndef _DOUBLE_DOUBLE +texture pos_tex; +texture ts1_tex; +texture ts2_tex; +texture ts3_tex; +texture ts4_tex; +texture ts5_tex; +texture ts6_tex; +#else +texture pos_tex; +texture ts1_tex; +texture ts2_tex; +texture ts3_tex; +texture ts4_tex; +texture ts5_tex; +texture ts6_tex; +#endif + +#else +#define pos_tex x_ +#define ts1_tex ts1 +#define ts2_tex ts2 +#define ts3_tex ts3 +#define ts4_tex ts4 +#define ts5_tex ts5 +#define ts6_tex ts6 +#endif + +//#define THREE_CONCURRENT + +#define TWOTHIRD (numtyp)0.66666666666666666667 + +#define zeta_idx(nbor_mem, packed_mem, nbor_pitch, n_stride, t_per_atom, \ + i, nbor_j, offset_j, idx) \ + if (nbor_mem==packed_mem) { \ + int jj = (nbor_j-offset_j-2*nbor_pitch)/n_stride; \ + idx = jj*n_stride + i*t_per_atom + offset_j; \ + } else { \ + idx = nbor_j; \ + } + +#if (ARCH < 300) + +#define store_answers_p(f, energy, virial, ii, inum, tid, t_per_atom, \ + offset, eflag, vflag, ans, engv) \ + if (t_per_atom>1) { \ + __local acctyp red_acc[6][BLOCK_PAIR]; \ + red_acc[0][tid]=f.x; \ + red_acc[1][tid]=f.y; \ + red_acc[2][tid]=f.z; \ + red_acc[3][tid]=energy; \ + for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \ + if (offset < s) { \ + for (int r=0; r<4; r++) \ + red_acc[r][tid] += red_acc[r][tid+s]; \ + } \ + } \ + f.x=red_acc[0][tid]; \ + f.y=red_acc[1][tid]; \ + f.z=red_acc[2][tid]; \ + energy=red_acc[3][tid]; \ + if (vflag>0) { \ + for (int r=0; r<6; r++) \ + red_acc[r][tid]=virial[r]; \ + for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \ + if (offset < s) { \ + for (int r=0; r<6; r++) \ + red_acc[r][tid] += red_acc[r][tid+s]; \ + } \ + } \ + for (int r=0; r<6; r++) \ + virial[r]=red_acc[r][tid]; \ + } \ + } \ + if (offset==0) { \ + int ei=ii; \ + if (eflag>0) { \ + engv[ei]+=energy*(acctyp)0.5; \ + ei+=inum; \ + } \ + if (vflag>0) { \ + for (int i=0; i<6; i++) { \ + engv[ei]+=virial[i]*(acctyp)0.5; \ + ei+=inum; \ + } \ + } \ + acctyp4 old=ans[ii]; \ + old.x+=f.x; \ + old.y+=f.y; \ + old.z+=f.z; \ + ans[ii]=old; \ + } + +#define store_zeta(z, tid, t_per_atom, offset) \ + if (t_per_atom>1) { \ + __local acctyp red_acc[BLOCK_PAIR]; \ + red_acc[tid]=z; \ + for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \ + if (offset < s) { \ + red_acc[tid] += red_acc[tid+s]; \ + } \ + } \ + z=red_acc[tid]; \ + } + +#else + +#define store_answers_p(f, energy, virial, ii, inum, tid, t_per_atom, \ + offset, eflag, vflag, ans, engv) \ + if (t_per_atom>1) { \ + for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \ + f.x += shfl_xor(f.x, s, t_per_atom); \ + f.y += shfl_xor(f.y, s, t_per_atom); \ + f.z += shfl_xor(f.z, s, t_per_atom); \ + energy += shfl_xor(energy, s, t_per_atom); \ + } \ + if (vflag>0) { \ + for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \ + for (int r=0; r<6; r++) \ + virial[r] += shfl_xor(virial[r], s, t_per_atom); \ + } \ + } \ + } \ + if (offset==0) { \ + int ei=ii; \ + if (eflag>0) { \ + engv[ei]+=energy*(acctyp)0.5; \ + ei+=inum; \ + } \ + if (vflag>0) { \ + for (int i=0; i<6; i++) { \ + engv[ei]+=virial[i]*(acctyp)0.5; \ + ei+=inum; \ + } \ + } \ + acctyp4 old=ans[ii]; \ + old.x+=f.x; \ + old.y+=f.y; \ + old.z+=f.z; \ + ans[ii]=old; \ + } + +#define store_zeta(z, tid, t_per_atom, offset) \ + if (t_per_atom>1) { \ + for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \ + z += shfl_xor(z, s, t_per_atom); \ + } \ + } + +#endif + +// Tersoff is currently used for 3 elements at most: 3*3*3 = 27 entries +// while the block size should never be less than 32. +// SHARED_SIZE = 32 for now to reduce the pressure on the shared memory per block +// must be increased if there will be more than 3 elements in the future. + +#define SHARED_SIZE 32 + +__kernel void k_tersoff_zbl_zeta(const __global numtyp4 *restrict x_, + const __global numtyp4 *restrict ts1_in, + const __global numtyp4 *restrict ts2_in, + const __global numtyp4 *restrict ts3_in, + const __global numtyp4 *restrict ts4_in, + const __global numtyp4 *restrict ts5_in, + const __global numtyp4 *restrict ts6_in, + const __global numtyp *restrict cutsq, + const __global int *restrict map, + const __global int *restrict elem2param, + const int nelements, const int nparams, + __global acctyp4 * zetaij, + const __global int * dev_nbor, + const __global int * dev_packed, + const int eflag, const int nall, const int inum, + const int nbor_pitch, const int t_per_atom) { + __local int tpa_sq,n_stride; + tpa_sq = fast_mul(t_per_atom,t_per_atom); + + int tid, ii, offset; + atom_info(tpa_sq,ii,tid,offset); + + // must be increased if there will be more than 3 elements in the future. + __local numtyp4 ts1[SHARED_SIZE]; + __local numtyp4 ts2[SHARED_SIZE]; + __local numtyp4 ts3[SHARED_SIZE]; + __local numtyp4 ts4[SHARED_SIZE]; + __local numtyp4 ts5[SHARED_SIZE]; + __local numtyp4 ts6[SHARED_SIZE]; + if (tid cutsq[ijparam]) continue; + + // compute zeta_ij + z = (acctyp)0; + + int nbor_k = nborj_start-offset_j+offset_k; + for ( ; nbor_k < nbor_end; nbor_k+=n_stride) { + int k=dev_packed[nbor_k]; + k &= NEIGHMASK; + + if (k == j) continue; + + numtyp4 kx; fetch4(kx,k,pos_tex); //x_[k]; + int ktype=kx.w; + ktype=map[ktype]; + int ijkparam=elem2param[itype*nelements*nelements+jtype*nelements+ktype]; + + // Compute rik + delr2.x = kx.x-ix.x; + delr2.y = kx.y-ix.y; + delr2.z = kx.z-ix.z; + numtyp rsq2 = delr2.x*delr2.x+delr2.y*delr2.y+delr2.z*delr2.z; + + if (rsq2 > cutsq[ijkparam]) continue; + + numtyp4 ts1_ijkparam = ts1[ijkparam]; //fetch4(ts1_ijkparam,ijkparam,ts1_tex); + numtyp ijkparam_lam3 = ts1_ijkparam.z; + numtyp ijkparam_powermint = ts1_ijkparam.w; + numtyp4 ts2_ijkparam = ts2[ijkparam]; //fetch4(ts2_ijkparam,ijkparam,ts2_tex); + numtyp ijkparam_bigr = ts2_ijkparam.z; + numtyp ijkparam_bigd = ts2_ijkparam.w; + numtyp4 ts4_ijkparam = ts4[ijkparam]; //fetch4(ts4_ijkparam,ijkparam,ts4_tex); + numtyp ijkparam_c = ts4_ijkparam.x; + numtyp ijkparam_d = ts4_ijkparam.y; + numtyp ijkparam_h = ts4_ijkparam.z; + numtyp ijkparam_gamma = ts4_ijkparam.w; + z += zeta(ijkparam_powermint, ijkparam_lam3, ijkparam_bigr, ijkparam_bigd, + ijkparam_c, ijkparam_d, ijkparam_h, ijkparam_gamma, + rsq1, rsq2, delr1, delr2); + } + + //int jj = (nbor_j-offset_j-2*nbor_pitch)/n_stride; + //int idx = jj*n_stride + i*t_per_atom + offset_j; + int idx; + zeta_idx(dev_nbor,dev_packed, nbor_pitch, n_stride, t_per_atom, + i, nbor_j, offset_j, idx); + store_zeta(z, tid, t_per_atom, offset_k); + + numtyp4 ts1_ijparam = ts1[ijparam]; //fetch4(ts1_ijparam,ijparam,ts1_tex); + numtyp ijparam_lam2 = ts1_ijparam.y; + numtyp4 ts2_ijparam = ts2[ijparam]; //fetch4(ts2_ijparam,ijparam,ts2_tex); + numtyp ijparam_bigb = ts2_ijparam.y; + numtyp ijparam_bigr = ts2_ijparam.z; + numtyp ijparam_bigd = ts2_ijparam.w; + numtyp4 ts3_ijparam = ts3[ijparam]; //fetch4(ts3_ijparam,ijparam,ts3_tex); + numtyp ijparam_c1 = ts3_ijparam.x; + numtyp ijparam_c2 = ts3_ijparam.y; + numtyp ijparam_c3 = ts3_ijparam.z; + numtyp ijparam_c4 = ts3_ijparam.w; + numtyp4 ts5_ijparam = ts5[ijparam]; //fetch4(ts5_ijparam,ijparam,ts5_tex); + numtyp ijparam_beta = ts5_ijparam.x; + numtyp ijparam_powern = ts5_ijparam.y; + numtyp4 ts6_ijparam = ts6[ijparam]; //fetch4(ts6_ijparam,ijparam,ts6_tex); + numtyp ijparam_ZBLcut = ts6_ijparam.z; + numtyp ijparam_ZBLexpscale = ts6_ijparam.w; + + if (offset_k == 0) { + numtyp fpfeng[4]; + force_zeta(ijparam_bigb, ijparam_bigr, ijparam_bigd, ijparam_lam2, + ijparam_beta, ijparam_powern, ijparam_c1, ijparam_c2, ijparam_c3, + ijparam_c4, ijparam_ZBLcut, ijparam_ZBLexpscale, rsq1, z, eflag, fpfeng); + acctyp4 zij; + zij.x = fpfeng[0]; + zij.y = fpfeng[1]; + zij.z = fpfeng[2]; + zij.w = z; + zetaij[idx] = zij; + } + + } // for nbor + } // if ii +} + +__kernel void k_tersoff_zbl_repulsive(const __global numtyp4 *restrict x_, + const __global numtyp4 *restrict ts1_in, + const __global numtyp4 *restrict ts2_in, + const __global numtyp4 *restrict ts6_in, + const numtyp global_e, const numtyp global_a_0, + const numtyp global_epsilon_0, + const __global numtyp *restrict cutsq, + const __global int *restrict map, + const __global int *restrict elem2param, + const int nelements, const int nparams, + const __global int * dev_nbor, + const __global int * dev_packed, + __global acctyp4 *restrict ans, + __global acctyp *restrict engv, + const int eflag, const int vflag, + const int inum, const int nbor_pitch, + const int t_per_atom) { + __local int n_stride; + int tid, ii, offset; + atom_info(t_per_atom,ii,tid,offset); + + __local numtyp4 ts1[SHARED_SIZE]; + __local numtyp4 ts2[SHARED_SIZE]; + __local numtyp4 ts6[SHARED_SIZE]; + if (tid0) + energy+=feng[1]; + if (vflag>0) { + virial[0] += delx*delx*force; + virial[1] += dely*dely*force; + virial[2] += delz*delz*force; + virial[3] += delx*dely*force; + virial[4] += delx*delz*force; + virial[5] += dely*delz*force; + } + } + } // for nbor + + store_answers(f,energy,virial,ii,inum,tid,t_per_atom,offset,eflag,vflag, + ans,engv); + } // if ii + +} + +__kernel void k_tersoff_zbl_three_center(const __global numtyp4 *restrict x_, + const __global numtyp4 *restrict ts1_in, + const __global numtyp4 *restrict ts2_in, + const __global numtyp4 *restrict ts4_in, + const __global numtyp *restrict cutsq, + const __global int *restrict map, + const __global int *restrict elem2param, + const int nelements, const int nparams, + const __global acctyp4 *restrict zetaij, + const __global int * dev_nbor, + const __global int * dev_packed, + __global acctyp4 *restrict ans, + __global acctyp *restrict engv, + const int eflag, const int vflag, + const int inum, const int nbor_pitch, + const int t_per_atom, const int evatom) { + __local int tpa_sq, n_stride; + tpa_sq=fast_mul(t_per_atom,t_per_atom); + numtyp lam3, powermint, bigr, bigd, c, d, h, gamma; + + int tid, ii, offset; + atom_info(tpa_sq,ii,tid,offset); // offset ranges from 0 to tpa_sq-1 + + __local numtyp4 ts1[SHARED_SIZE]; + __local numtyp4 ts2[SHARED_SIZE]; + __local numtyp4 ts4[SHARED_SIZE]; + if (tid cutsq[ijparam]) continue; + numtyp r1 = ucl_sqrt(rsq1); + numtyp r1inv = ucl_rsqrt(rsq1); + + // look up for zeta_ij + + //int jj = (nbor_j-offset_j-2*nbor_pitch) / n_stride; + //int idx = jj*n_stride + i*t_per_atom + offset_j; + int idx; + zeta_idx(dev_nbor,dev_packed, nbor_pitch, n_stride, t_per_atom, + i, nbor_j, offset_j, idx); + acctyp4 zeta_ij = zetaij[idx]; // fetch(zeta_ij,idx,zeta_tex); + numtyp force = zeta_ij.x*tpainv; + numtyp prefactor = zeta_ij.y; + f.x += delr1[0]*force; + f.y += delr1[1]*force; + f.z += delr1[2]*force; + + if (eflag>0) { + energy+=zeta_ij.z*tpainv; + } + if (vflag>0) { + numtyp mforce = -force; + virial[0] += delr1[0]*delr1[0]*mforce; + virial[1] += delr1[1]*delr1[1]*mforce; + virial[2] += delr1[2]*delr1[2]*mforce; + virial[3] += delr1[0]*delr1[1]*mforce; + virial[4] += delr1[0]*delr1[2]*mforce; + virial[5] += delr1[1]*delr1[2]*mforce; + } + + int nbor_k=nborj_start-offset_j+offset_k; + for ( ; nbor_k cutsq[ijkparam]) continue; + numtyp r2 = ucl_sqrt(rsq2); + numtyp r2inv = ucl_rsqrt(rsq2); + + numtyp fi[3], fj[3], fk[3]; + numtyp4 ts1_ijkparam = ts1[ijkparam]; //fetch4(ts1_ijkparam,ijkparam,ts1_tex); + lam3 = ts1_ijkparam.z; + powermint = ts1_ijkparam.w; + numtyp4 ts2_ijkparam = ts2[ijkparam]; //fetch4(ts2_ijkparam,ijkparam,ts2_tex); + bigr = ts2_ijkparam.z; + bigd = ts2_ijkparam.w; + numtyp4 ts4_ijkparam = ts4[ijkparam]; //fetch4(ts4_ijkparam,ijkparam,ts4_tex); + c = ts4_ijkparam.x; + d = ts4_ijkparam.y; + h = ts4_ijkparam.z; + gamma = ts4_ijkparam.w; + if (vflag>0) + attractive(bigr, bigd, powermint, lam3, c, d, h, gamma, + prefactor, r1, r1inv, r2, r2inv, delr1, delr2, fi, fj, fk); + else + attractive_fi(bigr, bigd, powermint, lam3, c, d, h, gamma, + prefactor, r1, r1inv, r2, r2inv, delr1, delr2, fi); + f.x += fi[0]; + f.y += fi[1]; + f.z += fi[2]; + + if (vflag>0) { + acctyp v[6]; + numtyp pre = (numtyp)2.0; + if (evatom==1) pre = TWOTHIRD; + v[0] = pre*(delr1[0]*fj[0] + delr2[0]*fk[0]); + v[1] = pre*(delr1[1]*fj[1] + delr2[1]*fk[1]); + v[2] = pre*(delr1[2]*fj[2] + delr2[2]*fk[2]); + v[3] = pre*(delr1[0]*fj[1] + delr2[0]*fk[1]); + v[4] = pre*(delr1[0]*fj[2] + delr2[0]*fk[2]); + v[5] = pre*(delr1[1]*fj[2] + delr2[1]*fk[2]); + + virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2]; + virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5]; + } + } // nbor_k + } // for nbor_j + + store_answers_p(f,energy,virial,ii,inum,tid,tpa_sq, + offset,eflag,vflag,ans,engv); + } // if ii +} + +__kernel void k_tersoff_zbl_three_end(const __global numtyp4 *restrict x_, + const __global numtyp4 *restrict ts1_in, + const __global numtyp4 *restrict ts2_in, + const __global numtyp4 *restrict ts4_in, + const __global numtyp *restrict cutsq, + const __global int *restrict map, + const __global int *restrict elem2param, + const int nelements, const int nparams, + const __global acctyp4 *restrict zetaij, + const __global int * dev_nbor, + const __global int * dev_packed, + __global acctyp4 *restrict ans, + __global acctyp *restrict engv, + const int eflag, const int vflag, + const int inum, const int nbor_pitch, + const int t_per_atom) { + __local int tpa_sq, n_stride; + tpa_sq=fast_mul(t_per_atom,t_per_atom); + numtyp lam3, powermint, bigr, bigd, c, d, h, gamma; + + int tid, ii, offset; + atom_info(tpa_sq,ii,tid,offset); + + __local numtyp4 ts1[SHARED_SIZE]; + __local numtyp4 ts2[SHARED_SIZE]; + __local numtyp4 ts4[SHARED_SIZE]; + if (tid cutsq[ijparam]) continue; + + numtyp mdelr1[3]; + mdelr1[0] = -delr1[0]; + mdelr1[1] = -delr1[1]; + mdelr1[2] = -delr1[2]; + + int nbor_k=j+nbor_pitch; + int numk=dev_nbor[nbor_k]; + if (dev_nbor==dev_packed) { + nbor_k+=nbor_pitch+fast_mul(j,t_per_atom-1); + k_end=nbor_k+fast_mul(numk/t_per_atom,n_stride)+(numk & (t_per_atom-1)); + nbor_k+=offset_k; + } else { + nbor_k+=nbor_pitch; + nbor_k=dev_nbor[nbor_k]; + k_end=nbor_k+numk; + nbor_k+=offset_k; + } + int nbork_start = nbor_k; + + // look up for zeta_ji: find i in the j's neighbor list + int m = tid / t_per_atom; + int ijnum = -1; + for ( ; nbor_k= 0) { + offset_kf = offset_k; + } else { + ijnum = red_acc[2*m+0]; + offset_kf = red_acc[2*m+1]; + } + + //int iix = (ijnum - offset_kf - 2*nbor_pitch) / n_stride; + //int idx = iix*n_stride + j*t_per_atom + offset_kf; + int idx; + zeta_idx(dev_nbor,dev_packed, nbor_pitch, n_stride, t_per_atom, + j, ijnum, offset_kf, idx); + acctyp4 zeta_ji = zetaij[idx]; // fetch(zeta_ji,idx,zeta_tex); + numtyp force = zeta_ji.x*tpainv; + numtyp prefactor_ji = zeta_ji.y; + f.x += delr1[0]*force; + f.y += delr1[1]*force; + f.z += delr1[2]*force; + + if (eflag>0) { + energy+=zeta_ji.z*tpainv; + } + if (vflag>0) { + numtyp mforce = -force; + virial[0] += mdelr1[0]*mdelr1[0]*mforce; + virial[1] += mdelr1[1]*mdelr1[1]*mforce; + virial[2] += mdelr1[2]*mdelr1[2]*mforce; + virial[3] += mdelr1[0]*mdelr1[1]*mforce; + virial[4] += mdelr1[0]*mdelr1[2]*mforce; + virial[5] += mdelr1[1]*mdelr1[2]*mforce; + } + + // attractive forces + for (nbor_k = nbork_start ; nbor_k cutsq[jikparam]) continue; + numtyp r2 = ucl_sqrt(rsq2); + numtyp r2inv = ucl_rsqrt(rsq2); + numtyp4 ts1_param, ts2_param, ts4_param; + numtyp fi[3]; + + ts1_param = ts1[jikparam]; //fetch4(ts1_jikparam,jikparam,ts1_tex); + lam3 = ts1_param.z; + powermint = ts1_param.w; + ts2_param = ts2[jikparam]; //fetch4(ts2_jikparam,jikparam,ts2_tex); + bigr = ts2_param.z; + bigd = ts2_param.w; + ts4_param = ts4[jikparam]; //fetch4(ts4_jikparam,jikparam,ts4_tex); + c = ts4_param.x; + d = ts4_param.y; + h = ts4_param.z; + gamma = ts4_param.w; + attractive_fj(bigr, bigd, powermint, lam3, c, d, h, gamma, + prefactor_ji, r1, r1inv, r2, r2inv, mdelr1, delr2, fi); + f.x += fi[0]; + f.y += fi[1]; + f.z += fi[2]; + + //int kk = (nbor_k - offset_k - 2*nbor_pitch) / n_stride; + //int idx = kk*n_stride + j*t_per_atom + offset_k; + int idx; + zeta_idx(dev_nbor,dev_packed, nbor_pitch, n_stride, t_per_atom, + j, nbor_k, offset_k, idx); + acctyp4 zeta_jk = zetaij[idx]; // fetch(zeta_jk,idx,zeta_tex); + numtyp prefactor_jk = zeta_jk.y; + int jkiparam=elem2param[jtype*nelements*nelements+ktype*nelements+itype]; + ts1_param = ts1[jkiparam]; //fetch4(ts1_jkiparam,jkiparam,ts1_tex); + lam3 = ts1_param.z; + powermint = ts1_param.w; + ts2_param = ts2[jkiparam]; //fetch4(ts2_jkiparam,jkiparam,ts2_tex); + bigr = ts2_param.z; + bigd = ts2_param.w; + ts4_param = ts4[jkiparam]; //fetch4(ts4_jkiparam,jkiparam,ts4_tex); + c = ts4_param.x; + d = ts4_param.y; + h = ts4_param.z; + gamma = ts4_param.w; + attractive_fk(bigr, bigd, powermint, lam3, c, d, h, gamma, + prefactor_jk, r2, r2inv, r1, r1inv, delr2, mdelr1, fi); + f.x += fi[0]; + f.y += fi[1]; + f.z += fi[2]; + } // for nbor_k + } // for nbor_j + + #ifdef THREE_CONCURRENT + store_answers(f,energy,virial,ii,inum,tid,tpa_sq,offset, + eflag,vflag,ans,engv); + #else + store_answers_p(f,energy,virial,ii,inum,tid,tpa_sq,offset, + eflag,vflag,ans,engv); + #endif + } // if ii +} + +__kernel void k_tersoff_zbl_three_end_vatom(const __global numtyp4 *restrict x_, + const __global numtyp4 *restrict ts1_in, + const __global numtyp4 *restrict ts2_in, + const __global numtyp4 *restrict ts4_in, + const __global numtyp *restrict cutsq, + const __global int *restrict map, + const __global int *restrict elem2param, + const int nelements, const int nparams, + const __global acctyp4 *restrict zetaij, + const __global int * dev_nbor, + const __global int * dev_packed, + __global acctyp4 *restrict ans, + __global acctyp *restrict engv, + const int eflag, const int vflag, + const int inum, const int nbor_pitch, + const int t_per_atom) { + __local int tpa_sq, n_stride; + tpa_sq=fast_mul(t_per_atom,t_per_atom); + numtyp lam3, powermint, bigr, bigd, c, d, h, gamma; + + int tid, ii, offset; + atom_info(tpa_sq,ii,tid,offset); + + __local numtyp4 ts1[SHARED_SIZE]; + __local numtyp4 ts2[SHARED_SIZE]; + __local numtyp4 ts4[SHARED_SIZE]; + if (tid cutsq[ijparam]) continue; + + numtyp mdelr1[3]; + mdelr1[0] = -delr1[0]; + mdelr1[1] = -delr1[1]; + mdelr1[2] = -delr1[2]; + + int nbor_k=j+nbor_pitch; + int numk=dev_nbor[nbor_k]; + if (dev_nbor==dev_packed) { + nbor_k+=nbor_pitch+fast_mul(j,t_per_atom-1); + k_end=nbor_k+fast_mul(numk/t_per_atom,n_stride)+(numk & (t_per_atom-1)); + nbor_k+=offset_k; + } else { + nbor_k+=nbor_pitch; + nbor_k=dev_nbor[nbor_k]; + k_end=nbor_k+numk; + nbor_k+=offset_k; + } + int nbork_start = nbor_k; + + // look up for zeta_ji + int m = tid / t_per_atom; + int ijnum = -1; + for ( ; nbor_k= 0) { + offset_kf = offset_k; + } else { + ijnum = red_acc[2*m+0]; + offset_kf = red_acc[2*m+1]; + } + + //int iix = (ijnum - offset_kf - 2*nbor_pitch) / n_stride; + //int idx = iix*n_stride + j*t_per_atom + offset_kf; + int idx; + zeta_idx(dev_nbor,dev_packed, nbor_pitch, n_stride, t_per_atom, + j, ijnum, offset_kf, idx); + acctyp4 zeta_ji = zetaij[idx]; // fetch(zeta_ji,idx,zeta_tex); + numtyp force = zeta_ji.x*tpainv; + numtyp prefactor_ji = zeta_ji.y; + f.x += delr1[0]*force; + f.y += delr1[1]*force; + f.z += delr1[2]*force; + + if (eflag>0) { + energy+=zeta_ji.z*tpainv; + } + if (vflag>0) { + numtyp mforce = -force; + virial[0] += mdelr1[0]*mdelr1[0]*mforce; + virial[1] += mdelr1[1]*mdelr1[1]*mforce; + virial[2] += mdelr1[2]*mdelr1[2]*mforce; + virial[3] += mdelr1[0]*mdelr1[1]*mforce; + virial[4] += mdelr1[0]*mdelr1[2]*mforce; + virial[5] += mdelr1[1]*mdelr1[2]*mforce; + } + + // attractive forces + for (nbor_k = nbork_start; nbor_k cutsq[jikparam]) continue; + numtyp r2 = ucl_sqrt(rsq2); + numtyp r2inv = ucl_rsqrt(rsq2); + + numtyp fi[3], fj[3], fk[3]; + numtyp4 ts1_param, ts2_param, ts4_param; + ts1_param = ts1[jikparam]; //fetch4(ts1_jikparam,jikparam,ts1_tex); + lam3 = ts1_param.z; + powermint = ts1_param.w; + ts2_param = ts2[jikparam]; //fetch4(ts2_jikparam,jikparam,ts2_tex); + bigr = ts2_param.z; + bigd = ts2_param.w; + ts4_param = ts4[jikparam]; //fetch4(ts4_jikparam,jikparam,ts4_tex); + c = ts4_param.x; + d = ts4_param.y; + h = ts4_param.z; + gamma = ts4_param.w; + attractive(bigr, bigd, powermint, lam3, c, d, h, gamma, + prefactor_ji, r1, r1inv, r2, r2inv, mdelr1, delr2, fi, fj, fk); + f.x += fj[0]; + f.y += fj[1]; + f.z += fj[2]; + + virial[0] += TWOTHIRD*(mdelr1[0]*fj[0] + delr2[0]*fk[0]); + virial[1] += TWOTHIRD*(mdelr1[1]*fj[1] + delr2[1]*fk[1]); + virial[2] += TWOTHIRD*(mdelr1[2]*fj[2] + delr2[2]*fk[2]); + virial[3] += TWOTHIRD*(mdelr1[0]*fj[1] + delr2[0]*fk[1]); + virial[4] += TWOTHIRD*(mdelr1[0]*fj[2] + delr2[0]*fk[2]); + virial[5] += TWOTHIRD*(mdelr1[1]*fj[2] + delr2[1]*fk[2]); + + //int kk = (nbor_k - offset_k - 2*nbor_pitch) / n_stride; + //int idx = kk*n_stride + j*t_per_atom + offset_k; + int idx; + zeta_idx(dev_nbor,dev_packed, nbor_pitch, n_stride, t_per_atom, + j, nbor_k, offset_k, idx); + acctyp4 zeta_jk = zetaij[idx]; // fetch(zeta_jk,idx,zeta_tex); + numtyp prefactor_jk = zeta_jk.y; + + int jkiparam=elem2param[jtype*nelements*nelements+ktype*nelements+itype]; + ts1_param = ts1[jkiparam]; //fetch4(ts1_jkiparam,jkiparam,ts1_tex); + lam3 = ts1_param.z; + powermint = ts1_param.w; + ts2_param = ts2[jkiparam]; //fetch4(ts2_jkiparam,jkiparam,ts2_tex); + bigr = ts2_param.z; + bigd = ts2_param.w; + ts4_param = ts4[jkiparam]; //fetch4(ts4_jkiparam,jkiparam,ts4_tex); + c = ts4_param.x; + d = ts4_param.y; + h = ts4_param.z; + gamma = ts4_param.w; + attractive(bigr, bigd, powermint, lam3, c, d, h, gamma, + prefactor_jk, r2, r2inv, r1, r1inv, delr2, mdelr1, fi, fj, fk); + f.x += fk[0]; + f.y += fk[1]; + f.z += fk[2]; + + virial[0] += TWOTHIRD*(delr2[0]*fj[0] + mdelr1[0]*fk[0]); + virial[1] += TWOTHIRD*(delr2[1]*fj[1] + mdelr1[1]*fk[1]); + virial[2] += TWOTHIRD*(delr2[2]*fj[2] + mdelr1[2]*fk[2]); + virial[3] += TWOTHIRD*(delr2[0]*fj[1] + mdelr1[0]*fk[1]); + virial[4] += TWOTHIRD*(delr2[0]*fj[2] + mdelr1[0]*fk[2]); + virial[5] += TWOTHIRD*(delr2[1]*fj[2] + mdelr1[1]*fk[2]); + } + } // for nbor + + #ifdef THREE_CONCURRENT + store_answers(f,energy,virial,ii,inum,tid,tpa_sq,offset, + eflag,vflag,ans,engv); + #else + store_answers_p(f,energy,virial,ii,inum,tid,tpa_sq,offset, + eflag,vflag,ans,engv); + #endif + } // if ii +} + diff --git a/lib/gpu/lal_tersoff_zbl.h b/lib/gpu/lal_tersoff_zbl.h new file mode 100644 index 0000000000..cc0b848845 --- /dev/null +++ b/lib/gpu/lal_tersoff_zbl.h @@ -0,0 +1,123 @@ +/*************************************************************************** + tersoff_zbl.h + ------------------- + Trung Dac Nguyen + + Class for acceleration of the tersoff pair style. + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : + email : ndactrung@gmail.com + ***************************************************************************/ + +#ifndef LAL_TERSOFF_ZBL_H +#define LAL_TERSOFF_ZBL_H + +#include "lal_base_three.h" + +namespace LAMMPS_AL { + +template +class TersoffZBL : public BaseThree { + public: + TersoffZBL(); + ~TersoffZBL(); + + /// 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* Z_i, const double* Z_j, + const double* ZBLcut, const double* ZBLexpscale, const double global_e, + const double global_a_0, const double global_epsilon_0, 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; + /// ts6.x = Z_i, ts6.y = Z_j, ts6.z = ZBLcut, ts6.w = ZBLexpscale + UCL_D_Vec ts6; + + UCL_D_Vec cutsq; + + UCL_D_Vec elem2param; + UCL_D_Vec map; + int _nparams,_nelements; + + /// Per-atom arrays: + /// zetaij.x = force, zetaij.y = prefactor, zetaij.z = evdwl, + /// zetaij.w = zetaij + UCL_D_Vec _zetaij; + + UCL_Kernel k_zeta; + UCL_Texture ts1_tex, ts2_tex, ts3_tex, ts4_tex, ts5_tex, ts6_tex; + + int _max_nbors; + numtyp _global_e,_global_a_0,_global_epsilon_0; + + private: + bool _allocated; + void loop(const bool _eflag, const bool _vflag, const int evatom); +}; + +} + +#endif + diff --git a/lib/gpu/lal_tersoff_zbl_ext.cpp b/lib/gpu/lal_tersoff_zbl_ext.cpp new file mode 100644 index 0000000000..fce240f8fe --- /dev/null +++ b/lib/gpu/lal_tersoff_zbl_ext.cpp @@ -0,0 +1,146 @@ +/*************************************************************************** + tersoff_zbl_ext.cpp + ------------------- + Trung Dac Nguyen + + Functions for LAMMPS access to tersoff/zbl acceleration routines. + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : + email : ndactrung@gmail.com + ***************************************************************************/ + +#include +#include +#include + +#include "lal_tersoff_zbl.h" + +using namespace std; +using namespace LAMMPS_AL; + +static TersoffZBL TSZMF; + +// --------------------------------------------------------------------------- +// Allocate memory on host and device and copy constants to device +// --------------------------------------------------------------------------- +int tersoff_zbl_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_Z_i, + const double* ts_Z_j, const double* ts_ZBLcut, + const double* ts_ZBLexpscale, const double global_e, + const double global_a_0, const double global_epsilon_0, + const double* ts_cutsq) { + TSZMF.clear(); + gpu_mode=TSZMF.device->gpu_mode(); + double gpu_split=TSZMF.device->particle_split(); + int first_gpu=TSZMF.device->first_device(); + int last_gpu=TSZMF.device->last_device(); + int world_me=TSZMF.device->world_me(); + int gpu_rank=TSZMF.device->gpu_rank(); + int procs_per_gpu=TSZMF.device->procs_per_gpu(); + + // disable host/device split for now + if (gpu_split != 1.0) + return -8; + + TSZMF.device->init_message(screen,"tersoff/zbl/gpu",first_gpu,last_gpu); + + bool message=false; + if (TSZMF.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=TSZMF.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_Z_i, ts_Z_j, + ts_ZBLcut, ts_ZBLexpscale, global_e, global_a_0, + global_epsilon_0, ts_cutsq); + + TSZMF.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) + TSZMF.estimate_gpu_overhead(); + return init_ok; +} + +void tersoff_zbl_gpu_clear() { + TSZMF.clear(); +} + +int ** tersoff_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 TSZMF.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_zbl_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) { + TSZMF.compute(ago,nlocal,nall,nlist,host_x,host_type,ilist,numj, + firstneigh,eflag,vflag,eatom,vatom,host_start,cpu_time,success); +} + +double tersoff_zbl_gpu_bytes() { + return TSZMF.host_memory_usage(); +} + + diff --git a/lib/gpu/lal_tersoff_zbl_extra.h b/lib/gpu/lal_tersoff_zbl_extra.h new file mode 100644 index 0000000000..79afb4de82 --- /dev/null +++ b/lib/gpu/lal_tersoff_zbl_extra.h @@ -0,0 +1,690 @@ +/// ************************************************************************** +// tersoff_zbl_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_ZBL_EXTRA_H +#define LAL_TERSOFF_ZBL_EXTRA_H + +#ifdef NV_KERNEL +#include "lal_aux_fun1.h" +#else +#endif + +#define MY_PI (numtyp)3.14159265358979323846 +#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 F_fermi(const numtyp r, const numtyp param_ZBLcut, + const numtyp param_ZBLexpscale) +{ + return ucl_recip((numtyp)1.0+ucl_exp(-param_ZBLexpscale*(r-param_ZBLcut))); +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline numtyp F_fermi_d(const numtyp r, const numtyp param_ZBLcut, + const numtyp param_ZBLexpscale) +{ + numtyp a = ucl_exp(-param_ZBLexpscale*(r-param_ZBLcut)); + numtyp b = (numtyp)1.0 + a; + return param_ZBLexpscale*a*ucl_recip(b*b); +} + +/* ---------------------------------------------------------------------- */ + +ucl_inline numtyp ters_fa(const numtyp r, + const numtyp param_bigb, + const numtyp param_bigr, + const numtyp param_bigd, + const numtyp param_lam2, + const numtyp param_ZBLcut, + const numtyp param_ZBLexpscale) +{ + 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)*F_fermi(r,param_ZBLcut,param_ZBLexpscale); +} + +/* ---------------------------------------------------------------------- */ + +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, + const numtyp param_ZBLcut, + const numtyp param_ZBLexpscale) +{ + if (r > param_bigr + param_bigd) return (numtyp)0.0; + numtyp f = F_fermi(r,param_ZBLcut,param_ZBLexpscale); + return param_bigb * ucl_exp(-param_lam2 * r) * + (param_lam2 * ters_fc(r,param_bigr,param_bigd) * f - + ters_fc_d(r,param_bigr,param_bigd) * f - + ters_fc(r,param_bigr,param_bigd) * F_fermi_d(r,param_ZBLcut,param_ZBLexpscale)); +} + +/* ---------------------------------------------------------------------- */ + +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) * + // error in negligible 2nd term fixed 9/30/2015 + // (1.0 - 0.5*(1.0 + 1.0/(2.0*param->powern)) * + ((numtyp)1.0 - ((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 = (numtyp)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 = (numtyp)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 = (numtyp)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 = (numtyp)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 param_Z_i, + const numtyp param_Z_j, + const numtyp param_ZBLcut, + const numtyp param_ZBLexpscale, + const numtyp global_e, + const numtyp global_a_0, + const numtyp global_epsilon_0, + const numtyp rsq, + const int eflag, + numtyp *ans) +{ + numtyp r,tmp_fc,tmp_fc_d,tmp_exp; + + // Tersoff repulsive portion + + 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); + + numtyp fforce_ters = param_biga * tmp_exp * (tmp_fc_d - tmp_fc*param_lam1); + numtyp eng_ters = tmp_fc * param_biga * tmp_exp; + + // ZBL repulsive portion + + numtyp esq = global_e*global_e; + numtyp a_ij = ((numtyp)0.8854*global_a_0) / + (ucl_powr(param_Z_i,(numtyp)0.23) + ucl_powr(param_Z_j,(numtyp)0.23)); + numtyp premult = (param_Z_i * param_Z_j * esq)/((numtyp)4.0*MY_PI*global_epsilon_0); + numtyp r_ov_a = r/a_ij; + numtyp t1 = (numtyp)0.1818*ucl_exp((numtyp)-3.2*r_ov_a); + numtyp t2 = (numtyp)0.5099*ucl_exp((numtyp)-0.9423*r_ov_a); + numtyp t3 = (numtyp)0.2802*ucl_exp((numtyp)-0.4029*r_ov_a); + numtyp t4 = (numtyp)0.02817*ucl_exp((numtyp)-0.2016*r_ov_a); + numtyp phi = t1 + t2 + t3 + t4; + numtyp dphi = (numtyp)-3.2*t1 - (numtyp)0.9423*t2 - (numtyp)0.4029*t3 - + (numtyp)0.2016*t4; + dphi *= ucl_recip(a_ij); +/* + numtyp phi = (numtyp)0.1818*ucl_exp((numtyp)-3.2*r_ov_a) + + (numtyp)0.5099*ucl_exp((numtyp)-0.9423*r_ov_a) + + (numtyp)0.2802*ucl_exp((numtyp)-0.4029*r_ov_a) + + (numtyp)0.02817*ucl_exp((numtyp)-0.2016*r_ov_a); + numtyp dphi = ucl_recip(a_ij) * ((numtyp)-3.2*(numtyp)0.1818*ucl_exp((numtyp)-3.2*r_ov_a) - + (numtyp)0.9423*(numtyp)0.5099*ucl_exp((numtyp)-0.9423*r_ov_a) - + (numtyp)0.4029*(numtyp)0.2802*ucl_exp((numtyp)-0.4029*r_ov_a) - + (numtyp)0.2016*(numtyp)0.02817*ucl_exp((numtyp)-0.2016*r_ov_a)); +*/ + numtyp rinv = ucl_recip(r); + numtyp fforce_ZBL = premult*(-phi)/rsq + premult*dphi*rinv; + numtyp eng_ZBL = premult*rinv*phi; + + // combine two parts with smoothing by Fermi-like function + // ans[0] = fforce + numtyp f = F_fermi(r,param_ZBLcut,param_ZBLexpscale); + numtyp f_d = F_fermi_d(r,param_ZBLcut,param_ZBLexpscale); + ans[0] = -(-f_d * eng_ZBL + ((numtyp)1.0 - f)*fforce_ZBL + f_d*eng_ters + + f*fforce_ters) * rinv; + + // ans[1] = eng + if (eflag) ans[1] = ((numtyp)1.0 - f)*eng_ZBL + f*eng_ters; +} + +/* ---------------------------------------------------------------------- */ + +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 param_ZBLcut, + const numtyp param_ZBLexpscale, + 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,param_ZBLcut,param_ZBLexpscale); + fa_d = ters_fa_d(r,param_bigb,param_bigr,param_bigd,param_lam2,param_ZBLcut,param_ZBLexpscale); + 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 + +