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

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13907 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2015-08-19 15:16:24 +00:00
parent 21f8abda24
commit 071ded95fc
10 changed files with 2048 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

159
lib/gpu/lal_lj_cubic.cpp Normal file
View File

@ -0,0 +1,159 @@
/***************************************************************************
lj_cubic.cpp
-------------------
Trung Dac Nguyen
Class for acceleration of the lj/cubic pair style.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
__________________________________________________________________________
begin :
email : ndactrung@gmail.com
***************************************************************************/
#if defined(USE_OPENCL)
#include "lj_cubic_cl.h"
#elif defined(USE_CUDART)
const char *lj_cubic=0;
#else
#include "lj_cubic_cubin.h"
#endif
#include "lal_lj_cubic.h"
#include <cassert>
using namespace LAMMPS_AL;
#define LJCubicT LJCubic<numtyp, acctyp>
extern Device<PRECISION,ACC_PRECISION> device;
template <class numtyp, class acctyp>
LJCubicT::LJCubic() : BaseAtomic<numtyp,acctyp>(), _allocated(false) {
}
template <class numtyp, class acctyp>
LJCubicT::~LJCubic() {
clear();
}
template <class numtyp, class acctyp>
int LJCubicT::bytes_per_atom(const int max_nbors) const {
return this->bytes_per_atom_atomic(max_nbors);
}
template <class numtyp, class acctyp>
int LJCubicT::init(const int ntypes,
double **host_cutsq, double **host_cut_inner_sq,
double **host_cut_inner, double **host_sigma,
double **host_epsilon, double **host_lj1,
double **host_lj2, double **host_lj3, double **host_lj4,
double *host_special_lj, const int nlocal,
const int nall, const int max_nbors,
const int maxspecial, const double cell_size,
const double gpu_split, FILE *_screen) {
int success;
success=this->init_atomic(nlocal,nall,max_nbors,maxspecial,cell_size,gpu_split,
_screen,lj_cubic,"k_lj_cubic");
if (success!=0)
return success;
// If atom type constants fit in shared memory use fast kernel
int lj_types=ntypes;
shared_types=false;
int max_shared_types=this->device->max_shared_types();
if (lj_types<=max_shared_types && this->_block_size>=max_shared_types) {
lj_types=max_shared_types;
shared_types=true;
}
_lj_types=lj_types;
// Allocate a host write buffer for data initialization
UCL_H_Vec<numtyp> host_write(lj_types*lj_types*32,*(this->ucl_device),
UCL_WRITE_ONLY);
for (int i=0; i<lj_types*lj_types; i++)
host_write[i]=0.0;
lj1.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
this->atom->type_pack4(ntypes,lj_types,lj1,host_write,host_lj1,host_lj2,
host_cutsq);
lj2.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
this->atom->type_pack4(ntypes,lj_types,lj2,host_write,host_cut_inner_sq,
host_cut_inner,host_sigma,host_epsilon);
lj3.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
this->atom->type_pack2(ntypes,lj_types,lj3,host_write,host_lj3,host_lj4);
UCL_H_Vec<double> dview;
sp_lj.alloc(4,*(this->ucl_device),UCL_READ_ONLY);
dview.view(host_special_lj,4,*(this->ucl_device));
ucl_copy(sp_lj,dview,false);
_allocated=true;
this->_max_bytes=lj1.row_bytes()+lj2.row_bytes()+lj3.row_bytes()
+sp_lj.row_bytes();
return 0;
}
template <class numtyp, class acctyp>
void LJCubicT::clear() {
if (!_allocated)
return;
_allocated=false;
lj1.clear();
lj2.clear();
lj3.clear();
sp_lj.clear();
this->clear_atomic();
}
template <class numtyp, class acctyp>
double LJCubicT::host_memory_usage() const {
return this->host_memory_usage_atomic()+sizeof(LJCubic<numtyp,acctyp>);
}
// ---------------------------------------------------------------------------
// Calculate energies, forces, and torques
// ---------------------------------------------------------------------------
template <class numtyp, class acctyp>
void LJCubicT::loop(const bool _eflag, const bool _vflag) {
// Compute the block size and grid size to keep all cores busy
const int BX=this->block_size();
int eflag, vflag;
if (_eflag)
eflag=1;
else
eflag=0;
if (_vflag)
vflag=1;
else
vflag=0;
int GX=static_cast<int>(ceil(static_cast<double>(this->ans->inum())/
(BX/this->_threads_per_atom)));
int ainum=this->ans->inum();
int nbor_pitch=this->nbor->nbor_pitch();
this->time_pair.start();
if (shared_types) {
this->k_pair_fast.set_size(GX,BX);
this->k_pair_fast.run(&this->atom->x, &lj1, &lj2, &lj3, &sp_lj,
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
&this->ans->force, &this->ans->engv, &eflag,
&vflag, &ainum, &nbor_pitch,
&this->_threads_per_atom);
} else {
this->k_pair.set_size(GX,BX);
this->k_pair.run(&this->atom->x, &lj1, &lj2, &lj3, &_lj_types, &sp_lj,
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
&this->ans->force, &this->ans->engv, &eflag, &vflag,
&ainum, &nbor_pitch, &this->_threads_per_atom);
}
this->time_pair.stop();
}
template class LJCubic<PRECISION,ACC_PRECISION>;

81
lib/gpu/lal_lj_cubic.h Normal file
View File

@ -0,0 +1,81 @@
/***************************************************************************
lj_cubic.h
-------------------
Trung Dac Nguyen
Class for acceleration of the lj/cubic pair style.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
__________________________________________________________________________
begin :
email : ndactrung@gmail.com
***************************************************************************/
#ifndef LAL_LJ_H
#define LAL_LJ_H
#include "lal_base_atomic.h"
namespace LAMMPS_AL {
template <class numtyp, class acctyp>
class LJCubic : public BaseAtomic<numtyp, acctyp> {
public:
LJCubic();
~LJCubic();
/// Clear any previous data and set up for a new LAMMPS run
/** \param max_nbors initial number of rows in the neighbor matrix
* \param cell_size cutoff + skin
* \param gpu_split fraction of particles handled by device
*
* Returns:
* - 0 if successfull
* - -1 if fix gpu not found
* - -3 if there is an out of memory error
* - -4 if the GPU library was not compiled for GPU
* - -5 Double precision is not supported on card **/
int init(const int ntypes, double **host_cutsq, double **host_cut_inner_sq,
double **host_cut_inner, double **host_sigma, double **host_epsilon,
double **host_lj1, double **host_lj2, double **host_lj3,
double **host_lj4, double *host_special_lj, const int nlocal,
const int nall, const int max_nbors, const int maxspecial,
const double cell_size, const double gpu_split, FILE *screen);
/// Clear all host and device data
/** \note This is called at the beginning of the init() routine **/
void clear();
/// Returns memory usage on device per atom
int bytes_per_atom(const int max_nbors) const;
/// Total host memory used by library for pair style
double host_memory_usage() const;
// --------------------------- TYPE DATA --------------------------
/// lj1.x = lj1, lj1.y = lj2, lj1.z = cutsq
UCL_D_Vec<numtyp4> lj1;
/// lj2.x = cut_inner_sq, lj2.y = cut_inner, lj2.z = sigma, lj2.w = epsilon
UCL_D_Vec<numtyp4> lj2;
/// lj3.x = lj3, lj3.y = lj4
UCL_D_Vec<numtyp2> lj3;
/// Special LJ values
UCL_D_Vec<numtyp> sp_lj;
/// If atom type constants fit in shared memory, use fast kernels
bool shared_types;
/// Number of atom types
int _lj_types;
private:
bool _allocated;
void loop(const bool _eflag, const bool _vflag);
};
}
#endif

124
lib/gpu/lal_lj_cubic_ext.cpp Normal file
View File

