Merge pull request #3585 from yury-lysogorskiy/feature/pace-extrapolation-kokkos

Feature/pace extrapolation kokkos

cmake/Modules/Packages/ML-PACE.cmake

@@ -1,6 +1,6 @@
-set(PACELIB_URL "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2022.10.15.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.01.3.tar.gz" CACHE STRING "URL for PACE evaluator library sources")
 
-set(PACELIB_MD5 "848ad6a6cc79fa82745927001fb1c9b5" CACHE STRING "MD5 checksum of PACE evaluator library tarball")
+set(PACELIB_MD5 "f418d32b60e531063ac4285bf702b468" CACHE STRING "MD5 checksum of PACE evaluator library tarball")
 mark_as_advanced(PACELIB_URL)
 mark_as_advanced(PACELIB_MD5)
 

doc/src/fix_pair.rst

@@ -1,7 +1,7 @@
 .. index:: fix pair
 
 fix pair command
-=======================
+================
 
 Syntax
 """"""
@@ -47,7 +47,12 @@ These are example use cases:
 The *N* argument determines how often the fix is invoked.
 
 The *pstyle* argument is the name of the pair style.  It can be a
-sub-style used in a :doc:`pair_style hybrid <pair_hybrid>` command.
+sub-style used in a :doc:`pair_style hybrid <pair_hybrid>` command.  If
+there are multiple sub-styles using the same pair style, then *pstyle*
+should be specified as "style:N", where *N* is the number of the
+instance of the pair style you wish monitor (e.g., the first or second).
+For example, *pstyle* could be specified as "pace/extrapolation" or
+"amoeba" or "eam:1" or "eam:2".
 
 One or more *name/flag* pairs of arguments follow.  Each *name* is a
 per-atom quantity which the pair style must recognize as an extraction

doc/src/pair_pace.rst

@@ -1,6 +1,7 @@
 .. index:: pair_style pace
 .. index:: pair_style pace/kk
 .. index:: pair_style pace/extrapolation
+.. index:: pair_style pace/extrapolation/kk
 
 pair_style pace command
 =======================
@@ -127,6 +128,9 @@ but not more often than every 20 steps.
 
 On all other steps `pair_style pace recursive` will be used.
 
+When using the pair style *pace/extrapolation* with the KOKKOS package on GPUs
+product B-basis evaluator is always used and only *linear* ASI is supported.
+
 ----------
 
 See the :doc:`pair_coeff <pair_coeff>` page for alternate ways
@@ -186,4 +190,4 @@ recursive, chunksize = 4096,
 
 .. _Lysogorskiy2022:
 
-**(Lysogorskiy2022)** Lysogorskiy, Bochkarev, Mrovec, Drautz, TBS (2022).
+**(Lysogorskiy2022)** Lysogorskiy, Bochkarev, Mrovec, Drautz, arXiv:2212.08716 (2022).

lib/pace/Install.py

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

src/KOKKOS/Install.sh

@@ -320,6 +320,8 @@ action pair_multi_lucy_rx_kokkos.cpp pair_multi_lucy_rx.cpp
 action pair_multi_lucy_rx_kokkos.h pair_multi_lucy_rx.h
 action pair_pace_kokkos.cpp pair_pace.cpp
 action pair_pace_kokkos.h pair_pace.h
+action pair_pace_extrapolation_kokkos.cpp pair_pace_extrapolation.cpp
+action pair_pace_extrapolation_kokkos.h pair_pace_extrapolation.h
 action pair_reaxff_kokkos.cpp pair_reaxff.cpp
 action pair_reaxff_kokkos.h pair_reaxff.h
 action pair_snap_kokkos.cpp pair_snap.cpp

src/KOKKOS/pair_pace_extrapolation_kokkos.h

