add per-atom array corerep_factor to all four pair styles: pace, pace/extrapolation, pace/kk, pace/extrapolation/kk

upd PACELIB_URL to v.2023.11.25 (+md5sum) in ML-PACE.cmake and Install.py

cmake/Modules/Packages/ML-PACE.cmake

@@ -1,6 +1,6 @@
-set(PACELIB_URL "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2023.11.22.tar.gz" CACHE STRING "URL for PACE evaluator library sources")
+set(PACELIB_URL "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2023.11.25.tar.gz" CACHE STRING "URL for PACE evaluator library sources")
 
-set(PACELIB_MD5 "c8e811f96d761ef8863f5b88a3fd36f4" CACHE STRING "MD5 checksum of PACE evaluator library tarball")
+set(PACELIB_MD5 "fda61c88a6c6a335d5fc10a86a5aa710" CACHE STRING "MD5 checksum of PACE evaluator library tarball")
 mark_as_advanced(PACELIB_URL)
 mark_as_advanced(PACELIB_MD5)
 GetFallbackURL(PACELIB_URL PACELIB_FALLBACK)

lib/pace/Install.py

@@ -18,11 +18,11 @@ from install_helpers import fullpath, geturl, checkmd5sum, getfallback
 # settings
 
 thisdir = fullpath('.')
-version ='v.2023.11.22'
+version ='v.2023.11.25'
 
 # known checksums for different PACE versions. used to validate the download.
 checksums = { \
-    'v.2023.11.22': 'c8e811f96d761ef8863f5b88a3fd36f4'
+    'v.2023.11.25': 'fda61c88a6c6a335d5fc10a86a5aa710'
 }
 
 parser = ArgumentParser(prog='Install.py', description="LAMMPS library build wrapper script")

src/KOKKOS/pair_pace_extrapolation_kokkos.cpp

@@ -126,6 +126,7 @@ void PairPACEExtrapolationKokkos<DeviceType>::grow(int natom, int maxneigh)
     MemKK::realloc_kokkos(dF_dfcut, "pace:dF_dfcut", natom);
     MemKK::realloc_kokkos(d_d_min, "pace:r_min_pair", natom);
     MemKK::realloc_kokkos(d_jj_min, "pace:j_min_pair", natom);
+    MemKK::realloc_kokkos(d_corerep, "pace:corerep", natom); // per-atom corerep
 
     MemKK::realloc_kokkos(dB_flatten, "pace:dB_flatten", natom, idx_ms_combs_max, basis_set->rankmax);
 
@@ -594,13 +595,20 @@ void PairPACEExtrapolationKokkos<DeviceType>::compute(int eflag_in, int vflag_in
     d_vatom = k_vatom.view<DeviceType>();
   }
 
