Fix broken charge history in fix_qeq_reax_kokkos

src/KOKKOS/fix_qeq_reax_kokkos.cpp

@@ -72,6 +72,9 @@ template<class DeviceType>
 FixQEqReaxKokkos<DeviceType>::~FixQEqReaxKokkos()
 {
   if (copymode) return;
+
+  memoryKK->destroy_kokkos(k_s_hist,s_hist);
+  memoryKK->destroy_kokkos(k_t_hist,t_hist);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -235,6 +238,8 @@ void FixQEqReaxKokkos<DeviceType>::pre_force(int vflag)
   Kokkos::parallel_scan(inum,computeH_functor);
 
   // init_matvec
+  k_s_hist.template sync<DeviceType>();
+  k_t_hist.template sync<DeviceType>();
   FixQEqReaxKokkosMatVecFunctor<DeviceType> matvec_functor(this);
   Kokkos::parallel_for(inum,matvec_functor);
 
@@ -262,6 +267,8 @@ void FixQEqReaxKokkos<DeviceType>::pre_force(int vflag)
 
   // calculate_Q();
   calculate_q();
+  k_s_hist.template modify<DeviceType>();
+  k_t_hist.template modify<DeviceType>();
 
   copymode = 0;
 
@@ -462,7 +469,7 @@ double FixQEqReaxKokkos<DeviceType>::calculate_H_k(const F_FLOAT &r, const F_FLO
 
 template<class DeviceType>
 KOKKOS_INLINE_FUNCTION
-void FixQEqReaxKokkos<DeviceType>::mat_vec_item(int ii) const
+void FixQEqReaxKokkos<DeviceType>::matvec_item(int ii) const
 {
   const int i = d_ilist[ii];
   const int itype = type(i);
@@ -1183,7 +1190,76 @@ double FixQEqReaxKokkos<DeviceType>::memory_usage()
   return bytes;
 }
 
-/* ---------------------------------------------------------------------- */\
+/* ----------------------------------------------------------------------
+   allocate fictitious charge arrays
+------------------------------------------------------------------------- */
+
+template<class DeviceType>
+void FixQEqReaxKokkos<DeviceType>::grow_arrays(int nmax)
+{
+  k_s_hist.template sync<LMPHostType>(); // force reallocation on host
+  k_t_hist.template sync<LMPHostType>();
+
+  memoryKK->grow_kokkos(k_s_hist,s_hist,nmax,nprev,"qeq:s_hist");
+  memoryKK->grow_kokkos(k_t_hist,t_hist,nmax,nprev,"qeq:t_hist");
+
+  k_s_hist.template modify<LMPHostType>();
+  k_t_hist.template modify<LMPHostType>();
+
+}
+
+/* ----------------------------------------------------------------------
+   copy values within fictitious charge arrays
+------------------------------------------------------------------------- */
+
+template<class DeviceType>
+void FixQEqReaxKokkos<DeviceType>::copy_arrays(int i, int j, int delflag)
+{
+  k_s_hist.template sync<LMPHostType>();
+  k_t_hist.template sync<LMPHostType>();
+
+  for (int m = 0; m < nprev; m++) {
+    s_hist[j][m] = s_hist[i][m];
+    t_hist[j][m] = t_hist[i][m];
+  }
+
+  k_s_hist.template modify<LMPHostType>();
+  k_t_hist.template modify<LMPHostType>();
+
+}
+
+/* ----------------------------------------------------------------------
+   pack values in local atom-based array for exchange with another proc
+------------------------------------------------------------------------- */
+
+template<class DeviceType>
+int FixQEqReaxKokkos<DeviceType>::pack_exchange(int i, double *buf)
+{
+  k_s_hist.template sync<LMPHostType>();
+  k_t_hist.template sync<LMPHostType>();
+
+  for (int m = 0; m < nprev; m++) buf[m] = s_hist[i][m];
+  for (int m = 0; m < nprev; m++) buf[nprev+m] = t_hist[i][m];
+  return nprev*2;
+}
+
+/* ----------------------------------------------------------------------
+   unpack values in local atom-based array from exchange with another proc
+------------------------------------------------------------------------- */
+
+template<class DeviceType>
+int FixQEqReaxKokkos<DeviceType>::unpack_exchange(int nlocal, double *buf)
+{
+  for (int m = 0; m < nprev; m++) s_hist[nlocal][m] = buf[m];
+  for (int m = 0; m < nprev; m++) t_hist[nlocal][m] = buf[nprev+m];
+
+  k_s_hist.template modify<LMPHostType>();
+  k_t_hist.template modify<LMPHostType>();
+
+  return nprev*2;
+}
+
+/* ---------------------------------------------------------------------- */
 
 namespace LAMMPS_NS {
 template class FixQEqReaxKokkos<LMPDeviceType>;

src/KOKKOS/fix_qeq_reax_kokkos.h

@@ -57,7 +57,7 @@ class FixQEqReaxKokkos : public FixQEqReax {
   void compute_h_item(int, int &, const bool &) const;
 
   KOKKOS_INLINE_FUNCTION
-  void mat_vec_item(int) const;
+  void matvec_item(int) const;
 
   KOKKOS_INLINE_FUNCTION
   void sparse12_item(int) const;
@@ -145,7 +145,7 @@ class FixQEqReaxKokkos : public FixQEqReax {
   void unpack_reverse_comm(int, int *, double *);
   double memory_usage();
 
- protected:
+ private:
   int inum;
   int allocated_flag;
 
@@ -210,6 +210,10 @@ class FixQEqReaxKokkos : public FixQEqReax {
   typename AT::t_int_2d d_sendlist;
   typename AT::t_xfloat_1d_um v_buf;
 
+  void grow_arrays(int);
+  void copy_arrays(int, int, int);
+  int pack_exchange(int, double *);
+  int unpack_exchange(int, double *);
 };
 
 template <class DeviceType>
@@ -235,7 +239,7 @@ struct FixQEqReaxKokkosMatVecFunctor  {
   };
   KOKKOS_INLINE_FUNCTION
   void operator()(const int ii) const {
-    c.mat_vec_item(ii);
+    c.matvec_item(ii);
   }
 };
 

src/USER-REAXC/fix_qeq_reax.h

@@ -122,15 +122,15 @@ class FixQEqReax : public Fix {
   //int GMRES(double*,double*);
   virtual void sparse_matvec(sparse_matrix*,double*,double*);
 
-  int pack_forward_comm(int, int *, double *, int, int *);
-  void unpack_forward_comm(int, int, double *);
-  int pack_reverse_comm(int, int, double *);
-  void unpack_reverse_comm(int, int *, double *);
-  double memory_usage();
-  void grow_arrays(int);
-  void copy_arrays(int, int, int);
-  int pack_exchange(int, double *);
-  int unpack_exchange(int, double *);
+  virtual int pack_forward_comm(int, int *, double *, int, int *);
+  virtual void unpack_forward_comm(int, int, double *);
+  virtual int pack_reverse_comm(int, int, double *);
+  virtual void unpack_reverse_comm(int, int *, double *);
+  virtual double memory_usage();
+  virtual void grow_arrays(int);
+  virtual void copy_arrays(int, int, int);
+  virtual int pack_exchange(int, double *);
+  virtual int unpack_exchange(int, double *);
 
   virtual double parallel_norm( double*, int );
   virtual double parallel_dot( double*, double*, int );