@@ -0,0 +1,339 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+// clang-format off
+PairStyle(pace/extrapolation/kk,PairPACEExtrapolationKokkos<LMPDeviceType>);
+PairStyle(pace/extrapolation/kk/device,PairPACEExtrapolationKokkos<LMPDeviceType>);
+PairStyle(pace/extrapolation/kk/host,PairPACEExtrapolationKokkos<LMPHostType>);
+// clang-format on
+#else
+
+// clang-format off
+#ifndef LMP_PAIR_PACE_EXTRAPOLATION_KOKKOS_H
+#define LMP_PAIR_PACE_EXTRAPOLATION_KOKKOS_H
+
+#include "pair_pace_extrapolation.h"
+#include "kokkos_type.h"
+#include "pair_kokkos.h"
+
+class SplineInterpolator;
+
+namespace LAMMPS_NS {
+
+template<class DeviceType>
+class PairPACEExtrapolationKokkos : public PairPACEExtrapolation {
+ public:
+  struct TagPairPACEComputeNeigh{};
+  struct TagPairPACEComputeRadial{};
+  struct TagPairPACEComputeYlm{};
+  struct TagPairPACEComputeAi{};
+  struct TagPairPACEConjugateAi{};
+  struct TagPairPACEComputeRho{};
+  struct TagPairPACEComputeFS{};
+  struct TagPairPACEComputeGamma{};
+  struct TagPairPACEComputeWeights{};
+  struct TagPairPACEComputeDerivative{};
+
+  template<int NEIGHFLAG, int EVFLAG>
+  struct TagPairPACEComputeForce{};
+
+  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
+  typedef EV_FLOAT value_type;
+  using complex = SNAComplex<double>;
+
+  PairPACEExtrapolationKokkos(class LAMMPS *);
+  ~PairPACEExtrapolationKokkos() override;
+
+  void compute(int, int) override;
+  void coeff(int, char **) override;
+  void init_style() override;
+  double init_one(int, int) override;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator() (TagPairPACEComputeNeigh,const typename Kokkos::TeamPolicy<DeviceType, TagPairPACEComputeNeigh>::member_type& team) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator() (TagPairPACEComputeRadial,const typename Kokkos::TeamPolicy<DeviceType, TagPairPACEComputeRadial>::member_type& team) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator() (TagPairPACEComputeYlm,const typename Kokkos::TeamPolicy<DeviceType, TagPairPACEComputeYlm>::member_type& team) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator() (TagPairPACEComputeAi,const typename Kokkos::TeamPolicy<DeviceType, TagPairPACEComputeAi>::member_type& team) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator() (TagPairPACEConjugateAi,const int& ii) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator() (TagPairPACEComputeRho,const int& iter) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator() (TagPairPACEComputeFS,const int& ii) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator() (TagPairPACEComputeGamma, const int& ii) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator() (TagPairPACEComputeWeights,const int& iter) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator() (TagPairPACEComputeDerivative,const typename Kokkos::TeamPolicy<DeviceType, TagPairPACEComputeDerivative>::member_type& team) const;
+
+  template<int NEIGHFLAG, int EVFLAG>
+  KOKKOS_INLINE_FUNCTION
+  void operator() (TagPairPACEComputeForce<NEIGHFLAG,EVFLAG>,const int& ii) const;
+
+  template<int NEIGHFLAG, int EVFLAG>
+  KOKKOS_INLINE_FUNCTION
+  void operator() (TagPairPACEComputeForce<NEIGHFLAG,EVFLAG>,const int& ii, EV_FLOAT&) const;
+
+
+  void *extract(const char *str, int &dim);
+  void *extract_peratom(const char *str, int &ncol);
+
+ protected:
+  int inum, maxneigh, chunk_size, chunk_offset, idx_ms_combs_max, total_num_functions_max;
+  int host_flag;
+  int gamma_flag;
+
+  int eflag, vflag;
+
+  int neighflag, max_ndensity;
+  int nelements, lmax, nradmax, nradbase;
+
+  typename AT::t_neighbors_2d d_neighbors;
+  typename AT::t_int_1d_randomread d_ilist;
+  typename AT::t_int_1d_randomread d_numneigh;
+
+  DAT::tdual_efloat_1d k_eatom;
+  DAT::tdual_virial_array k_vatom;
+  typename AT::t_efloat_1d d_eatom;
+  typename AT::t_virial_array d_vatom;
+
+  typename AT::t_x_array_randomread x;
+  typename AT::t_f_array f;
+  typename AT::t_int_1d_randomread type;
+
+  typedef Kokkos::DualView<F_FLOAT**, DeviceType> tdual_fparams;
+  tdual_fparams k_cutsq, k_scale;
+  typedef Kokkos::View<F_FLOAT**, DeviceType> t_fparams;
+  t_fparams d_cutsq, d_scale;
+
+  typename AT::t_int_1d d_map;
+
+  int need_dup;
+
+  using KKDeviceType = typename KKDevice<DeviceType>::value;
+
+  template<typename DataType, typename Layout>
+  using DupScatterView = KKScatterView<DataType, Layout, KKDeviceType, KKScatterSum, KKScatterDuplicated>;
+
+  template<typename DataType, typename Layout>
+  using NonDupScatterView = KKScatterView<DataType, Layout, KKDeviceType, KKScatterSum, KKScatterNonDuplicated>;
+
+  DupScatterView<F_FLOAT*[3], typename DAT::t_f_array::array_layout> dup_f;
+  DupScatterView<F_FLOAT*[6], typename DAT::t_virial_array::array_layout> dup_vatom;
+
+  NonDupScatterView<F_FLOAT*[3], typename DAT::t_f_array::array_layout> ndup_f;
+  NonDupScatterView<F_FLOAT*[6], typename DAT::t_virial_array::array_layout> ndup_vatom;
+
+  friend void pair_virial_fdotr_compute<PairPACEExtrapolationKokkos>(PairPACEExtrapolationKokkos*);
+
+  void grow(int, int);
+  void copy_pertype();
+  void copy_splines();
+  void copy_tilde();
+  void allocate() override;
+  void precompute_harmonics();
+  double memory_usage() override;
+
+  template<int NEIGHFLAG>
+  KOKKOS_INLINE_FUNCTION
+  void v_tally_xyz(EV_FLOAT &ev, const int &i, const int &j,
+      const F_FLOAT &fx, const F_FLOAT &fy, const F_FLOAT &fz,
+      const F_FLOAT &delx, const F_FLOAT &dely, const F_FLOAT &delz) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void compute_barplm(int, int, double, int) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void compute_ylm(int, int, double, double, double, int) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void cutoff_func_poly(const double, const double, const double, double &, double &) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void Fexp(const double, const double, double &, double &) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void FexpShiftedScaled(const double, const double, double &, double &) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void inner_cutoff(const double, const double, const double, double &, double &) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void FS_values_and_derivatives(const int, double&, const int) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void evaluate_splines(const int, const int, double, int, int, int, int) const;
+
+  template<class TagStyle>
+  void check_team_size_for(int, int&, int);
+
+  template<class TagStyle>
+  void check_team_size_reduce(int, int&, int);
+
+  // Utility routine which wraps computing per-team scratch size requirements for
+  // ComputeNeigh, ComputeUi, and ComputeFusedDeidrj
+  template <typename scratch_type>
+  int scratch_size_helper(int values_per_team);
+
+  typedef Kokkos::View<int*, DeviceType> t_ace_1i;
+  typedef Kokkos::View<int**, DeviceType> t_ace_2i;
+  typedef Kokkos::View<int***, DeviceType> t_ace_3i;
+  typedef Kokkos::View<int****, DeviceType> t_ace_4i;
+  typedef Kokkos::View<double*, DeviceType> t_ace_1d;
+  typedef Kokkos::View<double**, DeviceType> t_ace_2d;
+  typedef Kokkos::View<double*[3], DeviceType> t_ace_2d3;
+  typedef Kokkos::View<double***, DeviceType> t_ace_3d;
+  typedef Kokkos::View<const double***, DeviceType> tc_ace_3d;
+  typedef Kokkos::View<double**[3], DeviceType> t_ace_3d3;
+  typedef Kokkos::View<double**[4], DeviceType> t_ace_3d4;
+  typedef Kokkos::View<double****, DeviceType> t_ace_4d;
+  typedef Kokkos::View<complex*, DeviceType> t_ace_1c;
+  typedef Kokkos::View<complex**, DeviceType> t_ace_2c;
+  typedef Kokkos::View<complex***, DeviceType> t_ace_3c;
+  typedef Kokkos::View<complex**[3], DeviceType> t_ace_3c3;
+  typedef Kokkos::View<complex****, DeviceType> t_ace_4c;
+  typedef Kokkos::View<complex***[3], DeviceType> t_ace_4c3;
+
+  typedef Kokkos::View<double*>::HostMirror th_ace_1d;
+
+  t_ace_3d A_rank1;
+  t_ace_4c A;
+
+  t_ace_3c A_list;
+  t_ace_3c A_forward_prod;
+
+  t_ace_3d weights_rank1;
+  t_ace_4c weights;
+
+  t_ace_1d e_atom;
+  t_ace_2d rhos;
+  t_ace_2d dF_drho;
+
+  // 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;
+
+  // radial functions
+  t_ace_4d fr;
+  t_ace_4d dfr;
+  t_ace_3d gr;
+  t_ace_3d dgr;
+  t_ace_3d d_values;
+  t_ace_3d d_derivatives;
+
+  // inverted active set
+  tc_ace_3d d_ASI;
+  t_ace_2d projections;
+  t_ace_1d d_gamma;
+  th_ace_1d h_gamma;
+
+  // Spherical Harmonics
+  void pre_compute_harmonics(int);
+
+  KOKKOS_INLINE_FUNCTION
+  void compute_barplm(double rz, int lmaxi);
+
+  KOKKOS_INLINE_FUNCTION
+  void compute_ylm(double rx, double ry, double rz, int lmaxi);
+
+  t_ace_1d alm;
+  t_ace_1d blm;
+  t_ace_1d cl;
+  t_ace_1d dl;
+
+  t_ace_3d plm;
+  t_ace_3d dplm;
+
+  t_ace_3c ylm;
+  t_ace_4c3 dylm;
+
+  // short neigh list
+  t_ace_1i d_ncount;
+  t_ace_2d d_mu;
+  t_ace_2d d_rnorms;
+  t_ace_3d3 d_rhats;
+  t_ace_2i d_nearest;
+
+  // per-type
+  t_ace_1i d_ndensity;
+  t_ace_1i d_npoti;
+  t_ace_1d d_rho_core_cutoff;
+  t_ace_1d d_drho_core_cutoff;
+  t_ace_1d d_E0vals;
+  t_ace_2d d_wpre;
+  t_ace_2d d_mexp;
+
+  // tilde
+  t_ace_1i d_idx_ms_combs_count;
+  t_ace_1i d_total_basis_size;
+  t_ace_2i d_rank;
+  t_ace_2i d_num_ms_combs;
+  t_ace_2i d_func_inds;
+  t_ace_3i d_mus;
+  t_ace_3i d_ns;
+  t_ace_3i d_ls;
+  t_ace_3i d_ms_combs;
+//  t_ace_3d d_ctildes;
+  t_ace_2d d_gen_cgs;
+  t_ace_3d d_coeffs;
+
+  t_ace_3d3 f_ij;
+
+ public:
+  struct SplineInterpolatorKokkos {
+    int ntot, nlut, num_of_functions;
+    double cutoff, deltaSplineBins, invrscalelookup, rscalelookup;
+
+    t_ace_3d4 lookupTable;
+
+    void operator=(const SplineInterpolator &spline);
+
+    void deallocate() {
+      lookupTable = t_ace_3d4();
+    }
+
+    double memory_usage() {
+      return lookupTable.span() * sizeof(typename decltype(lookupTable)::value_type);
+    }
+
+    KOKKOS_INLINE_FUNCTION
+    void calcSplines(const int ii, const int jj, const double r, const t_ace_3d &d_values, const t_ace_3d &d_derivatives) const;
+  };
+
+  Kokkos::DualView<SplineInterpolatorKokkos**, DeviceType> k_splines_gk;
+  Kokkos::DualView<SplineInterpolatorKokkos**, DeviceType> k_splines_rnl;
+  Kokkos::DualView<SplineInterpolatorKokkos**, DeviceType> k_splines_hc;
+
+};
+}    // namespace LAMMPS_NS
+
+#endif
+#endif