-  if (gamma_flag && atom->nlocal > nmax) {
+  if (flag_compute_extrapolation_grade && atom->nlocal > nmax) {
         memory->destroy(extrapolation_grade_gamma);
         nmax = atom->nlocal;
         memory->create(extrapolation_grade_gamma, nmax, "pace/atom:gamma");
         //zeroify array
         memset(extrapolation_grade_gamma, 0, nmax * sizeof(*extrapolation_grade_gamma));
   }
+  if (flag_corerep_factor && atom->nlocal > nmax_corerep) {
+    memory->destroy(corerep_factor);
+    nmax_corerep = atom->nlocal;
+    memory->create(corerep_factor, nmax_corerep, "pace/atom:corerep");
+    //zeroify array
+    memset(corerep_factor, 0, nmax_corerep * sizeof(*corerep_factor));
+  }
 
   copymode = 1;
   if (!force->newton_pair)
@@ -658,6 +666,7 @@ void PairPACEExtrapolationKokkos<DeviceType>::compute(int eflag_in, int vflag_in
 
     Kokkos::deep_copy(projections, 0.0);
     Kokkos::deep_copy(d_gamma, 0.0);
+    Kokkos::deep_copy(d_corerep, 0.0);
 
     EV_FLOAT ev_tmp;
 
@@ -721,7 +730,7 @@ void PairPACEExtrapolationKokkos<DeviceType>::compute(int eflag_in, int vflag_in
     }
 
     //ComputeGamma
-    if (gamma_flag) {
+    if (flag_compute_extrapolation_grade) {
       typename Kokkos::RangePolicy<DeviceType,TagPairPACEComputeGamma> policy_gamma(0,chunk_size);
       Kokkos::parallel_for("ComputeGamma",policy_gamma,*this);
     }
@@ -763,12 +772,17 @@ void PairPACEExtrapolationKokkos<DeviceType>::compute(int eflag_in, int vflag_in
     }
     ev += ev_tmp;
 
-    //if gamma_flag - copy current d_gamma to extrapolation_grade_gamma
+    //if flag_compute_extrapolation_grade - copy current d_gamma to extrapolation_grade_gamma
-    if (gamma_flag){
+    if (flag_compute_extrapolation_grade){
         h_gamma = Kokkos::create_mirror_view(d_gamma);
         Kokkos::deep_copy(h_gamma, d_gamma);
         memcpy(extrapolation_grade_gamma+chunk_offset, (void *) h_gamma.data(), sizeof(double)*chunk_size);
     }
+    if (flag_corerep_factor) {
+        h_corerep = Kokkos::create_mirror_view(d_corerep);
+        Kokkos::deep_copy(h_corerep,d_corerep);
+        memcpy(corerep_factor+chunk_offset, (void *) h_corerep.data(), sizeof(double)*chunk_size);
+    }
 
     chunk_offset += chunk_size;
   } // end while
@@ -1050,7 +1064,7 @@ void PairPACEExtrapolationKokkos<DeviceType>::operator() (TagPairPACEComputeRho,
 
 
     //gamma_i
-    if (gamma_flag)
+    if (flag_compute_extrapolation_grade)
         Kokkos::atomic_add(&projections(ii, func_ind),  d_gen_cgs(mu_i, idx_ms_comb) * A_cur);
 
   } else { // rank > 1
@@ -1089,7 +1103,7 @@ void PairPACEExtrapolationKokkos<DeviceType>::operator() (TagPairPACEComputeRho,
       Kokkos::atomic_add(&rhos(ii, p), B.real_part_product(d_coeffs(mu_i, func_ind, p) * d_gen_cgs(mu_i, idx_ms_comb)));
     }
     //gamma_i
-    if (gamma_flag)
+    if (flag_compute_extrapolation_grade)
       Kokkos::atomic_add(&projections(ii, func_ind),  B.real_part_product(d_gen_cgs(mu_i, idx_ms_comb)));
   }
 }
@@ -1144,6 +1158,9 @@ void PairPACEExtrapolationKokkos<DeviceType>::operator() (TagPairPACEComputeFS,
     evdwl_cut += d_E0vals(mu_i);
     e_atom(ii) = evdwl_cut;
   }
+
+  if (flag_corerep_factor)
+    d_corerep(ii) = 1-fcut;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1858,6 +1875,7 @@ double PairPACEExtrapolationKokkos<DeviceType>::memory_usage()
   bytes += MemKK::memory_usage(rho_core);
   bytes += MemKK::memory_usage(dF_drho_core);
   bytes += MemKK::memory_usage(dF_dfcut);
+  bytes += MemKK::memory_usage(d_corerep);
   bytes += MemKK::memory_usage(dB_flatten);
   bytes += MemKK::memory_usage(fr);
   bytes += MemKK::memory_usage(dfr);
@@ -1925,9 +1943,10 @@ double PairPACEExtrapolationKokkos<DeviceType>::memory_usage()
 template<class DeviceType>
 void *PairPACEExtrapolationKokkos<DeviceType>::extract(const char *str, int &dim)
 {
-  //check if str=="gamma_flag" then compute extrapolation grades on this iteration
   dim = 0;
-  if (strcmp(str, "gamma_flag") == 0) return (void *) &gamma_flag;
+  //check if str=="flag_compute_extrapolation_grade" then compute extrapolation grades on this iteration
+  if (strcmp(str, "gamma_flag") == 0) return (void *) &flag_compute_extrapolation_grade;
+  if (strcmp(str, "corerep_flag") == 0) return (void *) &flag_corerep_factor;
 
   dim = 2;
   if (strcmp(str, "scale") == 0) return (void *) scale;
@@ -1950,6 +1969,10 @@ void *PairPACEExtrapolationKokkos<DeviceType>::extract_peratom(const char *str,
     ncol = 0;
     return (void *) extrapolation_grade_gamma;
   }
+  if (strcmp(str, "corerep") == 0) {
+    ncol = 0;
+    return (void *) corerep_factor;
+  }
 
   return nullptr;
 }

src/KOKKOS/pair_pace_extrapolation_kokkos.h

@@ -106,7 +106,7 @@ class PairPACEExtrapolationKokkos : public PairPACEExtrapolation {
  protected:
   int inum, maxneigh, chunk_size, chunk_offset, idx_ms_combs_max, total_num_functions_max;
   int host_flag;
-  int gamma_flag;
+  //int gamma_flag;
 
   int eflag, vflag;
 
@@ -235,13 +235,16 @@ class PairPACEExtrapolationKokkos : public PairPACEExtrapolation {
   t_ace_2d rhos;
   t_ace_2d dF_drho;
 
+  t_ace_3c dB_flatten;
+
   // hard-core repulsion
   t_ace_1d rho_core;
-  t_ace_3c dB_flatten;
   t_ace_2d cr;
   t_ace_2d dcr;
   t_ace_1d dF_drho_core;
   t_ace_1d dF_dfcut;
+  t_ace_1d d_corerep;
+  th_ace_1d h_corerep;
 
   // radial functions
   t_ace_4d fr;

src/KOKKOS/pair_pace_kokkos.cpp

@@ -124,6 +124,8 @@ void PairPACEKokkos<DeviceType>::grow(int natom, int maxneigh)
     MemKK::realloc_kokkos(dF_dfcut, "pace:dF_dfcut", natom);
     MemKK::realloc_kokkos(d_d_min, "pace:r_min_pair", natom);
     MemKK::realloc_kokkos(d_jj_min, "pace:j_min_pair", natom);
+    MemKK::realloc_kokkos(d_corerep, "pace:corerep", natom); // per-atom corerep
+
     MemKK::realloc_kokkos(dB_flatten, "pace:dB_flatten", natom, idx_rho_max, basis_set->rankmax);
   }
 
@@ -555,6 +557,13 @@ void PairPACEKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
     memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom");
     d_vatom = k_vatom.view<DeviceType>();
   }
+  if (flag_corerep_factor && atom->nlocal > nmax_corerep) {
+    memory->destroy(corerep_factor);
+    nmax_corerep = atom->nlocal;
+    memory->create(corerep_factor, nmax_corerep, "pace/atom:corerep");
+    //zeroify array
+    memset(corerep_factor, 0, nmax_corerep * sizeof(*corerep_factor));
+  }
 
   copymode = 1;
   if (!force->newton_pair)
@@ -610,6 +619,7 @@ void PairPACEKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
     Kokkos::deep_copy(rho_core, 0.0);
     Kokkos::deep_copy(d_d_min, PairPACE::aceimpl->basis_set->cutoffmax);
     Kokkos::deep_copy(d_jj_min, -1);
+    Kokkos::deep_copy(d_corerep, 0.0);
 
     EV_FLOAT ev_tmp;
 
@@ -709,6 +719,13 @@ void PairPACEKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
       }
     }
     ev += ev_tmp;