@ -0,0 +1,124 @@
/***************************************************************************
lj_cubic_ext.cpp
-------------------
Trung Dac Nguyen
Functions for LAMMPS access to lj/cubic acceleration routines.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
__________________________________________________________________________
begin :
email : ndactrung@gmail.com
***************************************************************************/
#include <iostream>
#include <cassert>
#include <math.h>
#include "lal_lj_cubic.h"
using namespace std;
using namespace LAMMPS_AL;
static LJCubic<PRECISION,ACC_PRECISION> LJCubicLMF;
// ---------------------------------------------------------------------------
// Allocate memory on host and device and copy constants to device
// ---------------------------------------------------------------------------
int ljcb_gpu_init(const int ntypes, double **cutsq, double **cut_inner_sq,
double **cut_inner, double **sigma, double **epsilon,
double **host_lj1, double **host_lj2, double **host_lj3,
double **host_lj4, double *special_lj,
const int inum, const int nall, const int max_nbors,
const int maxspecial, const double cell_size,
int &gpu_mode, FILE *screen) {
LJCubicLMF.clear();
gpu_mode=LJCubicLMF.device->gpu_mode();
double gpu_split=LJCubicLMF.device->particle_split();
int first_gpu=LJCubicLMF.device->first_device();
int last_gpu=LJCubicLMF.device->last_device();
int world_me=LJCubicLMF.device->world_me();
int gpu_rank=LJCubicLMF.device->gpu_rank();
int procs_per_gpu=LJCubicLMF.device->procs_per_gpu();
LJCubicLMF.device->init_message(screen,"lj/cubic",first_gpu,last_gpu);
bool message=false;
if (LJCubicLMF.device->replica_me()==0 && screen)
message=true;
if (message) {
fprintf(screen,"Initializing Device and compiling on process 0...");
fflush(screen);
}
int init_ok=0;
if (world_me==0)
init_ok=LJCubicLMF.init(ntypes, cutsq, cut_inner_sq, cut_inner, sigma,
epsilon, host_lj1, host_lj2, host_lj3, host_lj4,
special_lj, inum, nall, 300, maxspecial,
cell_size, gpu_split, screen);
LJCubicLMF.device->world_barrier();
if (message)
fprintf(screen,"Done.\n");
for (int i=0; i<procs_per_gpu; i++) {
if (message) {
if (last_gpu-first_gpu==0)
fprintf(screen,"Initializing Device %d on core %d...",first_gpu,i);
else
fprintf(screen,"Initializing Devices %d-%d on core %d...",first_gpu,
last_gpu,i);
fflush(screen);
}
if (gpu_rank==i && world_me!=0)
init_ok=LJCubicLMF.init(ntypes, cutsq, cut_inner_sq, cut_inner, sigma,
epsilon, host_lj1, host_lj2, host_lj3, host_lj4,
special_lj, inum, nall, 300, maxspecial,
cell_size, gpu_split, screen);
LJCubicLMF.device->gpu_barrier();
if (message)
fprintf(screen,"Done.\n");
}
if (message)
fprintf(screen,"\n");
if (init_ok==0)
LJCubicLMF.estimate_gpu_overhead();
return init_ok;
}
void ljcb_gpu_clear() {
LJCubicLMF.clear();
}
int ** ljcb_gpu_compute_n(const int ago, const int inum_full,
const int nall, double **host_x, int *host_type,
double *sublo, double *subhi, tagint *tag, int **nspecial,
tagint **special, const bool eflag, const bool vflag,
const bool eatom, const bool vatom, int &host_start,
int **ilist, int **jnum, const double cpu_time,
bool &success) {
return LJCubicLMF.compute(ago, inum_full, nall, host_x, host_type, sublo,
subhi, tag, nspecial, special, eflag, vflag, eatom,
vatom, host_start, ilist, jnum, cpu_time, success);
}
void ljcb_gpu_compute(const int ago, const int inum_full, const int nall,
double **host_x, int *host_type, int *ilist, int *numj,
int **firstneigh, const bool eflag, const bool vflag,
const bool eatom, const bool vatom, int &host_start,
const double cpu_time, bool &success) {
LJCubicLMF.compute(ago,inum_full,nall,host_x,host_type,ilist,numj,
firstneigh,eflag,vflag,eatom,vatom,host_start,cpu_time,success);
}
double ljcb_gpu_bytes() {
return LJCubicLMF.host_memory_usage();
}

451
lib/gpu/lal_tersoff.cpp Normal file
View File

