diff --git a/src/KOKKOS/pair_reax_c_kokkos.cpp b/src/KOKKOS/pair_reax_c_kokkos.cpp
index eba6fbb0b5a14abbe205992969b477cd3a997a33..eb4679a7ddff6ac134ae36ece6355f33c4ab28b6 100644
--- a/src/KOKKOS/pair_reax_c_kokkos.cpp
+++ b/src/KOKKOS/pair_reax_c_kokkos.cpp
@@ -660,6 +660,8 @@ void PairReaxCKokkos<DeviceType>::LR_vdW_Coulomb( int i, int j, double r_ij, LR_
template<class DeviceType>
void PairReaxCKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
{
+ copymode = 1;
+
bocnt = hbcnt = 0;
eflag = eflag_in;
@@ -703,8 +705,6 @@ void PairReaxCKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
pvector[i] = 0.0;
}
- copymode = 1;
-
EV_FLOAT_REAX ev;
EV_FLOAT_REAX ev_all;
@@ -1363,6 +1363,23 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxZero, const int &n) const {
d_dDeltap_self(n,j) = 0.0;
}
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void PairReaxCKokkos<DeviceType>::operator()(PairReaxZeroEAtom, const int &i) const {
+ v_eatom(i) = 0.0;
+}
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void PairReaxCKokkos<DeviceType>::operator()(PairReaxZeroVAtom, const int &i) const {
+ v_vatom(i,0) = 0.0;
+ v_vatom(i,1) = 0.0;
+ v_vatom(i,2) = 0.0;
+ v_vatom(i,3) = 0.0;
+ v_vatom(i,4) = 0.0;
+ v_vatom(i,5) = 0.0;
+}
+
/* ---------------------------------------------------------------------- */
template<class DeviceType>
@@ -2398,7 +2415,7 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeAngular<NEIGHFLAG,EV
F_FLOAT BOA_ij, BOA_ik, rij, bo_ij, bo_ik;
F_FLOAT dcos_theta_di[3], dcos_theta_dj[3], dcos_theta_dk[3];
F_FLOAT eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3];
- F_FLOAT delij[3], delik[3];
+ F_FLOAT delij[3], delik[3], delji[3], delki[3];
p_val6 = gp[14];
p_val8 = gp[33];
@@ -2648,8 +2665,10 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeAngular<NEIGHFLAG,EV
if (EVFLAG) {
eng_tmp = e_ang + e_pen + e_coa;
//if (eflag_atom) this->template ev_tally<NEIGHFLAG>(ev,i,j,eng_tmp,0.0,0.0,0.0,0.0);
+ for (int d = 0; d < 3; d++) delki[d] = -1.0 * delik[d];
+ for (int d = 0; d < 3; d++) delji[d] = -1.0 * delij[d];
if (eflag_atom) this->template e_tally<NEIGHFLAG>(ev,i,j,eng_tmp);
- if (vflag_either) this->template v_tally3<NEIGHFLAG>(ev,i,j,k,fj_tmp,fk_tmp,delij,delik);
+ if (vflag_either) this->template v_tally3<NEIGHFLAG>(ev,i,j,k,fj_tmp,fk_tmp,delji,delki);
}
}
@@ -3728,12 +3747,12 @@ void PairReaxCKokkos<DeviceType>::v_tally3(EV_FLOAT_REAX &ev, const int &i, cons
F_FLOAT v[6];
- v[0] = (drij[0]*fj[0] + drik[0]*fk[0]);
- v[1] = (drij[1]*fj[1] + drik[1]*fk[1]);
- v[2] = (drij[2]*fj[2] + drik[2]*fk[2]);
- v[3] = (drij[0]*fj[1] + drik[0]*fk[1]);
- v[4] = (drij[0]*fj[2] + drik[0]*fk[2]);
- v[5] = (drij[1]*fj[2] + drik[1]*fk[2]);
+ v[0] = drij[0]*fj[0] + drik[0]*fk[0];
+ v[1] = drij[1]*fj[1] + drik[1]*fk[1];
+ v[2] = drij[2]*fj[2] + drik[2]*fk[2];
+ v[3] = drij[0]*fj[1] + drik[0]*fk[1];
+ v[4] = drij[0]*fj[2] + drik[0]*fk[2];
+ v[5] = drij[1]*fj[2] + drik[1]*fk[2];
if (vflag_global) {
ev.v[0] += v[0];
@@ -3745,12 +3764,12 @@ void PairReaxCKokkos<DeviceType>::v_tally3(EV_FLOAT_REAX &ev, const int &i, cons
}
if (vflag_atom) {
- a_vatom(i,0) += THIRD*v[0]; a_vatom(i,1) += THIRD*v[1]; a_vatom(i,2) += THIRD*v[2];
- a_vatom(i,3) += THIRD*v[3]; a_vatom(i,4) += THIRD*v[4]; a_vatom(i,5) += THIRD*v[5];
- a_vatom(j,0) += THIRD*v[0]; a_vatom(j,1) += THIRD*v[1]; a_vatom(j,2) += THIRD*v[2];
- a_vatom(j,3) += THIRD*v[3]; a_vatom(j,4) += THIRD*v[4]; a_vatom(j,5) += THIRD*v[5];
- a_vatom(k,0) += THIRD*v[0]; a_vatom(k,1) += THIRD*v[1]; a_vatom(k,2) += THIRD*v[2];
- a_vatom(k,3) += THIRD*v[3]; a_vatom(k,4) += THIRD*v[4]; a_vatom(k,5) += THIRD*v[5];
+ a_vatom(i,0) += THIRD * v[0]; a_vatom(i,1) += THIRD * v[1]; a_vatom(i,2) += THIRD * v[2];
+ a_vatom(i,3) += THIRD * v[3]; a_vatom(i,4) += THIRD * v[4]; a_vatom(i,5) += THIRD * v[5];
+ a_vatom(j,0) += THIRD * v[0]; a_vatom(j,1) += THIRD * v[1]; a_vatom(j,2) += THIRD * v[2];
+ a_vatom(j,3) += THIRD * v[3]; a_vatom(j,4) += THIRD * v[4]; a_vatom(j,5) += THIRD * v[5];
+ a_vatom(k,0) += THIRD * v[0]; a_vatom(k,1) += THIRD * v[1]; a_vatom(k,2) += THIRD * v[2];
+ a_vatom(k,3) += THIRD * v[3]; a_vatom(k,4) += THIRD * v[4]; a_vatom(k,5) += THIRD * v[5];
}
}
@@ -3767,12 +3786,12 @@ void PairReaxCKokkos<DeviceType>::v_tally4(EV_FLOAT_REAX &ev, const int &i, cons
// The vatom array is atomic for Half/Thread neighbor style
F_FLOAT v[6];
- v[0] = 0.25 * (dril[0]*fi[0] + drjl[0]*fj[0] + drkl[0]*fk[0]);
- v[1] = 0.25 * (dril[1]*fi[1] + drjl[1]*fj[1] + drkl[1]*fk[1]);
- v[2] = 0.25 * (dril[2]*fi[2] + drjl[2]*fj[2] + drkl[2]*fk[2]);
- v[3] = 0.25 * (dril[0]*fi[1] + drjl[0]*fj[1] + drkl[0]*fk[1]);
- v[4] = 0.25 * (dril[0]*fi[2] + drjl[0]*fj[2] + drkl[0]*fk[2]);
- v[5] = 0.25 * (dril[1]*fi[2] + drjl[1]*fj[2] + drkl[1]*fk[2]);
+ v[0] = dril[0]*fi[0] + drjl[0]*fj[0] + drkl[0]*fk[0];
+ v[1] = dril[1]*fi[1] + drjl[1]*fj[1] + drkl[1]*fk[1];
+ v[2] = dril[2]*fi[2] + drjl[2]*fj[2] + drkl[2]*fk[2];
+ v[3] = dril[0]*fi[1] + drjl[0]*fj[1] + drkl[0]*fk[1];
+ v[4] = dril[0]*fi[2] + drjl[0]*fj[2] + drkl[0]*fk[2];
+ v[5] = dril[1]*fi[2] + drjl[1]*fj[2] + drkl[1]*fk[2];
if (vflag_global) {
ev.v[0] += v[0];
@@ -3785,14 +3804,14 @@ void PairReaxCKokkos<DeviceType>::v_tally4(EV_FLOAT_REAX &ev, const int &i, cons
if (vflag_atom) {
Kokkos::View<F_FLOAT*[6], typename DAT::t_virial_array::array_layout,DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > a_vatom = v_vatom;
- a_vatom(i,0) += v[0]; a_vatom(i,1) += v[1]; a_vatom(i,2) += v[2];
- a_vatom(i,3) += v[3]; a_vatom(i,4) += v[4]; a_vatom(i,5) += v[5];
- a_vatom(j,0) += v[0]; a_vatom(j,1) += v[1]; a_vatom(j,2) += v[2];
- a_vatom(j,3) += v[3]; a_vatom(j,4) += v[4]; a_vatom(j,5) += v[5];
- a_vatom(k,0) += v[0]; a_vatom(k,1) += v[1]; a_vatom(k,2) += v[2];
- a_vatom(k,3) += v[3]; a_vatom(k,4) += v[4]; a_vatom(k,5) += v[5];
- a_vatom(l,0) += v[0]; a_vatom(l,1) += v[1]; a_vatom(l,2) += v[2];
- a_vatom(l,3) += v[3]; a_vatom(l,4) += v[4]; a_vatom(l,5) += v[5];
+ a_vatom(i,0) += 0.25 * v[0]; a_vatom(i,1) += 0.25 * v[1]; a_vatom(i,2) += 0.25 * v[2];
+ a_vatom(i,3) += 0.25 * v[3]; a_vatom(i,4) += 0.25 * v[4]; a_vatom(i,5) += 0.25 * v[5];
+ a_vatom(j,0) += 0.25 * v[0]; a_vatom(j,1) += 0.25 * v[1]; a_vatom(j,2) += 0.25 * v[2];
+ a_vatom(j,3) += 0.25 * v[3]; a_vatom(j,4) += 0.25 * v[4]; a_vatom(j,5) += 0.25 * v[5];
+ a_vatom(k,0) += 0.25 * v[0]; a_vatom(k,1) += 0.25 * v[1]; a_vatom(k,2) += 0.25 * v[2];
+ a_vatom(k,3) += 0.25 * v[3]; a_vatom(k,4) += 0.25 * v[4]; a_vatom(k,5) += 0.25 * v[5];
+ a_vatom(l,0) += 0.25 * v[0]; a_vatom(l,1) += 0.25 * v[1]; a_vatom(l,2) += 0.25 * v[2];
+ a_vatom(l,3) += 0.25 * v[3]; a_vatom(l,4) += 0.25 * v[4]; a_vatom(l,5) += 0.25 * v[5];
}
}
@@ -3894,6 +3913,14 @@ void PairReaxCKokkos<DeviceType>::ev_setup(int eflag, int vflag)
if (eflag_global) eng_vdwl = eng_coul = 0.0;
if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0;
+ if (eflag_atom) {
+ Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, PairReaxZeroEAtom>(0,maxeatom),*this);
+ DeviceType::fence();
+ }
+ if (vflag_atom) {
+ Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, PairReaxZeroVAtom>(0,maxvatom),*this);
+ DeviceType::fence();
+ }
// if vflag_global = 2 and pair::compute() calls virial_fdotr_compute()
// compute global virial via (F dot r) instead of via pairwise summation
diff --git a/src/KOKKOS/pair_reax_c_kokkos.h b/src/KOKKOS/pair_reax_c_kokkos.h
index 8bba6c50d064a39ed593fd4faefaea020e60656a..28166df3bd5592dab245b9c74ca98feceb9d424d 100644
--- a/src/KOKKOS/pair_reax_c_kokkos.h
+++ b/src/KOKKOS/pair_reax_c_kokkos.h
@@ -84,6 +84,10 @@ struct PairReaxBuildListsHalf_LessAtomics{};
struct PairReaxZero{};
+struct PairReaxZeroEAtom{};
+
+struct PairReaxZeroVAtom{};
+
struct PairReaxBondOrder1{};
struct PairReaxBondOrder1_LessAtomics{};
@@ -172,6 +176,12 @@ class PairReaxCKokkos : public PairReaxC {
KOKKOS_INLINE_FUNCTION
void operator()(PairReaxZero, const int&) const;
+ KOKKOS_INLINE_FUNCTION
+ void operator()(PairReaxZeroEAtom, const int&) const;
+
+ KOKKOS_INLINE_FUNCTION
+ void operator()(PairReaxZeroVAtom, const int&) const;
+
KOKKOS_INLINE_FUNCTION
void operator()(PairReaxBondOrder1, const int&) const;