src/KOKKOS/pair_pace_kokkos.cpp

@@ -1574,7 +1574,25 @@ void PairPACEKokkos<DeviceType>::evaluate_splines(const int ii, const int jj, do
 }
 
 /* ---------------------------------------------------------------------- */
+template<class DeviceType>
+void PairPACEKokkos<DeviceType>::SplineInterpolatorKokkos::operator=(const SplineInterpolator &spline) {
+    cutoff = spline.cutoff;
+    deltaSplineBins = spline.deltaSplineBins;
+    ntot = spline.ntot;
+    nlut = spline.nlut;
+    invrscalelookup = spline.invrscalelookup;
+    rscalelookup = spline.rscalelookup;
+    num_of_functions = spline.num_of_functions;
 
+    lookupTable = t_ace_3d4("lookupTable", ntot+1, num_of_functions);
+    auto h_lookupTable = Kokkos::create_mirror_view(lookupTable);
+    for (int i = 0; i < ntot+1; i++)
+        for (int j = 0; j < num_of_functions; j++)
+            for (int k = 0; k < 4; k++)
+                h_lookupTable(i, j, k) = spline.lookupTable(i, j, k);
+    Kokkos::deep_copy(lookupTable, h_lookupTable);
+}
+/* ---------------------------------------------------------------------- */
 template<class DeviceType>
 KOKKOS_INLINE_FUNCTION
 void PairPACEKokkos<DeviceType>::SplineInterpolatorKokkos::calcSplines(const int ii, const int jj, const double r, const t_ace_3d &d_values, const t_ace_3d &d_derivatives) const