+
+    if (flag_corerep_factor) {
+      h_corerep = Kokkos::create_mirror_view(d_corerep);
+      Kokkos::deep_copy(h_corerep,d_corerep);
+      memcpy(corerep_factor+chunk_offset, (void *) h_corerep.data(), sizeof(double)*chunk_size);
+    }
+
     chunk_offset += chunk_size;
 
   } // end while
@@ -1073,6 +1090,9 @@ void PairPACEKokkos<DeviceType>::operator() (TagPairPACEComputeFS, const int& ii
       evdwl_cut += d_E0vals(mu_i);
       e_atom(ii) = evdwl_cut;
   }
+
+  if (flag_corerep_factor)
+    d_corerep(ii) = 1-fcut;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1759,6 +1779,7 @@ double PairPACEKokkos<DeviceType>::memory_usage()
   bytes += MemKK::memory_usage(rho_core);
   bytes += MemKK::memory_usage(dF_drho_core);
   bytes += MemKK::memory_usage(dF_dfcut);
+  bytes += MemKK::memory_usage(d_corerep);
   bytes += MemKK::memory_usage(dB_flatten);
   bytes += MemKK::memory_usage(fr);
   bytes += MemKK::memory_usage(dfr);
@@ -1814,6 +1835,42 @@ double PairPACEKokkos<DeviceType>::memory_usage()
   return bytes;
 }
 