@ -0,0 +1,451 @@
/***************************************************************************
tersoff.cpp
-------------------
Trung Dac Nguyen
Class for acceleration of the tersoff pair style.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
__________________________________________________________________________
begin : Thu April 17, 2014
email : ndactrung@gmail.com
***************************************************************************/
#if defined(USE_OPENCL)
#include "tersoff_cl.h"
#elif defined(USE_CUDART)
const char *tersoff=0;
#else
#include "tersoff_cubin.h"
#endif
#include "lal_tersoff.h"
#include <cassert>
using namespace LAMMPS_AL;
#define TersoffT Tersoff<numtyp, acctyp>
extern Device<PRECISION,ACC_PRECISION> device;
template <class numtyp, class acctyp>
TersoffT::Tersoff() : BaseThree<numtyp,acctyp>(), _allocated(false) {
}
template <class numtyp, class acctyp>
TersoffT::~Tersoff() {
clear();
}
template <class numtyp, class acctyp>
int TersoffT::bytes_per_atom(const int max_nbors) const {
return this->bytes_per_atom_atomic(max_nbors);
}
template <class numtyp, class acctyp>
int TersoffT::init(const int ntypes, const int nlocal, const int nall, const int max_nbors,
const double cell_size, const double gpu_split, FILE *_screen,
int* host_map, const int nelements, int*** host_elem2param, const int nparams,
const double* lam1, const double* lam2, const double* lam3,const double* powermint,
const double* biga, const double* bigb, const double* bigr, const double* bigd,
const double* c1, const double* c2, const double* c3, const double* c4,
const double* c, const double* d, const double* h, const double* gamma,
const double* beta, const double* powern, const double* host_cutsq)
{
int success;
success=this->init_three(nlocal,nall,max_nbors,0,cell_size,gpu_split,
_screen,tersoff,"k_tersoff_repulsive",
"k_tersoff_three_center", "k_tersoff_three_end");
if (success!=0)
return success;
int ef_nall=nall;
if (ef_nall==0)
ef_nall=2000;
_max_zij_size=static_cast<int>(static_cast<double>(ef_nall)*1.10);
_zetaij.alloc(_max_zij_size*max_nbors,*(this->ucl_device),UCL_READ_WRITE);
zeta_tex.get_texture(*(this->pair_program),"zeta_tex");
zeta_tex.bind_float(_zetaij,1);
k_zeta.set_function(*(this->pair_program),"k_tersoff_zeta");
// If atom type constants fit in shared memory use fast kernel
int lj_types=ntypes;
shared_types=false;
int max_shared_types=this->device->max_shared_types();
if (lj_types<=max_shared_types && this->_block_size>=max_shared_types) {
lj_types=max_shared_types;
shared_types=true;
}
_lj_types=lj_types;
_nparams = nparams;
_nelements = nelements;
UCL_H_Vec<numtyp4> dview(nparams,*(this->ucl_device),
UCL_WRITE_ONLY);
for (int i=0; i<nparams; i++) {
dview[i].x=(numtyp)0;
dview[i].y=(numtyp)0;
dview[i].z=(numtyp)0;
dview[i].w=(numtyp)0;
}
// pack coefficients into arrays
ts1.alloc(nparams,*(this->ucl_device),UCL_READ_ONLY);
for (int i=0; i<nparams; i++) {
dview[i].x=static_cast<numtyp>(lam1[i]);
dview[i].y=static_cast<numtyp>(lam2[i]);
dview[i].z=static_cast<numtyp>(lam3[i]);
dview[i].w=static_cast<numtyp>(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<nparams; i++) {
dview[i].x=static_cast<numtyp>(biga[i]);
dview[i].y=static_cast<numtyp>(bigb[i]);
dview[i].z=static_cast<numtyp>(bigr[i]);
dview[i].w=static_cast<numtyp>(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<nparams; i++) {
dview[i].x=static_cast<numtyp>(c1[i]);
dview[i].y=static_cast<numtyp>(c2[i]);
dview[i].z=static_cast<numtyp>(c3[i]);
dview[i].w=static_cast<numtyp>(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<nparams; i++) {
dview[i].x=static_cast<numtyp>(c[i]);
dview[i].y=static_cast<numtyp>(d[i]);
dview[i].z=static_cast<numtyp>(h[i]);
dview[i].w=static_cast<numtyp>(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<nparams; i++) {
dview[i].x=static_cast<numtyp>(beta[i]);
dview[i].y=static_cast<numtyp>(powern[i]);
dview[i].z=(numtyp)0;
dview[i].w=(numtyp)0;
}
ucl_copy(ts5,dview,false);
ts5_tex.get_texture(*(this->pair_program),"ts5_tex");
ts5_tex.bind_float(ts5,4);
UCL_H_Vec<numtyp> cutsq_view(nparams,*(this->ucl_device),
UCL_WRITE_ONLY);
for (int i=0; i<nparams; i++)
cutsq_view[i]=static_cast<numtyp>(host_cutsq[i]);
cutsq.alloc(nparams,*(this->ucl_device),UCL_READ_ONLY);
ucl_copy(cutsq,cutsq_view,false);
UCL_H_Vec<int> 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<int> dview_map(lj_types, *(this->ucl_device), UCL_WRITE_ONLY);
for (int i = 0; i < lj_types; i++)
dview_map[i] = host_map[i];
map.alloc(lj_types,*(this->ucl_device), UCL_READ_ONLY);
ucl_copy(map,dview_map,false);
_allocated=true;
this->_max_bytes=ts1.row_bytes()+ts2.row_bytes()+ts3.row_bytes()+
ts4.row_bytes()+ts5.row_bytes()+cutsq.row_bytes()+
map.row_bytes()+elem2param.row_bytes()+_zetaij.row_bytes();
return 0;
}
template <class numtyp, class acctyp>
void TersoffT::clear() {
if (!_allocated)
return;
_allocated=false;
ts1.clear();
ts2.clear();
ts3.clear();
ts4.clear();
ts5.clear();
cutsq.clear();
map.clear();
elem2param.clear();
_zetaij.clear();
k_zeta.clear();
this->clear_atomic();
}
template <class numtyp, class acctyp>
double TersoffT::host_memory_usage() const {
return this->host_memory_usage_atomic()+sizeof(Tersoff<numtyp,acctyp>);
}
#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 <class numtyp, class acctyp>
void TersoffT::compute(const int f_ago, const int nlocal, const int nall,
const int nlist, double **host_x, int *host_type,
int *ilist, int *numj, int **firstneigh,
const bool eflag, const bool vflag, const bool eatom,
const bool vatom, int &host_start,
const double cpu_time, bool &success) {
this->acc_timers();
if (nlist==0) {
host_start=0;
// Make sure textures are correct if realloc by a different hybrid style
this->resize_atom(0,nall,success);
this->zero_timers();
return;
}
int ago=this->hd_balancer.ago_first(f_ago);
int inum=this->hd_balancer.balance(ago,nlocal,cpu_time);
this->ans->inum(inum);
#ifdef THREE_CONCURRENT
this->ans2->inum(inum);
#endif
host_start=inum;
if (ago==0) {
this->reset_nbors(nall, inum, nlist, ilist, numj, firstneigh, success);
if (!success)
return;
}
this->atom->cast_x_data(host_x,host_type);
this->hd_balancer.start_timer();
this->atom->add_x_data(host_x,host_type);
// re-allocate zetaij if necessary
if (nall>_max_zij_size) {
_max_nbors = this->nbor->max_nbor_loop(nlist,numj,ilist);
_max_zij_size=static_cast<int>(static_cast<double>(nall)*1.10);
_zetaij.resize(_max_nbors*_max_zij_size);
zeta_tex.bind_float(_zetaij,1);
}
int ainum=nall;
int nbor_pitch=this->nbor->nbor_pitch();
int BX=this->block_pair();
int GX=static_cast<int>(ceil(static_cast<double>(ainum)/
(BX/this->_threads_per_atom)));
this->k_zeta.set_size(GX,BX);
this->k_zeta.run(&this->atom->x, &ts1, &ts2, &ts4, &cutsq,
&map, &elem2param, &_nelements, &_nparams, &_zetaij,
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
&nall, &ainum, &nbor_pitch, &this->_threads_per_atom);
int evatom=0;
if (eatom || vatom)
evatom=1;
#ifdef THREE_CONCURRENT
this->ucl_device->sync();
#endif
loop(eflag,vflag,evatom);
this->ans->copy_answers(eflag,vflag,eatom,vatom,ilist);
this->device->add_ans_object(this->ans);
#ifdef THREE_CONCURRENT
this->ans2->copy_answers(eflag,vflag,eatom,vatom,ilist);
this->device->add_ans_object(this->ans2);
#endif
this->hd_balancer.stop_timer();
}
// ---------------------------------------------------------------------------
// Reneighbor on GPU if necessary and then compute forces, virials, energies
// ---------------------------------------------------------------------------
template <class numtyp, class acctyp>
int ** TersoffT::compute(const int ago, const int inum_full,
const int nall, double **host_x, int *host_type,
double *sublo, double *subhi, tagint *tag,
int **nspecial, tagint **special, const bool eflag,
const bool vflag, const bool eatom,
const bool vatom, int &host_start,
int **ilist, int **jnum,
const double cpu_time, bool &success) {
this->acc_timers();
if (inum_full==0) {
host_start=0;
// Make sure textures are correct if realloc by a different hybrid style
this->resize_atom(0,nall,success);
this->zero_timers();
return NULL;
}
this->hd_balancer.balance(cpu_time);
int inum=this->hd_balancer.get_gpu_count(ago,inum_full);
this->ans->inum(inum);
#ifdef THREE_CONCURRENT
this->ans2->inum(inum);
#endif
host_start=inum;
// Build neighbor list on GPU if necessary
if (ago==0) {
_max_nbors = this->build_nbor_list(inum, inum_full-inum, nall, host_x, host_type,
sublo, subhi, tag, nspecial, special, success);
if (!success)
return NULL;
this->hd_balancer.start_timer();
} else {
this->atom->cast_x_data(host_x,host_type);
this->hd_balancer.start_timer();
this->atom->add_x_data(host_x,host_type);
}
*ilist=this->nbor->host_ilist.begin();
*jnum=this->nbor->host_acc.begin();
// re-allocate zetaij if necessary
if (nall>_max_zij_size) {
_max_zij_size=static_cast<int>(static_cast<double>(nall)*1.10);
_zetaij.resize(_max_nbors*_max_zij_size);
zeta_tex.bind_float(_zetaij,1);
}
int ainum=nall;
int nbor_pitch=this->nbor->nbor_pitch();
int BX=this->block_pair();
int GX=static_cast<int>(ceil(static_cast<double>(ainum)/
(BX/this->_threads_per_atom)));
this->k_zeta.set_size(GX,BX);
this->k_zeta.run(&this->atom->x, &ts1, &ts2, &ts4, &cutsq, &map,
&elem2param, &_nelements, &_nparams, &_zetaij,
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
&nall, &ainum, &nbor_pitch, &this->_threads_per_atom);
int evatom=0;
if (eatom || vatom)
evatom=1;
#ifdef THREE_CONCURRENT
this->ucl_device->sync();
#endif
loop(eflag,vflag,evatom);
this->ans->copy_answers(eflag,vflag,eatom,vatom);
this->device->add_ans_object(this->ans);
#ifdef THREE_CONCURRENT
this->ans2->copy_answers(eflag,vflag,eatom,vatom);
this->device->add_ans_object(this->ans2);
#endif
this->hd_balancer.stop_timer();
return this->nbor->host_jlist.begin()-host_start;
}
// ---------------------------------------------------------------------------
// Calculate energies, forces, and torques
// ---------------------------------------------------------------------------
template <class numtyp, class acctyp>
void TersoffT::loop(const bool _eflag, const bool _vflag, const int evatom) {
// Compute the block size and grid size to keep all cores busy
int BX=this->block_pair();
int eflag, vflag;
if (_eflag)
eflag=1;
else
eflag=0;
if (_vflag)
vflag=1;
else
vflag=0;
int ainum=this->ans->inum();
int nbor_pitch=this->nbor->nbor_pitch();
int GX=static_cast<int>(ceil(static_cast<double>(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<int>(ceil(static_cast<double>(this->ans->inum())/
(BX/(KTHREADS*JTHREADS))));
this->k_three_center.set_size(GX,BX);
this->k_three_center.run(&this->atom->x, &ts1, &ts2, &ts3, &ts4, &ts5, &cutsq,
&map, &elem2param, &_nelements, &_nparams, &_zetaij,
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
&this->ans->force, &this->ans->engv, &eflag, &vflag, &ainum,
&nbor_pitch, &this->_threads_per_atom, &evatom);
Answer<numtyp,acctyp> *end_ans;
#ifdef THREE_CONCURRENT
end_ans=this->ans2;
#else
end_ans=this->ans;
#endif
if (evatom!=0) {
this->k_three_end_vatom.set_size(GX,BX);
this->k_three_end_vatom.run(&this->atom->x, &ts1, &ts2, &ts3, &ts4, &ts5, &cutsq,
&map, &elem2param, &_nelements, &_nparams, &_zetaij,
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
&end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum,
&nbor_pitch, &this->_threads_per_atom);
} else {
this->k_three_end.set_size(GX,BX);
this->k_three_end.run(&this->atom->x, &ts1, &ts2, &ts3, &ts4, &ts5, &cutsq,
&map, &elem2param, &_nelements, &_nparams, &_zetaij,
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
&end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum,
&nbor_pitch, &this->_threads_per_atom);
}
this->time_pair.stop();
}
template class Tersoff<PRECISION,ACC_PRECISION>;

116
lib/gpu/lal_tersoff.h Normal file
View File

@ -0,0 +1,116 @@
/***************************************************************************
tersoff.h
-------------------
Trung Dac Nguyen
Class for acceleration of the tersoff pair style.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
__________________________________________________________________________
begin : Thu April 17, 2014
email : ndactrung@gmail.com
***************************************************************************/
#ifndef LAL_TERSOFF_H
#define LAL_TERSOFF_H
#include "lal_base_three.h"
namespace LAMMPS_AL {
template <class numtyp, class acctyp>
class Tersoff : public BaseThree<numtyp, acctyp> {
public:
Tersoff();
~Tersoff();
/// Clear any previous data and set up for a new LAMMPS run for generic systems
/** \param max_nbors initial number of rows in the neighbor matrix
* \param cell_size cutoff + skin
* \param gpu_split fraction of particles handled by device
*
* Returns:
* - 0 if successfull
* - -1 if fix gpu not found
* - -3 if there is an out of memory error
* - -4 if the GPU library was not compiled for GPU
* - -5 Double precision is not supported on card **/
int init(const int ntypes, const int nlocal, const int nall, const int max_nbors,
const double cell_size, const double gpu_split, FILE *screen,
int* host_map, const int nelements, int*** host_elem2param, const int nparams,
const double* lam1, const double* lam2, const double* lam3,
const double* powermint, const double* biga, const double* bigb,
const double* bigr, const double* bigd, const double* c1, const double* c2,
const double* c3, const double* c4, const double* c, const double* d,
const double* h, const double* gamma, const double* beta,
const double* powern, const double* cutsq);
/// Pair loop with host neighboring
void compute(const int f_ago, const int inum_full, const int nall,
const int nlist, double **host_x, int *host_type,
int *ilist, int *numj, int **firstneigh, const bool eflag,
const bool vflag, const bool eatom, const bool vatom,
int &host_start, const double cpu_time, bool &success);
/// Pair loop with device neighboring
int ** compute(const int ago, const int inum_full,
const int nall, double **host_x, int *host_type, double *sublo,
double *subhi, tagint *tag, int **nspecial,
tagint **special, const bool eflag, const bool vflag,
const bool eatom, const bool vatom, int &host_start,
int **ilist, int **numj, const double cpu_time, bool &success);
/// Clear all host and device data
/** \note This is called at the beginning of the init() routine **/
void clear();
/// Returns memory usage on device per atom
int bytes_per_atom(const int max_nbors) const;
/// Total host memory used by library for pair style
double host_memory_usage() const;
// --------------------------- TYPE DATA --------------------------
/// If atom type constants fit in shared memory, use fast kernels
bool shared_types;
/// Number of atom types
int _lj_types;
/// ts1.x = lam1, ts1.y = lam2, ts1.z = lam3, ts1.w = powermint
UCL_D_Vec<numtyp4> ts1;
/// ts2.x = biga, ts2.y = bigb, ts2.z = bigr, ts2.w = bigd
UCL_D_Vec<numtyp4> ts2;
/// ts3.x = c1, ts3.y = c2, ts3.z = c3, ts3.w = c4
UCL_D_Vec<numtyp4> ts3;
/// ts4.x = c, ts4.y = d, ts4.z = h, ts4.w = gamma
UCL_D_Vec<numtyp4> ts4;
/// ts5.x = beta, ts5.y = powern
UCL_D_Vec<numtyp4> ts5;
UCL_D_Vec<numtyp> cutsq;
UCL_D_Vec<int> elem2param;
UCL_D_Vec<int> map;
int _nparams,_nelements;
/// Per-atom arrays
UCL_D_Vec<numtyp> _zetaij;
UCL_Kernel k_zeta;
UCL_Texture ts1_tex, ts2_tex, ts3_tex, ts4_tex, ts5_tex, zeta_tex;
int _max_zij_size, _max_nbors;
private:
bool _allocated;
void loop(const bool _eflag, const bool _vflag, const int evatom);
};
}
#endif

139
lib/gpu/lal_tersoff_ext.cpp Normal file
View File

@ -0,0 +1,139 @@
/***************************************************************************
tersoff_ext.cpp
-------------------
Trung Dac Nguyen
Functions for LAMMPS access to tersoff acceleration routines.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
__________________________________________________________________________
begin : Thu April 17, 2014
email : ndactrung@gmail.com
***************************************************************************/
#include <iostream>
#include <cassert>
#include <math.h>
#include "lal_tersoff.h"
using namespace std;
using namespace LAMMPS_AL;
static Tersoff<PRECISION,ACC_PRECISION> TSMF;
// ---------------------------------------------------------------------------
// Allocate memory on host and device and copy constants to device
// ---------------------------------------------------------------------------
int tersoff_gpu_init(const int ntypes, const int inum, const int nall, const int max_nbors,
const double cell_size, int &gpu_mode, FILE *screen,
int* host_map, const int nelements, int*** host_elem2param, const int nparams,
const double* ts_lam1, const double* ts_lam2, const double* ts_lam3,
const double* ts_powermint, const double* ts_biga, const double* ts_bigb,
const double* ts_bigr, const double* ts_bigd,
const double* ts_c1, const double* ts_c2, const double* ts_c3, const double* ts_c4,
const double* ts_c, const double* ts_d, const double* ts_h,
const double* ts_gamma, const double* ts_beta,
const double* ts_powern, const double* ts_cutsq) {
TSMF.clear();
gpu_mode=TSMF.device->gpu_mode();
double gpu_split=TSMF.device->particle_split();
int first_gpu=TSMF.device->first_device();
int last_gpu=TSMF.device->last_device();
int world_me=TSMF.device->world_me();
int gpu_rank=TSMF.device->gpu_rank();
int procs_per_gpu=TSMF.device->procs_per_gpu();
// disable host/device split for now
if (gpu_split != 1.0)
return -8;
// disable multiple threads per atom for now
if (TSMF.device->threads_per_atom() != 1)
return -10;
TSMF.device->init_message(screen,"tersoff/gpu",first_gpu,last_gpu);
bool message=false;
if (TSMF.device->replica_me()==0 && screen)
message=true;
if (message) {
fprintf(screen,"Initializing Device and compiling on process 0...");
fflush(screen);
}
int init_ok=0;
if (world_me==0)
init_ok=TSMF.init(ntypes, inum, nall, 300, cell_size, gpu_split, screen,
host_map, nelements, host_elem2param, nparams,
ts_lam1, ts_lam2, ts_lam3, ts_powermint,
ts_biga, ts_bigb, ts_bigr, ts_bigd,
ts_c1, ts_c2, ts_c3, ts_c4, ts_c, ts_d, ts_h,
ts_gamma, ts_beta, ts_powern, ts_cutsq);
TSMF.device->world_barrier();
if (message)
fprintf(screen,"Done.\n");
for (int i=0; i<procs_per_gpu; i++) {
if (message) {
if (last_gpu-first_gpu==0)
fprintf(screen,"Initializing Device %d on core %d...",first_gpu,i);
else
fprintf(screen,"Initializing Devices %d-%d on core %d...",first_gpu,
last_gpu,i);
fflush(screen);
}
if (gpu_rank==i && world_me!=0)
init_ok=TSMF.init(ntypes, inum, nall, 300, cell_size, gpu_split, screen,
host_map, nelements, host_elem2param, nparams,
ts_lam1, ts_lam2, ts_lam3, ts_powermint,
ts_biga, ts_bigb, ts_bigr, ts_bigd,
ts_c1, ts_c2, ts_c3, ts_c4, ts_c, ts_d, ts_h,
ts_gamma, ts_beta, ts_powern, ts_cutsq);
TSMF.device->gpu_barrier();
if (message)
fprintf(screen,"Done.\n");
}
if (message)
fprintf(screen,"\n");
if (init_ok==0)
TSMF.estimate_gpu_overhead();
return init_ok;
}
void tersoff_gpu_clear() {
TSMF.clear();
}
int ** tersoff_gpu_compute_n(const int ago, const int inum_full,
const int nall, double **host_x, int *host_type,
double *sublo, double *subhi, tagint *tag, int **nspecial,
tagint **special, const bool eflag, const bool vflag,
const bool eatom, const bool vatom, int &host_start,
int **ilist, int **jnum, const double cpu_time,
bool &success) {
return TSMF.compute(ago, inum_full, nall, host_x, host_type, sublo,
subhi, tag, nspecial, special, eflag, vflag, eatom,
vatom, host_start, ilist, jnum, cpu_time, success);
}
void tersoff_gpu_compute(const int ago, const int nlocal, const int nall,
const int nlist, double **host_x, int *host_type,
int *ilist, int *numj, int **firstneigh, const bool eflag,
const bool vflag, const bool eatom, const bool vatom,
int &host_start, const double cpu_time, bool &success) {
TSMF.compute(ago,nlocal,nall,nlist,host_x,host_type,ilist,numj,
firstneigh,eflag,vflag,eatom,vatom,host_start,cpu_time,success);
}
double tersoff_gpu_bytes() {
return TSMF.host_memory_usage();
}

612
lib/gpu/lal_tersoff_extra.h Normal file
View File

@ -0,0 +1,612 @@
/// **************************************************************************
// tersoff_extra.h
// -------------------
// Trung Dac Nguyen
//
// Device code for Tersoff math routines
//
// __________________________________________________________________________
// This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
// __________________________________________________________________________
//
// begin :
// email : ndactrung@gmail.com
// ***************************************************************************/*
#ifndef LAL_TERSOFF_EXTRA_H
#define LAL_TERSOFF_EXTRA_H
#ifdef NV_KERNEL
#include "lal_aux_fun1.h"
#else
#endif
#define MY_PI2 (numtyp)1.57079632679489661923
#define MY_PI4 (numtyp)0.78539816339744830962
/* ---------------------------------------------------------------------- */
ucl_inline numtyp vec3_dot(const numtyp x[3], const numtyp y[3])
{
return (x[0]*y[0] + x[1]*y[1] + x[2]*y[2]);
}
ucl_inline void vec3_add(const numtyp x[3], const numtyp y[3], numtyp z[3])
{
z[0] = x[0]+y[0]; z[1] = x[1]+y[1]; z[2] = x[2]+y[2];
}
ucl_inline void vec3_scale(const numtyp k, const numtyp x[3], numtyp y[3])
{
y[0] = k*x[0]; y[1] = k*x[1]; y[2] = k*x[2];
}
ucl_inline void vec3_scaleadd(const numtyp k, const numtyp x[3],
const numtyp y[3], numtyp z[3])
{
z[0] = k*x[0]+y[0]; z[1] = k*x[1]+y[1]; z[2] = k*x[2]+y[2];
}
/* ---------------------------------------------------------------------- */
ucl_inline numtyp ters_gijk(const numtyp costheta,
const numtyp param_c,
const numtyp param_d,
const numtyp param_h,
const numtyp param_gamma)
{
const numtyp ters_c = param_c * param_c;
const numtyp ters_d = param_d * param_d;
const numtyp hcth = param_h - costheta;
return param_gamma*((numtyp)1.0 + ters_c*ucl_recip(ters_d) -
ters_c *ucl_recip(ters_d + hcth*hcth));
}
/* ---------------------------------------------------------------------- */
ucl_inline numtyp ters_gijk_d(const numtyp costheta,
const numtyp param_c,
const numtyp param_d,
const numtyp param_h,
const numtyp param_gamma)
{
const numtyp ters_c = param_c * param_c;
const numtyp ters_d = param_d * param_d;
const numtyp hcth = param_h - costheta;
const numtyp numerator = (numtyp)-2.0 * ters_c * hcth;
const numtyp denominator = ucl_recip(ters_d + hcth*hcth);
return param_gamma*numerator*denominator*denominator;
}
/* ---------------------------------------------------------------------- */
ucl_inline void costheta_d(const numtyp rij_hat[3],
const numtyp rij,
const numtyp rik_hat[3],
const numtyp rik,
numtyp *dri,
numtyp *drj,
numtyp *drk)
{
// first element is derivative wrt Ri, second wrt Rj, third wrt Rk
numtyp cos_theta = vec3_dot(rij_hat,rik_hat);
vec3_scaleadd(-cos_theta,rij_hat,rik_hat,drj);
vec3_scale(ucl_recip(rij),drj,drj);
vec3_scaleadd(-cos_theta,rik_hat,rij_hat,drk);
vec3_scale(ucl_recip(rik),drk,drk);
vec3_add(drj,drk,dri);
vec3_scale((numtyp)-1.0,dri,dri);
}
/* ---------------------------------------------------------------------- */
ucl_inline numtyp ters_fc(const numtyp r,
const numtyp param_bigr,
const numtyp param_bigd)
{
if (r < param_bigr-param_bigd) return (numtyp)1.0;
if (r > param_bigr+param_bigd) return (numtyp)0.0;
return (numtyp)0.5*((numtyp)1.0 - sin(MY_PI2*(r - param_bigr)/param_bigd));
}
/* ---------------------------------------------------------------------- */
ucl_inline numtyp ters_fc_d(const numtyp r,
const numtyp param_bigr,
const numtyp param_bigd)
{
if (r < param_bigr-param_bigd) return (numtyp)0.0;
if (r > param_bigr+param_bigd) return (numtyp)0.0;
return -(MY_PI4/param_bigd) * cos(MY_PI2*(r - param_bigr)/param_bigd);
}
/* ---------------------------------------------------------------------- */
ucl_inline numtyp ters_fa(const numtyp r,
const numtyp param_bigb,
const numtyp param_bigr,
const numtyp param_bigd,
const numtyp param_lam2)
{
if (r > param_bigr + param_bigd) return (numtyp)0.0;
return -param_bigb * ucl_exp(-param_lam2 * r) *
ters_fc(r,param_bigr,param_bigd);
}
/* ---------------------------------------------------------------------- */
ucl_inline numtyp ters_fa_d(const numtyp r,
const numtyp param_bigb,
const numtyp param_bigr,
const numtyp param_bigd,
const numtyp param_lam2)
{
if (r > param_bigr + param_bigd) return (numtyp)0.0;
return param_bigb * ucl_exp(-param_lam2 * r) * (param_lam2 *
ters_fc(r,param_bigr,param_bigd) - ters_fc_d(r,param_bigr,param_bigd));
}
/* ---------------------------------------------------------------------- */
ucl_inline numtyp ters_bij(const numtyp zeta,
const numtyp param_beta,
const numtyp param_powern,
const numtyp param_c1,
const numtyp param_c2,
const numtyp param_c3,
const numtyp param_c4)
{
numtyp tmp = param_beta * zeta;
if (tmp > param_c1) return ucl_rsqrt(tmp);
if (tmp > param_c2)
return ((numtyp)1.0 - ucl_powr(tmp,-param_powern) /
((numtyp)2.0*param_powern))*ucl_rsqrt(tmp);
if (tmp < param_c4) return (numtyp)1.0;
if (tmp < param_c3)
return (numtyp)1.0 - ucl_powr(tmp,param_powern)/((numtyp)2.0*param_powern);
return ucl_powr((numtyp)1.0 + ucl_powr(tmp,param_powern),
(numtyp)-1.0/((numtyp)2.0*param_powern));
}
/* ---------------------------------------------------------------------- */
ucl_inline numtyp ters_bij_d(const numtyp zeta,
const numtyp param_beta,
const numtyp param_powern,
const numtyp param_c1,
const numtyp param_c2,
const numtyp param_c3,
const numtyp param_c4)
{
numtyp tmp = param_beta * zeta;
if (tmp > param_c1)
return param_beta * (numtyp)-0.5*ucl_powr(tmp,(numtyp)-1.5);
if (tmp > param_c2)
return param_beta * ((numtyp)-0.5*ucl_powr(tmp,(numtyp)-1.5) *
((numtyp)1.0 - (numtyp)0.5 * ((numtyp)1.0 + (numtyp)1.0 /
((numtyp)2.0 * param_powern)) * ucl_powr(tmp,-param_powern)));
if (tmp < param_c4) return (numtyp)0.0;
if (tmp < param_c3)
return (numtyp)-0.5*param_beta * ucl_powr(tmp,param_powern-(numtyp)1.0);
numtyp tmp_n = ucl_powr(tmp,param_powern);
return (numtyp)-0.5 * ucl_powr((numtyp)1.0+tmp_n, (numtyp) -
(numtyp)1.0-((numtyp)1.0 / ((numtyp)2.0 * param_powern)))*tmp_n / zeta;
}
/* ---------------------------------------------------------------------- */
ucl_inline void ters_zetaterm_d(const numtyp prefactor,
const numtyp rij_hat[3],
const numtyp rij,
const numtyp rik_hat[3],
const numtyp rik,
const numtyp param_bigr,
const numtyp param_bigd,
const numtyp param_powermint,
const numtyp param_lam3,
const numtyp param_c,
const numtyp param_d,
const numtyp param_h,
const numtyp param_gamma,
numtyp dri[3],
numtyp drj[3],
numtyp drk[3])
{
numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp;
numtyp dcosdri[3],dcosdrj[3],dcosdrk[3];
fc = ters_fc(rik,param_bigr,param_bigd);
dfc = ters_fc_d(rik,param_bigr,param_bigd);
numtyp t = param_lam3*(rij-rik);
if ((int)param_powermint == 3) tmp = t*t*t;
else tmp = t;
if (tmp > (numtyp)69.0776) ex_delr = (acctyp)1.e30;
else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
else ex_delr = ucl_exp(tmp);
if ((int)param_powermint == 3)
ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr;
else ex_delr_d = param_lam3 * ex_delr;
cos_theta = vec3_dot(rij_hat,rik_hat);
gijk = ters_gijk(cos_theta,param_c,param_d,param_h,param_gamma);
gijk_d = ters_gijk_d(cos_theta,param_c,param_d,param_h,param_gamma);
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
// compute the derivative wrt Ri
// dri = -dfc*gijk*ex_delr*rik_hat;
// dri += fc*gijk_d*ex_delr*dcosdri;
// dri += fc*gijk*ex_delr_d*(rik_hat - rij_hat);
vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri);
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri);
vec3_scaleadd(fc*gijk*ex_delr_d,rik_hat,dri,dri);
vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri);
vec3_scale(prefactor,dri,dri);
// compute the derivative wrt Rj
// drj = fc*gijk_d*ex_delr*dcosdrj;
// drj += fc*gijk*ex_delr_d*rij_hat;
vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj);
vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj);
vec3_scale(prefactor,drj,drj);
// compute the derivative wrt Rk
// drk = dfc*gijk*ex_delr*rik_hat;
// drk += fc*gijk_d*ex_delr*dcosdrk;
// drk += -fc*gijk*ex_delr_d*rik_hat;
vec3_scale(dfc*gijk*ex_delr,rik_hat,drk);
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk);
vec3_scaleadd(-fc*gijk*ex_delr_d,rik_hat,drk,drk);
vec3_scale(prefactor,drk,drk);
}
ucl_inline void ters_zetaterm_d_fi(const numtyp prefactor,
const numtyp rij_hat[3],
const numtyp rij,
const numtyp rik_hat[3],
const numtyp rik,
const numtyp param_bigr,
const numtyp param_bigd,
const numtyp param_powermint,
const numtyp param_lam3,
const numtyp param_c,
const numtyp param_d,
const numtyp param_h,
const numtyp param_gamma,
numtyp dri[3])
{
numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp;
numtyp dcosdri[3],dcosdrj[3],dcosdrk[3];
fc = ters_fc(rik,param_bigr,param_bigd);
dfc = ters_fc_d(rik,param_bigr,param_bigd);
numtyp t = param_lam3*(rij-rik);
if ((int)param_powermint == 3) tmp = t*t*t;
else tmp = t;
if (tmp > (numtyp)69.0776) ex_delr = (acctyp)1.e30;
else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
else ex_delr = ucl_exp(tmp);
if ((int)param_powermint == 3)
ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr;
else ex_delr_d = param_lam3 * ex_delr;
cos_theta = vec3_dot(rij_hat,rik_hat);
gijk = ters_gijk(cos_theta,param_c,param_d,param_h,param_gamma);
gijk_d = ters_gijk_d(cos_theta,param_c,param_d,param_h,param_gamma);
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
// compute the derivative wrt Ri
// dri = -dfc*gijk*ex_delr*rik_hat;
// dri += fc*gijk_d*ex_delr*dcosdri;
// dri += fc*gijk*ex_delr_d*(rik_hat - rij_hat);
vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri);
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri);
vec3_scaleadd(fc*gijk*ex_delr_d,rik_hat,dri,dri);
vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri);
vec3_scale(prefactor,dri,dri);
}
ucl_inline void ters_zetaterm_d_fj(const numtyp prefactor,
const numtyp rij_hat[3],
const numtyp rij,
const numtyp rik_hat[3],
const numtyp rik,
const numtyp param_bigr,
const numtyp param_bigd,
const numtyp param_powermint,
const numtyp param_lam3,
const numtyp param_c,
const numtyp param_d,
const numtyp param_h,
const numtyp param_gamma,
numtyp drj[3])
{
numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,cos_theta,tmp;
numtyp dcosdri[3],dcosdrj[3],dcosdrk[3];
fc = ters_fc(rik,param_bigr,param_bigd);
numtyp t = param_lam3*(rij-rik);
if ((int)param_powermint == 3) tmp = t*t*t;
else tmp = t;
if (tmp > (numtyp)69.0776) ex_delr = (acctyp)1.e30;
else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
else ex_delr = ucl_exp(tmp);
if ((int)param_powermint == 3)
ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr;
else ex_delr_d = param_lam3 * ex_delr;
cos_theta = vec3_dot(rij_hat,rik_hat);
gijk = ters_gijk(cos_theta,param_c,param_d,param_h,param_gamma);
gijk_d = ters_gijk_d(cos_theta,param_c,param_d,param_h,param_gamma);
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
// compute the derivative wrt Rj
// drj = fc*gijk_d*ex_delr*dcosdrj;
// drj += fc*gijk*ex_delr_d*rij_hat;
vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj);
vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj);
vec3_scale(prefactor,drj,drj);
}
ucl_inline void ters_zetaterm_d_fk(const numtyp prefactor,
const numtyp rij_hat[3],
const numtyp rij,
const numtyp rik_hat[3],
const numtyp rik,
const numtyp param_bigr,
const numtyp param_bigd,
const numtyp param_powermint,
const numtyp param_lam3,
const numtyp param_c,
const numtyp param_d,
const numtyp param_h,
const numtyp param_gamma,
numtyp drk[3])
{
numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp;
numtyp dcosdri[3],dcosdrj[3],dcosdrk[3];
fc = ters_fc(rik,param_bigr,param_bigd);
dfc = ters_fc_d(rik,param_bigr,param_bigd);
numtyp t = param_lam3*(rij-rik);
if ((int)param_powermint == 3) tmp = t*t*t;
else tmp = t;
if (tmp > (numtyp)69.0776) ex_delr = (acctyp)1.e30;
else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
else ex_delr = ucl_exp(tmp);
if ((int)param_powermint == 3)
ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr;
else ex_delr_d = param_lam3 * ex_delr;
cos_theta = vec3_dot(rij_hat,rik_hat);
gijk = ters_gijk(cos_theta,param_c,param_d,param_h,param_gamma);
gijk_d = ters_gijk_d(cos_theta,param_c,param_d,param_h,param_gamma);
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
// compute the derivative wrt Rk
// drk = dfc*gijk*ex_delr*rik_hat;
// drk += fc*gijk_d*ex_delr*dcosdrk;
// drk += -fc*gijk*ex_delr_d*rik_hat;
vec3_scale(dfc*gijk*ex_delr,rik_hat,drk);
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk);
vec3_scaleadd(-fc*gijk*ex_delr_d,rik_hat,drk,drk);
vec3_scale(prefactor,drk,drk);
}
/* ---------------------------------------------------------------------- */
ucl_inline void repulsive(const numtyp param_bigr,
const numtyp param_bigd,
const numtyp param_lam1,
const numtyp param_biga,
const numtyp rsq,
const int eflag,
numtyp *ans)
{
numtyp r,tmp_fc,tmp_fc_d,tmp_exp;
r = ucl_sqrt(rsq);
tmp_fc = ters_fc(r,param_bigr,param_bigd);
tmp_fc_d = ters_fc_d(r,param_bigr,param_bigd);
tmp_exp = ucl_exp(-param_lam1 * r);
// fforce
ans[0] = -param_biga*tmp_exp*(tmp_fc_d - tmp_fc*param_lam1)*ucl_recip(r);
// eng
if (eflag) ans[1] = tmp_fc * param_biga * tmp_exp;
}
/* ---------------------------------------------------------------------- */
ucl_inline numtyp zeta(const numtyp param_powermint,
const numtyp param_lam3,
const numtyp param_bigr,
const numtyp param_bigd,
const numtyp param_c,
const numtyp param_d,
const numtyp param_h,
const numtyp param_gamma,
const numtyp rsqij,
const numtyp rsqik,
const numtyp4 delrij,
const numtyp4 delrik)
{
numtyp rij,rik,costheta,arg,ex_delr;
rij = ucl_sqrt(rsqij);
rik = ucl_sqrt(rsqik);
costheta = (delrij.x*delrik.x + delrij.y*delrik.y +
delrij.z*delrik.z) / (rij*rik);
numtyp t = param_lam3*(rij-rik);
if ((int)param_powermint == 3) arg = t*t*t;
else arg = t;
if (arg > (numtyp)69.0776) ex_delr = (numtyp)1.e30;
else if (arg < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
else ex_delr = ucl_exp(arg);
return ters_fc(rik,param_bigr,param_bigd) *
ters_gijk(costheta,param_c, param_d, param_h, param_gamma) * ex_delr;
}
/* ---------------------------------------------------------------------- */
ucl_inline void force_zeta(const numtyp param_bigb,
const numtyp param_bigr,
const numtyp param_bigd,
const numtyp param_lam2,
const numtyp param_beta,
const numtyp param_powern,
const numtyp param_c1,
const numtyp param_c2,
const numtyp param_c3,
const numtyp param_c4,
const numtyp rsq,
const numtyp zeta_ij,
const int eflag,
numtyp fpfeng[4])
{
numtyp r,fa,fa_d,bij;
r = ucl_sqrt(rsq);
fa = ters_fa(r,param_bigb,param_bigr,param_bigd,param_lam2);
fa_d = ters_fa_d(r,param_bigb,param_bigr,param_bigd,param_lam2);
bij = ters_bij(zeta_ij,param_beta,param_powern,
param_c1,param_c2, param_c3, param_c4);
fpfeng[0] = (numtyp)0.5*bij*fa_d * ucl_recip(r); // fforce
fpfeng[1] = (numtyp)-0.5*fa * ters_bij_d(zeta_ij,param_beta, param_powern,
param_c1,param_c2, param_c3, param_c4); // prefactor
if (eflag) fpfeng[2] = (numtyp)0.5*bij*fa; // eng
}
/* ----------------------------------------------------------------------
attractive term
use param_ij cutoff for rij test
use param_ijk cutoff for rik test
------------------------------------------------------------------------- */
ucl_inline void attractive(const numtyp param_bigr,
const numtyp param_bigd,
const numtyp param_powermint,
const numtyp param_lam3,
const numtyp param_c,
const numtyp param_d,
const numtyp param_h,
const numtyp param_gamma,
const numtyp prefactor,
const numtyp rij,
const numtyp rijinv,
const numtyp rik,
const numtyp rikinv,
const numtyp delrij[3],
const numtyp delrik[3],
numtyp fi[3],
numtyp fj[3],
numtyp fk[3])
{
numtyp rij_hat[3],rik_hat[3];
vec3_scale(rijinv,delrij,rij_hat);
vec3_scale(rikinv,delrik,rik_hat);
ters_zetaterm_d(prefactor,rij_hat,rij,rik_hat,rik,
param_bigr, param_bigd, param_powermint, param_lam3,
param_c, param_d, param_h, param_gamma, fi, fj, fk);
}
ucl_inline void attractive_fi(const numtyp param_bigr,
const numtyp param_bigd,
const numtyp param_powermint,
const numtyp param_lam3,
const numtyp param_c,
const numtyp param_d,
const numtyp param_h,
const numtyp param_gamma,
const numtyp prefactor,
const numtyp rij,
const numtyp rijinv,
const numtyp rik,
const numtyp rikinv,
const numtyp delrij[3],
const numtyp delrik[3],
numtyp fi[3])
{
numtyp rij_hat[3],rik_hat[3];
vec3_scale(rijinv,delrij,rij_hat);
vec3_scale(rikinv,delrik,rik_hat);
ters_zetaterm_d_fi(prefactor,rij_hat,rij,rik_hat,rik,
param_bigr, param_bigd, param_powermint, param_lam3,
param_c, param_d, param_h, param_gamma, fi);
}
ucl_inline void attractive_fj(const numtyp param_bigr,
const numtyp param_bigd,
const numtyp param_powermint,
const numtyp param_lam3,
const numtyp param_c,
const numtyp param_d,
const numtyp param_h,
const numtyp param_gamma,
const numtyp prefactor,
const numtyp rij,
const numtyp rijinv,
const numtyp rik,
const numtyp rikinv,
const numtyp delrij[3],
const numtyp delrik[3],
numtyp fj[3])
{
numtyp rij_hat[3],rik_hat[3];
vec3_scale(rijinv,delrij,rij_hat);
vec3_scale(rikinv,delrik,rik_hat);
ters_zetaterm_d_fj(prefactor,rij_hat,rij,rik_hat,rik,
param_bigr, param_bigd, param_powermint, param_lam3,
param_c, param_d, param_h, param_gamma, fj);
}
ucl_inline void attractive_fk(const numtyp param_bigr,
const numtyp param_bigd,
const numtyp param_powermint,
const numtyp param_lam3,
const numtyp param_c,
const numtyp param_d,
const numtyp param_h,
const numtyp param_gamma,
const numtyp prefactor,
const numtyp rij,
const numtyp rijinv,
const numtyp rik,
const numtyp rikinv,
const numtyp delrij[3],
const numtyp delrik[3],
numtyp fk[3])
{
numtyp rij_hat[3],rik_hat[3];
vec3_scale(rijinv,delrij,rij_hat);
vec3_scale(rikinv,delrik,rik_hat);
ters_zetaterm_d_fk(prefactor,rij_hat,rij,rik_hat,rik,
param_bigr, param_bigd, param_powermint, param_lam3,
param_c, param_d, param_h, param_gamma, fk);
}
#endif

159
lib/gpu/lal_zbl.cpp Normal file
View File

@ -0,0 +1,159 @@
/***************************************************************************
zbl.cpp
-------------------
Trung Dac Nguyen
Class for acceleration of the zbl pair style.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
__________________________________________________________________________
begin :
email : ndactrung@gmail.com
***************************************************************************/
#ifdef USE_OPENCL
#include "zbl_cl.h"
#elif defined(USE_CUDART)
const char *zbl=0;
#else
#include "zbl_cubin.h"
#endif
#include "lal_zbl.h"
#include <cassert>
using namespace LAMMPS_AL;
#define ZBLT ZBL<numtyp, acctyp>
extern Device<PRECISION,ACC_PRECISION> device;
template <class numtyp, class acctyp>
ZBLT::ZBL() : BaseAtomic<numtyp,acctyp>(), _allocated(false) {
}
template <class numtyp, class acctyp>
ZBLT::~ZBL() {
clear();
}
template <class numtyp, class acctyp>
int ZBLT::bytes_per_atom(const int max_nbors) const {
return this->bytes_per_atom_atomic(max_nbors);
}
template <class numtyp, class acctyp>
int ZBLT::init(const int ntypes, double **host_cutsq,
double **host_sw1, double **host_sw2,
double **host_sw3, double **host_sw4,
double **host_sw5,
double **host_d1a, double **host_d2a,
double **host_d3a, double **host_d4a,
double **host_zze, double cut_globalsq,
double cut_innersq, double cut_inner,
const int nlocal, const int nall, const int max_nbors,
const int maxspecial, const double cell_size,
const double gpu_split, FILE *_screen) {
int success;
success=this->init_atomic(nlocal,nall,max_nbors,maxspecial,cell_size,gpu_split,
_screen,zbl,"k_zbl");
if (success!=0)
return success;
// If atom type constants fit in shared memory use fast kernel
int lj_types=ntypes;
shared_types=false;
int max_shared_types=this->device->max_shared_types();
if (lj_types<=max_shared_types && this->_block_size>=max_shared_types) {
lj_types=max_shared_types;
shared_types=true;
}
_lj_types=lj_types;
// Allocate a host write buffer for data initialization
UCL_H_Vec<numtyp> host_write(lj_types*lj_types*32,*(this->ucl_device),
UCL_WRITE_ONLY);
for (int i=0; i<lj_types*lj_types; i++)
host_write[i]=0.0;
coeff1.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
this->atom->type_pack4(ntypes,lj_types,coeff1,host_write,host_sw1,host_sw2,
host_zze, host_cutsq);
coeff2.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
this->atom->type_pack4(ntypes,lj_types,coeff2,host_write,host_d1a,host_d2a,
host_d3a,host_d4a);
coeff3.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
this->atom->type_pack4(ntypes,lj_types,coeff3,host_write,host_sw3,host_sw4,host_sw5);
_cut_globalsq = cut_globalsq;
_cut_innersq = cut_innersq;
_cut_inner = cut_inner;
printf("params: %f %f %f\n", _cut_globalsq, _cut_innersq, _cut_inner);
_allocated=true;
this->_max_bytes=coeff1.row_bytes()+coeff2.row_bytes()+coeff3.row_bytes();
return 0;
}
template <class numtyp, class acctyp>
void ZBLT::clear() {
if (!_allocated)
return;
_allocated=false;
coeff1.clear();
coeff2.clear();
coeff3.clear();
this->clear_atomic();
}
template <class numtyp, class acctyp>
double ZBLT::host_memory_usage() const {
return this->host_memory_usage_atomic()+sizeof(ZBL<numtyp,acctyp>);
}
// ---------------------------------------------------------------------------
// Calculate energies, forces, and torques
// ---------------------------------------------------------------------------
template <class numtyp, class acctyp>
void ZBLT::loop(const bool _eflag, const bool _vflag) {
// Compute the block size and grid size to keep all cores busy
const int BX=this->block_size();
int eflag, vflag;
if (_eflag)
eflag=1;
else
eflag=0;
if (_vflag)
vflag=1;
else
vflag=0;
int GX=static_cast<int>(ceil(static_cast<double>(this->ans->inum())/
(BX/this->_threads_per_atom)));
int ainum=this->ans->inum();
int nbor_pitch=this->nbor->nbor_pitch();
this->time_pair.start();
if (shared_types) {
this->k_pair_fast.set_size(GX,BX);
this->k_pair_fast.run(&this->atom->x, &coeff1, &coeff2, &coeff3,
&_cut_globalsq, &_cut_innersq, &_cut_inner,
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
&this->ans->force, &this->ans->engv, &eflag, &vflag,
&ainum, &nbor_pitch, &this->_threads_per_atom);
} else {
this->k_pair.set_size(GX,BX);
this->k_pair.run(&this->atom->x, &coeff1, &coeff2, &coeff3,
&_cut_globalsq, &_cut_innersq, &_cut_inner, &_lj_types,
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
&this->ans->force, &this->ans->engv, &eflag, &vflag,
&ainum, &nbor_pitch, &this->_threads_per_atom);
}
this->time_pair.stop();
}
template class ZBL<PRECISION,ACC_PRECISION>;

84
lib/gpu/lal_zbl.h Normal file
View File

@ -0,0 +1,84 @@
/***************************************************************************
zbl.h
-------------------
Trung Dac Nguyen (ORNL)
Class for acceleration of the zbl pair style.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
__________________________________________________________________________
begin :
email : ndactrung@gmail.com
***************************************************************************/
#ifndef LAL_ZBL_H
#define LAL_ZBL_H
#include "lal_base_atomic.h"
namespace LAMMPS_AL {
template <class numtyp, class acctyp>
class ZBL : public BaseAtomic<numtyp, acctyp> {
public:
ZBL();
~ZBL();
/// Clear any previous data and set up for a new LAMMPS run
/** \param max_nbors initial number of rows in the neighbor matrix
* \param cell_size cutoff + skin
* \param gpu_split fraction of particles handled by device
*
* Returns:
* - 0 if successfull
* - -1 if fix gpu not found
* - -3 if there is an out of memory error
* - -4 if the GPU library was not compiled for GPU
* - -5 Double precision is not supported on card **/
int init(const int ntypes, double **host_cutsq, double **host_sw1,
double **host_sw2, double **host_sw3, double **host_sw4, double **host_sw5,
double **host_d1a, double **host_d2a, double **host_d3a, double **host_d4a,
double **host_zze, double cut_globalsq, double cut_innersq, double cut_inner,
const int nlocal, const int nall, const int max_nbors,
const int maxspecial, const double cell_size,
const double gpu_split, FILE *screen);
/// Clear all host and device data
/** \note This is called at the beginning of the init() routine **/
void clear();
/// Returns memory usage on device per atom
int bytes_per_atom(const int max_nbors) const;
/// Total host memory used by library for pair style
double host_memory_usage() const;
// --------------------------- TYPE DATA --------------------------
/// coeff1.x = sw1, coeff1.y = sw2, coeff1.z = zze, coeff1.w = cutsq
UCL_D_Vec<numtyp4> coeff1;
/// coeff2.x = d1a, coeff2.y = d2a, coeff2.z = d3a, coeff2.w = d4a
UCL_D_Vec<numtyp4> coeff2;
/// coeff3.x = sw3, coeff3.y = sw4, coeff3.z = sw5;
UCL_D_Vec<numtyp4> coeff3;
/// If atom type constants fit in shared memory, use fast kernels
bool shared_types;
double _cut_globalsq;
double _cut_innersq;
double _cut_inner;
/// Number of atom types
int _lj_types;
private:
bool _allocated;
void loop(const bool _eflag, const bool _vflag);
};
}
#endif

123
lib/gpu/lal_zbl_ext.cpp Normal file
View File

@ -0,0 +1,123 @@
/***************************************************************************
zbl_ext.cpp
-------------------
Trung Dac Nguyen
Functions for LAMMPS access to zbl acceleration routines.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
__________________________________________________________________________
begin :
email : ndactrung@gmail.com
***************************************************************************/
#include <iostream>
#include <cassert>
#include <math.h>
#include "lal_zbl.h"
using namespace std;
using namespace LAMMPS_AL;
static ZBL<PRECISION,ACC_PRECISION> ZBLMF;
// ---------------------------------------------------------------------------
// Allocate memory on host and device and copy constants to device
// ---------------------------------------------------------------------------
int zbl_gpu_init(const int ntypes, double **cutsq, double **host_sw1,
double **host_sw2, double **host_sw3, double **host_sw4, double **host_sw5,
double **host_d1a, double **host_d2a, double **host_d3a, double **host_d4a,
double **host_zze, double cut_globalsq, double cut_innersq, double cut_inner,
const int inum, const int nall, const int max_nbors,
const int maxspecial, const double cell_size, int &gpu_mode, FILE *screen) {
ZBLMF.clear();
gpu_mode=ZBLMF.device->gpu_mode();
double gpu_split=ZBLMF.device->particle_split();
int first_gpu=ZBLMF.device->first_device();
int last_gpu=ZBLMF.device->last_device();
int world_me=ZBLMF.device->world_me();
int gpu_rank=ZBLMF.device->gpu_rank();
int procs_per_gpu=ZBLMF.device->procs_per_gpu();
ZBLMF.device->init_message(screen,"zbl",first_gpu,last_gpu);
bool message=false;
if (ZBLMF.device->replica_me()==0 && screen)
message=true;
if (message) {
fprintf(screen,"Initializing Device and compiling on process 0...");
fflush(screen);
}
int init_ok=0;
if (world_me==0)
init_ok=ZBLMF.init(ntypes, cutsq, host_sw1, host_sw2, host_sw3, host_sw4,
host_sw5, host_d1a, host_d2a, host_d3a, host_d4a, host_zze,
cut_globalsq, cut_innersq, cut_inner,
inum, nall, 300, maxspecial, cell_size, gpu_split, screen);
ZBLMF.device->world_barrier();
if (message)
fprintf(screen,"Done.\n");
for (int i=0; i<procs_per_gpu; i++) {
if (message) {
if (last_gpu-first_gpu==0)
fprintf(screen,"Initializing Device %d on core %d...",first_gpu,i);
else
fprintf(screen,"Initializing Devices %d-%d on core %d...",first_gpu,
last_gpu,i);
fflush(screen);
}
if (gpu_rank==i && world_me!=0)
init_ok=ZBLMF.init(ntypes, cutsq, host_sw1, host_sw2, host_sw3, host_sw4,
host_sw5, host_d1a, host_d2a, host_d3a, host_d4a, host_zze,
cut_globalsq, cut_innersq, cut_inner,
inum, nall, 300, maxspecial, cell_size, gpu_split, screen);
ZBLMF.device->gpu_barrier();
if (message)
fprintf(screen,"Done.\n");
}
if (message)
fprintf(screen,"\n");
if (init_ok==0)
ZBLMF.estimate_gpu_overhead();
return init_ok;
}
void zbl_gpu_clear() {
ZBLMF.clear();
}
int ** zbl_gpu_compute_n(const int ago, const int inum_full,
const int nall, double **host_x, int *host_type,
double *sublo, double *subhi, tagint *tag, int **nspecial,
tagint **special, const bool eflag, const bool vflag,
const bool eatom, const bool vatom, int &host_start,
int **ilist, int **jnum, const double cpu_time,
bool &success) {
return ZBLMF.compute(ago, inum_full, nall, host_x, host_type, sublo,
subhi, tag, nspecial, special, eflag, vflag, eatom,
vatom, host_start, ilist, jnum, cpu_time, success);
}
void zbl_gpu_compute(const int ago, const int inum_full, const int nall,
double **host_x, int *host_type, int *ilist, int *numj,
int **firstneigh, const bool eflag, const bool vflag,
const bool eatom, const bool vatom, int &host_start,
const double cpu_time, bool &success) {
ZBLMF.compute(ago,inum_full,nall,host_x,host_type,ilist,numj,
firstneigh,eflag,vflag,eatom,vatom,host_start,cpu_time,success);
}
double zbl_gpu_bytes() {
return ZBLMF.host_memory_usage();
}