src/KOKKOS/pair_pace_kokkos.h

@@ -24,10 +24,11 @@ PairStyle(pace/kk/host,PairPACEKokkos<LMPHostType>);
 #define LMP_PAIR_PACE_KOKKOS_H
 
 #include "pair_pace.h"
-#include "ace-evaluator/ace_radial.h"
 #include "kokkos_type.h"
 #include "pair_kokkos.h"
 
+class SplineInterpolator;
+
 namespace LAMMPS_NS {
 
 template<class DeviceType>
@@ -293,23 +294,7 @@ class PairPACEKokkos : public PairPACE {
 
     t_ace_3d4 lookupTable;
 
-    void operator=(const SplineInterpolator &spline) {
+    void operator=(const SplineInterpolator &spline);
-      cutoff = spline.cutoff;
-      deltaSplineBins = spline.deltaSplineBins;
-      ntot = spline.ntot;
-      nlut = spline.nlut;
-      invrscalelookup = spline.invrscalelookup;
-      rscalelookup = spline.rscalelookup;
-      num_of_functions = spline.num_of_functions;
-
-      lookupTable = t_ace_3d4("lookupTable", ntot+1, num_of_functions);
-      auto h_lookupTable = Kokkos::create_mirror_view(lookupTable);
-      for (int i = 0; i < ntot+1; i++)
-        for (int j = 0; j < num_of_functions; j++)
-          for (int k = 0; k < 4; k++)
-            h_lookupTable(i, j, k) = spline.lookupTable(i, j, k);
-      Kokkos::deep_copy(lookupTable, h_lookupTable);
-    }
 
     void deallocate() {
       lookupTable = t_ace_3d4();

src/ML-PACE/pair_pace_extrapolation.cpp

@@ -98,6 +98,8 @@ PairPACEExtrapolation::PairPACEExtrapolation(LAMMPS *lmp) : Pair(lmp)
   scale = nullptr;
   flag_compute_extrapolation_grade = 0;
   extrapolation_grade_gamma = nullptr;
+
+  chunksize = 4096;
 }
 
 /* ----------------------------------------------------------------------
@@ -133,16 +135,12 @@ void PairPACEExtrapolation::compute(int eflag, int vflag)
 
   double **x = atom->x;
   double **f = atom->f;
-  tagint *tag = atom->tag;
   int *type = atom->type;
   // number of atoms in cell
   int nlocal = atom->nlocal;
 
   int newton_pair = force->newton_pair;
 
-  // number of atoms including ghost atoms
-  int nall = nlocal + atom->nghost;
-
   // inum: length of the neighborlists list
   inum = list->inum;
 
@@ -283,7 +281,20 @@ void PairPACEExtrapolation::allocate()
 
 void PairPACEExtrapolation::settings(int narg, char **arg)
 {
-  if (narg > 0) error->all(FLERR, "Pair style pace/extrapolation supports no keywords");
+//  if (narg > 2) error->all(FLERR, "Pair style pace/extrapolation supports no keywords");
+  if (narg > 2) utils::missing_cmd_args(FLERR, "pair_style pace/extrapolation", error);
+  // ACE potentials are parameterized in metal units
+  if (strcmp("metal", update->unit_style) != 0)
+    error->all(FLERR, "ACE potentials require 'metal' units");
+
+  int iarg = 0;
+  while (iarg < narg) {
+      if (strcmp(arg[iarg], "chunksize") == 0) {
+          chunksize = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
+          iarg += 2;
+      } else
+          error->all(FLERR, "Unknown pair_style pace keyword: {}", arg[iarg]);
+  }
 
   if (comm->me == 0)
     utils::logmesg(lmp, "ACE/AL version: {}.{}.{}\n", VERSION_YEAR, VERSION_MONTH, VERSION_DAY);
@@ -343,7 +354,6 @@ void PairPACEExtrapolation::coeff(int narg, char **arg)
   aceimpl->rec_ace->element_type_mapping.init(atom->ntypes + 1);
   aceimpl->rec_ace->element_type_mapping.fill(-1);    //-1 means atom not included into potential
 
-  FILE *species_type_file = nullptr;
 
   const int n = atom->ntypes;
   element_names.resize(n);

src/ML-PACE/pair_pace_extrapolation.h

@@ -49,13 +49,15 @@ class PairPACEExtrapolation : public Pair {
   struct ACEALImpl *aceimpl;
   int nmax;
 
-  void allocate();
+  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
 
   int flag_compute_extrapolation_grade;
 
   double **scale;
+
+  int chunksize;
 };
 
 }    // namespace LAMMPS_NS

src/fix_pair.cpp

@@ -21,6 +21,7 @@
 #include "memory.h"
 #include "pair.h"
 #include "update.h"
+#include "fmt/format.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
@@ -36,7 +37,7 @@ FixPair::FixPair(LAMMPS *lmp, int narg, char **arg) :
   if (nevery < 1) error->all(FLERR,"Illegal fix pair every value: {}", nevery);
 
   pairname = utils::strdup(arg[4]);
-  pstyle = force->pair_match(pairname,1,0);
+  query_pstyle(lmp);
   if (pstyle == nullptr) error->all(FLERR,"Pair style {} for fix pair not found", pairname);
 
   nfield = (narg-5) / 2;
@@ -130,6 +131,28 @@ FixPair::FixPair(LAMMPS *lmp, int narg, char **arg) :
   lasttime = -1;
 }
 
+/* ---------------------------------------------------------------------- */
+
+void FixPair::query_pstyle(LAMMPS *lmp) {
+    char *cptr=nullptr;
+    int nsub = 0;
+    if ((cptr = strchr(pairname, ':'))) {
+        *cptr = '\0';
+        nsub = utils::inumeric(FLERR,cptr+1,false,lmp);
+    }
+    pstyle = nullptr;
+    if (lmp->suffix_enable) {
+        if (lmp->suffix) {
+            pstyle = force->pair_match(fmt::format("{}/{}", pairname, lmp->suffix), 1, nsub);
+            if (pstyle == nullptr && (lmp->suffix2)) {
+                pstyle = force->pair_match(fmt::format("{}/{}", pairname, lmp->suffix2), 1, nsub);
+            }
+        }
+    }
+    if (pstyle == nullptr) pstyle = force->pair_match(pairname, 1, nsub);
+}
+
+
 /* ---------------------------------------------------------------------- */
 
 FixPair::~FixPair()
@@ -171,7 +194,7 @@ void FixPair::init()
 {
   // insure pair style still exists
 
-  pstyle = force->pair_match(pairname,1,0);
+  query_pstyle(lmp);
   if (pstyle == nullptr) error->all(FLERR,"Pair style {} for fix pair not found", pairname);
 }
 

src/fix_pair.h

@@ -56,6 +56,8 @@ class FixPair : public Fix {
   class Pair *pstyle;
   double *vector;
   double **array;
+
+  void query_pstyle(LAMMPS *lmp);
 };
 
 }    // namespace LAMMPS_NS