ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added and renamed API to make room for additional kernels (udirect2b only computes the field and fieldp, not accumulating forces, energies, nor virials)

Browse Source
This commit is contained in:
Trung Nguyen
2021-09-01 14:37:11 -05:00
parent 07b60827c4
commit 785a794d39
9 changed files with 448 additions and 111 deletions
Download Patch File Download Diff File Expand all files Collapse all files

35
lib/gpu/lal_amoeba.cpp
View File

@ -125,10 +125,10 @@ double AmoebaT::host_memory_usage() const {
}
// ---------------------------------------------------------------------------
// Calculate energies, forces, and torques
// Calculate the polar real-space term, returning tep
// ---------------------------------------------------------------------------
template <class numtyp, class acctyp>
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<int>(ceil(static_cast<double>(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 <class numtyp, class acctyp>
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<int>(ceil(static_cast<double>(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<PRECISION,ACC_PRECISION>;
}

20
lib/gpu/lal_amoeba.cu
View File

@ -715,11 +715,6 @@ __kernel void k_amoeba_udirect2b(const __global numtyp4 *restrict x_,
//numtyp4 xi__;
if (ii<inum) {
int itype,igroup;
numtyp bn[4],bcn[3];
numtyp fid[3],fip[3];
numtyp ci,uix,uiy,uiz,uixp,uiyp,uizp;
int numj, nbor, nbor_end;
nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
n_stride,nbor_end,nbor);
@ -728,6 +723,11 @@ __kernel void k_amoeba_udirect2b(const __global numtyp4 *restrict x_,
//numtyp qtmp; fetch(qtmp,i,q_tex);
//int itype=ix.w;
int itype,igroup;
numtyp bn[4],bcn[3];
numtyp fid[3],fip[3];
numtyp ci,uix,uiy,uiz,uixp,uiyp,uizp;
ci = polar1[i].x; // rpole[i][0];
dix = polar1[i].y; // rpole[i][1];
diy = polar1[i].z; // rpole[i][2];
@ -748,9 +748,9 @@ __kernel void k_amoeba_udirect2b(const __global numtyp4 *restrict x_,
numtyp pti = damping[itype].y;
numtyp ddi = damping[itype].z;
numtyp alsq2 = (numtyp)2.0 * aewald*aewald;
numtyp alsq2n = (numtyp)0.0;
if (aewald > (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<nbor_end; nbor+=n_stride) {
@ -841,7 +841,7 @@ __kernel void k_amoeba_udirect2b(const __global numtyp4 *restrict x_,
if (pgamma != (numtyp)0.0) {
damp = pgamma * ucl_powr(r/damp,(numtyp)1.5);
if (damp < (numtyp)50.0) {
expdamp = ucl_exp(-damp) ;
numtyp expdamp = ucl_exp(-damp) ;
scale3 = (numtyp)1.0 - expdamp ;
scale5 = (numtyp)1.0 - expdamp*((numtyp)1.0+(numtyp)0.5*damp);
scale7 = (numtyp)1.0 - expdamp*((numtyp)1.0+(numtyp)0.65*damp + (numtyp)0.15*damp*damp);
@ -850,7 +850,7 @@ __kernel void k_amoeba_udirect2b(const __global numtyp4 *restrict x_,
pgamma = MIN(pti,damping[jtype].y); // thole[jtype]
damp = pgamma * ucl_powr(r/damp,3.0);
if (damp < (numtyp)50.0) {
expdamp = ucl_exp(-damp);
numtyp expdamp = ucl_exp(-damp);
scale3 = (numtyp)1.0 - expdamp;
scale5 = (numtyp)1.0 - expdamp*((numtyp)1.0+damp);
scale7 = (numtyp)1.0 - expdamp*((numtyp)1.0+damp + (numtyp)0.6*damp*damp);

5
lib/gpu/lal_amoeba.h
View File

@ -78,9 +78,10 @@ class Amoeba : public BaseAmoeba<numtyp, acctyp> {
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);
};
}

27
lib/gpu/lal_amoeba_ext.cpp
View File

@ -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,
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,
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 **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() {

129
lib/gpu/lal_base_amoeba.cpp
View File

@ -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<int>(static_cast<double>(ef_nall)*1.10);
_tep.alloc(_max_tep_size*4,*(this->ucl_device),UCL_READ_WRITE,UCL_READ_WRITE);
_max_alloc_size=static_cast<int>(static_cast<double>(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<int>(static_cast<double>(nall)*1.10);
_tep.resize(_max_tep_size*4);
if (nall>_max_alloc_size) {
_max_alloc_size=static_cast<int>(static_cast<double>(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 <class numtyp, class acctyp>
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<int>(static_cast<double>(nall)*1.10);
_tep.resize(_max_tep_size*4);
if (nall>_max_alloc_size) {
_max_alloc_size=static_cast<int>(static_cast<double>(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 <class numtyp, class acctyp>
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<int>(static_cast<double>(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 <class numtyp, class acctyp>
double BaseAmoebaT::host_memory_usage_atomic() const {
return device->atom.host_memory_usage()+nbor->host_memory_usage()+
@ -446,6 +543,7 @@ void BaseAmoebaT::cast_extra_data(int* amtype, int* amgroup, double** rpole,
pextra[idx+3] = (numtyp)amgroup[i];
}
if (uind) {
n += nstride*_nall;
for (int i = 0; i < _nall; i++) {
int idx = n+i*nstride;
@ -453,7 +551,9 @@ void BaseAmoebaT::cast_extra_data(int* amtype, int* amgroup, double** rpole,
pextra[idx+1] = uind[i][1];
pextra[idx+2] = uind[i][2];
}
}
if (uinp) {
n += nstride*_nall;
for (int i = 0; i < _nall; i++) {
int idx = n+i*nstride;
@ -462,6 +562,7 @@ void BaseAmoebaT::cast_extra_data(int* amtype, int* amgroup, double** rpole,
pextra[idx+2] = uinp[i][2];
}
}
}
template <class numtyp, class acctyp>
void BaseAmoebaT::compile_kernels(UCL_Device &dev, const void *pair_str,

24
lib/gpu/lal_base_amoeba.h
View File

@ -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<numtyp,numtyp> _tep;
int _max_tep_size;
UCL_Vector<numtyp,numtyp> _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;
};
}

2
src/AMOEBA/pair_amoeba.h
View File

@ -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 *);

241
src/GPU/pair_amoeba_gpu.cpp
View File

@ -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 <cmath>
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,7 +110,7 @@ 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];
@ -116,7 +124,7 @@ void PairAmoebaGPU::polar_real()
}
inum = atom->nlocal;
firstneigh = amoeba_gpu_compute_n(neighbor->ago, inum, nall, atom->x,
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,
@ -126,18 +134,7 @@ void PairAmoebaGPU::polar_real()
success, atom->q, domain->boxlo,
domain->prd, &tep_pinned);
} 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);
}
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()

4
src/GPU/pair_amoeba_gpu.h
View File

@ -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<class numtyp>
void compute_force_from_tep(const numtyp*);
};