+
+/* ----------------------------------------------------------------------
+    extract method for extracting value of scale variable
+ ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+void *PairPACEKokkos<DeviceType>::extract(const char *str, int &dim)
+{
+  dim = 0;
+  if (strcmp(str, "corerep_flag") == 0) return (void *) &flag_corerep_factor;
+
+  dim = 2;
+  if (strcmp(str, "scale") == 0) return (void *) scale;
+  return nullptr;
+}
+
+/* ----------------------------------------------------------------------
+   peratom requests from FixPair
+   return ptr to requested data
+   also return ncol = # of quantites per atom
+     0 = per-atom vector
+     1 or more = # of columns in per-atom array
+   return NULL if str is not recognized
+---------------------------------------------------------------------- */
+
+template<class DeviceType>
+void *PairPACEKokkos<DeviceType>::extract_peratom(const char *str, int &ncol)
+{
+  if (strcmp(str, "corerep") == 0) {
+    ncol = 0;
+    return (void *) corerep_factor;
+  }
+
+  return nullptr;
+}
+
 /* ---------------------------------------------------------------------- */
 
 namespace LAMMPS_NS {

src/KOKKOS/pair_pace_kokkos.h

@@ -95,6 +95,9 @@ class PairPACEKokkos : public PairPACE {
   KOKKOS_INLINE_FUNCTION
   void operator() (TagPairPACEComputeForce<NEIGHFLAG,EVFLAG>,const int& ii, EV_FLOAT&) const;
 
+  void *extract(const char *str, int &dim) override;
+  void *extract_peratom(const char *str, int &ncol) override;
+
  protected:
   int inum, maxneigh, chunk_size, chunk_offset, idx_rho_max;
   int host_flag;
@@ -210,6 +213,8 @@ class PairPACEKokkos : public PairPACE {
   typedef Kokkos::View<complex****, DeviceType> t_ace_4c;
   typedef Kokkos::View<complex***[3], DeviceType> t_ace_4c3;
 
+  typedef typename Kokkos::View<double*, DeviceType>::HostMirror th_ace_1d;
+
   t_ace_3d A_rank1;
   t_ace_4c A;
 
@@ -223,13 +228,16 @@ class PairPACEKokkos : public PairPACE {
   t_ace_2d rhos;
   t_ace_2d dF_drho;
 
+  t_ace_3c dB_flatten;
+
   // hard-core repulsion
   t_ace_1d rho_core;
-  t_ace_3c dB_flatten;
   t_ace_2d cr;
   t_ace_2d dcr;
   t_ace_1d dF_drho_core;
   t_ace_1d dF_dfcut;
+  t_ace_1d d_corerep;
+  th_ace_1d h_corerep;
 
   // radial functions
   t_ace_4d fr;

src/ML-PACE/pair_pace.cpp

@@ -88,6 +88,10 @@ PairPACE::PairPACE(LAMMPS *lmp) : Pair(lmp)
   one_coeff = 1;
   manybody_flag = 1;
 
+  nmax_corerep = 0;
+  flag_corerep_factor = 0;
+  corerep_factor = nullptr;
+
   aceimpl = new ACEImpl;
   recursive = false;
 
@@ -110,6 +114,7 @@ PairPACE::~PairPACE()
     memory->destroy(setflag);
     memory->destroy(cutsq);
     memory->destroy(scale);
+    memory->destroy(corerep_factor);
   }
 }
 
@@ -144,6 +149,14 @@ void PairPACE::compute(int eflag, int vflag)
   // the pointer to the list of neighbors of "i"
   firstneigh = list->firstneigh;
 
+  if (flag_corerep_factor && atom->nlocal > nmax_corerep) {
+    memory->destroy(corerep_factor);
+    nmax_corerep = atom->nlocal;
+    memory->create(corerep_factor, nmax_corerep, "pace/atom:corerep_factor");
+    //zeroify array
+    memset(corerep_factor, 0, nmax_corerep * sizeof(*corerep_factor));
+  }
+
   //determine the maximum number of neighbours
   int max_jnum = 0;
   int nei = 0;
@@ -182,6 +195,9 @@ void PairPACE::compute(int eflag, int vflag)
       error->one(FLERR, e.what());
     }
 
+    if (flag_corerep_factor)
+      corerep_factor[i] = 1 - aceimpl->ace->ace_fcut;
+
     // 'compute_atom' will update the `aceimpl->ace->e_atom` and `aceimpl->ace->neighbours_forces(jj, alpha)` arrays
 
     for (jj = 0; jj < jnum; jj++) {
@@ -382,7 +398,29 @@ double PairPACE::init_one(int i, int j)
  ---------------------------------------------------------------------- */
 void *PairPACE::extract(const char *str, int &dim)
 {
+  dim = 0;
+  //check if str=="corerep_flag" then compute extrapolation grades on this iteration
+  if (strcmp(str, "corerep_flag") == 0) return (void *) &flag_corerep_factor;
+
   dim = 2;
   if (strcmp(str, "scale") == 0) return (void *) scale;
   return nullptr;
 }
+
+/* ----------------------------------------------------------------------
+   peratom requests from FixPair
+   return ptr to requested data
+   also return ncol = # of quantites per atom
+     0 = per-atom vector
+     1 or more = # of columns in per-atom array
+   return NULL if str is not recognized
+---------------------------------------------------------------------- */
+void *PairPACE::extract_peratom(const char *str, int &ncol)
+{
+  if (strcmp(str, "corerep") == 0) {
+    ncol = 0;
+    return (void *) corerep_factor;
+  }
+
+  return nullptr;
+}

src/ML-PACE/pair_pace.h

@@ -48,11 +48,15 @@ class PairPACE : public Pair {
   double init_one(int, int) override;
 
   void *extract(const char *, int &) override;
+  void *extract_peratom(const char *, int &) override;
 
  protected:
   struct ACEImpl *aceimpl;
+  int nmax_corerep = 0;
 
   virtual void allocate();
+  double *corerep_factor;                    //per-atom core-rep factor (= 1 - fcut)
+  int flag_corerep_factor;
 
   double **scale;
   bool recursive;    // "recursive" option for ACERecursiveEvaluator

src/ML-PACE/pair_pace_extrapolation.cpp

@@ -93,11 +93,14 @@ PairPACEExtrapolation::PairPACEExtrapolation(LAMMPS *lmp) : Pair(lmp)
   manybody_flag = 1;
 
   nmax = 0;
+  nmax_corerep = 0;
 
   aceimpl = new ACEALImpl;
   scale = nullptr;
   flag_compute_extrapolation_grade = 0;
   extrapolation_grade_gamma = nullptr;
+  flag_corerep_factor = 0;
+  corerep_factor = nullptr;
 
   chunksize = 4096;
 }
@@ -118,6 +121,7 @@ PairPACEExtrapolation::~PairPACEExtrapolation()
     memory->destroy(scale);
     memory->destroy(map);
     memory->destroy(extrapolation_grade_gamma);
+    memory->destroy(corerep_factor);
   }
 }
 
