From 785a794d3933c56d495d2970518ab653c8d1ba6c Mon Sep 17 00:00:00 2001 From: Trung Nguyen Date: Wed, 1 Sep 2021 14:37:11 -0500 Subject: [PATCH] Added and renamed API to make room for additional kernels (udirect2b only computes the field and fieldp, not accumulating forces, energies, nor virials) --- lib/gpu/lal_amoeba.cpp | 35 ++++- lib/gpu/lal_amoeba.cu | 20 +-- lib/gpu/lal_amoeba.h | 5 +- lib/gpu/lal_amoeba_ext.cpp | 31 ++-- lib/gpu/lal_base_amoeba.cpp | 155 ++++++++++++++++---- lib/gpu/lal_base_amoeba.h | 24 ++- src/AMOEBA/pair_amoeba.h | 2 +- src/GPU/pair_amoeba_gpu.cpp | 283 ++++++++++++++++++++++++++++++------ src/GPU/pair_amoeba_gpu.h | 4 + 9 files changed, 448 insertions(+), 111 deletions(-) diff --git a/lib/gpu/lal_amoeba.cpp b/lib/gpu/lal_amoeba.cpp index a3bd653efd..c7b4872db0 100644 --- a/lib/gpu/lal_amoeba.cpp +++ b/lib/gpu/lal_amoeba.cpp @@ -125,10 +125,10 @@ double AmoebaT::host_memory_usage() const { } // --------------------------------------------------------------------------- -// Calculate energies, forces, and torques +// Calculate the polar real-space term, returning tep // --------------------------------------------------------------------------- template -int AmoebaT::loop(const int eflag, const int vflag) { +int AmoebaT::polar_real(const int eflag, const int vflag) { // Compute the block size and grid size to keep all cores busy const int BX=this->block_size(); int GX=static_cast(ceil(static_cast(this->ans->inum())/ @@ -140,9 +140,7 @@ int AmoebaT::loop(const int eflag, const int vflag) { this->time_pair.start(); this->k_polar.set_size(GX,BX); - - this->k_polar.run(&this->atom->x, &this->atom->extra, - &damping, &sp_polar, + this->k_polar.run(&this->atom->x, &this->atom->extra, &damping, &sp_polar, &this->nbor->dev_nbor, &this->_nbor_data->begin(), &this->ans->force, &this->ans->engv, &this->_tep, &eflag, &vflag, &ainum, &_nall, &nbor_pitch, @@ -152,5 +150,32 @@ int AmoebaT::loop(const int eflag, const int vflag) { return GX; } +// --------------------------------------------------------------------------- +// Calculate the polar real-space term, returning tep +// --------------------------------------------------------------------------- +template +int AmoebaT::udirect2b(const int eflag, const int vflag) { + // Compute the block size and grid size to keep all cores busy + const int BX=this->block_size(); + int GX=static_cast(ceil(static_cast(this->ans->inum())/ + (BX/this->_threads_per_atom))); + + int _nall=this->atom->nall(); + int ainum=this->ans->inum(); + int nbor_pitch=this->nbor->nbor_pitch(); + this->time_pair.start(); +/* + this->k_polar.set_size(GX,BX); + this->k_polar.run(&this->atom->x, &this->atom->extra, &damping, &sp_polar, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &this->ans->force, &this->ans->engv, &this->_tep, + &eflag, &vflag, &ainum, &_nall, &nbor_pitch, + &this->_threads_per_atom, + &_aewald, &_felec, &_off2, &_polar_dscale, &_polar_uscale); +*/ + this->time_pair.stop(); + return GX; +} + template class Amoeba; } diff --git a/lib/gpu/lal_amoeba.cu b/lib/gpu/lal_amoeba.cu index 1f5fb42438..3d28939d42 100644 --- a/lib/gpu/lal_amoeba.cu +++ b/lib/gpu/lal_amoeba.cu @@ -715,11 +715,6 @@ __kernel void k_amoeba_udirect2b(const __global numtyp4 *restrict x_, //numtyp4 xi__; if (ii (numtyp)0.0) alsq2n = (numtyp)1.0 / (MY_PIS*aewald); + numtyp aesq2 = (numtyp)2.0 * aewald*aewald; + numtyp aesq2n = (numtyp)0.0; + if (aewald > (numtyp)0.0) aesq2n = (numtyp)1.0 / (MY_PIS*aewald); for ( ; nbor { numtyp _aewald, _felec, _off2, _polar_dscale, _polar_uscale; numtyp _qqrd2e; - private: + protected: bool _allocated; - int loop(const int eflag, const int vflag); + int polar_real(const int eflag, const int vflag); + int udirect2b(const int eflag, const int vflag); }; } diff --git a/lib/gpu/lal_amoeba_ext.cpp b/lib/gpu/lal_amoeba_ext.cpp index a7959ed93e..9fa3c7f75b 100644 --- a/lib/gpu/lal_amoeba_ext.cpp +++ b/lib/gpu/lal_amoeba_ext.cpp @@ -105,7 +105,7 @@ void amoeba_gpu_clear() { AMOEBAMF.clear(); } -int** amoeba_gpu_compute_n(const int ago, const int inum_full, +int** amoeba_gpu_compute_polar_real(const int ago, const int inum_full, const int nall, double **host_x, int *host_type, int *host_amtype, int *host_amgroup, double **host_rpole, double **host_uind, double **host_uinp, @@ -116,7 +116,7 @@ int** amoeba_gpu_compute_n(const int ago, const int inum_full, int **ilist, int **jnum, const double cpu_time, bool &success, double *host_q, double *boxlo, double *prd, void **tep_ptr) { - return AMOEBAMF.compute(ago, inum_full, nall, host_x, host_type, + return AMOEBAMF.compute_polar_real(ago, inum_full, nall, host_x, host_type, host_amtype, host_amgroup, host_rpole, host_uind, host_uinp, sublo, subhi, tag, nspecial, special, nspecial15, special15, eflag, vflag, eatom, @@ -124,18 +124,21 @@ int** amoeba_gpu_compute_n(const int ago, const int inum_full, host_q, boxlo, prd, tep_ptr); } -void amoeba_gpu_compute(const int ago, const int inum_full, const int nall, - double **host_x, int *host_type, int *host_amtype, int *host_amgroup, - double **host_rpole, double **host_uind, double **host_uinp, - 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, double *host_q, - const int nlocal, double *boxlo, double *prd, void **tep_ptr) { - AMOEBAMF.compute(ago,inum_full, nall, host_x, host_type, - host_amtype, host_amgroup, host_rpole, host_uind, host_uinp, - ilist, numj, firstneigh, eflag, vflag, eatom, vatom, - host_start, cpu_time, success, host_q, nlocal, boxlo, prd, tep_ptr); +int** amoeba_gpu_compute_udirect2b(const int ago, const int inum_full, + const int nall, double **host_x, int *host_type, + int *host_amtype, int *host_amgroup, double **host_rpole, + double *sublo, double *subhi, tagint *tag, int **nspecial, + tagint **special, int *nspecial15, tagint** special15, + 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, double *host_q, double *boxlo, + double *prd, void **fieldp_ptr) { + return AMOEBAMF.compute_udirect2b(ago, inum_full, nall, host_x, host_type, + host_amtype, host_amgroup, host_rpole, sublo, subhi, + tag, nspecial, special, nspecial15, special15, + eflag, vflag, eatom, vatom, host_start, ilist, jnum, + cpu_time, success, host_q, boxlo, prd, fieldp_ptr); } double amoeba_gpu_bytes() { diff --git a/lib/gpu/lal_base_amoeba.cpp b/lib/gpu/lal_base_amoeba.cpp index c5f4a01222..0c9a422cec 100644 --- a/lib/gpu/lal_base_amoeba.cpp +++ b/lib/gpu/lal_base_amoeba.cpp @@ -118,8 +118,9 @@ int BaseAmoebaT::init_atomic(const int nlocal, const int nall, if (ef_nall==0) ef_nall=2000; - _max_tep_size=static_cast(static_cast(ef_nall)*1.10); - _tep.alloc(_max_tep_size*4,*(this->ucl_device),UCL_READ_WRITE,UCL_READ_WRITE); + _max_alloc_size=static_cast(static_cast(ef_nall)*1.10); + _fieldp.alloc(_max_alloc_size*6,*(this->ucl_device),UCL_READ_WRITE,UCL_READ_WRITE); + _tep.alloc(_max_alloc_size*4,*(this->ucl_device),UCL_READ_WRITE,UCL_READ_WRITE); dev_nspecial15.alloc(nall,*(this->ucl_device),UCL_READ_ONLY); dev_special15.alloc(_maxspecial15*nall,*(this->ucl_device),UCL_READ_ONLY); dev_special15_t.alloc(nall*_maxspecial15,*(this->ucl_device),UCL_READ_ONLY); @@ -149,6 +150,7 @@ void BaseAmoebaT::clear_atomic() { ans->clear(); _tep.clear(); + _fieldp.clear(); dev_nspecial15.clear(); dev_special15.clear(); dev_special15_t.clear(); @@ -250,9 +252,9 @@ void BaseAmoebaT::compute(const int f_ago, const int inum_full, const int nall, // ------------------- Resize _tep array ------------------------ - if (nall>_max_tep_size) { - _max_tep_size=static_cast(static_cast(nall)*1.10); - _tep.resize(_max_tep_size*4); + if (nall>_max_alloc_size) { + _max_alloc_size=static_cast(static_cast(nall)*1.10); + _tep.resize(_max_alloc_size*4); dev_nspecial15.clear(); dev_special15.clear(); @@ -296,17 +298,17 @@ void BaseAmoebaT::compute(const int f_ago, const int inum_full, const int nall, device->precompute(f_ago,nlocal,nall,host_x,host_type,success,host_q, boxlo, prd); - const int red_blocks=loop(eflag,vflag); + const int red_blocks=polar_real(eflag,vflag); ans->copy_answers(eflag_in,vflag_in,eatom,vatom,ilist,red_blocks); device->add_ans_object(ans); hd_balancer.stop_timer(); } // --------------------------------------------------------------------------- -// Reneighbor on GPU if necessary and then compute forces, virials, energies +// Reneighbor on GPU if necessary, and then compute polar real-space // --------------------------------------------------------------------------- template -int** BaseAmoebaT::compute(const int ago, const int inum_full, const int nall, +int** BaseAmoebaT::compute_polar_real(const int ago, const int inum_full, const int nall, double **host_x, int *host_type, int *host_amtype, int *host_amgroup, double **host_rpole, double **host_uind, double **host_uinp, @@ -336,9 +338,9 @@ int** BaseAmoebaT::compute(const int ago, const int inum_full, const int nall, // ------------------- Resize _tep array ------------------------ - if (nall>_max_tep_size) { - _max_tep_size=static_cast(static_cast(nall)*1.10); - _tep.resize(_max_tep_size*4); + if (nall>_max_alloc_size) { + _max_alloc_size=static_cast(static_cast(nall)*1.10); + _tep.resize(_max_alloc_size*4); dev_nspecial15.clear(); dev_special15.clear(); @@ -388,16 +390,16 @@ int** BaseAmoebaT::compute(const int ago, const int inum_full, const int nall, device->precompute(ago,inum_full,nall,host_x,host_type,success,host_q, boxlo, prd); - const int red_blocks=loop(eflag,vflag); + const int red_blocks=polar_real(eflag,vflag); ans->copy_answers(eflag_in,vflag_in,eatom,vatom,red_blocks); device->add_ans_object(ans); hd_balancer.stop_timer(); // copy tep from device to host - _tep.update_host(_max_tep_size*4,false); + _tep.update_host(_max_alloc_size*4,false); /* - printf("GPU lib: tep size = %d: max tep size = %d\n", this->_tep.cols(), _max_tep_size); + printf("GPU lib: tep size = %d: max tep size = %d\n", this->_tep.cols(), _max_alloc_size); for (int i = 0; i < 10; i++) { numtyp4* p = (numtyp4*)(&this->_tep[4*i]); printf("i = %d; tep = %f %f %f\n", i, p->x, p->y, p->z); @@ -406,6 +408,101 @@ int** BaseAmoebaT::compute(const int ago, const int inum_full, const int nall, return nbor->host_jlist.begin()-host_start; } +// --------------------------------------------------------------------------- +// Reneighbor on GPU if necessary, and then compute the direct real space part +// of the permanent field +// --------------------------------------------------------------------------- +template +int** BaseAmoebaT::compute_udirect2b(const int ago, const int inum_full, const int nall, + double **host_x, int *host_type, int *host_amtype, + int *host_amgroup, double **host_rpole, + double *sublo, double *subhi, tagint *tag, + int **nspecial, tagint **special, + int *nspecial15, tagint **special15, + const bool eflag_in, const bool vflag_in, + const bool eatom, const bool vatom, int &host_start, + int **ilist, int **jnum, const double cpu_time, + bool &success, double *host_q, double *boxlo, + double *prd, void** fieldp_ptr) { + acc_timers(); + int eflag, vflag; + if (eatom) eflag=2; + else if (eflag_in) eflag=1; + else eflag=0; + if (vatom) vflag=2; + else if (vflag_in) vflag=1; + else vflag=0; + + #ifdef LAL_NO_BLOCK_REDUCE + if (eflag) eflag=2; + if (vflag) vflag=2; + #endif + + set_kernel(eflag,vflag); + + // ------------------- Resize _fieldp array ------------------------ + + if (nall>_max_alloc_size) { + _max_alloc_size=static_cast(static_cast(nall)*1.10); + _fieldp.resize(_max_alloc_size*8); + + dev_nspecial15.clear(); + dev_special15.clear(); + dev_special15_t.clear(); + dev_nspecial15.alloc(nall,*(this->ucl_device),UCL_READ_ONLY); + dev_special15.alloc(_maxspecial15*nall,*(this->ucl_device),UCL_READ_ONLY); + dev_special15_t.alloc(nall*_maxspecial15,*(this->ucl_device),UCL_READ_ONLY); + } + *fieldp_ptr=_fieldp.host.begin(); + + if (inum_full==0) { + host_start=0; + // Make sure textures are correct if realloc by a different hybrid style + resize_atom(0,nall,success); + zero_timers(); + return nullptr; + } + + hd_balancer.balance(cpu_time); + int inum=hd_balancer.get_gpu_count(ago,inum_full); + ans->inum(inum); + host_start=inum; + + // Build neighbor list on GPU if necessary + if (ago==0) { + build_nbor_list(inum, inum_full-inum, nall, host_x, host_type, + sublo, subhi, tag, nspecial, special, nspecial15, special15, + success); + if (!success) + return nullptr; + atom->cast_q_data(host_q); + cast_extra_data(host_amtype, host_amgroup, host_rpole, nullptr, nullptr); + hd_balancer.start_timer(); + } else { + atom->cast_x_data(host_x,host_type); + atom->cast_q_data(host_q); + cast_extra_data(host_amtype, host_amgroup, host_rpole, nullptr, nullptr); + hd_balancer.start_timer(); + atom->add_x_data(host_x,host_type); + } + atom->add_q_data(); + atom->add_extra_data(); + + *ilist=nbor->host_ilist.begin(); + *jnum=nbor->host_acc.begin(); + + const int red_blocks=udirect2b(eflag,vflag); + //ans->copy_answers(eflag_in,vflag_in,eatom,vatom,red_blocks); + //device->add_ans_object(ans); + hd_balancer.stop_timer(); + + // copy field and fieldp from device to host + + //_fieldp.update_host(_max_field_size*8,false); + + return nbor->host_jlist.begin()-host_start; +} + template double BaseAmoebaT::host_memory_usage_atomic() const { return device->atom.host_memory_usage()+nbor->host_memory_usage()+ @@ -446,20 +543,24 @@ void BaseAmoebaT::cast_extra_data(int* amtype, int* amgroup, double** rpole, pextra[idx+3] = (numtyp)amgroup[i]; } - n += nstride*_nall; - for (int i = 0; i < _nall; i++) { - int idx = n+i*nstride; - pextra[idx] = uind[i][0]; - pextra[idx+1] = uind[i][1]; - pextra[idx+2] = uind[i][2]; + if (uind) { + n += nstride*_nall; + for (int i = 0; i < _nall; i++) { + int idx = n+i*nstride; + pextra[idx] = uind[i][0]; + pextra[idx+1] = uind[i][1]; + pextra[idx+2] = uind[i][2]; + } } - - n += nstride*_nall; - for (int i = 0; i < _nall; i++) { - int idx = n+i*nstride; - pextra[idx] = uinp[i][0]; - pextra[idx+1] = uinp[i][1]; - pextra[idx+2] = uinp[i][2]; + + if (uinp) { + n += nstride*_nall; + for (int i = 0; i < _nall; i++) { + int idx = n+i*nstride; + pextra[idx] = uinp[i][0]; + pextra[idx+1] = uinp[i][1]; + pextra[idx+2] = uinp[i][2]; + } } } diff --git a/lib/gpu/lal_base_amoeba.h b/lib/gpu/lal_base_amoeba.h index ac9c23e8a9..7ef94c776e 100644 --- a/lib/gpu/lal_base_amoeba.h +++ b/lib/gpu/lal_base_amoeba.h @@ -128,7 +128,7 @@ class BaseAmoeba { tagint **special, int *nspecial15, tagint **special15, bool &success); - /// Pair loop with host neighboring + /// Compute polar real-space with host neighboring (not active for now) void compute(const int f_ago, const int inum_full, const int nall, double **host_x, int *host_type, int *host_amtype, int *host_amgroup, double **host_rpole, double **host_uind, @@ -138,8 +138,8 @@ class BaseAmoeba { const double cpu_time, bool &success, double *charge, const int nlocal, double *boxlo, double *prd, void **tep_ptr); - /// Pair loop with device neighboring - int** compute(const int ago, const int inum_full, const int nall, + /// Compute polar real-space with device neighboring + int** compute_polar_real(const int ago, const int inum_full, const int nall, double **host_x, int *host_type, int *host_amtype, int *host_amgroup, double **host_rpole, double **host_uind, double **host_uinp, double *sublo, double *subhi, @@ -150,6 +150,17 @@ class BaseAmoeba { int **ilist, int **numj, const double cpu_time, bool &success, double *charge, double *boxlo, double *prd, void **tep_ptr); + /// Compute the direct real space part of the permanent field (udirect2b) with device neighboring + int** compute_udirect2b(const int ago, const int inum_full, const int nall, + double **host_x, int *host_type, int *host_amtype, + int *host_amgroup, double **host_rpole, double *sublo, double *subhi, + tagint *tag, int **nspecial, tagint **special, + int *nspecial15, tagint **special15, + 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, + double *charge, double *boxlo, double *prd, void **fieldp_ptr); + // -------------------------- DEVICE DATA ------------------------- /// Device Properties and Atom and Neighbor storage @@ -179,8 +190,8 @@ class BaseAmoeba { double** uind, double** uinp); /// Per-atom arrays - UCL_Vector _tep; - int _max_tep_size; + UCL_Vector _tep,_fieldp; + int _max_alloc_size; // ------------------------ FORCE/ENERGY DATA ----------------------- @@ -217,7 +228,8 @@ class BaseAmoeba { void compile_kernels(UCL_Device &dev, const void *pair_string, const char *k); - virtual int loop(const int eflag, const int vflag) = 0; + virtual int polar_real(const int eflag, const int vflag) = 0; + virtual int udirect2b(const int eflag, const int vflag) = 0; }; } diff --git a/src/AMOEBA/pair_amoeba.h b/src/AMOEBA/pair_amoeba.h index 4644d4a137..9d23fccdd8 100644 --- a/src/AMOEBA/pair_amoeba.h +++ b/src/AMOEBA/pair_amoeba.h @@ -369,7 +369,7 @@ class PairAmoeba : public Pair { void umutual1(double **, double **); void umutual2b(double **, double **); void udirect1(double **); - void udirect2b(double **, double **); + virtual void udirect2b(double **, double **); void dampmut(double, double, double, double *); void dampdir(double, double, double, double *, double *); void cholesky(int, double *, double *); diff --git a/src/GPU/pair_amoeba_gpu.cpp b/src/GPU/pair_amoeba_gpu.cpp index 09ba100e4e..3f4e72c0af 100644 --- a/src/GPU/pair_amoeba_gpu.cpp +++ b/src/GPU/pair_amoeba_gpu.cpp @@ -24,19 +24,24 @@ #include "error.h" #include "force.h" #include "gpu_extra.h" +#include "math_const.h" +#include "my_page.h" #include "neigh_list.h" #include "neigh_request.h" #include "neighbor.h" #include "suffix.h" - #include using namespace LAMMPS_NS; +using namespace MathConst; + +enum{MUTUAL,OPT,TCG,DIRECT}; // External functions from cuda library for atom decomposition int amoeba_gpu_init(const int ntypes, const int max_amtype, const double *host_pdamp, const double *host_thole, + const double *host_dirdamp, const double *host_special_polar_wscale, const double *host_special_polar_piscale, const double *host_special_polar_pscale, @@ -48,7 +53,17 @@ int amoeba_gpu_init(const int ntypes, const int max_amtype, const double polar_uscale, int& tep_size); void amoeba_gpu_clear(); -int ** amoeba_gpu_compute_n(const int ago, const int inum, const int nall, +int ** amoeba_gpu_compute_udirect2b(const int ago, const int inum, const int nall, + double **host_x, int *host_type, int *host_amtype, int *host_amgroup, + double **host_rpole, double *sublo, double *subhi, tagint *tag, int **nspecial, + tagint **special, int* nspecial15, tagint** special15, + 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, double *host_q, double *boxlo, double *prd, + void **fieldp_ptr); + +int ** amoeba_gpu_compute_polar_real(const int ago, const int inum, const int nall, double **host_x, int *host_type, int *host_amtype, int *host_amgroup, double **host_rpole, double **host_uind, double **host_uinp, double *sublo, double *subhi, tagint *tag, int **nspecial, @@ -58,15 +73,6 @@ int ** amoeba_gpu_compute_n(const int ago, const int inum, const int nall, int **ilist, int **jnum, const double cpu_time, bool &success, double *host_q, double *boxlo, double *prd, void **tep_ptr); -void amoeba_gpu_compute(const int ago, const int inum, - const int nall, double **host_x, int *host_type, - int *host_amtype, int *host_amgroup, - double **host_rpole, double **host_uind, double **host_uinp, - 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, double *host_q, const int nlocal, - double *boxlo, double *prd, void **tep_ptr); double amoeba_gpu_bytes(); @@ -80,6 +86,8 @@ PairAmoebaGPU::PairAmoebaGPU(LAMMPS *lmp) : PairAmoeba(lmp), gpu_mode(GPU_FORCE) reinitflag = 0; cpu_time = 0.0; suffix_flag |= Suffix::GPU; + fieldp_pinned = nullptr; + tep_pinned = nullptr; GPU_EXTRA::gpu_ready(lmp->modify, lmp->error); } @@ -102,42 +110,31 @@ void PairAmoebaGPU::polar_real() bool success = true; int *ilist, *numneigh, **firstneigh; - if (gpu_mode != GPU_FORCE) { - double sublo[3],subhi[3]; - if (domain->triclinic == 0) { - sublo[0] = domain->sublo[0]; - sublo[1] = domain->sublo[1]; - sublo[2] = domain->sublo[2]; - subhi[0] = domain->subhi[0]; - subhi[1] = domain->subhi[1]; - subhi[2] = domain->subhi[2]; - } else { - domain->bbox(domain->sublo_lamda,domain->subhi_lamda,sublo,subhi); - } - inum = atom->nlocal; - - firstneigh = amoeba_gpu_compute_n(neighbor->ago, inum, nall, atom->x, - atom->type, amtype, amgroup, - rpole, uind, uinp, sublo, subhi, - atom->tag, atom->nspecial, atom->special, - atom->nspecial15, atom->special15, - eflag, vflag, eflag_atom, vflag_atom, - host_start, &ilist, &numneigh, cpu_time, - success, atom->q, domain->boxlo, - domain->prd, &tep_pinned); - + + double sublo[3],subhi[3]; + if (domain->triclinic == 0) { + sublo[0] = domain->sublo[0]; + sublo[1] = domain->sublo[1]; + sublo[2] = domain->sublo[2]; + subhi[0] = domain->subhi[0]; + subhi[1] = domain->subhi[1]; + subhi[2] = domain->subhi[2]; } else { - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - - amoeba_gpu_compute(neighbor->ago, inum, nall, atom->x, atom->type, - amtype, amgroup, rpole, uind, uinp, - ilist, numneigh, firstneigh, eflag, vflag, eflag_atom, - vflag_atom, host_start, cpu_time, success, atom->q, - atom->nlocal, domain->boxlo, domain->prd, &tep_pinned); + domain->bbox(domain->sublo_lamda,domain->subhi_lamda,sublo,subhi); } + inum = atom->nlocal; + + firstneigh = amoeba_gpu_compute_polar_real(neighbor->ago, inum, nall, atom->x, + atom->type, amtype, amgroup, + rpole, uind, uinp, sublo, subhi, + atom->tag, atom->nspecial, atom->special, + atom->nspecial15, atom->special15, + eflag, vflag, eflag_atom, vflag_atom, + host_start, &ilist, &numneigh, cpu_time, + success, atom->q, domain->boxlo, + domain->prd, &tep_pinned); + + if (!success) error->one(FLERR,"Insufficient memory on accelerator"); @@ -248,6 +245,7 @@ void PairAmoebaGPU::init_style() } // select the cutoff (off2) for neighbor list builds (the polar term for now) + // NOTE: induce and polar terms are using the same flags here if (use_ewald) choose(POLAR_LONG); else choose(POLAR); @@ -268,7 +266,7 @@ void PairAmoebaGPU::init_style() double felec = 0.5 * electric / am_dielectric; - int success = amoeba_gpu_init(atom->ntypes+1, max_amtype, pdamp, thole, + int success = amoeba_gpu_init(atom->ntypes+1, max_amtype, pdamp, thole, dirdamp, special_polar_wscale, special_polar_piscale, special_polar_pscale, atom->nlocal, atom->nlocal+atom->nghost, mnf, maxspecial, @@ -286,6 +284,199 @@ void PairAmoebaGPU::init_style() tep_single = true; } +/* ---------------------------------------------------------------------- + udirect2b = Ewald real direct field via list + udirect2b computes the real space contribution of the permanent + atomic multipole moments to the field via a neighbor list +------------------------------------------------------------------------- */ + +void PairAmoebaGPU::udirect2b(double **field, double **fieldp) +{ + int eflag=1, vflag=1; + int nall = atom->nlocal + atom->nghost; + int inum, host_start; + + bool success = true; + int *ilist, *numneigh, **firstneigh; + + double sublo[3],subhi[3]; + if (domain->triclinic == 0) { + sublo[0] = domain->sublo[0]; + sublo[1] = domain->sublo[1]; + sublo[2] = domain->sublo[2]; + subhi[0] = domain->subhi[0]; + subhi[1] = domain->subhi[1]; + subhi[2] = domain->subhi[2]; + } else { + domain->bbox(domain->sublo_lamda,domain->subhi_lamda,sublo,subhi); + } + inum = atom->nlocal; + + firstneigh = amoeba_gpu_compute_udirect2b(neighbor->ago, inum, nall, atom->x, + atom->type, amtype, amgroup, rpole, sublo, + subhi, atom->tag, atom->nspecial, atom->special, + atom->nspecial15, atom->special15, + eflag, vflag, eflag_atom, vflag_atom, + host_start, &ilist, &numneigh, cpu_time, + success, atom->q, domain->boxlo, + domain->prd, &fieldp_pinned); + if (!success) + error->one(FLERR,"Insufficient memory on accelerator"); + + // rebuild dipole-dipole pair list and store pairwise dipole matrices + // done one atom at a time in real-space double loop over atoms & neighs + + udirect2b_cpu(); +} + +/* ---------------------------------------------------------------------- + udirect2b = Ewald real direct field via list + udirect2b computes the real space contribution of the permanent + atomic multipole moments to the field via a neighbor list +------------------------------------------------------------------------- */ + +void PairAmoebaGPU::udirect2b_cpu() +{ + int i,j,k,m,n,ii,jj,kk,kkk,jextra,ndip,itype,jtype,igroup,jgroup; + double xr,yr,zr,r,r2; + double rr1,rr2,rr3,rr5; + double bfac,exp2a; + double ralpha,aefac; + double aesq2,aesq2n; + double pdi,pti,ddi; + double pgamma; + double damp,expdamp; + double scale3,scale5; + double scale7,scalek; + double bn[4],bcn[3]; + double factor_dscale,factor_pscale,factor_uscale,factor_wscale; + + int inum,jnum; + int *ilist,*jlist,*numneigh,**firstneigh; + + // launching the kernel to compute field and fieldp + + // amoeba_gpu_compute_field(...); + + double **x = atom->x; + + // neigh list + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // NOTE: doesn't this have a problem if aewald is tiny ?? + + aesq2 = 2.0 * aewald * aewald; + aesq2n = 0.0; + if (aewald > 0.0) aesq2n = 1.0 / (MY_PIS*aewald); + + // rebuild dipole-dipole pair list and store pairwise dipole matrices + // done one atom at a time in real-space double loop over atoms & neighs + + int *neighptr; + double *tdipdip; + + // compute the real space portion of the Ewald summation + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + itype = amtype[i]; + igroup = amgroup[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + n = ndip = 0; + neighptr = ipage_dipole->vget(); + tdipdip = dpage_dipdip->vget(); + + pdi = pdamp[itype]; + pti = thole[itype]; + ddi = dirdamp[itype]; + + // evaluate all sites within the cutoff distance + + for (jj = 0; jj < jnum; jj++) { + jextra = jlist[jj]; + j = jextra & NEIGHMASK15; + + xr = x[j][0] - x[i][0]; + yr = x[j][1] - x[i][1]; + zr = x[j][2] - x[i][2]; + r2 = xr*xr + yr* yr + zr*zr; + if (r2 > off2) continue; + + jtype = amtype[j]; + jgroup = amgroup[j]; + + factor_wscale = special_polar_wscale[sbmask15(jextra)]; + if (igroup == jgroup) { + factor_pscale = special_polar_piscale[sbmask15(jextra)]; + factor_dscale = polar_dscale; + factor_uscale = polar_uscale; + } else { + factor_pscale = special_polar_pscale[sbmask15(jextra)]; + factor_dscale = factor_uscale = 1.0; + } + + r = sqrt(r2); + rr1 = 1.0 / r; + rr2 = rr1 * rr1; + rr3 = rr2 * rr1; + rr5 = 3.0 * rr2 * rr3; + + // calculate the real space Ewald error function terms + + ralpha = aewald * r; + bn[0] = erfc(ralpha) * rr1; + exp2a = exp(-ralpha*ralpha); + aefac = aesq2n; + for (m = 1; m <= 3; m++) { + bfac = m+m-1; + aefac = aesq2 * aefac; + bn[m] = (bfac*bn[m-1]+aefac*exp2a) * rr2; + } + + // find terms needed later to compute mutual polarization + + if (poltyp != DIRECT) { + scale3 = 1.0; + scale5 = 1.0; + damp = pdi * pdamp[jtype]; + if (damp != 0.0) { + pgamma = MIN(pti,thole[jtype]); + damp = pgamma * pow(r/damp,3.0); + if (damp < 50.0) { + expdamp = exp(-damp); + scale3 = 1.0 - expdamp; + scale5 = 1.0 - expdamp*(1.0+damp); + } + } + scalek = factor_uscale; + bcn[0] = bn[1] - (1.0-scalek*scale3)*rr3; + bcn[1] = bn[2] - (1.0-scalek*scale5)*rr5; + + neighptr[n++] = j; + tdipdip[ndip++] = -bcn[0] + bcn[1]*xr*xr; + tdipdip[ndip++] = bcn[1]*xr*yr; + tdipdip[ndip++] = bcn[1]*xr*zr; + tdipdip[ndip++] = -bcn[0] + bcn[1]*yr*yr; + tdipdip[ndip++] = bcn[1]*yr*zr; + tdipdip[ndip++] = -bcn[0] + bcn[1]*zr*zr; + } + + } // jj + + firstneigh_dipole[i] = neighptr; + firstneigh_dipdip[i] = tdipdip; + numneigh_dipole[i] = n; + ipage_dipole->vgot(n); + dpage_dipdip->vgot(ndip); + } +} + /* ---------------------------------------------------------------------- */ double PairAmoebaGPU::memory_usage() diff --git a/src/GPU/pair_amoeba_gpu.h b/src/GPU/pair_amoeba_gpu.h index 4d29bfaf34..e5d4aab176 100644 --- a/src/GPU/pair_amoeba_gpu.h +++ b/src/GPU/pair_amoeba_gpu.h @@ -34,13 +34,17 @@ class PairAmoebaGPU : public PairAmoeba { enum { GPU_FORCE, GPU_NEIGH, GPU_HYB_NEIGH }; virtual void polar_real(); + virtual void udirect2b(double **, double **); private: int gpu_mode; double cpu_time; void *tep_pinned; + void *fieldp_pinned; bool tep_single; + void udirect2b_cpu(); + template void compute_force_from_tep(const numtyp*); };