forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14808 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
67a4004f23
commit
008896a77d
|
|
@ -0,0 +1,456 @@
|
|||
/***************************************************************************
|
||||
tersoff_mod.cpp
|
||||
-------------------
|
||||
Trung Dac Nguyen
|
||||
|
||||
Class for acceleration of the tersoff pair style.
|
||||
|
||||
__________________________________________________________________________
|
||||
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
|
||||
__________________________________________________________________________
|
||||
|
||||
begin :
|
||||
email : ndactrung@gmail.com
|
||||
***************************************************************************/
|
||||
|
||||
#if defined(USE_OPENCL)
|
||||
#include "tersoff_mod_cl.h"
|
||||
#elif defined(USE_CUDART)
|
||||
const char *tersoff_mod=0;
|
||||
#else
|
||||
#include "tersoff_mod_cubin.h"
|
||||
#endif
|
||||
|
||||
#include "lal_tersoff_mod.h"
|
||||
#include <cassert>
|
||||
using namespace LAMMPS_AL;
|
||||
#define TersoffMT TersoffMod<numtyp, acctyp>
|
||||
|
||||
extern Device<PRECISION,ACC_PRECISION> device;
|
||||
|
||||
template <class numtyp, class acctyp>
|
||||
TersoffMT::TersoffMod() : BaseThree<numtyp,acctyp>(), _allocated(false) {
|
||||
}
|
||||
|
||||
template <class numtyp, class acctyp>
|
||||
TersoffMT::~TersoffMod() {
|
||||
clear();
|
||||
}
|
||||
|
||||
template <class numtyp, class acctyp>
|
||||
int TersoffMT::bytes_per_atom(const int max_nbors) const {
|
||||
return this->bytes_per_atom_atomic(max_nbors);
|
||||
}
|
||||
|
||||
template <class numtyp, class acctyp>
|
||||
int TersoffMT::init(const int ntypes, const int nlocal, const int nall, const int max_nbors,
|
||||
const double cell_size, const double gpu_split, FILE *_screen,
|
||||
int* host_map, const int nelements, int*** host_elem2param, const int nparams,
|
||||
const double* lam1, const double* lam2, const double* lam3,const double* powermint,
|
||||
const double* biga, const double* bigb, const double* bigr, const double* bigd,
|
||||
const double* c1, const double* c2, const double* c3, const double* c4,
|
||||
const double* c5, const double* h, const double* beta, const double* powern,
|
||||
const double* powern_del, const double* ca1, const double* host_cutsq)
|
||||
{
|
||||
int success;
|
||||
success=this->init_three(nlocal,nall,max_nbors,0,cell_size,gpu_split,
|
||||
_screen,tersoff_mod,"k_tersoff_mod_repulsive",
|
||||
"k_tersoff_mod_three_center", "k_tersoff_mod_three_end");
|
||||
if (success!=0)
|
||||
return success;
|
||||
|
||||
int ef_nall=nall;
|
||||
if (ef_nall==0)
|
||||
ef_nall=2000;
|
||||
_zetaij.alloc(ef_nall*max_nbors,*(this->ucl_device),UCL_READ_WRITE);
|
||||
|
||||
k_zeta.set_function(*(this->pair_program),"k_tersoff_mod_zeta");
|
||||
|
||||
// If atom type constants fit in shared memory use fast kernel
|
||||
int lj_types=ntypes;
|
||||
shared_types=false;
|
||||
int max_shared_types=this->device->max_shared_types();
|
||||
if (lj_types<=max_shared_types && this->_block_size>=max_shared_types) {
|
||||
lj_types=max_shared_types;
|
||||
shared_types=true;
|
||||
}
|
||||
_lj_types=lj_types;
|
||||
|
||||
_nparams = nparams;
|
||||
_nelements = nelements;
|
||||
|
||||
UCL_H_Vec<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>(beta[i]);
|
||||
dview[i].y=static_cast<numtyp>(powern[i]);
|
||||
dview[i].z=static_cast<numtyp>(powern_del[i]);
|
||||
dview[i].w=static_cast<numtyp>(ca1[i]);
|
||||
}
|
||||
|
||||
ucl_copy(ts3,dview,false);
|
||||
ts3_tex.get_texture(*(this->pair_program),"ts3_tex");
|
||||
ts3_tex.bind_float(ts3,4);
|
||||
|
||||
ts4.alloc(nparams,*(this->ucl_device),UCL_READ_ONLY);
|
||||
|
||||
for (int i=0; i<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(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>(c5[i]);
|
||||
dview[i].y=static_cast<numtyp>(h[i]);
|
||||
dview[i].z=(numtyp)0;
|
||||
dview[i].w=(numtyp)0;
|
||||
}
|
||||
|
||||
ucl_copy(ts5,dview,false);
|
||||
ts5_tex.get_texture(*(this->pair_program),"ts5_tex");
|
||||
ts5_tex.bind_float(ts5,4);
|
||||
|
||||
UCL_H_Vec<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 < ntypes; i++)
|
||||
dview_map[i] = host_map[i];
|
||||
|
||||
map.alloc(lj_types,*(this->ucl_device), UCL_READ_ONLY);
|
||||
ucl_copy(map,dview_map,false);
|
||||
|
||||
_allocated=true;
|
||||
this->_max_bytes=ts1.row_bytes()+ts2.row_bytes()+ts3.row_bytes()+
|
||||
ts4.row_bytes()+cutsq.row_bytes()+
|
||||
map.row_bytes()+elem2param.row_bytes()+_zetaij.row_bytes();
|
||||
return 0;
|
||||
}
|
||||
|
||||
template <class numtyp, class acctyp>
|
||||
void TersoffMT::clear() {
|
||||
if (!_allocated)
|
||||
return;
|
||||
_allocated=false;
|
||||
|
||||
ts1.clear();
|
||||
ts2.clear();
|
||||
ts3.clear();
|
||||
ts4.clear();
|
||||
ts5.clear();
|
||||
cutsq.clear();
|
||||
map.clear();
|
||||
elem2param.clear();
|
||||
_zetaij.clear();
|
||||
|
||||
k_zeta.clear();
|
||||
|
||||
this->clear_atomic();
|
||||
}
|
||||
|
||||
template <class numtyp, class acctyp>
|
||||
double TersoffMT::host_memory_usage() const {
|
||||
return this->host_memory_usage_atomic()+sizeof(TersoffMod<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 TersoffMT::compute(const int f_ago, const int nlocal, const int nall,
|
||||
const int nlist, double **host_x, int *host_type,
|
||||
int *ilist, int *numj, int **firstneigh,
|
||||
const bool eflag, const bool vflag, const bool eatom,
|
||||
const bool vatom, int &host_start,
|
||||
const double cpu_time, bool &success) {
|
||||
this->acc_timers();
|
||||
if (nlist==0) {
|
||||
host_start=0;
|
||||
// Make sure textures are correct if realloc by a different hybrid style
|
||||
this->resize_atom(0,nall,success);
|
||||
this->zero_timers();
|
||||
return;
|
||||
}
|
||||
|
||||
int ago=this->hd_balancer.ago_first(f_ago);
|
||||
int inum=this->hd_balancer.balance(ago,nlocal,cpu_time);
|
||||
this->ans->inum(inum);
|
||||
#ifdef THREE_CONCURRENT
|
||||
this->ans2->inum(inum);
|
||||
#endif
|
||||
host_start=inum;
|
||||
|
||||
if (ago==0) {
|
||||
this->reset_nbors(nall, inum, nlist, ilist, numj, firstneigh, success);
|
||||
if (!success)
|
||||
return;
|
||||
_max_nbors = this->nbor->max_nbor_loop(nlist,numj,ilist);
|
||||
}
|
||||
|
||||
this->atom->cast_x_data(host_x,host_type);
|
||||
this->hd_balancer.start_timer();
|
||||
this->atom->add_x_data(host_x,host_type);
|
||||
|
||||
// re-allocate zetaij if necessary
|
||||
if (nall*_max_nbors > _zetaij.cols()) {
|
||||
int _nmax=static_cast<int>(static_cast<double>(nall)*1.10);
|
||||
_zetaij.resize(_max_nbors*_nmax);
|
||||
}
|
||||
|
||||
int _eflag;
|
||||
if (eflag)
|
||||
_eflag=1;
|
||||
else
|
||||
_eflag=0;
|
||||
|
||||
int ainum=nall;
|
||||
int nbor_pitch=this->nbor->nbor_pitch();
|
||||
int BX=this->block_pair();
|
||||
int GX=static_cast<int>(ceil(static_cast<double>(ainum)/
|
||||
(BX/(JTHREADS*KTHREADS))));
|
||||
|
||||
this->k_zeta.set_size(GX,BX);
|
||||
this->k_zeta.run(&this->atom->x, &ts1, &ts2, &ts3, &ts4, &ts5, &cutsq,
|
||||
&map, &elem2param, &_nelements, &_nparams, &_zetaij,
|
||||
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
|
||||
&_eflag, &nall, &ainum, &nbor_pitch, &this->_threads_per_atom);
|
||||
|
||||
int evatom=0;
|
||||
if (eatom || vatom)
|
||||
evatom=1;
|
||||
#ifdef THREE_CONCURRENT
|
||||
this->ucl_device->sync();
|
||||
#endif
|
||||
loop(eflag,vflag,evatom);
|
||||
this->ans->copy_answers(eflag,vflag,eatom,vatom,ilist);
|
||||
this->device->add_ans_object(this->ans);
|
||||
#ifdef THREE_CONCURRENT
|
||||
this->ans2->copy_answers(eflag,vflag,eatom,vatom,ilist);
|
||||
this->device->add_ans_object(this->ans2);
|
||||
#endif
|
||||
this->hd_balancer.stop_timer();
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// Reneighbor on GPU if necessary and then compute forces, virials, energies
|
||||
// ---------------------------------------------------------------------------
|
||||
template <class numtyp, class acctyp>
|
||||
int ** TersoffMT::compute(const int ago, const int inum_full,
|
||||
const int nall, double **host_x, int *host_type,
|
||||
double *sublo, double *subhi, tagint *tag,
|
||||
int **nspecial, tagint **special, const bool eflag,
|
||||
const bool vflag, const bool eatom,
|
||||
const bool vatom, int &host_start,
|
||||
int **ilist, int **jnum,
|
||||
const double cpu_time, bool &success) {
|
||||
this->acc_timers();
|
||||
|
||||
if (inum_full==0) {
|
||||
host_start=0;
|
||||
// Make sure textures are correct if realloc by a different hybrid style
|
||||
this->resize_atom(0,nall,success);
|
||||
this->zero_timers();
|
||||
return NULL;
|
||||
}
|
||||
|
||||
this->hd_balancer.balance(cpu_time);
|
||||
int inum=this->hd_balancer.get_gpu_count(ago,inum_full);
|
||||
this->ans->inum(inum);
|
||||
#ifdef THREE_CONCURRENT
|
||||
this->ans2->inum(inum);
|
||||
#endif
|
||||
host_start=inum;
|
||||
|
||||
// Build neighbor list on GPU if necessary
|
||||
if (ago==0) {
|
||||
_max_nbors = this->build_nbor_list(inum, inum_full-inum, nall, host_x, host_type,
|
||||
sublo, subhi, tag, nspecial, special, success);
|
||||
if (!success)
|
||||
return NULL;
|
||||
this->hd_balancer.start_timer();
|
||||
} else {
|
||||
this->atom->cast_x_data(host_x,host_type);
|
||||
this->hd_balancer.start_timer();
|
||||
this->atom->add_x_data(host_x,host_type);
|
||||
}
|
||||
*ilist=this->nbor->host_ilist.begin();
|
||||
*jnum=this->nbor->host_acc.begin();
|
||||
|
||||
// re-allocate zetaij if necessary
|
||||
if (nall*_max_nbors > _zetaij.cols()) {
|
||||
int _nmax=static_cast<int>(static_cast<double>(nall)*1.10);
|
||||
_zetaij.resize(_max_nbors*_nmax);
|
||||
}
|
||||
|
||||
int _eflag;
|
||||
if (eflag)
|
||||
_eflag=1;
|
||||
else
|
||||
_eflag=0;
|
||||
|
||||
int ainum=nall;
|
||||
int nbor_pitch=this->nbor->nbor_pitch();
|
||||
int BX=this->block_pair();
|
||||
int GX=static_cast<int>(ceil(static_cast<double>(ainum)/
|
||||
(BX/(JTHREADS*KTHREADS))));
|
||||
|
||||
this->k_zeta.set_size(GX,BX);
|
||||
this->k_zeta.run(&this->atom->x, &ts1, &ts2, &ts3, &ts4, &ts5, &cutsq,
|
||||
&map, &elem2param, &_nelements, &_nparams, &_zetaij,
|
||||
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
|
||||
&_eflag, &nall, &ainum, &nbor_pitch, &this->_threads_per_atom);
|
||||
|
||||
int evatom=0;
|
||||
if (eatom || vatom)
|
||||
evatom=1;
|
||||
#ifdef THREE_CONCURRENT
|
||||
this->ucl_device->sync();
|
||||
#endif
|
||||
loop(eflag,vflag,evatom);
|
||||
this->ans->copy_answers(eflag,vflag,eatom,vatom);
|
||||
this->device->add_ans_object(this->ans);
|
||||
#ifdef THREE_CONCURRENT
|
||||
this->ans2->copy_answers(eflag,vflag,eatom,vatom);
|
||||
this->device->add_ans_object(this->ans2);
|
||||
#endif
|
||||
this->hd_balancer.stop_timer();
|
||||
|
||||
return this->nbor->host_jlist.begin()-host_start;
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// Calculate energies, forces, and torques
|
||||
// ---------------------------------------------------------------------------
|
||||
template <class numtyp, class acctyp>
|
||||
void TersoffMT::loop(const bool _eflag, const bool _vflag, const int evatom) {
|
||||
// Compute the block size and grid size to keep all cores busy
|
||||
int BX=this->block_pair();
|
||||
int eflag, vflag;
|
||||
if (_eflag)
|
||||
eflag=1;
|
||||
else
|
||||
eflag=0;
|
||||
|
||||
if (_vflag)
|
||||
vflag=1;
|
||||
else
|
||||
vflag=0;
|
||||
|
||||
int ainum=this->ans->inum();
|
||||
int nbor_pitch=this->nbor->nbor_pitch();
|
||||
int GX=static_cast<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, &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, &ts4, &ts5, &cutsq,
|
||||
&map, &elem2param, &_nelements, &_nparams, &_zetaij,
|
||||
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
|
||||
&end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum,
|
||||
&nbor_pitch, &this->_threads_per_atom);
|
||||
|
||||
} else {
|
||||
this->k_three_end.set_size(GX,BX);
|
||||
this->k_three_end.run(&this->atom->x, &ts1, &ts2, &ts4, &ts5, &cutsq,
|
||||
&map, &elem2param, &_nelements, &_nparams, &_zetaij,
|
||||
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
|
||||
&end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum,
|
||||
&nbor_pitch, &this->_threads_per_atom);
|
||||
}
|
||||
|
||||
this->time_pair.stop();
|
||||
}
|
||||
|
||||
template class TersoffMod<PRECISION,ACC_PRECISION>;
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,118 @@
|
|||
/***************************************************************************
|
||||
tersoff_mod.h
|
||||
-------------------
|
||||
Trung Dac Nguyen
|
||||
|
||||
Class for acceleration of the tersoff/mod pair style.
|
||||
|
||||
__________________________________________________________________________
|
||||
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
|
||||
__________________________________________________________________________
|
||||
|
||||
begin :
|
||||
email : ndactrung@gmail.com
|
||||
***************************************************************************/
|
||||
|
||||
#ifndef LAL_TERSOFF_MOD_H
|
||||
#define LAL_TERSOFF_MOD_H
|
||||
|
||||
#include "lal_base_three.h"
|
||||
|
||||
namespace LAMMPS_AL {
|
||||
|
||||
template <class numtyp, class acctyp>
|
||||
class TersoffMod : public BaseThree<numtyp, acctyp> {
|
||||
public:
|
||||
TersoffMod();
|
||||
~TersoffMod();
|
||||
|
||||
/// Clear any previous data and set up for a new LAMMPS run for generic systems
|
||||
/** \param max_nbors initial number of rows in the neighbor matrix
|
||||
* \param cell_size cutoff + skin
|
||||
* \param gpu_split fraction of particles handled by device
|
||||
*
|
||||
* Returns:
|
||||
* - 0 if successfull
|
||||
* - -1 if fix gpu not found
|
||||
* - -3 if there is an out of memory error
|
||||
* - -4 if the GPU library was not compiled for GPU
|
||||
* - -5 Double precision is not supported on card **/
|
||||
int init(const int ntypes, const int nlocal, const int nall, const int max_nbors,
|
||||
const double cell_size, const double gpu_split, FILE *screen,
|
||||
int* host_map, const int nelements, int*** host_elem2param, const int nparams,
|
||||
const double* lam1, const double* lam2, const double* lam3,
|
||||
const double* powermint, const double* biga, const double* bigb,
|
||||
const double* bigr, const double* bigd, const double* c1, const double* c2,
|
||||
const double* c3, const double* c4, const double* c5,
|
||||
const double* h, const double* beta, const double* powern,
|
||||
const double* powern_del, const double* ca1, const double* cutsq);
|
||||
|
||||
/// Pair loop with host neighboring
|
||||
void compute(const int f_ago, const int inum_full, const int nall,
|
||||
const int nlist, double **host_x, int *host_type,
|
||||
int *ilist, int *numj, int **firstneigh, const bool eflag,
|
||||
const bool vflag, const bool eatom, const bool vatom,
|
||||
int &host_start, const double cpu_time, bool &success);
|
||||
|
||||
/// Pair loop with device neighboring
|
||||
int ** compute(const int ago, const int inum_full,
|
||||
const int nall, double **host_x, int *host_type, double *sublo,
|
||||
double *subhi, tagint *tag, int **nspecial,
|
||||
tagint **special, const bool eflag, const bool vflag,
|
||||
const bool eatom, const bool vatom, int &host_start,
|
||||
int **ilist, int **numj, const double cpu_time, bool &success);
|
||||
|
||||
/// Clear all host and device data
|
||||
/** \note This is called at the beginning of the init() routine **/
|
||||
void clear();
|
||||
|
||||
/// Returns memory usage on device per atom
|
||||
int bytes_per_atom(const int max_nbors) const;
|
||||
|
||||
/// Total host memory used by library for pair style
|
||||
double host_memory_usage() const;
|
||||
|
||||
// --------------------------- TYPE DATA --------------------------
|
||||
|
||||
/// If atom type constants fit in shared memory, use fast kernels
|
||||
bool shared_types;
|
||||
|
||||
/// Number of atom types
|
||||
int _lj_types;
|
||||
|
||||
/// ts1.x = lam1, ts1.y = lam2, ts1.z = lam3, ts1.w = powermint
|
||||
UCL_D_Vec<numtyp4> ts1;
|
||||
/// ts2.x = biga, ts2.y = bigb, ts2.z = bigr, ts2.w = bigd
|
||||
UCL_D_Vec<numtyp4> ts2;
|
||||
/// ts3.x = beta, ts3.y = powern, ts3.z = powern_del, ts3.w = ca1
|
||||
UCL_D_Vec<numtyp4> ts3;
|
||||
/// ts4.x = c1, ts4.y = c2, ts4.z = c3, ts4.w = c4
|
||||
UCL_D_Vec<numtyp4> ts4;
|
||||
/// ts5.x = c5, ts5.y = h
|
||||
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:
|
||||
/// zetaij.x = force, zetaij.y = prefactor, zetaij.z = evdwl,
|
||||
/// zetaij.w = zetaij
|
||||
UCL_D_Vec<acctyp4> _zetaij;
|
||||
|
||||
UCL_Kernel k_zeta;
|
||||
UCL_Texture ts1_tex, ts2_tex, ts3_tex, ts4_tex, ts5_tex;
|
||||
|
||||
int _max_nbors;
|
||||
|
||||
private:
|
||||
bool _allocated;
|
||||
void loop(const bool _eflag, const bool _vflag, const int evatom);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
|
|
@ -0,0 +1,135 @@
|
|||
/***************************************************************************
|
||||
tersoff_mod_ext.cpp
|
||||
-------------------
|
||||
Trung Dac Nguyen
|
||||
|
||||
Functions for LAMMPS access to tersoff acceleration routines.
|
||||
|
||||
__________________________________________________________________________
|
||||
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
|
||||
__________________________________________________________________________
|
||||
|
||||
begin : Thu April 17, 2014
|
||||
email : ndactrung@gmail.com
|
||||
***************************************************************************/
|
||||
|
||||
#include <iostream>
|
||||
#include <cassert>
|
||||
#include <math.h>
|
||||
|
||||
#include "lal_tersoff_mod.h"
|
||||
|
||||
using namespace std;
|
||||
using namespace LAMMPS_AL;
|
||||
|
||||
static TersoffMod<PRECISION,ACC_PRECISION> TSMMF;
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// Allocate memory on host and device and copy constants to device
|
||||
// ---------------------------------------------------------------------------
|
||||
int tersoff_mod_gpu_init(const int ntypes, const int inum, const int nall,
|
||||
const int max_nbors, const double cell_size, int &gpu_mode, FILE *screen,
|
||||
int* host_map, const int nelements, int*** host_elem2param, const int nparams,
|
||||
const double* ts_lam1, const double* ts_lam2, const double* ts_lam3,
|
||||
const double* ts_powermint, const double* ts_biga, const double* ts_bigb,
|
||||
const double* ts_bigr, const double* ts_bigd, const double* ts_c1,
|
||||
const double* ts_c2, const double* ts_c3, const double* ts_c4,
|
||||
const double* ts_c5, const double* ts_h, const double* ts_beta,
|
||||
const double* ts_powern, const double* ts_powern_del,
|
||||
const double* ts_ca1, const double* ts_cutsq) {
|
||||
TSMMF.clear();
|
||||
gpu_mode=TSMMF.device->gpu_mode();
|
||||
double gpu_split=TSMMF.device->particle_split();
|
||||
int first_gpu=TSMMF.device->first_device();
|
||||
int last_gpu=TSMMF.device->last_device();
|
||||
int world_me=TSMMF.device->world_me();
|
||||
int gpu_rank=TSMMF.device->gpu_rank();
|
||||
int procs_per_gpu=TSMMF.device->procs_per_gpu();
|
||||
|
||||
// disable host/device split for now
|
||||
if (gpu_split != 1.0)
|
||||
return -8;
|
||||
|
||||
TSMMF.device->init_message(screen,"tersoff/mod/gpu",first_gpu,last_gpu);
|
||||
|
||||
bool message=false;
|
||||
if (TSMMF.device->replica_me()==0 && screen)
|
||||
message=true;
|
||||
|
||||
if (message) {
|
||||
fprintf(screen,"Initializing Device and compiling on process 0...");
|
||||
fflush(screen);
|
||||
}
|
||||
|
||||
int init_ok=0;
|
||||
if (world_me==0)
|
||||
init_ok=TSMMF.init(ntypes, inum, nall, 300, cell_size, gpu_split, screen,
|
||||
host_map, nelements, host_elem2param, nparams,
|
||||
ts_lam1, ts_lam2, ts_lam3, ts_powermint,
|
||||
ts_biga, ts_bigb, ts_bigr, ts_bigd, ts_c1, ts_c2,
|
||||
ts_c3, ts_c4, ts_c5, ts_h, ts_beta, ts_powern,
|
||||
ts_powern_del, ts_ca1, ts_cutsq);
|
||||
|
||||
TSMMF.device->world_barrier();
|
||||
if (message)
|
||||
fprintf(screen,"Done.\n");
|
||||
|
||||
for (int i=0; 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=TSMMF.init(ntypes, inum, nall, 300, cell_size, gpu_split, screen,
|
||||
host_map, nelements, host_elem2param, nparams,
|
||||
ts_lam1, ts_lam2, ts_lam3, ts_powermint,
|
||||
ts_biga, ts_bigb, ts_bigr, ts_bigd, ts_c1, ts_c2,
|
||||
ts_c3, ts_c4, ts_c5, ts_h, ts_beta, ts_powern,
|
||||
ts_powern_del, ts_ca1, ts_cutsq);
|
||||
|
||||
TSMMF.device->gpu_barrier();
|
||||
if (message)
|
||||
fprintf(screen,"Done.\n");
|
||||
}
|
||||
if (message)
|
||||
fprintf(screen,"\n");
|
||||
|
||||
if (init_ok==0)
|
||||
TSMMF.estimate_gpu_overhead();
|
||||
return init_ok;
|
||||
}
|
||||
|
||||
void tersoff_mod_gpu_clear() {
|
||||
TSMMF.clear();
|
||||
}
|
||||
|
||||
int ** tersoff_mod_gpu_compute_n(const int ago, const int inum_full,
|
||||
const int nall, double **host_x, int *host_type,
|
||||
double *sublo, double *subhi, tagint *tag, int **nspecial,
|
||||
tagint **special, const bool eflag, const bool vflag,
|
||||
const bool eatom, const bool vatom, int &host_start,
|
||||
int **ilist, int **jnum, const double cpu_time,
|
||||
bool &success) {
|
||||
return TSMMF.compute(ago, inum_full, nall, host_x, host_type, sublo,
|
||||
subhi, tag, nspecial, special, eflag, vflag, eatom,
|
||||
vatom, host_start, ilist, jnum, cpu_time, success);
|
||||
}
|
||||
|
||||
void tersoff_mod_gpu_compute(const int ago, const int nlocal, const int nall,
|
||||
const int nlist, double **host_x, int *host_type,
|
||||
int *ilist, int *numj, int **firstneigh, const bool eflag,
|
||||
const bool vflag, const bool eatom, const bool vatom,
|
||||
int &host_start, const double cpu_time, bool &success) {
|
||||
TSMMF.compute(ago,nlocal,nall,nlist,host_x,host_type,ilist,numj,
|
||||
firstneigh,eflag,vflag,eatom,vatom,host_start,cpu_time,success);
|
||||
}
|
||||
|
||||
double tersoff_mod_gpu_bytes() {
|
||||
return TSMMF.host_memory_usage();
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,627 @@
|
|||
/// **************************************************************************
|
||||
// tersoff_mod_extra.h
|
||||
// -------------------
|
||||
// Trung Dac Nguyen
|
||||
//
|
||||
// Device code for Tersoff math routines
|
||||
//
|
||||
// __________________________________________________________________________
|
||||
// This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
|
||||
// __________________________________________________________________________
|
||||
//
|
||||
// begin :
|
||||
// email : ndactrung@gmail.com
|
||||
// ***************************************************************************/*
|
||||
|
||||
#ifndef LAL_TERSOFF_MOD_EXTRA_H
|
||||
#define LAL_TERSOFF_MOD_EXTRA_H
|
||||
|
||||
#ifdef NV_KERNEL
|
||||
#include "lal_aux_fun1.h"
|
||||
#else
|
||||
#endif
|
||||
|
||||
#define MY_PI2 (numtyp)1.57079632679489661923
|
||||
#define MY_PI4 (numtyp)0.78539816339744830962
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp vec3_dot(const numtyp x[3], const numtyp y[3])
|
||||
{
|
||||
return (x[0]*y[0] + x[1]*y[1] + x[2]*y[2]);
|
||||
}
|
||||
|
||||
ucl_inline void vec3_add(const numtyp x[3], const numtyp y[3], numtyp z[3])
|
||||
{
|
||||
z[0] = x[0]+y[0]; z[1] = x[1]+y[1]; z[2] = x[2]+y[2];
|
||||
}
|
||||
|
||||
ucl_inline void vec3_scale(const numtyp k, const numtyp x[3], numtyp y[3])
|
||||
{
|
||||
y[0] = k*x[0]; y[1] = k*x[1]; y[2] = k*x[2];
|
||||
}
|
||||
|
||||
ucl_inline void vec3_scaleadd(const numtyp k, const numtyp x[3],
|
||||
const numtyp y[3], numtyp z[3])
|
||||
{
|
||||
z[0] = k*x[0]+y[0]; z[1] = k*x[1]+y[1]; z[2] = k*x[2]+y[2];
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp ters_gijk_mod(const numtyp costheta,
|
||||
const numtyp param_c1,
|
||||
const numtyp param_c2,
|
||||
const numtyp param_c3,
|
||||
const numtyp param_c4,
|
||||
const numtyp param_c5,
|
||||
const numtyp param_h)
|
||||
{
|
||||
const numtyp tmp_h = (param_h - costheta)*(param_h - costheta);
|
||||
return param_c1 + (param_c2*tmp_h/(param_c3 + tmp_h)) *
|
||||
((numtyp)1.0 + param_c4*ucl_exp(-param_c5*tmp_h));
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp ters_gijk_d_mod(const numtyp costheta,
|
||||
const numtyp param_c2,
|
||||
const numtyp param_c3,
|
||||
const numtyp param_c4,
|
||||
const numtyp param_c5,
|
||||
const numtyp param_h)
|
||||
{
|
||||
const numtyp tmp_h = (param_h - costheta)*(param_h - costheta);
|
||||
const numtyp g1 = (param_h - costheta)/(param_c3 + tmp_h);
|
||||
const numtyp g2 = ucl_exp(-param_c5*tmp_h);
|
||||
return (numtyp)-2.0*param_c2*g1*((1 + param_c4*g2) *
|
||||
(1 + g1*(costheta - param_h)) - tmp_h*param_c4*param_c5*g2);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline void costheta_d(const numtyp rij_hat[3],
|
||||
const numtyp rij,
|
||||
const numtyp rik_hat[3],
|
||||
const numtyp rik,
|
||||
numtyp *dri,
|
||||
numtyp *drj,
|
||||
numtyp *drk)
|
||||
{
|
||||
// first element is derivative wrt Ri, second wrt Rj, third wrt Rk
|
||||
|
||||
numtyp cos_theta = vec3_dot(rij_hat,rik_hat);
|
||||
|
||||
vec3_scaleadd(-cos_theta,rij_hat,rik_hat,drj);
|
||||
vec3_scale(ucl_recip(rij),drj,drj);
|
||||
vec3_scaleadd(-cos_theta,rik_hat,rij_hat,drk);
|
||||
vec3_scale(ucl_recip(rik),drk,drk);
|
||||
vec3_add(drj,drk,dri);
|
||||
vec3_scale((numtyp)-1.0,dri,dri);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp ters_fc(const numtyp r,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd)
|
||||
{
|
||||
if (r < param_bigr-param_bigd) return (numtyp)1.0;
|
||||
if (r > param_bigr+param_bigd) return (numtyp)0.0;
|
||||
return (numtyp)0.5*((numtyp)1.0 -
|
||||
(numtyp)1.125*sin(MY_PI2*(r - param_bigr)/param_bigd) -
|
||||
(numtyp)0.125*sin(3*MY_PI2*(r - param_bigr)/param_bigd));
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp ters_fc_d(const numtyp r,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd)
|
||||
{
|
||||
if (r < param_bigr-param_bigd) return (numtyp)0.0;
|
||||
if (r > param_bigr+param_bigd) return (numtyp)0.0;
|
||||
return -((numtyp)0.375*MY_PI4/param_bigd) *
|
||||
((numtyp)3*cos(MY_PI2*(r - param_bigr)/param_bigd) +
|
||||
cos((numtyp)3*MY_PI2*(r - param_bigr)/param_bigd));
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp ters_fa(const numtyp r,
|
||||
const numtyp param_bigb,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_lam2)
|
||||
{
|
||||
if (r > param_bigr + param_bigd) return (numtyp)0.0;
|
||||
return -param_bigb * ucl_exp(-param_lam2 * r) *
|
||||
ters_fc(r,param_bigr,param_bigd);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp ters_fa_d(const numtyp r,
|
||||
const numtyp param_bigb,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_lam2)
|
||||
{
|
||||
if (r > param_bigr + param_bigd) return (numtyp)0.0;
|
||||
return param_bigb * ucl_exp(-param_lam2 * r) * (param_lam2 *
|
||||
ters_fc(r,param_bigr,param_bigd) - ters_fc_d(r,param_bigr,param_bigd));
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp ters_bij(const numtyp zeta,
|
||||
const numtyp param_beta,
|
||||
const numtyp param_powern,
|
||||
const numtyp param_powern_del,
|
||||
const numtyp param_ca1,
|
||||
const numtyp param_ca4)
|
||||
{
|
||||
numtyp tmp = param_beta * zeta;
|
||||
if (tmp > param_ca1)
|
||||
return ucl_powr(tmp, -param_powern/((numtyp)2.0*param_powern_del));
|
||||
if (tmp < param_ca4) return (numtyp)1.0;
|
||||
return ucl_powr((numtyp)1.0 + ucl_powr(tmp,param_powern),
|
||||
(numtyp)-1.0/((numtyp)2.0*param_powern_del));
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp ters_bij_d(const numtyp zeta,
|
||||
const numtyp param_beta,
|
||||
const numtyp param_powern,
|
||||
const numtyp param_powern_del,
|
||||
const numtyp param_ca1,
|
||||
const numtyp param_ca4)
|
||||
{
|
||||
numtyp tmp = param_beta * zeta;
|
||||
if (tmp > param_ca1) return (numtyp)-0.5*(param_powern/param_powern_del) *
|
||||
ucl_powr(tmp,(numtyp)-0.5*(param_powern/param_powern_del)) / zeta;
|
||||
if (tmp < param_ca4) return (numtyp)0.0;
|
||||
|
||||
numtyp tmp_n = ucl_powr(tmp,param_powern);
|
||||
return (numtyp)-0.5 *(param_powern/param_powern_del) *
|
||||
ucl_powr((numtyp)1.0+tmp_n, (numtyp)-1.0-((numtyp)1.0 /
|
||||
((numtyp)2.0*param_powern_del)))*tmp_n / zeta;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline void ters_zetaterm_d(const numtyp prefactor,
|
||||
const numtyp rij_hat[3],
|
||||
const numtyp rij,
|
||||
const numtyp rik_hat[3],
|
||||
const numtyp rik,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_h,
|
||||
const numtyp param_c1,
|
||||
const numtyp param_c2,
|
||||
const numtyp param_c3,
|
||||
const numtyp param_c4,
|
||||
const numtyp param_c5,
|
||||
numtyp dri[3],
|
||||
numtyp drj[3],
|
||||
numtyp drk[3])
|
||||
{
|
||||
numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp;
|
||||
numtyp dcosdri[3],dcosdrj[3],dcosdrk[3];
|
||||
|
||||
fc = ters_fc(rik,param_bigr,param_bigd);
|
||||
dfc = ters_fc_d(rik,param_bigr,param_bigd);
|
||||
|
||||
numtyp t = param_lam3*(rij-rik);
|
||||
if ((int)param_powermint == 3) tmp = t*t*t;
|
||||
else tmp = t;
|
||||
|
||||
if (tmp > (numtyp)69.0776) ex_delr = (numtyp)1.e30;
|
||||
else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
|
||||
else ex_delr = ucl_exp(tmp);
|
||||
|
||||
if ((int)param_powermint == 3)
|
||||
ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr;
|
||||
else ex_delr_d = param_lam3 * ex_delr;
|
||||
|
||||
cos_theta = vec3_dot(rij_hat,rik_hat);
|
||||
gijk = ters_gijk_mod(cos_theta,param_c1,param_c2,param_c3,param_c4,param_c5,param_h);
|
||||
gijk_d = ters_gijk_d_mod(cos_theta,param_c2,param_c3,param_c4,param_c5,param_h);
|
||||
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
|
||||
|
||||
// compute the derivative wrt Ri
|
||||
// dri = -dfc*gijk*ex_delr*rik_hat;
|
||||
// dri += fc*gijk_d*ex_delr*dcosdri;
|
||||
// dri += fc*gijk*ex_delr_d*(rik_hat - rij_hat);
|
||||
|
||||
vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri);
|
||||
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri);
|
||||
vec3_scaleadd(fc*gijk*ex_delr_d,rik_hat,dri,dri);
|
||||
vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri);
|
||||
vec3_scale(prefactor,dri,dri);
|
||||
|
||||
// compute the derivative wrt Rj
|
||||
// drj = fc*gijk_d*ex_delr*dcosdrj;
|
||||
// drj += fc*gijk*ex_delr_d*rij_hat;
|
||||
|
||||
vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj);
|
||||
vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj);
|
||||
vec3_scale(prefactor,drj,drj);
|
||||
|
||||
// compute the derivative wrt Rk
|
||||
// drk = dfc*gijk*ex_delr*rik_hat;
|
||||
// drk += fc*gijk_d*ex_delr*dcosdrk;
|
||||
// drk += -fc*gijk*ex_delr_d*rik_hat;
|
||||
|
||||
vec3_scale(dfc*gijk*ex_delr,rik_hat,drk);
|
||||
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk);
|
||||
vec3_scaleadd(-fc*gijk*ex_delr_d,rik_hat,drk,drk);
|
||||
vec3_scale(prefactor,drk,drk);
|
||||
}
|
||||
|
||||
ucl_inline void ters_zetaterm_d_fi(const numtyp prefactor,
|
||||
const numtyp rij_hat[3],
|
||||
const numtyp rij,
|
||||
const numtyp rik_hat[3],
|
||||
const numtyp rik,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_h,
|
||||
const numtyp param_c1,
|
||||
const numtyp param_c2,
|
||||
const numtyp param_c3,
|
||||
const numtyp param_c4,
|
||||
const numtyp param_c5,
|
||||
numtyp dri[3])
|
||||
{
|
||||
numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp;
|
||||
numtyp dcosdri[3],dcosdrj[3],dcosdrk[3];
|
||||
|
||||
fc = ters_fc(rik,param_bigr,param_bigd);
|
||||
dfc = ters_fc_d(rik,param_bigr,param_bigd);
|
||||
|
||||
numtyp t = param_lam3*(rij-rik);
|
||||
if ((int)param_powermint == 3) tmp = t*t*t;
|
||||
else tmp = t;
|
||||
|
||||
if (tmp > (numtyp)69.0776) ex_delr = (numtyp)1.e30;
|
||||
else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
|
||||
else ex_delr = ucl_exp(tmp);
|
||||
|
||||
if ((int)param_powermint == 3)
|
||||
ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr;
|
||||
else ex_delr_d = param_lam3 * ex_delr;
|
||||
|
||||
cos_theta = vec3_dot(rij_hat,rik_hat);
|
||||
gijk = ters_gijk_mod(cos_theta,param_c1,param_c2,param_c3,param_c4,param_c5,param_h);
|
||||
gijk_d = ters_gijk_d_mod(cos_theta,param_c2,param_c3,param_c4,param_c5,param_h);
|
||||
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
|
||||
|
||||
// compute the derivative wrt Ri
|
||||
// dri = -dfc*gijk*ex_delr*rik_hat;
|
||||
// dri += fc*gijk_d*ex_delr*dcosdri;
|
||||
// dri += fc*gijk*ex_delr_d*(rik_hat - rij_hat);
|
||||
|
||||
vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri);
|
||||
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri);
|
||||
vec3_scaleadd(fc*gijk*ex_delr_d,rik_hat,dri,dri);
|
||||
vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri);
|
||||
vec3_scale(prefactor,dri,dri);
|
||||
}
|
||||
|
||||
ucl_inline void ters_zetaterm_d_fj(const numtyp prefactor,
|
||||
const numtyp rij_hat[3],
|
||||
const numtyp rij,
|
||||
const numtyp rik_hat[3],
|
||||
const numtyp rik,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_h,
|
||||
const numtyp param_c1,
|
||||
const numtyp param_c2,
|
||||
const numtyp param_c3,
|
||||
const numtyp param_c4,
|
||||
const numtyp param_c5,
|
||||
numtyp drj[3])
|
||||
{
|
||||
numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,cos_theta,tmp;
|
||||
numtyp dcosdri[3],dcosdrj[3],dcosdrk[3];
|
||||
|
||||
fc = ters_fc(rik,param_bigr,param_bigd);
|
||||
|
||||
numtyp t = param_lam3*(rij-rik);
|
||||
if ((int)param_powermint == 3) tmp = t*t*t;
|
||||
else tmp = t;
|
||||
|
||||
if (tmp > (numtyp)69.0776) ex_delr = (numtyp)1.e30;
|
||||
else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
|
||||
else ex_delr = ucl_exp(tmp);
|
||||
|
||||
if ((int)param_powermint == 3)
|
||||
ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr;
|
||||
else ex_delr_d = param_lam3 * ex_delr;
|
||||
|
||||
cos_theta = vec3_dot(rij_hat,rik_hat);
|
||||
gijk = ters_gijk_mod(cos_theta,param_c1,param_c2,param_c3,param_c4,param_c5,param_h);
|
||||
gijk_d = ters_gijk_d_mod(cos_theta,param_c2,param_c3,param_c4,param_c5,param_h);
|
||||
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
|
||||
|
||||
// compute the derivative wrt Rj
|
||||
// drj = fc*gijk_d*ex_delr*dcosdrj;
|
||||
// drj += fc*gijk*ex_delr_d*rij_hat;
|
||||
|
||||
vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj);
|
||||
vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj);
|
||||
vec3_scale(prefactor,drj,drj);
|
||||
}
|
||||
|
||||
ucl_inline void ters_zetaterm_d_fk(const numtyp prefactor,
|
||||
const numtyp rij_hat[3],
|
||||
const numtyp rij,
|
||||
const numtyp rik_hat[3],
|
||||
const numtyp rik,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_h,
|
||||
const numtyp param_c1,
|
||||
const numtyp param_c2,
|
||||
const numtyp param_c3,
|
||||
const numtyp param_c4,
|
||||
const numtyp param_c5,
|
||||
numtyp drk[3])
|
||||
{
|
||||
numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp;
|
||||
numtyp dcosdri[3],dcosdrj[3],dcosdrk[3];
|
||||
|
||||
fc = ters_fc(rik,param_bigr,param_bigd);
|
||||
dfc = ters_fc_d(rik,param_bigr,param_bigd);
|
||||
|
||||
numtyp t = param_lam3*(rij-rik);
|
||||
if ((int)param_powermint == 3) tmp = t*t*t;
|
||||
else tmp = t;
|
||||
|
||||
if (tmp > (numtyp)69.0776) ex_delr = (numtyp)1.e30;
|
||||
else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
|
||||
else ex_delr = ucl_exp(tmp);
|
||||
|
||||
if ((int)param_powermint == 3)
|
||||
ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr;
|
||||
else ex_delr_d = param_lam3 * ex_delr;
|
||||
|
||||
cos_theta = vec3_dot(rij_hat,rik_hat);
|
||||
gijk = ters_gijk_mod(cos_theta,param_c1,param_c2,param_c3,param_c4,param_c5,param_h);
|
||||
gijk_d = ters_gijk_d_mod(cos_theta,param_c2,param_c3,param_c4,param_c5,param_h);
|
||||
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
|
||||
|
||||
// compute the derivative wrt Rk
|
||||
// drk = dfc*gijk*ex_delr*rik_hat;
|
||||
// drk += fc*gijk_d*ex_delr*dcosdrk;
|
||||
// drk += -fc*gijk*ex_delr_d*rik_hat;
|
||||
|
||||
vec3_scale(dfc*gijk*ex_delr,rik_hat,drk);
|
||||
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk);
|
||||
vec3_scaleadd(-fc*gijk*ex_delr_d,rik_hat,drk,drk);
|
||||
vec3_scale(prefactor,drk,drk);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline void repulsive(const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_lam1,
|
||||
const numtyp param_biga,
|
||||
const numtyp rsq,
|
||||
const int eflag,
|
||||
numtyp *ans)
|
||||
{
|
||||
numtyp r,tmp_fc,tmp_fc_d,tmp_exp;
|
||||
r = ucl_sqrt(rsq);
|
||||
tmp_fc = ters_fc(r,param_bigr,param_bigd);
|
||||
tmp_fc_d = ters_fc_d(r,param_bigr,param_bigd);
|
||||
tmp_exp = ucl_exp(-param_lam1 * r);
|
||||
// fforce
|
||||
ans[0] = -param_biga*tmp_exp*(tmp_fc_d - tmp_fc*param_lam1)*ucl_recip(r);
|
||||
// eng
|
||||
if (eflag) ans[1] = tmp_fc * param_biga * tmp_exp;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp zeta(const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_h,
|
||||
const numtyp param_c1,
|
||||
const numtyp param_c2,
|
||||
const numtyp param_c3,
|
||||
const numtyp param_c4,
|
||||
const numtyp param_c5,
|
||||
const numtyp rsqij,
|
||||
const numtyp rsqik,
|
||||
const numtyp4 delrij,
|
||||
const numtyp4 delrik)
|
||||
{
|
||||
numtyp rij,rik,costheta,arg,ex_delr;
|
||||
|
||||
rij = ucl_sqrt(rsqij);
|
||||
rik = ucl_sqrt(rsqik);
|
||||
costheta = (delrij.x*delrik.x + delrij.y*delrik.y +
|
||||
delrij.z*delrik.z) / (rij*rik);
|
||||
|
||||
numtyp t = param_lam3*(rij-rik);
|
||||
if ((int)param_powermint == 3) arg = t*t*t;
|
||||
else arg = t;
|
||||
|
||||
if (arg > (numtyp)69.0776) ex_delr = (numtyp)1.e30;
|
||||
else if (arg < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
|
||||
else ex_delr = ucl_exp(arg);
|
||||
|
||||
return ters_fc(rik,param_bigr,param_bigd) *
|
||||
ters_gijk_mod(costheta,param_c1,param_c2,param_c3,param_c4,param_c5,
|
||||
param_h) * ex_delr;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline void force_zeta(const numtyp param_bigb,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_lam2,
|
||||
const numtyp param_beta,
|
||||
const numtyp param_powern,
|
||||
const numtyp param_powern_del,
|
||||
const numtyp param_ca1,
|
||||
const numtyp param_ca4,
|
||||
const numtyp rsq,
|
||||
const numtyp zeta_ij,
|
||||
const int eflag,
|
||||
numtyp fpfeng[4])
|
||||
{
|
||||
numtyp r,fa,fa_d,bij;
|
||||
|
||||
r = ucl_sqrt(rsq);
|
||||
fa = ters_fa(r,param_bigb,param_bigr,param_bigd,param_lam2);
|
||||
fa_d = ters_fa_d(r,param_bigb,param_bigr,param_bigd,param_lam2);
|
||||
bij = ters_bij(zeta_ij,param_beta,param_powern,
|
||||
param_powern_del,param_ca1,param_ca4);
|
||||
fpfeng[0] = (numtyp)0.5*bij*fa_d * ucl_recip(r); // fforce
|
||||
fpfeng[1] = (numtyp)-0.5*fa * ters_bij_d(zeta_ij,param_beta, param_powern,
|
||||
param_powern_del,param_ca1,param_ca4); // prefactor
|
||||
if (eflag) fpfeng[2] = (numtyp)0.5*bij*fa; // eng
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
attractive term
|
||||
use param_ij cutoff for rij test
|
||||
use param_ijk cutoff for rik test
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline void attractive(const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_h,
|
||||
const numtyp param_c1,
|
||||
const numtyp param_c2,
|
||||
const numtyp param_c3,
|
||||
const numtyp param_c4,
|
||||
const numtyp param_c5,
|
||||
const numtyp prefactor,
|
||||
const numtyp rij,
|
||||
const numtyp rijinv,
|
||||
const numtyp rik,
|
||||
const numtyp rikinv,
|
||||
const numtyp delrij[3],
|
||||
const numtyp delrik[3],
|
||||
numtyp fi[3],
|
||||
numtyp fj[3],
|
||||
numtyp fk[3])
|
||||
{
|
||||
numtyp rij_hat[3],rik_hat[3];
|
||||
vec3_scale(rijinv,delrij,rij_hat);
|
||||
vec3_scale(rikinv,delrik,rik_hat);
|
||||
ters_zetaterm_d(prefactor,rij_hat,rij,rik_hat,rik,
|
||||
param_bigr, param_bigd, param_powermint, param_lam3,
|
||||
param_h, param_c1, param_c2, param_c3, param_c4, param_c5,
|
||||
fi, fj, fk);
|
||||
}
|
||||
|
||||
ucl_inline void attractive_fi(const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_h,
|
||||
const numtyp param_c1,
|
||||
const numtyp param_c2,
|
||||
const numtyp param_c3,
|
||||
const numtyp param_c4,
|
||||
const numtyp param_c5,
|
||||
const numtyp prefactor,
|
||||
const numtyp rij,
|
||||
const numtyp rijinv,
|
||||
const numtyp rik,
|
||||
const numtyp rikinv,
|
||||
const numtyp delrij[3],
|
||||
const numtyp delrik[3],
|
||||
numtyp fi[3])
|
||||
{
|
||||
numtyp rij_hat[3],rik_hat[3];
|
||||
vec3_scale(rijinv,delrij,rij_hat);
|
||||
vec3_scale(rikinv,delrik,rik_hat);
|
||||
ters_zetaterm_d_fi(prefactor,rij_hat,rij,rik_hat,rik,
|
||||
param_bigr, param_bigd, param_powermint, param_lam3,
|
||||
param_h, param_c1, param_c2, param_c3, param_c4, param_c5,
|
||||
fi);
|
||||
}
|
||||
|
||||
ucl_inline void attractive_fj(const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_h,
|
||||
const numtyp param_c1,
|
||||
const numtyp param_c2,
|
||||
const numtyp param_c3,
|
||||
const numtyp param_c4,
|
||||
const numtyp param_c5,
|
||||
const numtyp prefactor,
|
||||
const numtyp rij,
|
||||
const numtyp rijinv,
|
||||
const numtyp rik,
|
||||
const numtyp rikinv,
|
||||
const numtyp delrij[3],
|
||||
const numtyp delrik[3],
|
||||
numtyp fj[3])
|
||||
{
|
||||
numtyp rij_hat[3],rik_hat[3];
|
||||
vec3_scale(rijinv,delrij,rij_hat);
|
||||
vec3_scale(rikinv,delrik,rik_hat);
|
||||
ters_zetaterm_d_fj(prefactor,rij_hat,rij,rik_hat,rik,
|
||||
param_bigr, param_bigd, param_powermint, param_lam3,
|
||||
param_h, param_c1, param_c2, param_c3, param_c4, param_c5,
|
||||
fj);
|
||||
}
|
||||
|
||||
ucl_inline void attractive_fk(const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_h,
|
||||
const numtyp param_c1,
|
||||
const numtyp param_c2,
|
||||
const numtyp param_c3,
|
||||
const numtyp param_c4,
|
||||
const numtyp param_c5,
|
||||
const numtyp prefactor,
|
||||
const numtyp rij,
|
||||
const numtyp rijinv,
|
||||
const numtyp rik,
|
||||
const numtyp rikinv,
|
||||
const numtyp delrij[3],
|
||||
const numtyp delrik[3],
|
||||
numtyp fk[3])
|
||||
{
|
||||
numtyp rij_hat[3],rik_hat[3];
|
||||
vec3_scale(rijinv,delrij,rij_hat);
|
||||
vec3_scale(rikinv,delrik,rik_hat);
|
||||
ters_zetaterm_d_fk(prefactor,rij_hat,rij,rik_hat,rik,
|
||||
param_bigr, param_bigd, param_powermint, param_lam3,
|
||||
param_h, param_c1, param_c2, param_c3, param_c4, param_c5,
|
||||
fk);
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,482 @@
|
|||
/***************************************************************************
|
||||
tersoff_zbl.cpp
|
||||
-------------------
|
||||
Trung Dac Nguyen
|
||||
|
||||
Class for acceleration of the tersoff/zbl pair style.
|
||||
|
||||
__________________________________________________________________________
|
||||
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
|
||||
__________________________________________________________________________
|
||||
|
||||
begin :
|
||||
email : ndactrung@gmail.com
|
||||
***************************************************************************/
|
||||
|
||||
#if defined(USE_OPENCL)
|
||||
#include "tersoff_zbl_cl.h"
|
||||
#elif defined(USE_CUDART)
|
||||
const char *tersoff_zbl=0;
|
||||
#else
|
||||
#include "tersoff_zbl_cubin.h"
|
||||
#endif
|
||||
|
||||
#include "lal_tersoff_zbl.h"
|
||||
#include <cassert>
|
||||
using namespace LAMMPS_AL;
|
||||
#define TersoffZT TersoffZBL<numtyp, acctyp>
|
||||
|
||||
extern Device<PRECISION,ACC_PRECISION> device;
|
||||
|
||||
template <class numtyp, class acctyp>
|
||||
TersoffZT::TersoffZBL() : BaseThree<numtyp,acctyp>(), _allocated(false) {
|
||||
}
|
||||
|
||||
template <class numtyp, class acctyp>
|
||||
TersoffZT::~TersoffZBL() {
|
||||
clear();
|
||||
}
|
||||
|
||||
template <class numtyp, class acctyp>
|
||||
int TersoffZT::bytes_per_atom(const int max_nbors) const {
|
||||
return this->bytes_per_atom_atomic(max_nbors);
|
||||
}
|
||||
|
||||
template <class numtyp, class acctyp>
|
||||
int TersoffZT::init(const int ntypes, const int nlocal, const int nall,
|
||||
const int max_nbors, const double cell_size,
|
||||
const double gpu_split, FILE *_screen, int* host_map,
|
||||
const int nelements, int*** host_elem2param,
|
||||
const int nparams, const double* lam1, const double* lam2,
|
||||
const double* lam3, const double* powermint,
|
||||
const double* biga, const double* bigb, const double* bigr,
|
||||
const double* bigd, const double* c1, const double* c2,
|
||||
const double* c3, const double* c4, const double* c,
|
||||
const double* d, const double* h, const double* gamma,
|
||||
const double* beta, const double* powern, const double* Z_i,
|
||||
const double* Z_j, const double* ZBLcut,
|
||||
const double* ZBLexpscale, const double global_e,
|
||||
const double global_a_0, const double global_epsilon_0,
|
||||
const double* host_cutsq)
|
||||
{
|
||||
int success;
|
||||
success=this->init_three(nlocal,nall,max_nbors,0,cell_size,gpu_split,
|
||||
_screen,tersoff_zbl,"k_tersoff_zbl_repulsive",
|
||||
"k_tersoff_zbl_three_center", "k_tersoff_zbl_three_end");
|
||||
if (success!=0)
|
||||
return success;
|
||||
|
||||
int ef_nall=nall;
|
||||
if (ef_nall==0)
|
||||
ef_nall=2000;
|
||||
_zetaij.alloc(ef_nall*max_nbors,*(this->ucl_device),UCL_READ_WRITE);
|
||||
|
||||
k_zeta.set_function(*(this->pair_program),"k_tersoff_zbl_zeta");
|
||||
|
||||
// If atom type constants fit in shared memory use fast kernel
|
||||
int lj_types=ntypes;
|
||||
shared_types=false;
|
||||
int max_shared_types=this->device->max_shared_types();
|
||||
if (lj_types<=max_shared_types && this->_block_size>=max_shared_types) {
|
||||
lj_types=max_shared_types;
|
||||
shared_types=true;
|
||||
}
|
||||
_lj_types=lj_types;
|
||||
|
||||
_nparams = nparams;
|
||||
_nelements = nelements;
|
||||
|
||||
UCL_H_Vec<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);
|
||||
|
||||
ts6.alloc(nparams,*(this->ucl_device),UCL_READ_ONLY);
|
||||
|
||||
for (int i=0; i<nparams; i++) {
|
||||
dview[i].x=static_cast<numtyp>(Z_i[i]);
|
||||
dview[i].y=static_cast<numtyp>(Z_j[i]);
|
||||
dview[i].z=static_cast<numtyp>(ZBLcut[i]);
|
||||
dview[i].w=static_cast<numtyp>(ZBLexpscale[i]);
|
||||
}
|
||||
|
||||
ucl_copy(ts6,dview,false);
|
||||
ts6_tex.get_texture(*(this->pair_program),"ts6_tex");
|
||||
ts6_tex.bind_float(ts6,4);
|
||||
|
||||
UCL_H_Vec<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 < ntypes; i++)
|
||||
dview_map[i] = host_map[i];
|
||||
|
||||
map.alloc(lj_types,*(this->ucl_device), UCL_READ_ONLY);
|
||||
ucl_copy(map,dview_map,false);
|
||||
|
||||
_global_e = global_e;
|
||||
_global_a_0 = global_a_0;
|
||||
_global_epsilon_0 = global_epsilon_0;
|
||||
|
||||
_allocated=true;
|
||||
this->_max_bytes=ts1.row_bytes()+ts2.row_bytes()+ts3.row_bytes()+
|
||||
ts4.row_bytes()+ts5.row_bytes()+cutsq.row_bytes()+
|
||||
map.row_bytes()+elem2param.row_bytes()+_zetaij.row_bytes();
|
||||
return 0;
|
||||
}
|
||||
|
||||
template <class numtyp, class acctyp>
|
||||
void TersoffZT::clear() {
|
||||
if (!_allocated)
|
||||
return;
|
||||
_allocated=false;
|
||||
|
||||
ts1.clear();
|
||||
ts2.clear();
|
||||
ts3.clear();
|
||||
ts4.clear();
|
||||
ts5.clear();
|
||||
ts6.clear();
|
||||
cutsq.clear();
|
||||
map.clear();
|
||||
elem2param.clear();
|
||||
_zetaij.clear();
|
||||
|
||||
k_zeta.clear();
|
||||
|
||||
this->clear_atomic();
|
||||
}
|
||||
|
||||
template <class numtyp, class acctyp>
|
||||
double TersoffZT::host_memory_usage() const {
|
||||
return this->host_memory_usage_atomic()+sizeof(TersoffZBL<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 TersoffZT::compute(const int f_ago, const int nlocal, const int nall,
|
||||
const int nlist, double **host_x, int *host_type,
|
||||
int *ilist, int *numj, int **firstneigh,
|
||||
const bool eflag, const bool vflag, const bool eatom,
|
||||
const bool vatom, int &host_start,
|
||||
const double cpu_time, bool &success) {
|
||||
this->acc_timers();
|
||||
if (nlist==0) {
|
||||
host_start=0;
|
||||
// Make sure textures are correct if realloc by a different hybrid style
|
||||
this->resize_atom(0,nall,success);
|
||||
this->zero_timers();
|
||||
return;
|
||||
}
|
||||
|
||||
int ago=this->hd_balancer.ago_first(f_ago);
|
||||
int inum=this->hd_balancer.balance(ago,nlocal,cpu_time);
|
||||
this->ans->inum(inum);
|
||||
#ifdef THREE_CONCURRENT
|
||||
this->ans2->inum(inum);
|
||||
#endif
|
||||
host_start=inum;
|
||||
|
||||
if (ago==0) {
|
||||
this->reset_nbors(nall, inum, nlist, ilist, numj, firstneigh, success);
|
||||
if (!success)
|
||||
return;
|
||||
_max_nbors = this->nbor->max_nbor_loop(nlist,numj,ilist);
|
||||
}
|
||||
|
||||
this->atom->cast_x_data(host_x,host_type);
|
||||
this->hd_balancer.start_timer();
|
||||
this->atom->add_x_data(host_x,host_type);
|
||||
|
||||
// re-allocate zetaij if necessary
|
||||
if (nall*_max_nbors > _zetaij.cols()) {
|
||||
int _nmax=static_cast<int>(static_cast<double>(nall)*1.10);
|
||||
_zetaij.resize(_max_nbors*_nmax);
|
||||
}
|
||||
|
||||
int _eflag;
|
||||
if (eflag)
|
||||
_eflag=1;
|
||||
else
|
||||
_eflag=0;
|
||||
|
||||
int ainum=nall;
|
||||
int nbor_pitch=this->nbor->nbor_pitch();
|
||||
int BX=this->block_pair();
|
||||
int GX=static_cast<int>(ceil(static_cast<double>(ainum)/
|
||||
(BX/(JTHREADS*KTHREADS))));
|
||||
|
||||
this->k_zeta.set_size(GX,BX);
|
||||
this->k_zeta.run(&this->atom->x, &ts1, &ts2, &ts3, &ts4, &ts5, &ts6, &cutsq,
|
||||
&map, &elem2param, &_nelements, &_nparams, &_zetaij,
|
||||
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
|
||||
&_eflag, &nall, &ainum, &nbor_pitch, &this->_threads_per_atom);
|
||||
|
||||
int evatom=0;
|
||||
if (eatom || vatom)
|
||||
evatom=1;
|
||||
#ifdef THREE_CONCURRENT
|
||||
this->ucl_device->sync();
|
||||
#endif
|
||||
loop(eflag,vflag,evatom);
|
||||
this->ans->copy_answers(eflag,vflag,eatom,vatom,ilist);
|
||||
this->device->add_ans_object(this->ans);
|
||||
#ifdef THREE_CONCURRENT
|
||||
this->ans2->copy_answers(eflag,vflag,eatom,vatom,ilist);
|
||||
this->device->add_ans_object(this->ans2);
|
||||
#endif
|
||||
this->hd_balancer.stop_timer();
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// Reneighbor on GPU if necessary and then compute forces, virials, energies
|
||||
// ---------------------------------------------------------------------------
|
||||
template <class numtyp, class acctyp>
|
||||
int ** TersoffZT::compute(const int ago, const int inum_full,
|
||||
const int nall, double **host_x, int *host_type,
|
||||
double *sublo, double *subhi, tagint *tag,
|
||||
int **nspecial, tagint **special, const bool eflag,
|
||||
const bool vflag, const bool eatom,
|
||||
const bool vatom, int &host_start,
|
||||
int **ilist, int **jnum,
|
||||
const double cpu_time, bool &success) {
|
||||
this->acc_timers();
|
||||
|
||||
if (inum_full==0) {
|
||||
host_start=0;
|
||||
// Make sure textures are correct if realloc by a different hybrid style
|
||||
this->resize_atom(0,nall,success);
|
||||
this->zero_timers();
|
||||
return NULL;
|
||||
}
|
||||
|
||||
this->hd_balancer.balance(cpu_time);
|
||||
int inum=this->hd_balancer.get_gpu_count(ago,inum_full);
|
||||
this->ans->inum(inum);
|
||||
#ifdef THREE_CONCURRENT
|
||||
this->ans2->inum(inum);
|
||||
#endif
|
||||
host_start=inum;
|
||||
|
||||
// Build neighbor list on GPU if necessary
|
||||
if (ago==0) {
|
||||
_max_nbors = this->build_nbor_list(inum, inum_full-inum, nall, host_x, host_type,
|
||||
sublo, subhi, tag, nspecial, special, success);
|
||||
if (!success)
|
||||
return NULL;
|
||||
this->hd_balancer.start_timer();
|
||||
} else {
|
||||
this->atom->cast_x_data(host_x,host_type);
|
||||
this->hd_balancer.start_timer();
|
||||
this->atom->add_x_data(host_x,host_type);
|
||||
}
|
||||
*ilist=this->nbor->host_ilist.begin();
|
||||
*jnum=this->nbor->host_acc.begin();
|
||||
|
||||
// re-allocate zetaij if necessary
|
||||
if (nall*_max_nbors > _zetaij.cols()) {
|
||||
int _nmax=static_cast<int>(static_cast<double>(nall)*1.10);
|
||||
_zetaij.resize(_max_nbors*_nmax);
|
||||
}
|
||||
|
||||
int _eflag;
|
||||
if (eflag)
|
||||
_eflag=1;
|
||||
else
|
||||
_eflag=0;
|
||||
|
||||
int ainum=nall;
|
||||
int nbor_pitch=this->nbor->nbor_pitch();
|
||||
int BX=this->block_pair();
|
||||
int GX=static_cast<int>(ceil(static_cast<double>(ainum)/
|
||||
(BX/(JTHREADS*KTHREADS))));
|
||||
|
||||
this->k_zeta.set_size(GX,BX);
|
||||
this->k_zeta.run(&this->atom->x, &ts1, &ts2, &ts3, &ts4, &ts5, &ts6, &cutsq,
|
||||
&map, &elem2param, &_nelements, &_nparams, &_zetaij,
|
||||
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
|
||||
&_eflag, &nall, &ainum, &nbor_pitch, &this->_threads_per_atom);
|
||||
|
||||
int evatom=0;
|
||||
if (eatom || vatom)
|
||||
evatom=1;
|
||||
#ifdef THREE_CONCURRENT
|
||||
this->ucl_device->sync();
|
||||
#endif
|
||||
loop(eflag,vflag,evatom);
|
||||
this->ans->copy_answers(eflag,vflag,eatom,vatom);
|
||||
this->device->add_ans_object(this->ans);
|
||||
#ifdef THREE_CONCURRENT
|
||||
this->ans2->copy_answers(eflag,vflag,eatom,vatom);
|
||||
this->device->add_ans_object(this->ans2);
|
||||
#endif
|
||||
this->hd_balancer.stop_timer();
|
||||
|
||||
return this->nbor->host_jlist.begin()-host_start;
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// Calculate energies, forces, and torques
|
||||
// ---------------------------------------------------------------------------
|
||||
template <class numtyp, class acctyp>
|
||||
void TersoffZT::loop(const bool _eflag, const bool _vflag, const int evatom) {
|
||||
// Compute the block size and grid size to keep all cores busy
|
||||
int BX=this->block_pair();
|
||||
int eflag, vflag;
|
||||
if (_eflag)
|
||||
eflag=1;
|
||||
else
|
||||
eflag=0;
|
||||
|
||||
if (_vflag)
|
||||
vflag=1;
|
||||
else
|
||||
vflag=0;
|
||||
|
||||
int ainum=this->ans->inum();
|
||||
int nbor_pitch=this->nbor->nbor_pitch();
|
||||
int GX=static_cast<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, &ts6,
|
||||
&_global_e, &_global_a_0, &_global_epsilon_0, &cutsq,
|
||||
&map, &elem2param, &_nelements, &_nparams,
|
||||
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
|
||||
&this->ans->force, &this->ans->engv,
|
||||
&eflag, &vflag, &ainum, &nbor_pitch,
|
||||
&this->_threads_per_atom);
|
||||
|
||||
BX=this->block_size();
|
||||
GX=static_cast<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, &ts4, &cutsq,
|
||||
&map, &elem2param, &_nelements, &_nparams, &_zetaij,
|
||||
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
|
||||
&this->ans->force, &this->ans->engv, &eflag, &vflag, &ainum,
|
||||
&nbor_pitch, &this->_threads_per_atom, &evatom);
|
||||
|
||||
Answer<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, &ts4, &cutsq,
|
||||
&map, &elem2param, &_nelements, &_nparams, &_zetaij,
|
||||
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
|
||||
&end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum,
|
||||
&nbor_pitch, &this->_threads_per_atom);
|
||||
|
||||
} else {
|
||||
this->k_three_end.set_size(GX,BX);
|
||||
this->k_three_end.run(&this->atom->x, &ts1, &ts2, &ts4, &cutsq,
|
||||
&map, &elem2param, &_nelements, &_nparams, &_zetaij,
|
||||
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
|
||||
&end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum,
|
||||
&nbor_pitch, &this->_threads_per_atom);
|
||||
}
|
||||
|
||||
this->time_pair.stop();
|
||||
}
|
||||
|
||||
template class TersoffZBL<PRECISION,ACC_PRECISION>;
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,123 @@
|
|||
/***************************************************************************
|
||||
tersoff_zbl.h
|
||||
-------------------
|
||||
Trung Dac Nguyen
|
||||
|
||||
Class for acceleration of the tersoff pair style.
|
||||
|
||||
__________________________________________________________________________
|
||||
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
|
||||
__________________________________________________________________________
|
||||
|
||||
begin :
|
||||
email : ndactrung@gmail.com
|
||||
***************************************************************************/
|
||||
|
||||
#ifndef LAL_TERSOFF_ZBL_H
|
||||
#define LAL_TERSOFF_ZBL_H
|
||||
|
||||
#include "lal_base_three.h"
|
||||
|
||||
namespace LAMMPS_AL {
|
||||
|
||||
template <class numtyp, class acctyp>
|
||||
class TersoffZBL : public BaseThree<numtyp, acctyp> {
|
||||
public:
|
||||
TersoffZBL();
|
||||
~TersoffZBL();
|
||||
|
||||
/// Clear any previous data and set up for a new LAMMPS run for generic systems
|
||||
/** \param max_nbors initial number of rows in the neighbor matrix
|
||||
* \param cell_size cutoff + skin
|
||||
* \param gpu_split fraction of particles handled by device
|
||||
*
|
||||
* Returns:
|
||||
* - 0 if successfull
|
||||
* - -1 if fix gpu not found
|
||||
* - -3 if there is an out of memory error
|
||||
* - -4 if the GPU library was not compiled for GPU
|
||||
* - -5 Double precision is not supported on card **/
|
||||
int init(const int ntypes, const int nlocal, const int nall, const int max_nbors,
|
||||
const double cell_size, const double gpu_split, FILE *screen,
|
||||
int* host_map, const int nelements, int*** host_elem2param, const int nparams,
|
||||
const double* lam1, const double* lam2, const double* lam3,
|
||||
const double* powermint, const double* biga, const double* bigb,
|
||||
const double* bigr, const double* bigd, const double* c1, const double* c2,
|
||||
const double* c3, const double* c4, const double* c, const double* d,
|
||||
const double* h, const double* gamma, const double* beta,
|
||||
const double* powern, const double* Z_i, const double* Z_j,
|
||||
const double* ZBLcut, const double* ZBLexpscale, const double global_e,
|
||||
const double global_a_0, const double global_epsilon_0, const double* cutsq);
|
||||
|
||||
/// Pair loop with host neighboring
|
||||
void compute(const int f_ago, const int inum_full, const int nall,
|
||||
const int nlist, double **host_x, int *host_type,
|
||||
int *ilist, int *numj, int **firstneigh, const bool eflag,
|
||||
const bool vflag, const bool eatom, const bool vatom,
|
||||
int &host_start, const double cpu_time, bool &success);
|
||||
|
||||
/// Pair loop with device neighboring
|
||||
int ** compute(const int ago, const int inum_full,
|
||||
const int nall, double **host_x, int *host_type, double *sublo,
|
||||
double *subhi, tagint *tag, int **nspecial,
|
||||
tagint **special, const bool eflag, const bool vflag,
|
||||
const bool eatom, const bool vatom, int &host_start,
|
||||
int **ilist, int **numj, const double cpu_time, bool &success);
|
||||
|
||||
/// Clear all host and device data
|
||||
/** \note This is called at the beginning of the init() routine **/
|
||||
void clear();
|
||||
|
||||
/// Returns memory usage on device per atom
|
||||
int bytes_per_atom(const int max_nbors) const;
|
||||
|
||||
/// Total host memory used by library for pair style
|
||||
double host_memory_usage() const;
|
||||
|
||||
// --------------------------- TYPE DATA --------------------------
|
||||
|
||||
/// If atom type constants fit in shared memory, use fast kernels
|
||||
bool shared_types;
|
||||
|
||||
/// Number of atom types
|
||||
int _lj_types;
|
||||
|
||||
/// ts1.x = lam1, ts1.y = lam2, ts1.z = lam3, ts1.w = powermint
|
||||
UCL_D_Vec<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;
|
||||
/// ts6.x = Z_i, ts6.y = Z_j, ts6.z = ZBLcut, ts6.w = ZBLexpscale
|
||||
UCL_D_Vec<numtyp4> ts6;
|
||||
|
||||
UCL_D_Vec<numtyp> cutsq;
|
||||
|
||||
UCL_D_Vec<int> elem2param;
|
||||
UCL_D_Vec<int> map;
|
||||
int _nparams,_nelements;
|
||||
|
||||
/// Per-atom arrays:
|
||||
/// zetaij.x = force, zetaij.y = prefactor, zetaij.z = evdwl,
|
||||
/// zetaij.w = zetaij
|
||||
UCL_D_Vec<acctyp4> _zetaij;
|
||||
|
||||
UCL_Kernel k_zeta;
|
||||
UCL_Texture ts1_tex, ts2_tex, ts3_tex, ts4_tex, ts5_tex, ts6_tex;
|
||||
|
||||
int _max_nbors;
|
||||
numtyp _global_e,_global_a_0,_global_epsilon_0;
|
||||
|
||||
private:
|
||||
bool _allocated;
|
||||
void loop(const bool _eflag, const bool _vflag, const int evatom);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
|
|
@ -0,0 +1,146 @@
|
|||
/***************************************************************************
|
||||
tersoff_zbl_ext.cpp
|
||||
-------------------
|
||||
Trung Dac Nguyen
|
||||
|
||||
Functions for LAMMPS access to tersoff/zbl acceleration routines.
|
||||
|
||||
__________________________________________________________________________
|
||||
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
|
||||
__________________________________________________________________________
|
||||
|
||||
begin :
|
||||
email : ndactrung@gmail.com
|
||||
***************************************************************************/
|
||||
|
||||
#include <iostream>
|
||||
#include <cassert>
|
||||
#include <math.h>
|
||||
|
||||
#include "lal_tersoff_zbl.h"
|
||||
|
||||
using namespace std;
|
||||
using namespace LAMMPS_AL;
|
||||
|
||||
static TersoffZBL<PRECISION,ACC_PRECISION> TSZMF;
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// Allocate memory on host and device and copy constants to device
|
||||
// ---------------------------------------------------------------------------
|
||||
int tersoff_zbl_gpu_init(const int ntypes, const int inum, const int nall,
|
||||
const int max_nbors, const double cell_size, int &gpu_mode,
|
||||
FILE *screen, int* host_map, const int nelements,
|
||||
int*** host_elem2param, const int nparams,
|
||||
const double* ts_lam1, const double* ts_lam2,
|
||||
const double* ts_lam3, const double* ts_powermint,
|
||||
const double* ts_biga, const double* ts_bigb,
|
||||
const double* ts_bigr, const double* ts_bigd,
|
||||
const double* ts_c1, const double* ts_c2,
|
||||
const double* ts_c3, const double* ts_c4,
|
||||
const double* ts_c, const double* ts_d, const double* ts_h,
|
||||
const double* ts_gamma, const double* ts_beta,
|
||||
const double* ts_powern, const double* ts_Z_i,
|
||||
const double* ts_Z_j, const double* ts_ZBLcut,
|
||||
const double* ts_ZBLexpscale, const double global_e,
|
||||
const double global_a_0, const double global_epsilon_0,
|
||||
const double* ts_cutsq) {
|
||||
TSZMF.clear();
|
||||
gpu_mode=TSZMF.device->gpu_mode();
|
||||
double gpu_split=TSZMF.device->particle_split();
|
||||
int first_gpu=TSZMF.device->first_device();
|
||||
int last_gpu=TSZMF.device->last_device();
|
||||
int world_me=TSZMF.device->world_me();
|
||||
int gpu_rank=TSZMF.device->gpu_rank();
|
||||
int procs_per_gpu=TSZMF.device->procs_per_gpu();
|
||||
|
||||
// disable host/device split for now
|
||||
if (gpu_split != 1.0)
|
||||
return -8;
|
||||
|
||||
TSZMF.device->init_message(screen,"tersoff/zbl/gpu",first_gpu,last_gpu);
|
||||
|
||||
bool message=false;
|
||||
if (TSZMF.device->replica_me()==0 && screen)
|
||||
message=true;
|
||||
|
||||
if (message) {
|
||||
fprintf(screen,"Initializing Device and compiling on process 0...");
|
||||
fflush(screen);
|
||||
}
|
||||
|
||||
int init_ok=0;
|
||||
if (world_me==0)
|
||||
init_ok=TSZMF.init(ntypes, inum, nall, 300, cell_size, gpu_split, screen,
|
||||
host_map, nelements, host_elem2param, nparams,
|
||||
ts_lam1, ts_lam2, ts_lam3, ts_powermint,
|
||||
ts_biga, ts_bigb, ts_bigr, ts_bigd,
|
||||
ts_c1, ts_c2, ts_c3, ts_c4, ts_c, ts_d, ts_h,
|
||||
ts_gamma, ts_beta, ts_powern, ts_Z_i, ts_Z_j,
|
||||
ts_ZBLcut, ts_ZBLexpscale, global_e, global_a_0,
|
||||
global_epsilon_0, ts_cutsq);
|
||||
|
||||
TSZMF.device->world_barrier();
|
||||
if (message)
|
||||
fprintf(screen,"Done.\n");
|
||||
|
||||
for (int i=0; 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=TSZMF.init(ntypes, inum, nall, 300, cell_size, gpu_split, screen,
|
||||
host_map, nelements, host_elem2param, nparams,
|
||||
ts_lam1, ts_lam2, ts_lam3, ts_powermint,
|
||||
ts_biga, ts_bigb, ts_bigr, ts_bigd,
|
||||
ts_c1, ts_c2, ts_c3, ts_c4, ts_c, ts_d, ts_h,
|
||||
ts_gamma, ts_beta, ts_powern, ts_Z_i, ts_Z_j,
|
||||
ts_ZBLcut, ts_ZBLexpscale, global_e, global_a_0,
|
||||
global_epsilon_0, ts_cutsq);
|
||||
|
||||
TSZMF.device->gpu_barrier();
|
||||
if (message)
|
||||
fprintf(screen,"Done.\n");
|
||||
}
|
||||
if (message)
|
||||
fprintf(screen,"\n");
|
||||
|
||||
if (init_ok==0)
|
||||
TSZMF.estimate_gpu_overhead();
|
||||
return init_ok;
|
||||
}
|
||||
|
||||
void tersoff_zbl_gpu_clear() {
|
||||
TSZMF.clear();
|
||||
}
|
||||
|
||||
int ** tersoff_zbl_gpu_compute_n(const int ago, const int inum_full,
|
||||
const int nall, double **host_x, int *host_type,
|
||||
double *sublo, double *subhi, tagint *tag, int **nspecial,
|
||||
tagint **special, const bool eflag, const bool vflag,
|
||||
const bool eatom, const bool vatom, int &host_start,
|
||||
int **ilist, int **jnum, const double cpu_time,
|
||||
bool &success) {
|
||||
return TSZMF.compute(ago, inum_full, nall, host_x, host_type, sublo,
|
||||
subhi, tag, nspecial, special, eflag, vflag, eatom,
|
||||
vatom, host_start, ilist, jnum, cpu_time, success);
|
||||
}
|
||||
|
||||
void tersoff_zbl_gpu_compute(const int ago, const int nlocal, const int nall,
|
||||
const int nlist, double **host_x, int *host_type,
|
||||
int *ilist, int *numj, int **firstneigh, const bool eflag,
|
||||
const bool vflag, const bool eatom, const bool vatom,
|
||||
int &host_start, const double cpu_time, bool &success) {
|
||||
TSZMF.compute(ago,nlocal,nall,nlist,host_x,host_type,ilist,numj,
|
||||
firstneigh,eflag,vflag,eatom,vatom,host_start,cpu_time,success);
|
||||
}
|
||||
|
||||
double tersoff_zbl_gpu_bytes() {
|
||||
return TSZMF.host_memory_usage();
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,690 @@
|
|||
/// **************************************************************************
|
||||
// tersoff_zbl_extra.h
|
||||
// -------------------
|
||||
// Trung Dac Nguyen
|
||||
//
|
||||
// Device code for Tersoff math routines
|
||||
//
|
||||
// __________________________________________________________________________
|
||||
// This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
|
||||
// __________________________________________________________________________
|
||||
//
|
||||
// begin :
|
||||
// email : ndactrung@gmail.com
|
||||
// ***************************************************************************/*
|
||||
|
||||
#ifndef LAL_TERSOFF_ZBL_EXTRA_H
|
||||
#define LAL_TERSOFF_ZBL_EXTRA_H
|
||||
|
||||
#ifdef NV_KERNEL
|
||||
#include "lal_aux_fun1.h"
|
||||
#else
|
||||
#endif
|
||||
|
||||
#define MY_PI (numtyp)3.14159265358979323846
|
||||
#define MY_PI2 (numtyp)1.57079632679489661923
|
||||
#define MY_PI4 (numtyp)0.78539816339744830962
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp vec3_dot(const numtyp x[3], const numtyp y[3])
|
||||
{
|
||||
return (x[0]*y[0] + x[1]*y[1] + x[2]*y[2]);
|
||||
}
|
||||
|
||||
ucl_inline void vec3_add(const numtyp x[3], const numtyp y[3], numtyp z[3])
|
||||
{
|
||||
z[0] = x[0]+y[0]; z[1] = x[1]+y[1]; z[2] = x[2]+y[2];
|
||||
}
|
||||
|
||||
ucl_inline void vec3_scale(const numtyp k, const numtyp x[3], numtyp y[3])
|
||||
{
|
||||
y[0] = k*x[0]; y[1] = k*x[1]; y[2] = k*x[2];
|
||||
}
|
||||
|
||||
ucl_inline void vec3_scaleadd(const numtyp k, const numtyp x[3],
|
||||
const numtyp y[3], numtyp z[3])
|
||||
{
|
||||
z[0] = k*x[0]+y[0]; z[1] = k*x[1]+y[1]; z[2] = k*x[2]+y[2];
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp ters_gijk(const numtyp costheta,
|
||||
const numtyp param_c,
|
||||
const numtyp param_d,
|
||||
const numtyp param_h,
|
||||
const numtyp param_gamma)
|
||||
{
|
||||
const numtyp ters_c = param_c * param_c;
|
||||
const numtyp ters_d = param_d * param_d;
|
||||
const numtyp hcth = param_h - costheta;
|
||||
return param_gamma*((numtyp)1.0 + ters_c*ucl_recip(ters_d) -
|
||||
ters_c *ucl_recip(ters_d + hcth*hcth));
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp ters_gijk_d(const numtyp costheta,
|
||||
const numtyp param_c,
|
||||
const numtyp param_d,
|
||||
const numtyp param_h,
|
||||
const numtyp param_gamma)
|
||||
{
|
||||
const numtyp ters_c = param_c * param_c;
|
||||
const numtyp ters_d = param_d * param_d;
|
||||
const numtyp hcth = param_h - costheta;
|
||||
const numtyp numerator = (numtyp)-2.0 * ters_c * hcth;
|
||||
const numtyp denominator = ucl_recip(ters_d + hcth*hcth);
|
||||
return param_gamma*numerator*denominator*denominator;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline void costheta_d(const numtyp rij_hat[3],
|
||||
const numtyp rij,
|
||||
const numtyp rik_hat[3],
|
||||
const numtyp rik,
|
||||
numtyp *dri,
|
||||
numtyp *drj,
|
||||
numtyp *drk)
|
||||
{
|
||||
// first element is derivative wrt Ri, second wrt Rj, third wrt Rk
|
||||
|
||||
numtyp cos_theta = vec3_dot(rij_hat,rik_hat);
|
||||
|
||||
vec3_scaleadd(-cos_theta,rij_hat,rik_hat,drj);
|
||||
vec3_scale(ucl_recip(rij),drj,drj);
|
||||
vec3_scaleadd(-cos_theta,rik_hat,rij_hat,drk);
|
||||
vec3_scale(ucl_recip(rik),drk,drk);
|
||||
vec3_add(drj,drk,dri);
|
||||
vec3_scale((numtyp)-1.0,dri,dri);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp ters_fc(const numtyp r,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd)
|
||||
{
|
||||
if (r < param_bigr-param_bigd) return (numtyp)1.0;
|
||||
if (r > param_bigr+param_bigd) return (numtyp)0.0;
|
||||
return (numtyp)0.5*((numtyp)1.0 - sin(MY_PI2*(r - param_bigr)/param_bigd));
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp ters_fc_d(const numtyp r,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd)
|
||||
{
|
||||
if (r < param_bigr-param_bigd) return (numtyp)0.0;
|
||||
if (r > param_bigr+param_bigd) return (numtyp)0.0;
|
||||
return -(MY_PI4/param_bigd) * cos(MY_PI2*(r - param_bigr)/param_bigd);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp F_fermi(const numtyp r, const numtyp param_ZBLcut,
|
||||
const numtyp param_ZBLexpscale)
|
||||
{
|
||||
return ucl_recip((numtyp)1.0+ucl_exp(-param_ZBLexpscale*(r-param_ZBLcut)));
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp F_fermi_d(const numtyp r, const numtyp param_ZBLcut,
|
||||
const numtyp param_ZBLexpscale)
|
||||
{
|
||||
numtyp a = ucl_exp(-param_ZBLexpscale*(r-param_ZBLcut));
|
||||
numtyp b = (numtyp)1.0 + a;
|
||||
return param_ZBLexpscale*a*ucl_recip(b*b);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp ters_fa(const numtyp r,
|
||||
const numtyp param_bigb,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_lam2,
|
||||
const numtyp param_ZBLcut,
|
||||
const numtyp param_ZBLexpscale)
|
||||
{
|
||||
if (r > param_bigr + param_bigd) return (numtyp)0.0;
|
||||
return -param_bigb * ucl_exp(-param_lam2 * r) *
|
||||
ters_fc(r,param_bigr,param_bigd)*F_fermi(r,param_ZBLcut,param_ZBLexpscale);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp ters_fa_d(const numtyp r,
|
||||
const numtyp param_bigb,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_lam2,
|
||||
const numtyp param_ZBLcut,
|
||||
const numtyp param_ZBLexpscale)
|
||||
{
|
||||
if (r > param_bigr + param_bigd) return (numtyp)0.0;
|
||||
numtyp f = F_fermi(r,param_ZBLcut,param_ZBLexpscale);
|
||||
return param_bigb * ucl_exp(-param_lam2 * r) *
|
||||
(param_lam2 * ters_fc(r,param_bigr,param_bigd) * f -
|
||||
ters_fc_d(r,param_bigr,param_bigd) * f -
|
||||
ters_fc(r,param_bigr,param_bigd) * F_fermi_d(r,param_ZBLcut,param_ZBLexpscale));
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp ters_bij(const numtyp zeta,
|
||||
const numtyp param_beta,
|
||||
const numtyp param_powern,
|
||||
const numtyp param_c1,
|
||||
const numtyp param_c2,
|
||||
const numtyp param_c3,
|
||||
const numtyp param_c4)
|
||||
{
|
||||
numtyp tmp = param_beta * zeta;
|
||||
if (tmp > param_c1) return ucl_rsqrt(tmp);
|
||||
if (tmp > param_c2)
|
||||
return ((numtyp)1.0 - ucl_powr(tmp,-param_powern) /
|
||||
((numtyp)2.0*param_powern))*ucl_rsqrt(tmp);
|
||||
if (tmp < param_c4) return (numtyp)1.0;
|
||||
if (tmp < param_c3)
|
||||
return (numtyp)1.0 - ucl_powr(tmp,param_powern)/((numtyp)2.0*param_powern);
|
||||
return ucl_powr((numtyp)1.0 + ucl_powr(tmp,param_powern),
|
||||
(numtyp)-1.0/((numtyp)2.0*param_powern));
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp ters_bij_d(const numtyp zeta,
|
||||
const numtyp param_beta,
|
||||
const numtyp param_powern,
|
||||
const numtyp param_c1,
|
||||
const numtyp param_c2,
|
||||
const numtyp param_c3,
|
||||
const numtyp param_c4)
|
||||
{
|
||||
numtyp tmp = param_beta * zeta;
|
||||
if (tmp > param_c1)
|
||||
return param_beta * (numtyp)-0.5*ucl_powr(tmp,(numtyp)-1.5);
|
||||
if (tmp > param_c2)
|
||||
return param_beta * ((numtyp)-0.5*ucl_powr(tmp,(numtyp)-1.5) *
|
||||
// error in negligible 2nd term fixed 9/30/2015
|
||||
// (1.0 - 0.5*(1.0 + 1.0/(2.0*param->powern)) *
|
||||
((numtyp)1.0 - ((numtyp)1.0 + (numtyp)1.0 /((numtyp)2.0 * param_powern)) *
|
||||
ucl_powr(tmp,-param_powern)));
|
||||
if (tmp < param_c4) return (numtyp)0.0;
|
||||
if (tmp < param_c3)
|
||||
return (numtyp)-0.5*param_beta * ucl_powr(tmp,param_powern-(numtyp)1.0);
|
||||
|
||||
numtyp tmp_n = ucl_powr(tmp,param_powern);
|
||||
return (numtyp)-0.5 * ucl_powr((numtyp)1.0+tmp_n, (numtyp) -
|
||||
(numtyp)1.0-((numtyp)1.0 / ((numtyp)2.0 * param_powern)))*tmp_n / zeta;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline void ters_zetaterm_d(const numtyp prefactor,
|
||||
const numtyp rij_hat[3],
|
||||
const numtyp rij,
|
||||
const numtyp rik_hat[3],
|
||||
const numtyp rik,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_c,
|
||||
const numtyp param_d,
|
||||
const numtyp param_h,
|
||||
const numtyp param_gamma,
|
||||
numtyp dri[3],
|
||||
numtyp drj[3],
|
||||
numtyp drk[3])
|
||||
{
|
||||
numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp;
|
||||
numtyp dcosdri[3],dcosdrj[3],dcosdrk[3];
|
||||
|
||||
fc = ters_fc(rik,param_bigr,param_bigd);
|
||||
dfc = ters_fc_d(rik,param_bigr,param_bigd);
|
||||
|
||||
numtyp t = param_lam3*(rij-rik);
|
||||
if ((int)param_powermint == 3) tmp = t*t*t;
|
||||
else tmp = t;
|
||||
|
||||
if (tmp > (numtyp)69.0776) ex_delr = (numtyp)1.e30;
|
||||
else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
|
||||
else ex_delr = ucl_exp(tmp);
|
||||
|
||||
if ((int)param_powermint == 3)
|
||||
ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr;
|
||||
else ex_delr_d = param_lam3 * ex_delr;
|
||||
|
||||
cos_theta = vec3_dot(rij_hat,rik_hat);
|
||||
gijk = ters_gijk(cos_theta,param_c,param_d,param_h,param_gamma);
|
||||
gijk_d = ters_gijk_d(cos_theta,param_c,param_d,param_h,param_gamma);
|
||||
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
|
||||
|
||||
// compute the derivative wrt Ri
|
||||
// dri = -dfc*gijk*ex_delr*rik_hat;
|
||||
// dri += fc*gijk_d*ex_delr*dcosdri;
|
||||
// dri += fc*gijk*ex_delr_d*(rik_hat - rij_hat);
|
||||
|
||||
vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri);
|
||||
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri);
|
||||
vec3_scaleadd(fc*gijk*ex_delr_d,rik_hat,dri,dri);
|
||||
vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri);
|
||||
vec3_scale(prefactor,dri,dri);
|
||||
|
||||
// compute the derivative wrt Rj
|
||||
// drj = fc*gijk_d*ex_delr*dcosdrj;
|
||||
// drj += fc*gijk*ex_delr_d*rij_hat;
|
||||
|
||||
vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj);
|
||||
vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj);
|
||||
vec3_scale(prefactor,drj,drj);
|
||||
|
||||
// compute the derivative wrt Rk
|
||||
// drk = dfc*gijk*ex_delr*rik_hat;
|
||||
// drk += fc*gijk_d*ex_delr*dcosdrk;
|
||||
// drk += -fc*gijk*ex_delr_d*rik_hat;
|
||||
|
||||
vec3_scale(dfc*gijk*ex_delr,rik_hat,drk);
|
||||
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk);
|
||||
vec3_scaleadd(-fc*gijk*ex_delr_d,rik_hat,drk,drk);
|
||||
vec3_scale(prefactor,drk,drk);
|
||||
}
|
||||
|
||||
ucl_inline void ters_zetaterm_d_fi(const numtyp prefactor,
|
||||
const numtyp rij_hat[3],
|
||||
const numtyp rij,
|
||||
const numtyp rik_hat[3],
|
||||
const numtyp rik,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_c,
|
||||
const numtyp param_d,
|
||||
const numtyp param_h,
|
||||
const numtyp param_gamma,
|
||||
numtyp dri[3])
|
||||
{
|
||||
numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp;
|
||||
numtyp dcosdri[3],dcosdrj[3],dcosdrk[3];
|
||||
|
||||
fc = ters_fc(rik,param_bigr,param_bigd);
|
||||
dfc = ters_fc_d(rik,param_bigr,param_bigd);
|
||||
|
||||
numtyp t = param_lam3*(rij-rik);
|
||||
if ((int)param_powermint == 3) tmp = t*t*t;
|
||||
else tmp = t;
|
||||
|
||||
if (tmp > (numtyp)69.0776) ex_delr = (numtyp)1.e30;
|
||||
else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
|
||||
else ex_delr = ucl_exp(tmp);
|
||||
|
||||
if ((int)param_powermint == 3)
|
||||
ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr;
|
||||
else ex_delr_d = param_lam3 * ex_delr;
|
||||
|
||||
cos_theta = vec3_dot(rij_hat,rik_hat);
|
||||
gijk = ters_gijk(cos_theta,param_c,param_d,param_h,param_gamma);
|
||||
gijk_d = ters_gijk_d(cos_theta,param_c,param_d,param_h,param_gamma);
|
||||
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
|
||||
|
||||
// compute the derivative wrt Ri
|
||||
// dri = -dfc*gijk*ex_delr*rik_hat;
|
||||
// dri += fc*gijk_d*ex_delr*dcosdri;
|
||||
// dri += fc*gijk*ex_delr_d*(rik_hat - rij_hat);
|
||||
|
||||
vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri);
|
||||
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri);
|
||||
vec3_scaleadd(fc*gijk*ex_delr_d,rik_hat,dri,dri);
|
||||
vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri);
|
||||
vec3_scale(prefactor,dri,dri);
|
||||
}
|
||||
|
||||
ucl_inline void ters_zetaterm_d_fj(const numtyp prefactor,
|
||||
const numtyp rij_hat[3],
|
||||
const numtyp rij,
|
||||
const numtyp rik_hat[3],
|
||||
const numtyp rik,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_c,
|
||||
const numtyp param_d,
|
||||
const numtyp param_h,
|
||||
const numtyp param_gamma,
|
||||
numtyp drj[3])
|
||||
{
|
||||
numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,cos_theta,tmp;
|
||||
numtyp dcosdri[3],dcosdrj[3],dcosdrk[3];
|
||||
|
||||
fc = ters_fc(rik,param_bigr,param_bigd);
|
||||
|
||||
numtyp t = param_lam3*(rij-rik);
|
||||
if ((int)param_powermint == 3) tmp = t*t*t;
|
||||
else tmp = t;
|
||||
|
||||
if (tmp > (numtyp)69.0776) ex_delr = (numtyp)1.e30;
|
||||
else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
|
||||
else ex_delr = ucl_exp(tmp);
|
||||
|
||||
if ((int)param_powermint == 3)
|
||||
ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr;
|
||||
else ex_delr_d = param_lam3 * ex_delr;
|
||||
|
||||
cos_theta = vec3_dot(rij_hat,rik_hat);
|
||||
gijk = ters_gijk(cos_theta,param_c,param_d,param_h,param_gamma);
|
||||
gijk_d = ters_gijk_d(cos_theta,param_c,param_d,param_h,param_gamma);
|
||||
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
|
||||
|
||||
// compute the derivative wrt Rj
|
||||
// drj = fc*gijk_d*ex_delr*dcosdrj;
|
||||
// drj += fc*gijk*ex_delr_d*rij_hat;
|
||||
|
||||
vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj);
|
||||
vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj);
|
||||
vec3_scale(prefactor,drj,drj);
|
||||
}
|
||||
|
||||
ucl_inline void ters_zetaterm_d_fk(const numtyp prefactor,
|
||||
const numtyp rij_hat[3],
|
||||
const numtyp rij,
|
||||
const numtyp rik_hat[3],
|
||||
const numtyp rik,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_c,
|
||||
const numtyp param_d,
|
||||
const numtyp param_h,
|
||||
const numtyp param_gamma,
|
||||
numtyp drk[3])
|
||||
{
|
||||
numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp;
|
||||
numtyp dcosdri[3],dcosdrj[3],dcosdrk[3];
|
||||
|
||||
fc = ters_fc(rik,param_bigr,param_bigd);
|
||||
dfc = ters_fc_d(rik,param_bigr,param_bigd);
|
||||
|
||||
numtyp t = param_lam3*(rij-rik);
|
||||
if ((int)param_powermint == 3) tmp = t*t*t;
|
||||
else tmp = t;
|
||||
|
||||
if (tmp > (numtyp)69.0776) ex_delr = (numtyp)1.e30;
|
||||
else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
|
||||
else ex_delr = ucl_exp(tmp);
|
||||
|
||||
if ((int)param_powermint == 3)
|
||||
ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr;
|
||||
else ex_delr_d = param_lam3 * ex_delr;
|
||||
|
||||
cos_theta = vec3_dot(rij_hat,rik_hat);
|
||||
gijk = ters_gijk(cos_theta,param_c,param_d,param_h,param_gamma);
|
||||
gijk_d = ters_gijk_d(cos_theta,param_c,param_d,param_h,param_gamma);
|
||||
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
|
||||
|
||||
// compute the derivative wrt Rk
|
||||
// drk = dfc*gijk*ex_delr*rik_hat;
|
||||
// drk += fc*gijk_d*ex_delr*dcosdrk;
|
||||
// drk += -fc*gijk*ex_delr_d*rik_hat;
|
||||
|
||||
vec3_scale(dfc*gijk*ex_delr,rik_hat,drk);
|
||||
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk);
|
||||
vec3_scaleadd(-fc*gijk*ex_delr_d,rik_hat,drk,drk);
|
||||
vec3_scale(prefactor,drk,drk);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline void repulsive(const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_lam1,
|
||||
const numtyp param_biga,
|
||||
const numtyp param_Z_i,
|
||||
const numtyp param_Z_j,
|
||||
const numtyp param_ZBLcut,
|
||||
const numtyp param_ZBLexpscale,
|
||||
const numtyp global_e,
|
||||
const numtyp global_a_0,
|
||||
const numtyp global_epsilon_0,
|
||||
const numtyp rsq,
|
||||
const int eflag,
|
||||
numtyp *ans)
|
||||
{
|
||||
numtyp r,tmp_fc,tmp_fc_d,tmp_exp;
|
||||
|
||||
// Tersoff repulsive portion
|
||||
|
||||
r = ucl_sqrt(rsq);
|
||||
tmp_fc = ters_fc(r,param_bigr,param_bigd);
|
||||
tmp_fc_d = ters_fc_d(r,param_bigr,param_bigd);
|
||||
tmp_exp = ucl_exp(-param_lam1 * r);
|
||||
|
||||
numtyp fforce_ters = param_biga * tmp_exp * (tmp_fc_d - tmp_fc*param_lam1);
|
||||
numtyp eng_ters = tmp_fc * param_biga * tmp_exp;
|
||||
|
||||
// ZBL repulsive portion
|
||||
|
||||
numtyp esq = global_e*global_e;
|
||||
numtyp a_ij = ((numtyp)0.8854*global_a_0) /
|
||||
(ucl_powr(param_Z_i,(numtyp)0.23) + ucl_powr(param_Z_j,(numtyp)0.23));
|
||||
numtyp premult = (param_Z_i * param_Z_j * esq)/((numtyp)4.0*MY_PI*global_epsilon_0);
|
||||
numtyp r_ov_a = r/a_ij;
|
||||
numtyp t1 = (numtyp)0.1818*ucl_exp((numtyp)-3.2*r_ov_a);
|
||||
numtyp t2 = (numtyp)0.5099*ucl_exp((numtyp)-0.9423*r_ov_a);
|
||||
numtyp t3 = (numtyp)0.2802*ucl_exp((numtyp)-0.4029*r_ov_a);
|
||||
numtyp t4 = (numtyp)0.02817*ucl_exp((numtyp)-0.2016*r_ov_a);
|
||||
numtyp phi = t1 + t2 + t3 + t4;
|
||||
numtyp dphi = (numtyp)-3.2*t1 - (numtyp)0.9423*t2 - (numtyp)0.4029*t3 -
|
||||
(numtyp)0.2016*t4;
|
||||
dphi *= ucl_recip(a_ij);
|
||||
/*
|
||||
numtyp phi = (numtyp)0.1818*ucl_exp((numtyp)-3.2*r_ov_a) +
|
||||
(numtyp)0.5099*ucl_exp((numtyp)-0.9423*r_ov_a) +
|
||||
(numtyp)0.2802*ucl_exp((numtyp)-0.4029*r_ov_a) +
|
||||
(numtyp)0.02817*ucl_exp((numtyp)-0.2016*r_ov_a);
|
||||
numtyp dphi = ucl_recip(a_ij) * ((numtyp)-3.2*(numtyp)0.1818*ucl_exp((numtyp)-3.2*r_ov_a) -
|
||||
(numtyp)0.9423*(numtyp)0.5099*ucl_exp((numtyp)-0.9423*r_ov_a) -
|
||||
(numtyp)0.4029*(numtyp)0.2802*ucl_exp((numtyp)-0.4029*r_ov_a) -
|
||||
(numtyp)0.2016*(numtyp)0.02817*ucl_exp((numtyp)-0.2016*r_ov_a));
|
||||
*/
|
||||
numtyp rinv = ucl_recip(r);
|
||||
numtyp fforce_ZBL = premult*(-phi)/rsq + premult*dphi*rinv;
|
||||
numtyp eng_ZBL = premult*rinv*phi;
|
||||
|
||||
// combine two parts with smoothing by Fermi-like function
|
||||
// ans[0] = fforce
|
||||
numtyp f = F_fermi(r,param_ZBLcut,param_ZBLexpscale);
|
||||
numtyp f_d = F_fermi_d(r,param_ZBLcut,param_ZBLexpscale);
|
||||
ans[0] = -(-f_d * eng_ZBL + ((numtyp)1.0 - f)*fforce_ZBL + f_d*eng_ters +
|
||||
f*fforce_ters) * rinv;
|
||||
|
||||
// ans[1] = eng
|
||||
if (eflag) ans[1] = ((numtyp)1.0 - f)*eng_ZBL + f*eng_ters;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline numtyp zeta(const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_c,
|
||||
const numtyp param_d,
|
||||
const numtyp param_h,
|
||||
const numtyp param_gamma,
|
||||
const numtyp rsqij,
|
||||
const numtyp rsqik,
|
||||
const numtyp4 delrij,
|
||||
const numtyp4 delrik)
|
||||
{
|
||||
numtyp rij,rik,costheta,arg,ex_delr;
|
||||
|
||||
rij = ucl_sqrt(rsqij);
|
||||
rik = ucl_sqrt(rsqik);
|
||||
costheta = (delrij.x*delrik.x + delrij.y*delrik.y +
|
||||
delrij.z*delrik.z) / (rij*rik);
|
||||
|
||||
numtyp t = param_lam3*(rij-rik);
|
||||
if ((int)param_powermint == 3) arg = t*t*t;
|
||||
else arg = t;
|
||||
|
||||
if (arg > (numtyp)69.0776) ex_delr = (numtyp)1.e30;
|
||||
else if (arg < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
|
||||
else ex_delr = ucl_exp(arg);
|
||||
|
||||
return ters_fc(rik,param_bigr,param_bigd) *
|
||||
ters_gijk(costheta,param_c, param_d, param_h, param_gamma) * ex_delr;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline void force_zeta(const numtyp param_bigb,
|
||||
const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_lam2,
|
||||
const numtyp param_beta,
|
||||
const numtyp param_powern,
|
||||
const numtyp param_c1,
|
||||
const numtyp param_c2,
|
||||
const numtyp param_c3,
|
||||
const numtyp param_c4,
|
||||
const numtyp param_ZBLcut,
|
||||
const numtyp param_ZBLexpscale,
|
||||
const numtyp rsq,
|
||||
const numtyp zeta_ij,
|
||||
const int eflag,
|
||||
numtyp fpfeng[4])
|
||||
{
|
||||
numtyp r,fa,fa_d,bij;
|
||||
|
||||
r = ucl_sqrt(rsq);
|
||||
fa = ters_fa(r,param_bigb,param_bigr,param_bigd,param_lam2,param_ZBLcut,param_ZBLexpscale);
|
||||
fa_d = ters_fa_d(r,param_bigb,param_bigr,param_bigd,param_lam2,param_ZBLcut,param_ZBLexpscale);
|
||||
bij = ters_bij(zeta_ij,param_beta,param_powern,
|
||||
param_c1,param_c2, param_c3, param_c4);
|
||||
fpfeng[0] = (numtyp)0.5*bij*fa_d * ucl_recip(r); // fforce
|
||||
fpfeng[1] = (numtyp)-0.5*fa * ters_bij_d(zeta_ij,param_beta, param_powern,
|
||||
param_c1,param_c2, param_c3, param_c4); // prefactor
|
||||
if (eflag) fpfeng[2] = (numtyp)0.5*bij*fa; // eng
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
attractive term
|
||||
use param_ij cutoff for rij test
|
||||
use param_ijk cutoff for rik test
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
ucl_inline void attractive(const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_c,
|
||||
const numtyp param_d,
|
||||
const numtyp param_h,
|
||||
const numtyp param_gamma,
|
||||
const numtyp prefactor,
|
||||
const numtyp rij,
|
||||
const numtyp rijinv,
|
||||
const numtyp rik,
|
||||
const numtyp rikinv,
|
||||
const numtyp delrij[3],
|
||||
const numtyp delrik[3],
|
||||
numtyp fi[3],
|
||||
numtyp fj[3],
|
||||
numtyp fk[3])
|
||||
{
|
||||
numtyp rij_hat[3],rik_hat[3];
|
||||
vec3_scale(rijinv,delrij,rij_hat);
|
||||
vec3_scale(rikinv,delrik,rik_hat);
|
||||
ters_zetaterm_d(prefactor,rij_hat,rij,rik_hat,rik,
|
||||
param_bigr, param_bigd, param_powermint, param_lam3,
|
||||
param_c, param_d, param_h, param_gamma, fi, fj, fk);
|
||||
}
|
||||
|
||||
ucl_inline void attractive_fi(const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_c,
|
||||
const numtyp param_d,
|
||||
const numtyp param_h,
|
||||
const numtyp param_gamma,
|
||||
const numtyp prefactor,
|
||||
const numtyp rij,
|
||||
const numtyp rijinv,
|
||||
const numtyp rik,
|
||||
const numtyp rikinv,
|
||||
const numtyp delrij[3],
|
||||
const numtyp delrik[3],
|
||||
numtyp fi[3])
|
||||
{
|
||||
numtyp rij_hat[3],rik_hat[3];
|
||||
vec3_scale(rijinv,delrij,rij_hat);
|
||||
vec3_scale(rikinv,delrik,rik_hat);
|
||||
ters_zetaterm_d_fi(prefactor,rij_hat,rij,rik_hat,rik,
|
||||
param_bigr, param_bigd, param_powermint, param_lam3,
|
||||
param_c, param_d, param_h, param_gamma, fi);
|
||||
}
|
||||
|
||||
ucl_inline void attractive_fj(const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_c,
|
||||
const numtyp param_d,
|
||||
const numtyp param_h,
|
||||
const numtyp param_gamma,
|
||||
const numtyp prefactor,
|
||||
const numtyp rij,
|
||||
const numtyp rijinv,
|
||||
const numtyp rik,
|
||||
const numtyp rikinv,
|
||||
const numtyp delrij[3],
|
||||
const numtyp delrik[3],
|
||||
numtyp fj[3])
|
||||
{
|
||||
numtyp rij_hat[3],rik_hat[3];
|
||||
vec3_scale(rijinv,delrij,rij_hat);
|
||||
vec3_scale(rikinv,delrik,rik_hat);
|
||||
ters_zetaterm_d_fj(prefactor,rij_hat,rij,rik_hat,rik,
|
||||
param_bigr, param_bigd, param_powermint, param_lam3,
|
||||
param_c, param_d, param_h, param_gamma, fj);
|
||||
}
|
||||
|
||||
ucl_inline void attractive_fk(const numtyp param_bigr,
|
||||
const numtyp param_bigd,
|
||||
const numtyp param_powermint,
|
||||
const numtyp param_lam3,
|
||||
const numtyp param_c,
|
||||
const numtyp param_d,
|
||||
const numtyp param_h,
|
||||
const numtyp param_gamma,
|
||||
const numtyp prefactor,
|
||||
const numtyp rij,
|
||||
const numtyp rijinv,
|
||||
const numtyp rik,
|
||||
const numtyp rikinv,
|
||||
const numtyp delrij[3],
|
||||
const numtyp delrik[3],
|
||||
numtyp fk[3])
|
||||
{
|
||||
numtyp rij_hat[3],rik_hat[3];
|
||||
vec3_scale(rijinv,delrij,rij_hat);
|
||||
vec3_scale(rikinv,delrik,rik_hat);
|
||||
ters_zetaterm_d_fk(prefactor,rij_hat,rij,rik_hat,rik,
|
||||
param_bigr, param_bigd, param_powermint, param_lam3,
|
||||
param_c, param_d, param_h, param_gamma, fk);
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
Loading…
Reference in New Issue