@@ -166,6 +170,13 @@ void PairPACEExtrapolation::compute(int eflag, int vflag)
     //zeroify array
     memset(extrapolation_grade_gamma, 0, nmax * sizeof(*extrapolation_grade_gamma));
   }
+  if (flag_corerep_factor && atom->nlocal > nmax_corerep) {
+    memory->destroy(corerep_factor);
+    nmax_corerep = atom->nlocal;
+    memory->create(corerep_factor, nmax_corerep, "pace/atom:corerep_factor");
+    //zeroify array
+    memset(corerep_factor, 0, nmax_corerep * sizeof(*corerep_factor));
+  }
 
   //determine the maximum number of neighbours
   int max_jnum = 0;
@@ -216,6 +227,11 @@ void PairPACEExtrapolation::compute(int eflag, int vflag)
     if (flag_compute_extrapolation_grade)
       extrapolation_grade_gamma[i] = aceimpl->ace->max_gamma_grade;
 
+    if (flag_corerep_factor) {
+      corerep_factor[i] = 1 - (flag_compute_extrapolation_grade ? aceimpl->ace->ace_fcut
+                              : aceimpl->rec_ace->ace_fcut);
+    }
+
     Array2D<DOUBLE_TYPE> &neighbours_forces =
         (flag_compute_extrapolation_grade ? aceimpl->ace->neighbours_forces
                                           : aceimpl->rec_ace->neighbours_forces);
