tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #4318 from weinbe2/bugfix/reaxff-bond-int-overflow

Browse Source 
Fixes to ReaxFF/Kokkos integer overflow issues for bond tables + cleanup of unused compute
This commit is contained in:
Axel Kohlmeyer
2024-09-11 00:30:25 -04:00
committed by GitHub
parent fde43d5079 2be54aa454
commit 0f4aedaf39
2 changed files with 150 additions and 248 deletions
Download Patch File Download Diff File Expand all files Collapse all files

383
src/KOKKOS/pair_reaxff_kokkos.cpp
View File

@ -1502,21 +1502,11 @@ void PairReaxFFKokkos<DeviceType>::allocate_array()
}
if (cut_hbsq > 0.0) {
MemKK::realloc_kokkos(d_hb_first,"reaxff/kk:hb_first",nmax);
MemKK::realloc_kokkos(d_hb_num,"reaxff/kk:hb_num",nmax);
if (((bigint) nmax*maxhb) > MAXSMALLINT)
error->one(FLERR,"Too many hydrogen bonds in pair reaxff");
MemKK::realloc_kokkos(d_hb_list,"reaxff/kk:hb_list",nmax*maxhb);
MemKK::realloc_kokkos(d_hb_list,"reaxff/kk:hb_list", nmax, maxhb);
}
MemKK::realloc_kokkos(d_bo_first,"reaxff/kk:bo_first",nmax);
MemKK::realloc_kokkos(d_bo_num,"reaxff/kk:bo_num",nmax);
if (((bigint) nmax*maxbo) > MAXSMALLINT)
error->one(FLERR,"Too many bonds in pair reaxff");
MemKK::realloc_kokkos(d_bo_list,"reaxff/kk:bo_list",nmax*maxbo);
MemKK::realloc_kokkos(d_bo_list,"reaxff/kk:bo_list", nmax, maxbo);
MemKK::realloc_kokkos(d_BO,"reaxff/kk:BO",nmax,maxbo);
MemKK::realloc_kokkos(d_BO_s,"reaxff/kk:BO",nmax,maxbo);
@ -1547,9 +1537,9 @@ void PairReaxFFKokkos<DeviceType>::allocate_array()
MemKK::realloc_kokkos(d_Deltap,"reaxff/kk:Deltap",nmax);
MemKK::realloc_kokkos(d_total_bo,"reaxff/kk:total_bo",nmax);
MemKK::realloc_kokkos(d_Cdbo,"reaxff/kk:Cdbo",nmax,3*maxbo);
MemKK::realloc_kokkos(d_Cdbopi,"reaxff/kk:Cdbopi",nmax,3*maxbo);
MemKK::realloc_kokkos(d_Cdbopi2,"reaxff/kk:Cdbopi2",nmax,3*maxbo);
MemKK::realloc_kokkos(d_Cdbo,"reaxff/kk:Cdbo",nmax,maxbo);
MemKK::realloc_kokkos(d_Cdbopi,"reaxff/kk:Cdbopi",nmax,maxbo);
MemKK::realloc_kokkos(d_Cdbopi2,"reaxff/kk:Cdbopi2",nmax,maxbo);
MemKK::realloc_kokkos(d_Delta,"reaxff/kk:Delta",nmax);
MemKK::realloc_kokkos(d_Delta_boc,"reaxff/kk:Delta_boc",nmax);
@ -1606,19 +1596,10 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxBuildListsHalfBlocking<
F_FLOAT dDeltap_self_i[3] = {0.0,0.0,0.0};
F_FLOAT total_bo_i = 0.0;
d_bo_first[i] = i*maxbo;
const int bo_first_i = d_bo_first[i];
int ihb = -1;
int hb_first_i;
if (cut_hbsq > 0.0) {
if (cut_hbsq > 0.0)
ihb = paramssing(itype).p_hbond;
if (ihb == 1) {
d_hb_first[i] = i*maxhb;
hb_first_i = d_hb_first[i];
}
}
int nnz;
blocking_t selected_jj[blocksize];
@ -1632,9 +1613,6 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxBuildListsHalfBlocking<
int j = d_neighbors(i,jj);
j &= NEIGHMASK;
d_bo_first[j] = j*maxbo;
d_hb_first[j] = j*maxhb;
delij[0] = x(j,0) - xtmp;
delij[1] = x(j,1) - ytmp;
delij[2] = x(j,2) - ztmp;
@ -1663,7 +1641,7 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxBuildListsHalfBlocking<
const F_FLOAT rsq = delij[0]*delij[0] + delij[1]*delij[1] + delij[2]*delij[2];
// hbond list
build_hb_list<NEIGHFLAG>(rsq, i, hb_first_i, ihb, j, jtype);
build_hb_list<NEIGHFLAG>(rsq, i, ihb, j, jtype);
if (rsq > cut_bosq) continue;
@ -1680,23 +1658,23 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxBuildListsHalfBlocking<
BO = BO_s + BO_pi + BO_pi2;
if (BO < bo_cut) continue;
int ii_index = -1;
int jj_index = -1;
if (build_bo_list<NEIGHFLAG>(bo_first_i, i, j, ii_index, jj_index)) {
int i_index = -1;
int j_index = -1;
if (build_bo_list<NEIGHFLAG>(i, j, i_index, j_index)) {
// from BondOrder1
d_BO(i,jj_index) = BO;
d_BO_s(i,jj_index) = BO_s;
d_BO(i,j_index) = BO;
d_BO_s(i,j_index) = BO_s;
d_BO(j,ii_index) = BO;
d_BO_s(j,ii_index) = BO_s;
d_BO(j,i_index) = BO;
d_BO_s(j,i_index) = BO_s;
d_BO_pi(j,ii_index) = BO_pi;
d_BO_pi2(j,ii_index) = BO_pi2;
d_BO_pi(j,i_index) = BO_pi;
d_BO_pi2(j,i_index) = BO_pi2;
d_BO_pi(i,jj_index) = BO_pi;
d_BO_pi2(i,jj_index) = BO_pi2;
d_BO_pi(i,j_index) = BO_pi;
d_BO_pi2(i,j_index) = BO_pi2;
F_FLOAT Cln_BOp_s = p_bo2 * C12 / rij / rij;
F_FLOAT Cln_BOp_pi = p_bo4 * C34 / rij / rij;
@ -1709,18 +1687,18 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxBuildListsHalfBlocking<
for (int d = 0; d < 3; d++) dDeltap_self_i[d] += dBOp_i[d];
for (int d = 0; d < 3; d++) a_dDeltap_self(j,d) += -dBOp_i[d];
d_dln_BOp_pi(i,jj_index) = -(BO_pi*Cln_BOp_pi);
d_dln_BOp_pi(j,ii_index) = -(BO_pi*Cln_BOp_pi);
d_dln_BOp_pi(i,j_index) = -(BO_pi*Cln_BOp_pi);
d_dln_BOp_pi(j,i_index) = -(BO_pi*Cln_BOp_pi);
d_dln_BOp_pi2(i,jj_index) = -(BO_pi2*Cln_BOp_pi2);
d_dln_BOp_pi2(j,ii_index) = -(BO_pi2*Cln_BOp_pi2);
d_dln_BOp_pi2(i,j_index) = -(BO_pi2*Cln_BOp_pi2);
d_dln_BOp_pi2(j,i_index) = -(BO_pi2*Cln_BOp_pi2);
d_dBOp(i,jj_index) = -(BO_s*Cln_BOp_s+BO_pi*Cln_BOp_pi+BO_pi2*Cln_BOp_pi2);
d_dBOp(j,ii_index) = -(BO_s*Cln_BOp_s+BO_pi*Cln_BOp_pi+BO_pi2*Cln_BOp_pi2);
d_BO(i,jj_index) = BO - bo_cut;
d_BO(j,ii_index) = BO - bo_cut;
d_BO_s(i,jj_index) = BO_s - bo_cut;
d_BO_s(j,ii_index) = BO_s - bo_cut;
d_dBOp(i,j_index) = -(BO_s*Cln_BOp_s+BO_pi*Cln_BOp_pi+BO_pi2*Cln_BOp_pi2);
d_dBOp(j,i_index) = -(BO_s*Cln_BOp_s+BO_pi*Cln_BOp_pi+BO_pi2*Cln_BOp_pi2);
d_BO(i,j_index) = BO - bo_cut;
d_BO(j,i_index) = BO - bo_cut;
d_BO_s(i,j_index) = BO_s - bo_cut;
d_BO_s(j,i_index) = BO_s - bo_cut;
total_bo_i += (BO - bo_cut);
a_total_bo[j] += (BO - bo_cut);
}
@ -1750,19 +1728,10 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxBuildListsHalfBlockingP
F_FLOAT C12, C34, C56, BO_s, BO_pi, BO_pi2, BO, delij[3];
d_bo_first[i] = i*maxbo;
const int bo_first_i = d_bo_first[i];
int ihb = -1;
int hb_first_i;
if (cut_hbsq > 0.0) {
if (cut_hbsq > 0.0)
ihb = paramssing(itype).p_hbond;
if (ihb == 1) {
d_hb_first[i] = i*maxhb;
hb_first_i = d_hb_first[i];
}
}
int nnz;
blocking_t selected_jj[blocksize];
@ -1780,9 +1749,6 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxBuildListsHalfBlockingP
int j = d_neighbors(i,jj);
j &= NEIGHMASK;
d_bo_first[j] = j*maxbo;
d_hb_first[j] = j*maxhb;
delij[0] = x(j,0) - xtmp;
delij[1] = x(j,1) - ytmp;
delij[2] = x(j,2) - ztmp;
@ -1808,7 +1774,7 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxBuildListsHalfBlockingP
const F_FLOAT rsq = delij[0]*delij[0] + delij[1]*delij[1] + delij[2]*delij[2];
// hbond list
build_hb_list<NEIGHFLAG>(rsq, i, hb_first_i, ihb, j, jtype);
build_hb_list<NEIGHFLAG>(rsq, i, ihb, j, jtype);
if (rsq > cut_bosq) continue;
@ -1825,9 +1791,9 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxBuildListsHalfBlockingP
BO = BO_s + BO_pi + BO_pi2;
if (BO < bo_cut) continue;
int ii_index = -1;
int jj_index = -1;
build_bo_list<NEIGHFLAG>(bo_first_i, i, j, ii_index, jj_index);
int i_index = -1;
int j_index = -1;
build_bo_list<NEIGHFLAG>(i, j, i_index, j_index);
}
}
}
@ -1848,26 +1814,15 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxBuildListsHalfPreview<N
F_FLOAT C12, C34, C56, BO_s, BO_pi, BO_pi2, BO, delij[3];
d_bo_first[i] = i*maxbo;
const int bo_first_i = d_bo_first[i];
int ihb = -1;
int hb_first_i;
if (cut_hbsq > 0.0) {
if (cut_hbsq > 0.0)
ihb = paramssing(itype).p_hbond;
if (ihb == 1) {
d_hb_first[i] = i*maxhb;
hb_first_i = d_hb_first[i];
}
}
for (int jj = 0; jj < jnum; jj++) {
int j = d_neighbors(i,jj);
j &= NEIGHMASK;
d_bo_first[j] = j*maxbo;
d_hb_first[j] = j*maxhb;
const int jtype = type(j);
delij[0] = x(j,0) - xtmp;
@ -1876,7 +1831,7 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxBuildListsHalfPreview<N
const F_FLOAT rsq = delij[0]*delij[0] + delij[1]*delij[1] + delij[2]*delij[2];
// hbond list
build_hb_list<NEIGHFLAG>(rsq, i, hb_first_i, ihb, j, jtype);
build_hb_list<NEIGHFLAG>(rsq, i, ihb, j, jtype);
if (rsq > cut_bosq) continue;
@ -1893,10 +1848,10 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxBuildListsHalfPreview<N
BO = BO_s + BO_pi + BO_pi2;
if (BO < bo_cut) continue;
int ii_index = -1;
int jj_index = -1;
int i_index = -1;
int j_index = -1;
build_bo_list<NEIGHFLAG>(bo_first_i, i, j, ii_index, jj_index);
build_bo_list<NEIGHFLAG>(i, j, i_index, j_index);
}
}
@ -1905,7 +1860,7 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxBuildListsHalfPreview<N
template<class DeviceType>
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void PairReaxFFKokkos<DeviceType>::build_hb_list(F_FLOAT rsq, int i, int hb_first_i, int ihb, int j, int jtype) const {
void PairReaxFFKokkos<DeviceType>::build_hb_list(F_FLOAT rsq, int i, int ihb, int j, int jtype) const {
int i_index, j_index;
int jhb = -1;
@ -1913,30 +1868,26 @@ void PairReaxFFKokkos<DeviceType>::build_hb_list(F_FLOAT rsq, int i, int hb_firs
jhb = paramssing(jtype).p_hbond;
if (ihb == 1 && jhb == 2) {
if (NEIGHFLAG == HALF) {
j_index = hb_first_i + d_hb_num[i];
j_index = d_hb_num[i];
d_hb_num[i]++;
} else
j_index = hb_first_i + Kokkos::atomic_fetch_add(&d_hb_num[i],1);
j_index = Kokkos::atomic_fetch_add(&d_hb_num[i],1);
const int jj_index = j_index - hb_first_i;
if (jj_index >= maxhb)
d_resize_hb() = MAX(d_resize_hb(),jj_index+1);
if (j_index >= maxhb)
d_resize_hb() = MAX(d_resize_hb(), j_index+1);
else
d_hb_list[j_index] = j;
d_hb_list(i, j_index) = j;
} else if (j < nlocal && ihb == 2 && jhb == 1) {
if (NEIGHFLAG == HALF) {
i_index = d_hb_first[j] + d_hb_num[j];
i_index = d_hb_num[j];
d_hb_num[j]++;
} else
i_index = d_hb_first[j] + Kokkos::atomic_fetch_add(&d_hb_num[j],1);
i_index = Kokkos::atomic_fetch_add(&d_hb_num[j],1);
const int ii_index = i_index - d_hb_first[j];
if (ii_index >= maxhb)
d_resize_hb() = MAX(d_resize_hb(),ii_index+1);
if (i_index >= maxhb)
d_resize_hb() = MAX(d_resize_hb(), i_index+1);
else
d_hb_list[i_index] = i;
d_hb_list(j, i_index) = i;
}
}
@ -1947,31 +1898,27 @@ void PairReaxFFKokkos<DeviceType>::build_hb_list(F_FLOAT rsq, int i, int hb_firs
template<class DeviceType>
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
bool PairReaxFFKokkos<DeviceType>::build_bo_list(int bo_first_i, int i, int j, int& ii_index, int& jj_index) const {
int i_index, j_index;
bool PairReaxFFKokkos<DeviceType>::build_bo_list(int i, int j, int& i_index, int& j_index) const {
if (NEIGHFLAG == HALF) {
j_index = bo_first_i + d_bo_num[i];
i_index = d_bo_first[j] + d_bo_num[j];
j_index = d_bo_num[i];
i_index = d_bo_num[j];
d_bo_num[i]++;
d_bo_num[j]++;
} else {
j_index = bo_first_i + Kokkos::atomic_fetch_add(&d_bo_num[i],1);
i_index = d_bo_first[j] + Kokkos::atomic_fetch_add(&d_bo_num[j],1);
j_index = Kokkos::atomic_fetch_add(&d_bo_num[i],1);
i_index = Kokkos::atomic_fetch_add(&d_bo_num[j],1);
}
jj_index = j_index - bo_first_i;
ii_index = i_index - d_bo_first[j];
bool set_dB_flag = true;
if (jj_index >= maxbo || ii_index >= maxbo) {
const int max_val = MAX(ii_index+1,jj_index+1);
if (j_index >= maxbo || i_index >= maxbo) {
const int max_val = MAX(i_index + 1, j_index + 1);
d_resize_bo() = MAX(d_resize_bo(),max_val);
set_dB_flag = false;
} else {
d_bo_list[j_index] = j;
d_bo_list[i_index] = i;
d_bo_list(i, j_index) = j;
d_bo_list(j, i_index) = i;
set_dB_flag = true;
}
@ -1995,13 +1942,11 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxBuildListsFull, const i
F_FLOAT dDeltap_self_i[3] = {0.0,0.0,0.0};
F_FLOAT total_bo_i = 0.0;
const int j_start = d_bo_first[i];
const int j_end = j_start + d_bo_num[i];
for (int jj = j_start; jj < j_end; jj++) {
int j = d_bo_list[jj];
const int jnum = d_bo_num[i];
for (int j_index = 0; j_index < jnum; j_index++) {
int j = d_bo_list(i, j_index);
j &= NEIGHMASK;
const int jtype = type(j);
const int j_index = jj - j_start;
delij[0] = x(j,0) - xtmp;
delij[1] = x(j,1) - ytmp;
delij[2] = x(j,2) - ztmp;
@ -2110,23 +2055,19 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxBondOrder2, const int &
const int i = d_ilist[ii];
const int itype = type(i);
const int j_start = d_bo_first[i];
const int j_end = j_start + d_bo_num[i];
const int jnum = d_bo_num[i];
const F_FLOAT val_i = paramssing(itype).valency;
d_total_bo[i] = 0.0;
F_FLOAT total_bo = 0.0;
for (int jj = j_start; jj < j_end; jj++) {
int j = d_bo_list[jj];
for (int j_index = 0; j_index < jnum; j_index++) {
int j = d_bo_list(i, j_index);
j &= NEIGHMASK;
const int jtype = type(j);
const int j_index = jj - j_start;
const int i_index = maxbo+j_index;
// calculate corrected BO and total bond order
const F_FLOAT val_j = paramssing(jtype).valency;
const F_FLOAT ovc = paramstwbp(itype,jtype).ovc;
const F_FLOAT v13cor = paramstwbp(itype,jtype).v13cor;
@ -2221,10 +2162,6 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxBondOrder2, const int &
d_Cdbo(i,j_index) = 0.0;
d_Cdbopi(i,j_index) = 0.0;
d_Cdbopi2(i,j_index) = 0.0;
d_Cdbo(j,i_index) = 0.0;
d_Cdbopi(j,i_index) = 0.0;
d_Cdbopi2(j,i_index) = 0.0;
d_CdDelta[j] = 0.0;
}
d_CdDelta[i] = 0.0;
@ -2279,20 +2216,18 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxComputeMulti1, const in
if (imass > 21.0) dfvl = 0.0;
else dfvl = 1.0;
const int j_start = d_bo_first[i];
const int j_end = j_start + d_bo_num[i];
const int jnum = d_bo_num[i];
F_FLOAT sum_ovun1 = 0.0;
F_FLOAT sum_ovun2 = 0.0;
for (int jj = j_start; jj < j_end; jj++) {
int j = d_bo_list[jj];
for (int j_index = 0; j_index < jnum; j_index++) {
int j = d_bo_list(i, j_index);
j &= NEIGHMASK;
const int jtype = type(j);
const int j_index = jj - j_start;
sum_ovun1 += paramstwbp(itype,jtype).p_ovun1 * paramstwbp(itype,jtype).De_s * d_BO(i,j_index);
sum_ovun2 += (d_Delta[j] - dfvl * d_Delta_lp_temp[j]) * (d_BO_pi(i,j_index) + d_BO_pi2(i,j_index));
sum_ovun2 += (d_Delta[j] - dfvl * d_Delta_lp_temp[j]) * (d_BO_pi(i, j_index) + d_BO_pi2(i,j_index));
}
d_sum_ovun(i,1) += sum_ovun1;
d_sum_ovun(i,2) += sum_ovun2;
@ -2402,16 +2337,14 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxComputeMulti2<NEIGHFLAG
if (numbonds > 0 || enobondsflag)
a_CdDelta[i] += CEunder3;
const int j_start = d_bo_first[i];
const int j_end = j_start + d_bo_num[i];
const int jnum = d_bo_num[i];
F_FLOAT CdDelta_i = 0.0;
for (int jj = j_start; jj < j_end; jj++) {
int j = d_bo_list[jj];
for (int j_index = 0; j_index < jnum; j_index++) {
int j = d_bo_list(i, j_index);
j &= NEIGHMASK;
const int jtype = type(j);
const F_FLOAT jmass = paramssing(jtype).mass;
const int j_index = jj - j_start;
const F_FLOAT De_s = paramstwbp(itype,jtype).De_s;
// multibody lone pair: correction for C2
@ -2461,24 +2394,23 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxCountAngularTorsion<POP
const int i = d_ilist[ii];
const int itype = type(i);
const int j_start = d_bo_first[i];
const int j_end = j_start + d_bo_num[i];
const int jnum = d_bo_num[i];
if (POPULATE) {
// Computes and stores SBO2, CSBO2, dSBO1, dSBO2
compute_angular_sbo(i, itype, j_start, j_end);
compute_angular_sbo(i, itype, jnum);
}
// Angular
// Count buffer size for `i`
int location_angular = 0; // dummy declaration
int count_angular = preprocess_angular<false>(i, itype, j_start, j_end, location_angular);
int count_angular = preprocess_angular<false>(i, itype, jnum, location_angular);
location_angular = Kokkos::atomic_fetch_add(&d_count_angular_torsion(0), count_angular);
if (POPULATE) {
// Fill buffer for `i`
preprocess_angular<true>(i, itype, j_start, j_end, location_angular);
preprocess_angular<true>(i, itype, jnum, location_angular);
}
// Torsion
@ -2490,12 +2422,12 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxCountAngularTorsion<POP
// Count buffer size for `i`
int location_torsion = 0; // dummy declaration
int count_torsion = preprocess_torsion<false>(i, itype, itag, xtmp, ytmp, ztmp, j_start, j_end, location_torsion);
int count_torsion = preprocess_torsion<false>(i, itype, itag, xtmp, ytmp, ztmp, jnum, location_torsion);
location_torsion = Kokkos::atomic_fetch_add(&d_count_angular_torsion(1), count_torsion);
if (POPULATE) {
// Fill buffer for `i`
preprocess_torsion<true>(i, itype, itag, xtmp, ytmp, ztmp, j_start, j_end, location_torsion);
preprocess_torsion<true>(i, itype, itag, xtmp, ytmp, ztmp, jnum, location_torsion);
}
}
@ -2504,7 +2436,7 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxCountAngularTorsion<POP
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairReaxFFKokkos<DeviceType>::compute_angular_sbo(int i, int itype, int j_start, int j_end) const {
void PairReaxFFKokkos<DeviceType>::compute_angular_sbo(int i, int itype, int jnum) const {
F_FLOAT SBO2, CSBO2, dSBO1, dSBO2;
@ -2514,8 +2446,7 @@ void PairReaxFFKokkos<DeviceType>::compute_angular_sbo(int i, int itype, int j_s
F_FLOAT SBOp = 0.0;
F_FLOAT prod_SBO = 1.0;
for (int jj = j_start; jj < j_end; jj++) {
const int j_index = jj - j_start;
for (int j_index = 0; j_index < jnum; j_index++) {
const F_FLOAT bo_ij = d_BO(i,j_index);
SBOp += (d_BO_pi(i,j_index) + d_BO_pi2(i,j_index));
@ -2568,29 +2499,26 @@ void PairReaxFFKokkos<DeviceType>::compute_angular_sbo(int i, int itype, int j_s
template<class DeviceType>
template<bool POPULATE>
KOKKOS_INLINE_FUNCTION
int PairReaxFFKokkos<DeviceType>::preprocess_angular(int i, int itype, int j_start, int j_end, int location_angular) const {
int PairReaxFFKokkos<DeviceType>::preprocess_angular(int i, int itype, int jnum, int location_angular) const {
int count_angular = 0;
for (int jj = j_start; jj < j_end; jj++) {
int j = d_bo_list[jj];
for (int j_index = 0; j_index < jnum; j_index++) {
int j = d_bo_list(i, j_index);
j &= NEIGHMASK;
const int j_index = jj - j_start;
const F_FLOAT bo_ij = d_BO(i,j_index);
if (bo_ij <= thb_cut) continue;
if (i >= nlocal && j >= nlocal) continue;
const int i_index = maxbo + j_index;
const int jtype = type(j);
for (int kk = jj+1; kk < j_end; kk++) {
for (int k_index = j_index + 1; k_index < jnum; k_index++) {
//for (int kk = j_start; kk < j_end; kk++) {
int k = d_bo_list[kk];
int k = d_bo_list(i, k_index);
k &= NEIGHMASK;
if (k == j) continue;
const int k_index = kk - j_start;
const F_FLOAT bo_ik = d_BO(i,k_index);
if (bo_ij <= thb_cut || bo_ik <= thb_cut || bo_ij * bo_ik <= thb_cutsq) continue;
@ -2608,14 +2536,14 @@ int PairReaxFFKokkos<DeviceType>::preprocess_angular(int i, int itype, int j_sta
pack.i0 = i;
pack.i1 = j;
pack.i2 = k;
pack.i3 = j_start;
pack.i3 = jnum;
d_angular_pack(location_angular, 0) = pack;
// Second pack stores i_index, j_index, k_index, and j_end
pack.i0 = i_index;
// Second pack stores j_index and k_index
// i0 is unused because there's no i_index
pack.i1 = j_index;
pack.i2 = k_index;
pack.i3 = j_end;
// i3 is unused
d_angular_pack(location_angular, 1) = pack;
location_angular++;
@ -2634,17 +2562,16 @@ template<class DeviceType>
template<bool POPULATE>
KOKKOS_INLINE_FUNCTION
int PairReaxFFKokkos<DeviceType>::preprocess_torsion(int i, int /*itype*/, tagint itag,
F_FLOAT xtmp, F_FLOAT ytmp, F_FLOAT ztmp, int j_start, int j_end, int location_torsion) const {
F_FLOAT xtmp, F_FLOAT ytmp, F_FLOAT ztmp, int jknum, int location_torsion) const {
// in reaxff_torsion_angles: j = i, k = j, i = k;
int count_torsion = 0;
for (int jj = j_start; jj < j_end; jj++) {
int j = d_bo_list[jj];
for (int j_index = 0; j_index < jknum; j_index++) {
int j = d_bo_list(i, j_index);
j &= NEIGHMASK;
const tagint jtag = tag(j);
const int j_index = jj - j_start;
// skip half of the interactions
if (itag > jtag) {
@ -2660,23 +2587,20 @@ int PairReaxFFKokkos<DeviceType>::preprocess_torsion(int i, int /*itype*/, tagin
const F_FLOAT bo_ij = d_BO(i,j_index);
if (bo_ij < thb_cut) continue;
const int l_start = d_bo_first[j];
const int l_end = l_start + d_bo_num[j];
const int lnum = d_bo_num[j];
for (int kk = j_start; kk < j_end; kk++) {
int k = d_bo_list[kk];
for (int k_index = 0; k_index < jknum; k_index++) {
int k = d_bo_list(i, k_index);
k &= NEIGHMASK;
if (k == j) continue;
const int k_index = kk - j_start;
const F_FLOAT bo_ik = d_BO(i,k_index);
if (bo_ik < thb_cut) continue;
for (int ll = l_start; ll < l_end; ll++) {
int l = d_bo_list[ll];
for (int l_index = 0; l_index < lnum; l_index++) {
int l = d_bo_list(j, l_index);
l &= NEIGHMASK;
if (l == i) continue;
const int l_index = ll - l_start;
const F_FLOAT bo_jl = d_BO(j,l_index);
if (l == k || bo_jl < thb_cut || bo_ij*bo_ik*bo_jl < thb_cut) continue;
@ -2716,9 +2640,9 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxComputeAngularPreproces
auto v_f = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_f),decltype(ndup_f)>::get(dup_f,ndup_f);
auto a_f = v_f.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
Kokkos::View<F_FLOAT**, typename DAT::t_ffloat_2d_dl::array_layout,typename KKDevice<DeviceType>::value,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value>> a_Cdbo = d_Cdbo;
Kokkos::View<F_FLOAT**, typename DAT::t_ffloat_2d_dl::array_layout,typename KKDevice<DeviceType>::value,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value>> a_Cdbopi = d_Cdbopi;
Kokkos::View<F_FLOAT**, typename DAT::t_ffloat_2d_dl::array_layout,typename KKDevice<DeviceType>::value,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value>> a_Cdbopi2 = d_Cdbopi2;
Kokkos::View<F_FLOAT**, typename decltype(d_Cdbo)::array_layout,KKDeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value>> a_Cdbo = d_Cdbo;
Kokkos::View<F_FLOAT**, typename decltype(d_Cdbopi)::array_layout,KKDeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value>> a_Cdbopi = d_Cdbopi;
Kokkos::View<F_FLOAT**, typename decltype(d_Cdbopi2)::array_layout,KKDeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value>> a_Cdbopi2 = d_Cdbopi2;
auto v_CdDelta = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_CdDelta),decltype(ndup_CdDelta)>::get(dup_CdDelta,ndup_CdDelta);
auto a_CdDelta = v_CdDelta.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
@ -2758,13 +2682,13 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxComputeAngularPreproces
const int i = pack.i0;
const int j = pack.i1;
const int k = pack.i2;
const int j_start = pack.i3;
const int jnum = pack.i3;
pack = d_angular_pack(apack, 1);
const int i_index = pack.i0;
// i0 is unused
const int j_index = pack.i1;
const int k_index = pack.i2;
const int j_end = pack.i3;
// i3 is unused
const int itype = type(i);
const X_FLOAT xtmp = x(i,0);
@ -2914,17 +2838,13 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxComputeAngularPreproces
// Forces
a_Cdbo(i,j_index) += (CEval1 + CEpen2 + (CEcoa1 - CEcoa4));
a_Cdbo(j,i_index) += (CEval1 + CEpen2 + (CEcoa1 - CEcoa4));
a_Cdbo(i,k_index) += (CEval2 + CEpen3 + (CEcoa2 - CEcoa5));
a_Cdbo(k,i_index) += (CEval2 + CEpen3 + (CEcoa2 - CEcoa5));
CdDelta_i += ((CEval3 + CEval7) + CEpen1 + CEcoa3);
CdDelta_j += CEcoa4;
a_CdDelta[k] += CEcoa5;
for (int ll = j_start; ll < j_end; ll++) {
const int l_index = ll - j_start;
for (int l_index = 0; l_index < jnum; l_index++) {
temp_bo_jt = d_BO(i,l_index);
temp = temp_bo_jt * temp_bo_jt * temp_bo_jt;
pBOjt7 = temp * temp * temp_bo_jt;
@ -2978,8 +2898,8 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxComputeTorsionPreproces
auto v_CdDelta = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_CdDelta),decltype(ndup_CdDelta)>::get(dup_CdDelta,ndup_CdDelta);
auto a_CdDelta = v_CdDelta.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
Kokkos::View<F_FLOAT**, typename DAT::t_ffloat_2d_dl::array_layout,typename KKDevice<DeviceType>::value,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value>> a_Cdbo = d_Cdbo;
Kokkos::View<F_FLOAT**, typename DAT::t_ffloat_2d_dl::array_layout,typename KKDevice<DeviceType>::value,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value>> a_Cdbopi = d_Cdbopi;
Kokkos::View<F_FLOAT**, typename decltype(d_Cdbo)::array_layout,KKDeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value>> a_Cdbo = d_Cdbo;
Kokkos::View<F_FLOAT**, typename decltype(d_Cdbopi)::array_layout,KKDeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value>> a_Cdbopi = d_Cdbopi;
//auto a_Cdbo = dup_Cdbo.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
// in reaxff_torsion_angles: j = i, k = j, i = k;
@ -3340,21 +3260,18 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxComputeHydrogen<NEIGHFL
const X_FLOAT ytmp = x(i,1);
const X_FLOAT ztmp = x(i,2);
const int j_start = d_bo_first[i];
const int j_end = j_start + d_bo_num[i];
const int k_start = d_hb_first[i];
const int k_end = k_start + d_hb_num[i];
const int jnum = d_bo_num[i];
const int knum = d_hb_num[i];
int top = 0;
for (int jj = j_start; jj < j_end; jj++) {
int j = d_bo_list[jj];
for (int j_index = 0; j_index < jnum; j_index++) {
int j = d_bo_list(i, j_index);
j &= NEIGHMASK;
const int jtype = type(j);
const int j_index = jj - j_start;
const F_FLOAT bo_ij = d_BO(i,j_index);
if (paramssing(jtype).p_hbond == 2 && bo_ij >= HB_THRESHOLD) {
hblist[top] = jj;
hblist[top] = j_index;
top ++;
}
}
@ -3362,8 +3279,8 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxComputeHydrogen<NEIGHFL
F_FLOAT fitmp[3];
for (int d = 0; d < 3; d++) fitmp[d] = 0.0;
for (int kk = k_start; kk < k_end; kk++) {
int k = d_hb_list[kk];
for (int k_index = 0; k_index < knum; k_index++) {
int k = d_hb_list(i, k_index);
k &= NEIGHMASK;
const tagint ktag = tag(k);
const int ktype = type(k);
@ -3375,14 +3292,13 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxComputeHydrogen<NEIGHFL
const F_FLOAT rik = sqrt(rsqik);
for (int itr = 0; itr < top; itr++) {
const int jj = hblist[itr];
int j = d_bo_list[jj];
const int j_index = hblist[itr];
int j = d_bo_list(i, j_index);
j &= NEIGHMASK;
const tagint jtag = tag(j);
if (jtag == ktag) continue;
const int jtype = type(j);
const int j_index = jj - j_start;
const F_FLOAT bo_ij = d_BO(i,j_index);
delij[0] = x(j,0) - xtmp;
@ -3470,20 +3386,19 @@ template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxUpdateBond<NEIGHFLAG>, const int &ii) const {
Kokkos::View<F_FLOAT**, typename DAT::t_ffloat_2d_dl::array_layout,KKDeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value>> a_Cdbo = d_Cdbo;
Kokkos::View<F_FLOAT**, typename DAT::t_ffloat_2d_dl::array_layout,KKDeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value>> a_Cdbopi = d_Cdbopi;
Kokkos::View<F_FLOAT**, typename DAT::t_ffloat_2d_dl::array_layout,KKDeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value>> a_Cdbopi2 = d_Cdbopi2;
Kokkos::View<F_FLOAT**, typename decltype(d_Cdbo)::array_layout,KKDeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value>> a_Cdbo = d_Cdbo;
Kokkos::View<F_FLOAT**, typename decltype(d_Cdbopi)::array_layout,KKDeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value>> a_Cdbopi = d_Cdbopi;
Kokkos::View<F_FLOAT**, typename decltype(d_Cdbopi2)::array_layout,KKDeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value>> a_Cdbopi2 = d_Cdbopi2;
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 tagint itag = tag(i);
const int j_start = d_bo_first[i];
const int j_end = j_start + d_bo_num[i];
const int jnum = d_bo_num[i];
for (int jj = j_start; jj < j_end; jj++) {
int j = d_bo_list[jj];
for (int j_index = 0; j_index < jnum; j_index++) {
int j = d_bo_list(i, j_index);
j &= NEIGHMASK;
const tagint jtag = tag(j);
@ -3501,19 +3416,16 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxUpdateBond<NEIGHFLAG>,
if (!flag) continue;
const int j_index = jj - j_start;
const F_FLOAT Cdbo_i = d_Cdbo(i,j_index);
const F_FLOAT Cdbopi_i = d_Cdbopi(i,j_index);
const F_FLOAT Cdbopi2_i = d_Cdbopi2(i,j_index);
const int k_start = d_bo_first[j];
const int k_end = k_start + d_bo_num[j];
const int knum = d_bo_num[j];
for (int kk = k_start; kk < k_end; kk++) {
int k = d_bo_list[kk];
for (int k_index = 0; k_index < knum; k_index++) {
int k = d_bo_list(j, k_index);
k &= NEIGHMASK;
if (k != i) continue;
const int k_index = kk - k_start;
a_Cdbo(j,k_index) += Cdbo_i;
a_Cdbopi(j,k_index) += Cdbopi_i;
@ -3541,13 +3453,12 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxComputeBond1<NEIGHFLAG,
const int itype = type(i);
const tagint itag = tag(i);
const F_FLOAT imass = paramssing(itype).mass;
const int j_start = d_bo_first[i];
const int j_end = j_start + d_bo_num[i];
const int jnum = d_bo_num[i];
F_FLOAT CdDelta_i = 0.0;
for (int jj = j_start; jj < j_end; jj++) {
int j = d_bo_list[jj];
for (int j_index = 0; j_index < jnum; j_index++) {
int j = d_bo_list(i, j_index);
j &= NEIGHMASK;
const tagint jtag = tag(j);
@ -3562,7 +3473,6 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxComputeBond1<NEIGHFLAG,
}
const int jtype = type(j);
const int j_index = jj - j_start;
const F_FLOAT jmass = paramssing(jtype).mass;
// bond energy (nlocal only)
@ -3655,15 +3565,14 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxComputeBond2<NEIGHFLAG,
const X_FLOAT ytmp = x(i,1);
const X_FLOAT ztmp = x(i,2);
const tagint itag = tag(i);
const int j_start = d_bo_first[i];
const int j_end = j_start + d_bo_num[i];
const int jknum = d_bo_num[i];
F_FLOAT CdDelta_i = d_CdDelta[i];
F_FLOAT fitmp[3];
for (int j = 0; j < 3; j++) fitmp[j] = 0.0;
for (int jj = j_start; jj < j_end; jj++) {
int j = d_bo_list[jj];
for (int j_index = 0; j_index < jknum; j_index++) {
int j = d_bo_list(i, j_index);
j &= NEIGHMASK;
const tagint jtag = tag(j);
@ -3677,15 +3586,13 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxComputeBond2<NEIGHFLAG,
if (x(j,2) == ztmp && x(j,1) == ytmp && x(j,0) < xtmp) continue;
}
const int j_index = jj - j_start;
F_FLOAT CdDelta_j = d_CdDelta[j];
delij[0] = x(j,0) - xtmp;
delij[1] = x(j,1) - ytmp;
delij[2] = x(j,2) - ztmp;
const int k_start = d_bo_first[j];
const int k_end = k_start + d_bo_num[j];
const int knum = d_bo_num[j];
F_FLOAT coef_C1dbo, coef_C2dbo, coef_C3dbo, coef_C1dbopi, coef_C2dbopi, coef_C3dbopi, coef_C4dbopi;
F_FLOAT coef_C1dbopi2, coef_C2dbopi2, coef_C3dbopi2, coef_C4dbopi2, coef_C1dDelta, coef_C2dDelta, coef_C3dDelta;
@ -3769,10 +3676,9 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxComputeBond2<NEIGHFLAG,
}
// forces on k: i neighbor
for (int kk = j_start; kk < j_end; kk++) {
int k = d_bo_list[kk];
for (int k_index = 0; k_index < jknum; k_index++) {
int k = d_bo_list(i, k_index);
k &= NEIGHMASK;
const int k_index = kk - j_start;
delik[0] = x(k,0) - xtmp;
delik[1] = x(k,1) - ytmp;
@ -3799,10 +3705,9 @@ void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxComputeBond2<NEIGHFLAG,
}
// forces on k: j neighbor
for (int kk = k_start; kk < k_end; kk++) {
int k = d_bo_list[kk];
for (int k_index = 0; k_index < knum; k_index++) {
int k = d_bo_list(j, k_index);
k &= NEIGHMASK;
const int k_index = kk - k_start;
for (int d = 0; d < 3; d++) deljk[d] = x(k,d) - x(j,d);
@ -4203,15 +4108,13 @@ void PairReaxFFKokkos<DeviceType>::calculate_find_bond_item(int ii, int &numbond
int nj = 0;
if (mask[i] & groupbit) {
const int j_start = d_bo_first[i];
const int j_end = j_start + d_bo_num[i];
for (int jj = j_start; jj < j_end; jj++) {
int j = d_bo_list[jj];
const int jnum = d_bo_num[i];
for (int j_index = 0; j_index < jnum; j_index++) {
int j = d_bo_list(i, j_index);
j &= NEIGHMASK;
if (mask[j] & groupbit) {
const tagint jtag = tag[j];
const int j_index = jj - j_start;
double bo_tmp = d_BO(i,j_index);
double bo_tmp = d_BO(i, j_index);
if (bo_tmp > bo_cut_bond) {
d_neighid(i,nj) = jtag;
@ -4409,15 +4312,13 @@ KOKKOS_INLINE_FUNCTION
void PairReaxFFKokkos<DeviceType>::operator()(TagPairReaxFindBondSpecies, const int &i) const {
int nj = 0;
const int j_start = d_bo_first[i];
const int j_end = j_start + d_bo_num[i];
for (int jj = j_start; jj < j_end; jj++) {
int j = d_bo_list[jj];
const int jnum = d_bo_num[i];
for (int j_index = 0; j_index < jnum; j_index++) {
int j = d_bo_list(i, j_index);
j &= NEIGHMASK;
if (j < i) continue;
const int j_index = jj - j_start;
double bo_tmp = d_BO(i,j_index);
double bo_tmp = d_BO(i, j_index);
if (bo_tmp >= 0.10) { // Why is this a hardcoded value?
k_tmpid.view<DeviceType>()(i,nj) = j;

15
src/KOKKOS/pair_reaxff_kokkos.h
View File

@ -178,14 +178,14 @@ class PairReaxFFKokkos : public PairReaxFF {
// TagPairReaxBuildListsHalfBlocking, HalfBlockingPreview, HalfPreview
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void build_hb_list(F_FLOAT, int, int, int, int, int) const;
void build_hb_list(F_FLOAT, int, int, int, int) const;
// Isolated function that builds the bond order list, reused across
// TagPairReaxBuildListsHalfBlocking, HalfBlockingPreview, HalfPreview
// Returns if we need to populate d_d* functions or not
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
bool build_bo_list(int, int, int, int&, int&) const;
bool build_bo_list(int, int, int&, int&) const;
KOKKOS_INLINE_FUNCTION
void operator()(TagPairReaxBuildListsFull, const int&) const;
@ -245,17 +245,17 @@ class PairReaxFFKokkos : public PairReaxFF {
// Abstraction for computing SBSO2, CSBO2, dSBO1, dsBO2
KOKKOS_INLINE_FUNCTION
void compute_angular_sbo(int, int, int, int) const;
void compute_angular_sbo(int, int, int) const;
// Abstraction for counting and populating angular intermediates
template<bool POPULATE>
KOKKOS_INLINE_FUNCTION
int preprocess_angular(int, int, int, int, int) const;
int preprocess_angular(int, int, int, int) const;
// Abstraction for counting and populating torsion intermediated
template<bool POPULATE>
KOKKOS_INLINE_FUNCTION
int preprocess_torsion(int, int, tagint, F_FLOAT, F_FLOAT, F_FLOAT, int, int, int) const;
int preprocess_torsion(int, int, tagint, F_FLOAT, F_FLOAT, F_FLOAT, int, int) const;
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
@ -436,7 +436,7 @@ class PairReaxFFKokkos : public PairReaxFF {
typename AT::t_ffloat_2d_dl d_C1dbo, d_C2dbo, d_C3dbo;
typename AT::t_ffloat_2d_dl d_C1dbopi, d_C2dbopi, d_C3dbopi, d_C4dbopi;
typename AT::t_ffloat_2d_dl d_C1dbopi2, d_C2dbopi2, d_C3dbopi2, d_C4dbopi2;
typename AT::t_ffloat_2d_dl d_Cdbo, d_Cdbopi, d_Cdbopi2, d_dDeltap_self;
typename AT::t_ffloat_2d_dl d_dDeltap_self, d_Cdbo, d_Cdbopi, d_Cdbopi2;
int need_dup;
@ -470,7 +470,8 @@ class PairReaxFFKokkos : public PairReaxFF {
typename AT::t_int_1d_randomread d_ilist;
typename AT::t_int_1d_randomread d_numneigh;
typename AT::t_int_1d d_bo_first, d_bo_num, d_bo_list, d_hb_first, d_hb_num, d_hb_list;
typename AT::t_int_1d d_bo_num, d_hb_num;
typename AT::t_int_2d d_bo_list, d_hb_list;
DAT::tdual_int_scalar k_resize_bo, k_resize_hb;
typename AT::t_int_scalar d_resize_bo, d_resize_hb;