Homec4science
Diffusion lammps (lammm-devel) Browse

lammps/src/KOKKOS/pair_tersoff_zbl_kokkos.cpp
3c8e75ad590alammm-devel

Last Change
Raw File

pair_tersoff_zbl_kokkos.cpp
View Options

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Ray Shan (SNL)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_tersoff_zbl_kokkos.h"
#include "kokkos.h"
#include "atom_kokkos.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_request.h"
#include "neigh_list_kokkos.h"
#include "update.h"
#include "integrate.h"
#include "respa.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "atom_masks.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define KOKKOS_CUDA_MAX_THREADS 256
#define KOKKOS_CUDA_MIN_BLOCKS 8
/* ---------------------------------------------------------------------- */
template<class DeviceType>
PairTersoffZBLKokkos<DeviceType>::PairTersoffZBLKokkos(LAMMPS *lmp) : PairTersoffZBL(lmp)
{
respa_enable = 0;
atomKK = (AtomKokkos *) atom;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = X_MASK | F_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK;
datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK;
if (strcmp(update->unit_style,"metal") == 0) {
global_a_0 = 0.529;
global_epsilon_0 = 0.00552635;
global_e = 1.0;
} else if (strcmp(update->unit_style,"real") == 0) {
global_a_0 = 0.529;
global_epsilon_0 = 0.00552635 * 0.043365121;
global_e = 1.0;
} else error->all(FLERR,"Pair tersoff/zbl/kk requires metal or real units");
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
PairTersoffZBLKokkos<DeviceType>::~PairTersoffZBLKokkos()
{
if (!copymode) {
memory->destroy_kokkos(k_eatom,eatom);
memory->destroy_kokkos(k_vatom,vatom);
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairTersoffZBLKokkos<DeviceType>::allocate()
{
PairTersoffZBL::allocate();
int n = atom->ntypes;
k_params = Kokkos::DualView<params_ters***,Kokkos::LayoutRight,DeviceType>
("PairTersoffZBL::paramskk",n+1,n+1,n+1);
paramskk = k_params.d_view;
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
template<class DeviceType>
void PairTersoffZBLKokkos<DeviceType>::init_style()
{
PairTersoffZBL::init_style();
// irequest = neigh request made by parent class
neighflag = lmp->kokkos->neighflag;
int irequest = neighbor->nrequest - 1;
neighbor->requests[irequest]->
kokkos_host = Kokkos::Impl::is_same<DeviceType,LMPHostType>::value &&
!Kokkos::Impl::is_same<DeviceType,LMPDeviceType>::value;
neighbor->requests[irequest]->
kokkos_device = Kokkos::Impl::is_same<DeviceType,LMPDeviceType>::value;
if (neighflag == FULL)
error->all(FLERR,"Cannot (yet) use full neighbor list style with tersoff/zbl/kk");
if (neighflag == FULL || neighflag == HALF || neighflag == HALFTHREAD) {
neighbor->requests[irequest]->full = 1;
neighbor->requests[irequest]->half = 0;
if (neighflag == FULL)
neighbor->requests[irequest]->ghost = 1;
else
neighbor->requests[irequest]->ghost = 0;
} else {
error->all(FLERR,"Cannot use chosen neighbor list style with tersoff/zbl/kk");
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairTersoffZBLKokkos<DeviceType>::setup_params()
{
PairTersoffZBL::setup_params();
int i,j,k,m;
int n = atom->ntypes;
for (i = 1; i <= n; i++)
for (j = 1; j <= n; j++)
for (k = 1; k <= n; k++) {
m = elem2param[map[i]][map[j]][map[k]];
k_params.h_view(i,j,k).powerm = params[m].powerm;
k_params.h_view(i,j,k).gamma = params[m].gamma;
k_params.h_view(i,j,k).lam3 = params[m].lam3;
k_params.h_view(i,j,k).c = params[m].c;
k_params.h_view(i,j,k).d = params[m].d;
k_params.h_view(i,j,k).h = params[m].h;
k_params.h_view(i,j,k).powern = params[m].powern;
k_params.h_view(i,j,k).beta = params[m].beta;
k_params.h_view(i,j,k).lam2 = params[m].lam2;
k_params.h_view(i,j,k).bigb = params[m].bigb;
k_params.h_view(i,j,k).bigr = params[m].bigr;
k_params.h_view(i,j,k).bigd = params[m].bigd;
k_params.h_view(i,j,k).lam1 = params[m].lam1;
k_params.h_view(i,j,k).biga = params[m].biga;
k_params.h_view(i,j,k).cutsq = params[m].cutsq;
k_params.h_view(i,j,k).c1 = params[m].c1;
k_params.h_view(i,j,k).c2 = params[m].c2;
k_params.h_view(i,j,k).c3 = params[m].c3;
k_params.h_view(i,j,k).c4 = params[m].c4;
k_params.h_view(i,j,k).Z_i = params[m].Z_i;
k_params.h_view(i,j,k).Z_j = params[m].Z_j;
k_params.h_view(i,j,k).ZBLcut = params[m].ZBLcut;
k_params.h_view(i,j,k).ZBLexpscale = params[m].ZBLexpscale;
}
k_params.template modify<LMPHostType>();
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairTersoffZBLKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
{
eflag = eflag_in;
vflag = vflag_in;
if (neighflag == FULL) no_virial_fdotr_compute = 1;
if (eflag || vflag) ev_setup(eflag,vflag,0);
else evflag = vflag_fdotr = 0;
// reallocate per-atom arrays if necessary
if (eflag_atom) {
memory->destroy_kokkos(k_eatom,eatom);
memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
d_eatom = k_eatom.view<DeviceType>();
}
if (vflag_atom) {
memory->destroy_kokkos(k_vatom,vatom);
memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
d_vatom = k_vatom.view<DeviceType>();
}
atomKK->sync(execution_space,datamask_read);
k_params.template sync<DeviceType>();
if (eflag || vflag) atomKK->modified(execution_space,datamask_modify);
else atomKK->modified(execution_space,F_MASK);
x = atomKK->k_x.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
tag = atomKK->k_tag.view<DeviceType>();
nlocal = atom->nlocal;
nall = atom->nlocal + atom->nghost;
newton_pair = force->newton_pair;
inum = list->inum;
const int ignum = inum + list->gnum;
NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(list);
d_numneigh = k_list->d_numneigh;
d_neighbors = k_list->d_neighbors;
d_ilist = k_list->d_ilist;
copymode = 1;
EV_FLOAT ev;
EV_FLOAT ev_all;
// build short neighbor list
int max_neighs = d_neighbors.dimension_1();
if ((d_neighbors_short.dimension_1() != max_neighs) ||
(d_neighbors_short.dimension_0() != ignum)) {
d_neighbors_short = Kokkos::View<int**,DeviceType>("Tersoff::neighbors_short",ignum,max_neighs);
}
if (d_numneigh_short.dimension_0()!=ignum)
d_numneigh_short = Kokkos::View<int*,DeviceType>("Tersoff::numneighs_short",ignum);
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagPairTersoffZBLComputeShortNeigh>(0,neighflag==FULL?ignum:inum), *this);
if (neighflag == HALF) {
if (evflag)
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagPairTersoffZBLComputeHalf<HALF,1> >(0,inum),*this,ev);
else
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairTersoffZBLComputeHalf<HALF,0> >(0,inum),*this);
ev_all += ev;
} else if (neighflag == HALFTHREAD) {
if (evflag)
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagPairTersoffZBLComputeHalf<HALFTHREAD,1> >(0,inum),*this,ev);
else
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairTersoffZBLComputeHalf<HALFTHREAD,0> >(0,inum),*this);
ev_all += ev;
} else if (neighflag == FULL) {
if (evflag)
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagPairTersoffZBLComputeFullA<FULL,1> >(0,inum),*this,ev);
else
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairTersoffZBLComputeFullA<FULL,0> >(0,inum),*this);
ev_all += ev;
if (evflag)
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagPairTersoffZBLComputeFullB<FULL,1> >(0,ignum),*this,ev);
else
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairTersoffZBLComputeFullB<FULL,0> >(0,ignum),*this);
ev_all += ev;
}
if (eflag_global) eng_vdwl += ev_all.evdwl;
if (vflag_global) {
virial[0] += ev_all.v[0];
virial[1] += ev_all.v[1];
virial[2] += ev_all.v[2];
virial[3] += ev_all.v[3];
virial[4] += ev_all.v[4];
virial[5] += ev_all.v[5];
}
if (eflag_atom) {
k_eatom.template modify<DeviceType>();
k_eatom.template sync<LMPHostType>();
}
if (vflag_atom) {
k_vatom.template modify<DeviceType>();
k_vatom.template sync<LMPHostType>();
}
if (vflag_fdotr) pair_virial_fdotr_compute(this);
copymode = 0;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeShortNeigh, const int& ii) const {
const int i = d_ilist[ii];
const X_FLOAT xtmp = x(i,0);
const X_FLOAT ytmp = x(i,1);
const X_FLOAT ztmp = x(i,2);
const int jnum = d_numneigh[i];
int inside = 0;
for (int jj = 0; jj < jnum; jj++) {
int j = d_neighbors(i,jj);
j &= NEIGHMASK;
const X_FLOAT delx = xtmp - x(j,0);
const X_FLOAT dely = ytmp - x(j,1);
const X_FLOAT delz = ztmp - x(j,2);
const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cutmax*cutmax) {
d_neighbors_short(i,inside) = j;
inside++;
}
}
d_numneigh_short(i) = inside;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeHalf<NEIGHFLAG,EVFLAG>, const int &ii, EV_FLOAT& ev) const {
// The f array is atomic for Half/Thread neighbor style
Kokkos::View<F_FLOAT*[3], typename DAT::t_f_array::array_layout,DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > a_f = f;
const int i = d_ilist[ii];
if (i >= nlocal) return;
const X_FLOAT xtmp = x(i,0);
const X_FLOAT ytmp = x(i,1);
const X_FLOAT ztmp = x(i,2);
const int itype = type(i);
const tagint itag = tag(i);
F_FLOAT fi[3], fj[3], fk[3];
//const AtomNeighborsConst d_neighbors_i = k_list.get_neighbors_const(i);
const int jnum = d_numneigh_short[i];
// repulsive
F_FLOAT f_x = 0.0;
F_FLOAT f_y = 0.0;
F_FLOAT f_z = 0.0;
for (int jj = 0; jj < jnum; jj++) {
int j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
const int jtype = type(j);
const tagint jtag = tag(j);
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x(j,2) < ztmp) continue;
if (x(j,2) == ztmp && x(j,1) < ytmp) continue;
if (x(j,2) == ztmp && x(j,1) == ytmp && x(j,0) < xtmp) continue;
}
const X_FLOAT delx = xtmp - x(j,0);
const X_FLOAT dely = ytmp - x(j,1);
const X_FLOAT delz = ztmp - x(j,2);
const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
const F_FLOAT cutsq = paramskk(itype,jtype,jtype).cutsq;
if (rsq > cutsq) continue;
// Tersoff repulsive portion
const F_FLOAT r = sqrt(rsq);
const F_FLOAT tmp_fce = ters_fc_k(itype,jtype,jtype,r);
const F_FLOAT tmp_fcd = ters_dfc(itype,jtype,jtype,r);
const F_FLOAT tmp_exp = exp(-paramskk(itype,jtype,jtype).lam1 * r);
const F_FLOAT frep_t = paramskk(itype,jtype,jtype).biga * tmp_exp *
(tmp_fcd - tmp_fce*paramskk(itype,jtype,jtype).lam1);
const F_FLOAT eng_t = tmp_fce * paramskk(itype,jtype,jtype).biga * tmp_exp;
// ZBL repulsive portion
const F_FLOAT esq = pow(global_e,2.0);
const F_FLOAT a_ij = (0.8854*global_a_0) /
(pow(paramskk(itype,jtype,jtype).Z_i,0.23) + pow(paramskk(itype,jtype,jtype).Z_j,0.23));
const F_FLOAT premult = (paramskk(itype,jtype,jtype).Z_i * paramskk(itype,jtype,jtype).Z_j * esq)/
(4.0*MY_PI*global_epsilon_0);
const F_FLOAT r_ov_a = r/a_ij;
const F_FLOAT phi = 0.1818*exp(-3.2*r_ov_a) + 0.5099*exp(-0.9423*r_ov_a) +
0.2802*exp(-0.4029*r_ov_a) + 0.02817*exp(-0.2016*r_ov_a);
const F_FLOAT dphi = (1.0/a_ij) * (-3.2*0.1818*exp(-3.2*r_ov_a) -
0.9423*0.5099*exp(-0.9423*r_ov_a) -
0.4029*0.2802*exp(-0.4029*r_ov_a) -
0.2016*0.02817*exp(-0.2016*r_ov_a));
const F_FLOAT frep_z = premult*-phi/rsq + premult*dphi/r;
const F_FLOAT eng_z = premult*(1.0/r)*phi;
// combine two parts with smoothing by Fermi-like function
F_FLOAT frep, eng;
frep = -(-fermi_d_k(itype,jtype,jtype,r) * eng_z +
(1.0 - fermi_k(itype,jtype,jtype,r))*frep_z +
fermi_d_k(itype,jtype,jtype,r)*eng_t + fermi_k(itype,jtype,jtype,r)*frep_t) / r;
if (eflag)
eng = (1.0 - fermi_k(itype,jtype,jtype,r)) * eng_z +
fermi_k(itype,jtype,jtype,r) * eng_t;
f_x += delx*frep;
f_y += dely*frep;
f_z += delz*frep;
a_f(j,0) -= delx*frep;
a_f(j,1) -= dely*frep;
a_f(j,2) -= delz*frep;
if (EVFLAG) {
if (eflag) ev.evdwl += eng;
if (vflag_either || eflag_atom) this->template ev_tally<NEIGHFLAG>(ev,i,j,eng,frep,delx,dely,delz);
}
}
// attractive: bond order
for (int jj = 0; jj < jnum; jj++) {
int j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
const int jtype = type(j);
const F_FLOAT delx1 = xtmp - x(j,0);
const F_FLOAT dely1 = ytmp - x(j,1);
const F_FLOAT delz1 = ztmp - x(j,2);
const F_FLOAT rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
const F_FLOAT cutsq1 = paramskk(itype,jtype,jtype).cutsq;
F_FLOAT bo_ij = 0.0;
if (rsq1 > cutsq1) continue;
const F_FLOAT rij = sqrt(rsq1);
for (int kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
int k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
const int ktype = type(k);
const F_FLOAT delx2 = xtmp - x(k,0);
const F_FLOAT dely2 = ytmp - x(k,1);
const F_FLOAT delz2 = ztmp - x(k,2);
const F_FLOAT rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
const F_FLOAT cutsq2 = paramskk(itype,jtype,ktype).cutsq;
if (rsq2 > cutsq2) continue;
const F_FLOAT rik = sqrt(rsq2);
bo_ij += bondorder(itype,jtype,ktype,rij,delx1,dely1,delz1,rik,delx2,dely2,delz2);
}
// attractive: pairwise potential and force
const F_FLOAT fa = ters_fa_k(itype,jtype,jtype,rij);
const F_FLOAT dfa = ters_dfa(itype,jtype,jtype,rij);
const F_FLOAT bij = ters_bij_k(itype,jtype,jtype,bo_ij);
const F_FLOAT fatt = -0.5*bij * dfa / rij;
const F_FLOAT prefactor = 0.5*fa * ters_dbij(itype,jtype,jtype,bo_ij);
f_x += delx1*fatt;
f_y += dely1*fatt;
f_z += delz1*fatt;
F_FLOAT fj_x = -delx1*fatt;
F_FLOAT fj_y = -dely1*fatt;
F_FLOAT fj_z = -delz1*fatt;
if (EVFLAG) {
const F_FLOAT eng = 0.5*bij * fa;
if (eflag) ev.evdwl += eng;
if (vflag_either || eflag_atom)
this->template ev_tally<NEIGHFLAG>(ev,i,j,eng,fatt,delx1,dely1,delz1);
}
// attractive: three-body force
for (int kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
int k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
const int ktype = type(k);
const F_FLOAT delx2 = xtmp - x(k,0);
const F_FLOAT dely2 = ytmp - x(k,1);
const F_FLOAT delz2 = ztmp - x(k,2);
const F_FLOAT rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
const F_FLOAT cutsq2 = paramskk(itype,jtype,ktype).cutsq;
if (rsq2 > cutsq2) continue;
const F_FLOAT rik = sqrt(rsq2);
ters_dthb(itype,jtype,ktype,prefactor,rij,delx1,dely1,delz1,
rik,delx2,dely2,delz2,fi,fj,fk);
f_x += fi[0];
f_y += fi[1];
f_z += fi[2];
fj_x += fj[0];
fj_y += fj[1];
fj_z += fj[2];
a_f(k,0) += fk[0];
a_f(k,1) += fk[1];
a_f(k,2) += fk[2];
if (vflag_atom) {
F_FLOAT delrij[3], delrik[3];
delrij[0] = -delx1; delrij[1] = -dely1; delrij[2] = -delz1;
delrik[0] = -delx2; delrik[1] = -dely2; delrik[2] = -delz2;
if (vflag_either) this->template v_tally3<NEIGHFLAG>(ev,i,j,k,fj,fk,delrij,delrik);
}
}
a_f(j,0) += fj_x;
a_f(j,1) += fj_y;
a_f(j,2) += fj_z;
}
a_f(i,0) += f_x;
a_f(i,1) += f_y;
a_f(i,2) += f_z;
}
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeHalf<NEIGHFLAG,EVFLAG>, const int &ii) const {
EV_FLOAT ev;
this->template operator()<NEIGHFLAG,EVFLAG>(TagPairTersoffZBLComputeHalf<NEIGHFLAG,EVFLAG>(), ii, ev);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullA<NEIGHFLAG,EVFLAG>, const int &ii, EV_FLOAT& ev) const {
const int i = d_ilist[ii];
const X_FLOAT xtmp = x(i,0);
const X_FLOAT ytmp = x(i,1);
const X_FLOAT ztmp = x(i,2);
const int itype = type(i);
int j,k,jj,kk,jtype,ktype;
F_FLOAT rsq1, cutsq1, rsq2, cutsq2, rij, rik, bo_ij;
F_FLOAT fi[3], fj[3], fk[3];
X_FLOAT delx1, dely1, delz1, delx2, dely2, delz2;
//const AtomNeighborsConst d_neighbors_i = k_list.get_neighbors_const(i);
const int jnum = d_numneigh[i];
// repulsive
F_FLOAT f_x = 0.0;
F_FLOAT f_y = 0.0;
F_FLOAT f_z = 0.0;
for (jj = 0; jj < jnum; jj++) {
j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
const int jtype = type(j);
const X_FLOAT delx = xtmp - x(j,0);
const X_FLOAT dely = ytmp - x(j,1);
const X_FLOAT delz = ztmp - x(j,2);
const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
const F_FLOAT cutsq = paramskk(itype,jtype,jtype).cutsq;
if (rsq > cutsq) continue;
// Tersoff repulsive portion
const F_FLOAT r = sqrt(rsq);
const F_FLOAT tmp_fce = ters_fc_k(itype,jtype,jtype,r);
const F_FLOAT tmp_fcd = ters_dfc(itype,jtype,jtype,r);
const F_FLOAT tmp_exp = exp(-paramskk(itype,jtype,jtype).lam1 * r);
const F_FLOAT frep_t = paramskk(itype,jtype,jtype).biga * tmp_exp *
(tmp_fcd - tmp_fce*paramskk(itype,jtype,jtype).lam1);
const F_FLOAT eng_t = tmp_fce * paramskk(itype,jtype,jtype).biga * tmp_exp;
// ZBL repulsive portion
const F_FLOAT esq = pow(global_e,2.0);
const F_FLOAT a_ij = (0.8854*global_a_0) /
(pow(paramskk(itype,jtype,jtype).Z_i,0.23) + pow(paramskk(itype,jtype,jtype).Z_j,0.23));
const F_FLOAT premult = (paramskk(itype,jtype,jtype).Z_i * paramskk(itype,jtype,jtype).Z_j * esq)/
(4.0*MY_PI*global_epsilon_0);
const F_FLOAT r_ov_a = r/a_ij;
const F_FLOAT phi = 0.1818*exp(-3.2*r_ov_a) + 0.5099*exp(-0.9423*r_ov_a) +
0.2802*exp(-0.4029*r_ov_a) + 0.02817*exp(-0.2016*r_ov_a);
const F_FLOAT dphi = (1.0/a_ij) * (-3.2*0.1818*exp(-3.2*r_ov_a) -
0.9423*0.5099*exp(-0.9423*r_ov_a) -
0.4029*0.2802*exp(-0.4029*r_ov_a) -
0.2016*0.02817*exp(-0.2016*r_ov_a));
const F_FLOAT frep_z = premult*-phi/rsq + premult*dphi/r;
const F_FLOAT eng_z = premult*(1.0/r)*phi;
// combine two parts with smoothing by Fermi-like function
F_FLOAT frep, eng;
frep = -(-fermi_d_k(itype,jtype,jtype,r) * eng_z +
(1.0 - fermi_k(itype,jtype,jtype,r))*frep_z +
fermi_d_k(itype,jtype,jtype,r)*eng_t + fermi_k(itype,jtype,jtype,r)*frep_t) / r;
if (eflag)
eng = (1.0 - fermi_k(itype,jtype,jtype,r)) * eng_z +
fermi_k(itype,jtype,jtype,r) * eng_t;
f_x += delx*frep;
f_y += dely*frep;
f_z += delz*frep;
if (EVFLAG) {
if (eflag)
ev.evdwl += 0.5*eng;
if (vflag_either || eflag_atom)
this->template ev_tally<NEIGHFLAG>(ev,i,j,eng,frep,delx,dely,delz);
}
}
// attractive: bond order
for (jj = 0; jj < jnum; jj++) {
j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
jtype = type(j);
delx1 = xtmp - x(j,0);
dely1 = ytmp - x(j,1);
delz1 = ztmp - x(j,2);
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
cutsq1 = paramskk(itype,jtype,jtype).cutsq;
bo_ij = 0.0;
if (rsq1 > cutsq1) continue;
rij = sqrt(rsq1);
for (kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
ktype = type(k);
delx2 = xtmp - x(k,0);
dely2 = ytmp - x(k,1);
delz2 = ztmp - x(k,2);
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
cutsq2 = paramskk(itype,jtype,ktype).cutsq;
if (rsq2 > cutsq2) continue;
rik = sqrt(rsq2);
bo_ij += bondorder(itype,jtype,ktype,rij,delx1,dely1,delz1,rik,delx2,dely2,delz2);
}
// attractive: pairwise potential and force
const F_FLOAT fa = ters_fa_k(itype,jtype,jtype,rij);
const F_FLOAT dfa = ters_dfa(itype,jtype,jtype,rij);
const F_FLOAT bij = ters_bij_k(itype,jtype,jtype,bo_ij);
const F_FLOAT fatt = -0.5*bij * dfa / rij;
const F_FLOAT prefactor = 0.5*fa * ters_dbij(itype,jtype,jtype,bo_ij);
const F_FLOAT eng = 0.5*bij * fa;
f_x += delx1*fatt;
f_y += dely1*fatt;
f_z += delz1*fatt;
if (EVFLAG) {
if (eflag) ev.evdwl += 0.5*eng;
if (vflag_either || eflag_atom)
this->template ev_tally<NEIGHFLAG>(ev,i,j,eng,fatt,delx1,dely1,delz1);
}
// attractive: three-body force
for (kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
ktype = type(k);
delx2 = xtmp - x(k,0);
dely2 = ytmp - x(k,1);
delz2 = ztmp - x(k,2);
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
cutsq2 = paramskk(itype,jtype,ktype).cutsq;
if (rsq2 > cutsq2) continue;
rik = sqrt(rsq2);
ters_dthb(itype,jtype,ktype,prefactor,rij,delx1,dely1,delz1,
rik,delx2,dely2,delz2,fi,fj,fk);
f_x += fi[0];
f_y += fi[1];
f_z += fi[2];
if (vflag_atom) {
F_FLOAT delrij[3], delrik[3];
delrij[0] = -delx1; delrij[1] = -dely1; delrij[2] = -delz1;
delrik[0] = -delx2; delrik[1] = -dely2; delrik[2] = -delz2;
if (vflag_either) this->template v_tally3<NEIGHFLAG>(ev,i,j,k,fj,fk,delrij,delrik);
}
}
}
f(i,0) += f_x;
f(i,1) += f_y;
f(i,2) += f_z;
}
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullA<NEIGHFLAG,EVFLAG>, const int &ii) const {
EV_FLOAT ev;
this->template operator()<NEIGHFLAG,EVFLAG>(TagPairTersoffZBLComputeFullA<NEIGHFLAG,EVFLAG>(), ii, ev);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullB<NEIGHFLAG,EVFLAG>, const int &ii, EV_FLOAT& ev) const {
const int i = d_ilist[ii];
const X_FLOAT xtmp = x(i,0);
const X_FLOAT ytmp = x(i,1);
const X_FLOAT ztmp = x(i,2);
const int itype = type(i);
int j,k,jj,kk,jtype,ktype,j_jnum;
F_FLOAT rsq1, cutsq1, rsq2, cutsq2, rij, rik, bo_ij;
F_FLOAT fj[3], fk[3];
X_FLOAT delx1, dely1, delz1, delx2, dely2, delz2;
const int jnum = d_numneigh_short[i];
F_FLOAT f_x = 0.0;
F_FLOAT f_y = 0.0;
F_FLOAT f_z = 0.0;
// attractive: bond order
for (jj = 0; jj < jnum; jj++) {
j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
if (j >= nlocal) continue;
jtype = type(j);
delx1 = x(j,0) - xtmp;
dely1 = x(j,1) - ytmp;
delz1 = x(j,2) - ztmp;
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
cutsq1 = paramskk(jtype,itype,itype).cutsq;
bo_ij = 0.0;
if (rsq1 > cutsq1) continue;
rij = sqrt(rsq1);
j_jnum = d_numneigh_short[j];
for (kk = 0; kk < j_jnum; kk++) {
k = d_neighbors_short(j,kk);
if (k == i) continue;
k &= NEIGHMASK;
ktype = type(k);
delx2 = x(j,0) - x(k,0);
dely2 = x(j,1) - x(k,1);
delz2 = x(j,2) - x(k,2);
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
cutsq2 = paramskk(jtype,itype,ktype).cutsq;
if (rsq2 > cutsq2) continue;
rik = sqrt(rsq2);
bo_ij += bondorder(jtype,itype,ktype,rij,delx1,dely1,delz1,rik,delx2,dely2,delz2);
}
// attractive: pairwise potential and force
const F_FLOAT fa = ters_fa_k(jtype,itype,itype,rij);
const F_FLOAT dfa = ters_dfa(jtype,itype,itype,rij);
const F_FLOAT bij = ters_bij_k(jtype,itype,itype,bo_ij);
const F_FLOAT fatt = -0.5*bij * dfa / rij;
const F_FLOAT prefactor = 0.5*fa * ters_dbij(jtype,itype,itype,bo_ij);
const F_FLOAT eng = 0.5*bij * fa;
f_x -= delx1*fatt;
f_y -= dely1*fatt;
f_z -= delz1*fatt;
if (EVFLAG) {
if (eflag)
ev.evdwl += 0.5 * eng;
if (vflag_either || eflag_atom)
this->template ev_tally<NEIGHFLAG>(ev,i,j,eng,fatt,delx1,dely1,delz1);
}
// attractive: three-body force
for (kk = 0; kk < j_jnum; kk++) {
k = d_neighbors_short(j,kk);
if (k == i) continue;
k &= NEIGHMASK;
ktype = type(k);
delx2 = x(j,0) - x(k,0);
dely2 = x(j,1) - x(k,1);
delz2 = x(j,2) - x(k,2);
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
cutsq2 = paramskk(jtype,itype,ktype).cutsq;
if (rsq2 > cutsq2) continue;
rik = sqrt(rsq2);
ters_dthbj(jtype,itype,ktype,prefactor,rij,delx1,dely1,delz1,
rik,delx2,dely2,delz2,fj,fk);
f_x += fj[0];
f_y += fj[1];
f_z += fj[2];
if (vflag_atom) {
F_FLOAT delrji[3], delrjk[3];
delrji[0] = -delx1; delrji[1] = -dely1; delrji[2] = -delz1;
delrjk[0] = -delx2; delrjk[1] = -dely2; delrjk[2] = -delz2;
if (vflag_either) v_tally3_atom(ev,i,j,k,fj,fk,delrji,delrjk);
}
const F_FLOAT fa_jk = ters_fa_k(jtype,ktype,itype,rik);
const F_FLOAT prefactor_jk = 0.5*fa_jk * ters_dbij(jtype,ktype,itype,bo_ij);
ters_dthbk(jtype,ktype,itype,prefactor_jk,rik,delx2,dely2,delz2,
rij,delx1,dely1,delz1,fk);
f_x += fk[0];
f_y += fk[1];
f_z += fk[2];
}
}
f(i,0) += f_x;
f(i,1) += f_y;
f(i,2) += f_z;
}
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullB<NEIGHFLAG,EVFLAG>, const int &ii) const {
EV_FLOAT ev;
this->template operator()<NEIGHFLAG,EVFLAG>(TagPairTersoffZBLComputeFullB<NEIGHFLAG,EVFLAG>(), ii, ev);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
double PairTersoffZBLKokkos<DeviceType>::ters_fc_k(const int &i, const int &j,
const int &k, const F_FLOAT &r) const
{
const F_FLOAT ters_R = paramskk(i,j,k).bigr;
const F_FLOAT ters_D = paramskk(i,j,k).bigd;
if (r < ters_R-ters_D) return 1.0;
if (r > ters_R+ters_D) return 0.0;
return 0.5*(1.0 - sin(MY_PI2*(r - ters_R)/ters_D));
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
double PairTersoffZBLKokkos<DeviceType>::ters_dfc(const int &i, const int &j,
const int &k, const F_FLOAT &r) const
{
const F_FLOAT ters_R = paramskk(i,j,k).bigr;
const F_FLOAT ters_D = paramskk(i,j,k).bigd;
if (r < ters_R-ters_D) return 0.0;
if (r > ters_R+ters_D) return 0.0;
return -(MY_PI4/ters_D) * cos(MY_PI2*(r - ters_R)/ters_D);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
double PairTersoffZBLKokkos<DeviceType>::bondorder(const int &i, const int &j, const int &k,
const F_FLOAT &rij, const F_FLOAT &dx1, const F_FLOAT &dy1, const F_FLOAT &dz1,
const F_FLOAT &rik, const F_FLOAT &dx2, const F_FLOAT &dy2, const F_FLOAT &dz2) const
{
F_FLOAT arg, ex_delr;
const F_FLOAT costheta = (dx1*dx2 + dy1*dy2 + dz1*dz2)/(rij*rik);
if (int(paramskk(i,j,k).powerm) == 3) arg = pow(paramskk(i,j,k).lam3 * (rij-rik),3.0);
else arg = paramskk(i,j,k).lam3 * (rij-rik);
if (arg > 69.0776) ex_delr = 1.e30;
else if (arg < -69.0776) ex_delr = 0.0;
else ex_delr = exp(arg);
return ters_fc_k(i,j,k,rik) * ters_gijk(i,j,k,costheta) * ex_delr;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
double PairTersoffZBLKokkos<DeviceType>::
ters_gijk(const int &i, const int &j, const int &k, const F_FLOAT &cos) const
{
const F_FLOAT ters_c = paramskk(i,j,k).c * paramskk(i,j,k).c;
const F_FLOAT ters_d = paramskk(i,j,k).d * paramskk(i,j,k).d;
const F_FLOAT hcth = paramskk(i,j,k).h - cos;
return paramskk(i,j,k).gamma*(1.0 + ters_c/ters_d - ters_c/(ters_d+hcth*hcth));
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
double PairTersoffZBLKokkos<DeviceType>::
ters_dgijk(const int &i, const int &j, const int &k, const F_FLOAT &cos) const
{
const F_FLOAT ters_c = paramskk(i,j,k).c * paramskk(i,j,k).c;
const F_FLOAT ters_d = paramskk(i,j,k).d * paramskk(i,j,k).d;
const F_FLOAT hcth = paramskk(i,j,k).h - cos;
const F_FLOAT numerator = -2.0 * ters_c * hcth;
const F_FLOAT denominator = 1.0/(ters_d + hcth*hcth);
return paramskk(i,j,k).gamma * numerator * denominator * denominator;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
double PairTersoffZBLKokkos<DeviceType>::ters_fa_k(const int &i, const int &j,
const int &k, const F_FLOAT &r) const
{
if (r > paramskk(i,j,k).bigr + paramskk(i,j,k).bigd) return 0.0;
return -paramskk(i,j,k).bigb * exp(-paramskk(i,j,k).lam2 * r)
* ters_fc_k(i,j,k,r) * fermi_k(i,j,k,r);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
double PairTersoffZBLKokkos<DeviceType>::ters_dfa(const int &i, const int &j,
const int &k, const F_FLOAT &r) const
{
if (r > paramskk(i,j,k).bigr + paramskk(i,j,k).bigd) return 0.0;
return paramskk(i,j,k).bigb * exp(-paramskk(i,j,k).lam2 * r) *
(paramskk(i,j,k).lam2 * ters_fc_k(i,j,k,r) * fermi_k(i,j,k,r) -
ters_dfc(i,j,k,r) * fermi_k(i,j,k,r) - ters_fc_k(i,j,k,r) *
fermi_d_k(i,j,k,r));
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
double PairTersoffZBLKokkos<DeviceType>::ters_bij_k(const int &i, const int &j,
const int &k, const F_FLOAT &bo) const
{
const F_FLOAT tmp = paramskk(i,j,k).beta * bo;
if (tmp > paramskk(i,j,k).c1) return 1.0/sqrt(tmp);
if (tmp > paramskk(i,j,k).c2)
return (1.0 - pow(tmp,-paramskk(i,j,k).powern) / (2.0*paramskk(i,j,k).powern))/sqrt(tmp);
if (tmp < paramskk(i,j,k).c4) return 1.0;
if (tmp < paramskk(i,j,k).c3)
return 1.0 - pow(tmp,paramskk(i,j,k).powern)/(2.0*paramskk(i,j,k).powern);
return pow(1.0 + pow(tmp,paramskk(i,j,k).powern), -1.0/(2.0*paramskk(i,j,k).powern));
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
double PairTersoffZBLKokkos<DeviceType>::ters_dbij(const int &i, const int &j,
const int &k, const F_FLOAT &bo) const
{
const F_FLOAT tmp = paramskk(i,j,k).beta * bo;
if (tmp > paramskk(i,j,k).c1) return paramskk(i,j,k).beta * -0.5*pow(tmp,-1.5);
if (tmp > paramskk(i,j,k).c2)
return paramskk(i,j,k).beta * (-0.5*pow(tmp,-1.5) *
(1.0 - 0.5*(1.0 + 1.0/(2.0*paramskk(i,j,k).powern)) *
pow(tmp,-paramskk(i,j,k).powern)));
if (tmp < paramskk(i,j,k).c4) return 0.0;
if (tmp < paramskk(i,j,k).c3)
return -0.5*paramskk(i,j,k).beta * pow(tmp,paramskk(i,j,k).powern-1.0);
const F_FLOAT tmp_n = pow(tmp,paramskk(i,j,k).powern);
return -0.5 * pow(1.0+tmp_n, -1.0-(1.0/(2.0*paramskk(i,j,k).powern)))*tmp_n / bo;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairTersoffZBLKokkos<DeviceType>::ters_dthb(
const int &i, const int &j, const int &k, const F_FLOAT &prefactor,
const F_FLOAT &rij, const F_FLOAT &dx1, const F_FLOAT &dy1, const F_FLOAT &dz1,
const F_FLOAT &rik, const F_FLOAT &dx2, const F_FLOAT &dy2, const F_FLOAT &dz2,
F_FLOAT *fi, F_FLOAT *fj, F_FLOAT *fk) const
{
// from PairTersoffZBL::attractive
F_FLOAT rij_hat[3],rik_hat[3];
F_FLOAT rijinv,rikinv;
F_FLOAT delrij[3], delrik[3];
delrij[0] = dx1; delrij[1] = dy1; delrij[2] = dz1;
delrik[0] = dx2; delrik[1] = dy2; delrik[2] = dz2;
//rij = sqrt(rsq1);
rijinv = 1.0/rij;
vec3_scale(rijinv,delrij,rij_hat);
//rik = sqrt(rsq2);
rikinv = 1.0/rik;
vec3_scale(rikinv,delrik,rik_hat);
// from PairTersoffZBL::ters_zetaterm_d
F_FLOAT gijk,dgijk,ex_delr,dex_delr,fc,dfc,cos,tmp;
F_FLOAT dcosfi[3],dcosfj[3],dcosfk[3];
fc = ters_fc_k(i,j,k,rik);
dfc = ters_dfc(i,j,k,rik);
if (int(paramskk(i,j,k).powerm) == 3) tmp = pow(paramskk(i,j,k).lam3 * (rij-rik),3.0);
else tmp = paramskk(i,j,k).lam3 * (rij-rik);
if (tmp > 69.0776) ex_delr = 1.e30;
else if (tmp < -69.0776) ex_delr = 0.0;
else ex_delr = exp(tmp);
if (int(paramskk(i,j,k).powerm) == 3)
dex_delr = 3.0*pow(paramskk(i,j,k).lam3,3.0) * pow(rij-rik,2.0)*ex_delr;
else dex_delr = paramskk(i,j,k).lam3 * ex_delr;
cos = vec3_dot(rij_hat,rik_hat);
gijk = ters_gijk(i,j,k,cos);
dgijk = ters_dgijk(i,j,k,cos);
// from PairTersoffZBL::costheta_d
vec3_scaleadd(-cos,rij_hat,rik_hat,dcosfj);
vec3_scale(rijinv,dcosfj,dcosfj);
vec3_scaleadd(-cos,rik_hat,rij_hat,dcosfk);
vec3_scale(rikinv,dcosfk,dcosfk);
vec3_add(dcosfj,dcosfk,dcosfi);
vec3_scale(-1.0,dcosfi,dcosfi);
vec3_scale(-dfc*gijk*ex_delr,rik_hat,fi);
vec3_scaleadd(fc*dgijk*ex_delr,dcosfi,fi,fi);
vec3_scaleadd(fc*gijk*dex_delr,rik_hat,fi,fi);
vec3_scaleadd(-fc*gijk*dex_delr,rij_hat,fi,fi);
vec3_scale(prefactor,fi,fi);
vec3_scale(fc*dgijk*ex_delr,dcosfj,fj);
vec3_scaleadd(fc*gijk*dex_delr,rij_hat,fj,fj);
vec3_scale(prefactor,fj,fj);
vec3_scale(dfc*gijk*ex_delr,rik_hat,fk);
vec3_scaleadd(fc*dgijk*ex_delr,dcosfk,fk,fk);
vec3_scaleadd(-fc*gijk*dex_delr,rik_hat,fk,fk);
vec3_scale(prefactor,fk,fk);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairTersoffZBLKokkos<DeviceType>::ters_dthbj(
const int &i, const int &j, const int &k, const F_FLOAT &prefactor,
const F_FLOAT &rij, const F_FLOAT &dx1, const F_FLOAT &dy1, const F_FLOAT &dz1,
const F_FLOAT &rik, const F_FLOAT &dx2, const F_FLOAT &dy2, const F_FLOAT &dz2,
F_FLOAT *fj, F_FLOAT *fk) const
{
F_FLOAT rij_hat[3],rik_hat[3];
F_FLOAT rijinv,rikinv;
F_FLOAT delrij[3], delrik[3];
delrij[0] = dx1; delrij[1] = dy1; delrij[2] = dz1;
delrik[0] = dx2; delrik[1] = dy2; delrik[2] = dz2;
rijinv = 1.0/rij;
vec3_scale(rijinv,delrij,rij_hat);
rikinv = 1.0/rik;
vec3_scale(rikinv,delrik,rik_hat);
F_FLOAT gijk,dgijk,ex_delr,dex_delr,fc,dfc,cos,tmp;
F_FLOAT dcosfi[3],dcosfj[3],dcosfk[3];
fc = ters_fc_k(i,j,k,rik);
dfc = ters_dfc(i,j,k,rik);
if (int(paramskk(i,j,k).powerm) == 3) tmp = pow(paramskk(i,j,k).lam3 * (rij-rik),3.0);
else tmp = paramskk(i,j,k).lam3 * (rij-rik);
if (tmp > 69.0776) ex_delr = 1.e30;
else if (tmp < -69.0776) ex_delr = 0.0;
else ex_delr = exp(tmp);
if (int(paramskk(i,j,k).powerm) == 3)
dex_delr = 3.0*pow(paramskk(i,j,k).lam3,3.0) * pow(rij-rik,2.0)*ex_delr;
else dex_delr = paramskk(i,j,k).lam3 * ex_delr;
cos = vec3_dot(rij_hat,rik_hat);
gijk = ters_gijk(i,j,k,cos);
dgijk = ters_dgijk(i,j,k,cos);
vec3_scaleadd(-cos,rij_hat,rik_hat,dcosfj);
vec3_scale(rijinv,dcosfj,dcosfj);
vec3_scaleadd(-cos,rik_hat,rij_hat,dcosfk);
vec3_scale(rikinv,dcosfk,dcosfk);
vec3_add(dcosfj,dcosfk,dcosfi);
vec3_scale(-1.0,dcosfi,dcosfi);
vec3_scale(fc*dgijk*ex_delr,dcosfj,fj);
vec3_scaleadd(fc*gijk*dex_delr,rij_hat,fj,fj);
vec3_scale(prefactor,fj,fj);
vec3_scale(dfc*gijk*ex_delr,rik_hat,fk);
vec3_scaleadd(fc*dgijk*ex_delr,dcosfk,fk,fk);
vec3_scaleadd(-fc*gijk*dex_delr,rik_hat,fk,fk);
vec3_scale(prefactor,fk,fk);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairTersoffZBLKokkos<DeviceType>::ters_dthbk(
const int &i, const int &j, const int &k, const F_FLOAT &prefactor,
const F_FLOAT &rij, const F_FLOAT &dx1, const F_FLOAT &dy1, const F_FLOAT &dz1,
const F_FLOAT &rik, const F_FLOAT &dx2, const F_FLOAT &dy2, const F_FLOAT &dz2,
F_FLOAT *fk) const
{
F_FLOAT rij_hat[3],rik_hat[3];
F_FLOAT rijinv,rikinv;
F_FLOAT delrij[3], delrik[3];
delrij[0] = dx1; delrij[1] = dy1; delrij[2] = dz1;
delrik[0] = dx2; delrik[1] = dy2; delrik[2] = dz2;
rijinv = 1.0/rij;
vec3_scale(rijinv,delrij,rij_hat);
rikinv = 1.0/rik;
vec3_scale(rikinv,delrik,rik_hat);
F_FLOAT gijk,dgijk,ex_delr,dex_delr,fc,dfc,cos,tmp;
F_FLOAT dcosfi[3],dcosfj[3],dcosfk[3];
fc = ters_fc_k(i,j,k,rik);
dfc = ters_dfc(i,j,k,rik);
if (int(paramskk(i,j,k).powerm) == 3) tmp = pow(paramskk(i,j,k).lam3 * (rij-rik),3.0);
else tmp = paramskk(i,j,k).lam3 * (rij-rik);
if (tmp > 69.0776) ex_delr = 1.e30;
else if (tmp < -69.0776) ex_delr = 0.0;
else ex_delr = exp(tmp);
if (int(paramskk(i,j,k).powerm) == 3)
dex_delr = 3.0*pow(paramskk(i,j,k).lam3,3.0) * pow(rij-rik,2.0)*ex_delr;
else dex_delr = paramskk(i,j,k).lam3 * ex_delr;
cos = vec3_dot(rij_hat,rik_hat);
gijk = ters_gijk(i,j,k,cos);
dgijk = ters_dgijk(i,j,k,cos);
vec3_scaleadd(-cos,rij_hat,rik_hat,dcosfj);
vec3_scale(rijinv,dcosfj,dcosfj);
vec3_scaleadd(-cos,rik_hat,rij_hat,dcosfk);
vec3_scale(rikinv,dcosfk,dcosfk);
vec3_add(dcosfj,dcosfk,dcosfi);
vec3_scale(-1.0,dcosfi,dcosfi);
vec3_scale(dfc*gijk*ex_delr,rik_hat,fk);
vec3_scaleadd(fc*dgijk*ex_delr,dcosfk,fk,fk);
vec3_scaleadd(-fc*gijk*dex_delr,rik_hat,fk,fk);
vec3_scale(prefactor,fk,fk);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
double PairTersoffZBLKokkos<DeviceType>::fermi_k(const int &i, const int &j,
const int &k, const F_FLOAT &r) const
{
return 1.0 / (1.0 + exp(-paramskk(i,j,k).ZBLexpscale *
(r - paramskk(i,j,k).ZBLcut)));
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
double PairTersoffZBLKokkos<DeviceType>::fermi_d_k(const int &i, const int &j,
const int &k, const F_FLOAT &r) const
{
return paramskk(i,j,k).ZBLexpscale * exp(-paramskk(i,j,k).ZBLexpscale *
(r - paramskk(i,j,k).ZBLcut)) /
pow(1.0 + exp(-paramskk(i,j,k).ZBLexpscale *
(r - paramskk(i,j,k).ZBLcut)),2.0);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void PairTersoffZBLKokkos<DeviceType>::ev_tally(EV_FLOAT &ev, const int &i, const int &j,
const F_FLOAT &epair, const F_FLOAT &fpair, const F_FLOAT &delx,
const F_FLOAT &dely, const F_FLOAT &delz) const
{
const int VFLAG = vflag_either;
// The eatom and vatom arrays are atomic for Half/Thread neighbor style
Kokkos::View<E_FLOAT*, typename DAT::t_efloat_1d::array_layout,DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > v_eatom = k_eatom.view<DeviceType>();
Kokkos::View<F_FLOAT*[6], typename DAT::t_virial_array::array_layout,DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > v_vatom = k_vatom.view<DeviceType>();
if (eflag_atom) {
const E_FLOAT epairhalf = 0.5 * epair;
v_eatom[i] += epairhalf;
if (NEIGHFLAG != FULL) v_eatom[j] += epairhalf;
}
if (VFLAG) {
const E_FLOAT v0 = delx*delx*fpair;
const E_FLOAT v1 = dely*dely*fpair;
const E_FLOAT v2 = delz*delz*fpair;
const E_FLOAT v3 = delx*dely*fpair;
const E_FLOAT v4 = delx*delz*fpair;
const E_FLOAT v5 = dely*delz*fpair;
if (vflag_global) {
if (NEIGHFLAG != FULL) {
ev.v[0] += v0;
ev.v[1] += v1;
ev.v[2] += v2;
ev.v[3] += v3;
ev.v[4] += v4;
ev.v[5] += v5;
} else {
ev.v[0] += 0.5*v0;
ev.v[1] += 0.5*v1;
ev.v[2] += 0.5*v2;
ev.v[3] += 0.5*v3;
ev.v[4] += 0.5*v4;
ev.v[5] += 0.5*v5;
}
}
if (vflag_atom) {
v_vatom(i,0) += 0.5*v0;
v_vatom(i,1) += 0.5*v1;
v_vatom(i,2) += 0.5*v2;
v_vatom(i,3) += 0.5*v3;
v_vatom(i,4) += 0.5*v4;
v_vatom(i,5) += 0.5*v5;
if (NEIGHFLAG != FULL) {
v_vatom(j,0) += 0.5*v0;
v_vatom(j,1) += 0.5*v1;
v_vatom(j,2) += 0.5*v2;
v_vatom(j,3) += 0.5*v3;
v_vatom(j,4) += 0.5*v4;
v_vatom(j,5) += 0.5*v5;
}
}
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void PairTersoffZBLKokkos<DeviceType>::v_tally3(EV_FLOAT &ev, const int &i, const int &j, const int &k,
F_FLOAT *fj, F_FLOAT *fk, F_FLOAT *drij, F_FLOAT *drik) const
{
// The eatom and vatom arrays are atomic for Half/Thread neighbor style
Kokkos::View<F_FLOAT*[6], typename DAT::t_virial_array::array_layout,DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > v_vatom = k_vatom.view<DeviceType>();
F_FLOAT v[6];
v[0] = THIRD * (drij[0]*fj[0] + drik[0]*fk[0]);
v[1] = THIRD * (drij[1]*fj[1] + drik[1]*fk[1]);
v[2] = THIRD * (drij[2]*fj[2] + drik[2]*fk[2]);
v[3] = THIRD * (drij[0]*fj[1] + drik[0]*fk[1]);
v[4] = THIRD * (drij[0]*fj[2] + drik[0]*fk[2]);
v[5] = THIRD * (drij[1]*fj[2] + drik[1]*fk[2]);
if (vflag_global) {
ev.v[0] += v[0];
ev.v[1] += v[1];
ev.v[2] += v[2];
ev.v[3] += v[3];
ev.v[4] += v[4];
ev.v[5] += v[5];
}
if (vflag_atom) {
v_vatom(i,0) += v[0]; v_vatom(i,1) += v[1]; v_vatom(i,2) += v[2];
v_vatom(i,3) += v[3]; v_vatom(i,4) += v[4]; v_vatom(i,5) += v[5];
if (NEIGHFLAG != FULL) {
v_vatom(j,0) += v[0]; v_vatom(j,1) += v[1]; v_vatom(j,2) += v[2];
v_vatom(j,3) += v[3]; v_vatom(j,4) += v[4]; v_vatom(j,5) += v[5];
v_vatom(k,0) += v[0]; v_vatom(k,1) += v[1]; v_vatom(k,2) += v[2];
v_vatom(k,3) += v[3]; v_vatom(k,4) += v[4]; v_vatom(k,5) += v[5];
}
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairTersoffZBLKokkos<DeviceType>::v_tally3_atom(EV_FLOAT &ev, const int &i, const int &j, const int &k,
F_FLOAT *fj, F_FLOAT *fk, F_FLOAT *drji, F_FLOAT *drjk) const
{
F_FLOAT v[6];
v[0] = THIRD * (drji[0]*fj[0] + drjk[0]*fk[0]);
v[1] = THIRD * (drji[1]*fj[1] + drjk[1]*fk[1]);
v[2] = THIRD * (drji[2]*fj[2] + drjk[2]*fk[2]);
v[3] = THIRD * (drji[0]*fj[1] + drjk[0]*fk[1]);
v[4] = THIRD * (drji[0]*fj[2] + drjk[0]*fk[2]);
v[5] = THIRD * (drji[1]*fj[2] + drjk[1]*fk[2]);
if (vflag_global) {
ev.v[0] += v[0];
ev.v[1] += v[1];
ev.v[2] += v[2];
ev.v[3] += v[3];
ev.v[4] += v[4];
ev.v[5] += v[5];
}
if (vflag_atom) {
d_vatom(i,0) += v[0]; d_vatom(i,1) += v[1]; d_vatom(i,2) += v[2];
d_vatom(i,3) += v[3]; d_vatom(i,4) += v[4]; d_vatom(i,5) += v[5];
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
int PairTersoffZBLKokkos<DeviceType>::sbmask(const int& j) const {
return j >> SBBITS & 3;
}
namespace LAMMPS_NS {
template class PairTersoffZBLKokkos<LMPDeviceType>;
#ifdef KOKKOS_HAVE_CUDA
template class PairTersoffZBLKokkos<LMPHostType>;
#endif
}
c4science · Help