@@ -437,9 +453,11 @@ double PairPACEExtrapolation::init_one(int i, int j)
  ---------------------------------------------------------------------- */
 void *PairPACEExtrapolation::extract(const char *str, int &dim)
 {
-  //check if str=="gamma_flag" then compute extrapolation grades on this iteration
   dim = 0;
+  //check if str=="gamma_flag" then compute extrapolation grades on this iteration
   if (strcmp(str, "gamma_flag") == 0) return (void *) &flag_compute_extrapolation_grade;
+  //check if str=="corerep_flag" then compute extrapolation grades on this iteration
+  if (strcmp(str, "corerep_flag") == 0) return (void *) &flag_corerep_factor;
 
   dim = 2;
   if (strcmp(str, "scale") == 0) return (void *) scale;
@@ -461,5 +479,10 @@ void *PairPACEExtrapolation::extract_peratom(const char *str, int &ncol)
     return (void *) extrapolation_grade_gamma;
   }
 
+  if (strcmp(str, "corerep") == 0) {
+    ncol = 0;
+    return (void *) corerep_factor;
+  }
+
   return nullptr;
 }

src/ML-PACE/pair_pace_extrapolation.h

@@ -47,13 +47,15 @@ class PairPACEExtrapolation : public Pair {
 
  protected:
   struct ACEALImpl *aceimpl;
-  int nmax;
+  int nmax = 0, nmax_corerep = 0;
 
   virtual void allocate();
   std::vector<std::string> element_names;    // list of elements (used by dump pace/extrapolation)
-  double *extrapolation_grade_gamma;         //per-atom gamma value
+  double *extrapolation_grade_gamma = nullptr;         //per-atom gamma value
+  double *corerep_factor = nullptr;                    //per-atom core-rep factor (= 1 - fcut)
 
-  int flag_compute_extrapolation_grade;
+  int flag_compute_extrapolation_grade = 0;
+  int flag_corerep_factor = 0;
 
   double **scale;