From 001dcb6bb129e48cd1b52330919fadbb48145b4e Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 14 Jul 2020 09:02:57 -0600 Subject: [PATCH] Fix issue with Kokkos ReaxFF and pair hybrid --- src/KOKKOS/fix_qeq_reax_kokkos.cpp | 11 ++++---- src/KOKKOS/fix_qeq_reax_kokkos.h | 2 +- src/KOKKOS/pair_reaxc_kokkos.cpp | 42 ++++++++++-------------------- src/KOKKOS/pair_reaxc_kokkos.h | 1 - src/USER-REAXC/fix_qeq_reax.cpp | 25 +++++++++++------- src/USER-REAXC/pair_reaxc.h | 2 +- 6 files changed, 37 insertions(+), 46 deletions(-) diff --git a/src/KOKKOS/fix_qeq_reax_kokkos.cpp b/src/KOKKOS/fix_qeq_reax_kokkos.cpp index 10d33ead45..c1b652de82 100644 --- a/src/KOKKOS/fix_qeq_reax_kokkos.cpp +++ b/src/KOKKOS/fix_qeq_reax_kokkos.cpp @@ -195,7 +195,6 @@ void FixQEqReaxKokkos::pre_force(int vflag) type = atomKK->k_type.view(); mask = atomKK->k_mask.view(); nlocal = atomKK->nlocal; - nall = atom->nlocal + atom->nghost; newton_pair = force->newton_pair; k_params.template sync(); @@ -207,6 +206,7 @@ void FixQEqReaxKokkos::pre_force(int vflag) d_neighbors = k_list->d_neighbors; d_ilist = k_list->d_ilist; inum = list->inum; + ignum = inum + list->gnum; copymode = 1; @@ -373,7 +373,6 @@ void FixQEqReaxKokkos::allocate_array() } // init_storage - const int ignum = atom->nlocal + atom->nghost; FixQEqReaxKokkosZeroFunctor zero_functor(this); Kokkos::parallel_for(ignum,zero_functor); @@ -732,7 +731,7 @@ void FixQEqReaxKokkos::cg_solve1() FixQEqReaxKokkosSparse12Functor sparse12_functor(this); Kokkos::parallel_for(inum,sparse12_functor); if (neighflag != FULL) { - Kokkos::parallel_for(Kokkos::RangePolicy(nlocal,nlocal+atom->nghost),*this); + Kokkos::parallel_for(Kokkos::RangePolicy(inum,ignum),*this); if (neighflag == HALF) { FixQEqReaxKokkosSparse13Functor sparse13_functor(this); Kokkos::parallel_for(inum,sparse13_functor); @@ -788,7 +787,7 @@ void FixQEqReaxKokkos::cg_solve1() FixQEqReaxKokkosSparse22Functor sparse22_functor(this); Kokkos::parallel_for(inum,sparse22_functor); if (neighflag != FULL) { - Kokkos::parallel_for(Kokkos::RangePolicy(nlocal,nlocal+atom->nghost),*this); + Kokkos::parallel_for(Kokkos::RangePolicy(inum,ignum),*this); if (need_dup) dup_o.reset_except(d_o); if (neighflag == HALF) { @@ -869,7 +868,7 @@ void FixQEqReaxKokkos::cg_solve2() FixQEqReaxKokkosSparse32Functor sparse32_functor(this); Kokkos::parallel_for(inum,sparse32_functor); if (neighflag != FULL) { - Kokkos::parallel_for(Kokkos::RangePolicy(nlocal,nlocal+atom->nghost),*this); + Kokkos::parallel_for(Kokkos::RangePolicy(inum,ignum),*this); if (need_dup) dup_o.reset_except(d_o); if (neighflag == HALF) { @@ -927,7 +926,7 @@ void FixQEqReaxKokkos::cg_solve2() FixQEqReaxKokkosSparse22Functor sparse22_functor(this); Kokkos::parallel_for(inum,sparse22_functor); if (neighflag != FULL) { - Kokkos::parallel_for(Kokkos::RangePolicy(nlocal,nlocal+atom->nghost),*this); + Kokkos::parallel_for(Kokkos::RangePolicy(inum,ignum),*this); if (need_dup) dup_o.reset_except(d_o); if (neighflag == HALF) { diff --git a/src/KOKKOS/fix_qeq_reax_kokkos.h b/src/KOKKOS/fix_qeq_reax_kokkos.h index 55dec64d33..bd533a0542 100644 --- a/src/KOKKOS/fix_qeq_reax_kokkos.h +++ b/src/KOKKOS/fix_qeq_reax_kokkos.h @@ -151,7 +151,7 @@ class FixQEqReaxKokkos : public FixQEqReax { double memory_usage(); private: - int inum; + int inum,ignum; int allocated_flag; int need_dup; diff --git a/src/KOKKOS/pair_reaxc_kokkos.cpp b/src/KOKKOS/pair_reaxc_kokkos.cpp index 5326e926e5..2fdaddcf8d 100644 --- a/src/KOKKOS/pair_reaxc_kokkos.cpp +++ b/src/KOKKOS/pair_reaxc_kokkos.cpp @@ -199,6 +199,8 @@ void PairReaxCKokkos::setup() hbond_parameters *hbp; for (i = 1; i <= n; i++) { + if (map[i] == -1) continue; + // general k_params_sing.h_view(i).mass = system->reax_param.sbp[map[i]].mass; @@ -229,6 +231,8 @@ void PairReaxCKokkos::setup() k_params_sing.h_view(i).p_hbond = system->reax_param.sbp[map[i]].p_hbond; for (j = 1; j <= n; j++) { + if (map[j] == -1) continue; + twbp = &(system->reax_param.tbp[map[i]][map[j]]); // vdW @@ -270,6 +274,8 @@ void PairReaxCKokkos::setup() k_params_twbp.h_view(i,j).p_ovun1 = twbp->p_ovun1; for (k = 1; k <= n; k++) { + if (map[k] == -1) continue; + // Angular thbh = &(system->reax_param.thbp[map[i]][map[j]][map[k]]); thbp = &(thbh->prm[0]); @@ -290,6 +296,8 @@ void PairReaxCKokkos::setup() k_params_hbp.h_view(i,j,k).r0_hb = hbp->r0_hb; for (m = 1; m <= n; m++) { + if (map[m] == -1) continue; + // Torsion fbh = &(system->reax_param.fbp[map[i]][map[j]][map[k]][map[m]]); fbp = &(fbh->prm[0]); @@ -377,7 +385,9 @@ void PairReaxCKokkos::init_md() d_LR = k_LR.template view(); for (int i = 1; i <= ntypes; ++i) { + if (map[i] == -1) continue; for (int j = i; j <= ntypes; ++j) { + if (map[j] == -1) continue; int n = LR[i][j].n; if (n == 0) continue; k_LR.h_view(i,j).dx = LR[i][j].dx; @@ -539,7 +549,9 @@ void PairReaxCKokkos::Deallocate_Lookup_Tables() ntypes = atom->ntypes; for( i = 0; i <= ntypes; ++i ) { - for( j = i; j <= ntypes; ++j ) + if (map[i] == -1) continue; + for( j = i; j <= ntypes; ++j ) { + if (map[i] == -1) continue; if (LR[i][j].n) { sfree( control->error_ptr, LR[i][j].y, "LR[i,j].y" ); sfree( control->error_ptr, LR[i][j].H, "LR[i,j].H" ); @@ -548,6 +560,7 @@ void PairReaxCKokkos::Deallocate_Lookup_Tables() sfree( control->error_ptr, LR[i][j].ele, "LR[i,j].ele" ); sfree( control->error_ptr, LR[i][j].CEclmb, "LR[i,j].CEclmb" ); } + } sfree( control->error_ptr, LR[i], "LR[i]" ); } sfree( control->error_ptr, LR, "LR" ); @@ -3659,33 +3672,6 @@ void PairReaxCKokkos::v_tally3_atom(EV_FLOAT_REAX &ev, const int &i, } } -/* ---------------------------------------------------------------------- */ - -template -void *PairReaxCKokkos::extract(const char *str, int &dim) -{ - dim = 1; - if (strcmp(str,"chi") == 0 && chi) { - for (int i = 1; i <= atom->ntypes; i++) - if (map[i] >= 0) chi[i] = system->reax_param.sbp[map[i]].chi; - else chi[i] = 0.0; - return (void *) chi; - } - if (strcmp(str,"eta") == 0 && eta) { - for (int i = 1; i <= atom->ntypes; i++) - if (map[i] >= 0) eta[i] = system->reax_param.sbp[map[i]].eta; - else eta[i] = 0.0; - return (void *) eta; - } - if (strcmp(str,"gamma") == 0 && gamma) { - for (int i = 1; i <= atom->ntypes; i++) - if (map[i] >= 0) gamma[i] = system->reax_param.sbp[map[i]].gamma; - else gamma[i] = 0.0; - return (void *) gamma; - } - return NULL; -} - /* ---------------------------------------------------------------------- setup for energy, virial computation see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) diff --git a/src/KOKKOS/pair_reaxc_kokkos.h b/src/KOKKOS/pair_reaxc_kokkos.h index 24a56b7076..60865bfbe2 100644 --- a/src/KOKKOS/pair_reaxc_kokkos.h +++ b/src/KOKKOS/pair_reaxc_kokkos.h @@ -121,7 +121,6 @@ class PairReaxCKokkos : public PairReaxC { void ev_setup(int, int, int alloc = 1); void compute(int, int); - void *extract(const char *, int &); void init_style(); double memory_usage(); void FindBond(int &); diff --git a/src/USER-REAXC/fix_qeq_reax.cpp b/src/USER-REAXC/fix_qeq_reax.cpp index e3b9903cbf..b9cff75d65 100644 --- a/src/USER-REAXC/fix_qeq_reax.cpp +++ b/src/USER-REAXC/fix_qeq_reax.cpp @@ -463,19 +463,26 @@ void FixQEqReax::min_setup_pre_force(int vflag) void FixQEqReax::init_storage() { int NN; + int *ilist; - if (reaxc) + if (reaxc) { NN = reaxc->list->inum + reaxc->list->gnum; - else + ilist = reaxc->list->ilist; + } else { NN = list->inum + list->gnum; + ilist = list->ilist; + } - for (int i = 0; i < NN; i++) { - Hdia_inv[i] = 1. / eta[atom->type[i]]; - b_s[i] = -chi[atom->type[i]]; - b_t[i] = -1.0; - b_prc[i] = 0; - b_prm[i] = 0; - s[i] = t[i] = 0; + for (int ii = 0; ii < NN; ii++) { + int i = ilist[ii]; + if (atom->mask[i] & groupbit) { + Hdia_inv[i] = 1. / eta[atom->type[i]]; + b_s[i] = -chi[atom->type[i]]; + b_t[i] = -1.0; + b_prc[i] = 0; + b_prm[i] = 0; + s[i] = t[i] = 0; + } } } diff --git a/src/USER-REAXC/pair_reaxc.h b/src/USER-REAXC/pair_reaxc.h index a833ddbcf0..df04875027 100644 --- a/src/USER-REAXC/pair_reaxc.h +++ b/src/USER-REAXC/pair_reaxc.h @@ -44,7 +44,7 @@ class PairReaxC : public Pair { void coeff(int, char **); virtual void init_style(); double init_one(int, int); - void *extract(const char *, int &); + virtual void *extract(const char *, int &); int fixbond_flag, fixspecies_flag; int **tmpid; double **tmpbo,**tmpr;