diff --git a/src/colvarcomp_coordnums.cpp b/src/colvarcomp_coordnums.cpp
index 16dd3b22f..a2372a42b 100644
--- a/src/colvarcomp_coordnums.cpp
+++ b/src/colvarcomp_coordnums.cpp
@@ -14,6 +14,196 @@
 #include "colvarcomp.h"
 #include "colvarcomp_coordnums.h"
 
+template <bool use_group1_com, bool use_group2_com, int flags, int n, int m>
+void inline colvar::coordnum::main_loop()
+{
+  constexpr bool use_pairlist = bool(flags & ef_use_pairlist);
+  constexpr bool rebuild_pairlist = bool(flags & ef_rebuild_pairlist);
+  if constexpr (rebuild_pairlist) {
+    static_assert(use_pairlist, "rebuild_pairlist requires use_pairlist.");
+  }
+  size_t const group1_num_coords = use_group1_com ? 1 : group1->size();
+  size_t const group2_num_coords = use_group2_com ? 1 : group2->size();
+
+  cvm::atom_pos const group1_com = group1->center_of_mass();
+  cvm::atom_pos const group2_com = group2->center_of_mass();
+  cvm::rvector group1_com_grad, group2_com_grad;
+
+  bool *pairlist_elem = pairlist.get();
+
+  /// Bravais lattice vectors
+  cvm::rvector unit_cell_x, unit_cell_y, unit_cell_z;
+  /// Reciprocal lattice vectors
+  cvm::rvector reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z;
+
+  if (flags & ef_use_internal_pbc) {
+    cvmodule->proxy->get_lattice(
+      unit_cell_x, unit_cell_y, unit_cell_z,
+      reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z);
+  }
+  colvarproxy_system::Boundaries_type boundaries_type = cvmodule->proxy->get_boundary_type();
+  if (flags & ef_use_internal_pbc) {
+    if (boundaries_type == colvarproxy_system::Boundaries_type::boundaries_unsupported) {
+      cvmodule->error("Error: unsupported boundary conditions.\n", COLVARS_INPUT_ERROR);
+    }
+  }
+
+  for (size_t i = 0; i < group1_num_coords; ++i) {
+
+    cvm::real const x1 = use_group1_com ? group1_com.x : group1->pos_x(i);
+    cvm::real const y1 = use_group1_com ? group1_com.y : group1->pos_y(i);
+    cvm::real const z1 = use_group1_com ? group1_com.z : group1->pos_z(i);
+    cvm::rvector g1{0, 0, 0};
+
+    for (size_t j = 0; j < group2_num_coords; ++j) {
+
+      cvm::real const x2 = use_group2_com ? group2_com.x : group2->pos_x(j);
+      cvm::real const y2 = use_group2_com ? group2_com.y : group2->pos_y(j);
+      cvm::real const z2 = use_group2_com ? group2_com.z : group2->pos_z(j);
+
+      cvm::real &gx2 = use_group2_com ? group2_com_grad.x : group2->grad_x(j);
+      cvm::real &gy2 = use_group2_com ? group2_com_grad.y : group2->grad_y(j);
+      cvm::real &gz2 = use_group2_com ? group2_com_grad.z : group2->grad_z(j);
+
+      cvm::real partial = 0;
+
+      bool const within =
+          ((flags & ef_use_pairlist) && (*pairlist_elem || (flags & ef_rebuild_pairlist))) ||
+          !(flags & ef_use_pairlist);
+
+      if (within) {
+        cvm::rvector diff{0, 0, 0};
+        if (flags & ef_use_internal_pbc) {
+          if (boundaries_type == colvarproxy_system::Boundaries_type::boundaries_non_periodic) {
+            diff = cvm::rvector{x2, y2, z2} - cvm::rvector{x1, y1, z1};
+          } else {
+            diff = colvarproxy_system::position_distance_kernel(
+              cvm::rvector{x1, y1, z1}, cvm::rvector{x2, y2, z2},
+              unit_cell_x, unit_cell_y, unit_cell_z,
+              reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z,
+              true, true, true);
+          }
+        } else {
+          diff = cvmodule->proxy->position_distance(cvm::rvector{x1, y1, z1}, cvm::rvector{x2, y2, z2});
+        }
+        partial = compute_pair_coordnum<flags, n, m>(
+          inv_r0_vec, inv_r0sq_vec, diff, en, ed,
+          g1.x, g1.y, g1.z, gx2, gy2, gz2,
+          tolerance, tolerance_l2_max);
+      }
+
+      if ((flags & ef_use_pairlist) && (flags & ef_rebuild_pairlist)) {
+        *pairlist_elem = partial > 0.0 ? true : false;
+      }
+
+      x.real_value += partial;
+
+      if (flags & ef_use_pairlist) {
+        pairlist_elem++;
+      }
+    }
+    if (flags & ef_gradients) {
+      if (!use_group1_com) {
+        group1->grad_x(i) += g1.x;
+        group1->grad_y(i) += g1.y;
+        group1->grad_z(i) += g1.z;
+      } else {
+        group1_com_grad += g1;
+      }
+    }
+  }
+
+  if (use_group1_com) {
+    group1->set_weighted_gradient(group1_com_grad);
+  }
+  if (use_group2_com) {
+    group2->set_weighted_gradient(group2_com_grad);
+  }
+}
+
+template <int t_num_tables, typename T>
+class static_function_table_impl_base {
+public:
+  inline constexpr static const int max_n = 10;
+  inline constexpr static const int max_m = 20;
+  inline constexpr static const int num_ef_combinations = 16;
+  inline constexpr static const int num_tables = t_num_tables;
+  typedef void (T::*compute_pair_coordnum_type)();
+
+  template <int flags>
+  inline constexpr const int select() const {
+    return int(bool(flags & colvar::coordnum::ef_gradients)) + \
+          (int(bool(flags & colvar::coordnum::ef_use_internal_pbc)) << 1) + \
+          (int(bool(flags & colvar::coordnum::ef_use_pairlist)) << 2) + \
+          (int(bool(flags & colvar::coordnum::ef_rebuild_pairlist)) << 3);
+  }
+
+  template <int table_index, int flags>
+  inline constexpr compute_pair_coordnum_type get_func(int n, int m) {
+    if (n <= 0) return nullptr;
+    else if (m <= 0) return nullptr;
+    else if (n > max_n) return nullptr;
+    else if (m > max_m) return nullptr;
+    else return select<flags>() < num_ef_combinations ? funcs_[table_index][select<flags>()][n-1][m-1] : nullptr;
+  }
+
+  static_function_table_impl_base() {
+    init();
+  }
+
+  void init() {
+    std::memset(funcs_, 0, sizeof(compute_pair_coordnum_type)*t_num_tables*num_ef_combinations*max_n*max_m);
+  }
+protected:
+  compute_pair_coordnum_type funcs_[t_num_tables][num_ef_combinations][max_n][max_m];
+};
+
+class colvar::coordnum::static_function_table_impl:
+  public static_function_table_impl_base<4, colvar::coordnum> {
+private:
+  template <int table_index, int n, int m>
+  void set_func() {
+    constexpr bool use_group1_com = table_index & 1;
+    constexpr bool use_group2_com = (table_index >> 1) & 1;
+#define SET_FUNC(FLAG)                                                                \
+    do {                                                                              \
+      funcs_[table_index][select<FLAG>()][n-1][m-1] =                                 \
+        &colvar::coordnum::main_loop<use_group1_com, use_group2_com, FLAG, n, m>;     \
+    } while (0)
+
+    SET_FUNC(colvar::coordnum::ef_null);
+    SET_FUNC(colvar::coordnum::ef_gradients);
+    SET_FUNC(colvar::coordnum::ef_use_internal_pbc);
+    SET_FUNC(colvar::coordnum::ef_use_pairlist);
+    // SET_FUNC(colvar::coordnum::ef_rebuild_pairlist);
+    SET_FUNC(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_internal_pbc);
+    SET_FUNC(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_pairlist);
+    // SET_FUNC(colvar::coordnum::ef_gradients | colvar::coordnum::ef_rebuild_pairlist);
+    SET_FUNC(colvar::coordnum::ef_use_internal_pbc | colvar::coordnum::ef_use_pairlist);
+    // SET_FUNC(colvar::coordnum::ef_use_internal_pbc | colvar::coordnum::ef_rebuild_pairlist);
+    SET_FUNC(colvar::coordnum::ef_use_pairlist | colvar::coordnum::ef_rebuild_pairlist);
+    SET_FUNC(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_internal_pbc | colvar::coordnum::ef_use_pairlist);
+    // SET_FUNC(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_internal_pbc | colvar::coordnum::ef_rebuild_pairlist);
+    SET_FUNC(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_pairlist | colvar::coordnum::ef_rebuild_pairlist);
+    SET_FUNC(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_internal_pbc | colvar::coordnum::ef_use_pairlist | colvar::coordnum::ef_rebuild_pairlist);
+#undef SET_FUNC
+  }
+public:
+  template <bool use_group1_com, bool use_group2_com>
+  static inline constexpr int get_table_index() {
+    return int(use_group1_com) + (int(use_group2_com) << 1);
+  }
+
+  template <int n, int m>
+  void set_func() {
+    static_assert(n <= max_n, "n is larger than max_n!");
+    static_assert(m <= max_m, "m is larger than max_m!");
+    set_func<0, n, m>();
+    set_func<1, n, m>();
+    set_func<2, n, m>();
+    set_func<3, n, m>();
+  }
+};
 
 colvar::coordnum::coordnum()
 {
@@ -21,8 +211,17 @@ colvar::coordnum::coordnum()
   x.type(colvarvalue::type_scalar);
   cvm::real const r0 = cvmodule->proxy->angstrom_to_internal(4.0);
   update_cutoffs({r0, r0, r0});
-  b_use_internal_pbc = cvm::main()->proxy->use_internal_pbc();
+  b_use_internal_pbc = cvmodule->proxy->use_internal_pbc();
   // Boundaries will be set later, when the number of pairs is known
+  static_function_table = std::make_unique<static_function_table_impl>();
+#ifdef COLVARS_COORDNUM_ENABLE_MORE_STATIC_TEMPLATE
+  // Enable some commonly used combinations of (n,m)
+  static_function_table->set_func<2, 4>();
+  static_function_table->set_func<4, 8>();
+  static_function_table->set_func<8, 16>();
+  static_function_table->set_func<10, 20>();
+#endif // COLVARS_COORDNUM_ENABLE_MORE_STATIC_TEMPLATE
+  static_function_table->set_func<6, 12>();
 }
 
 
@@ -177,7 +376,7 @@ void colvar::coordnum::compute_tolerance_l2_max()
   size_t i;
   // Find the value of l2 such that F(l2) = 0 using the Newton method
   for (i = 0; i < num_iters_max; i++) {
-    F = switching_function<ef_use_pairlist | ef_gradients>(l2, dFdl2, en, ed, tolerance);
+    F = switching_function<ef_use_pairlist | ef_gradients, 0, 0>(l2, dFdl2, en, ed, tolerance);
     if ((std::fabs(F) < result_tol) || (std::fabs(dFdl2) < dF_tol)) {
       break;
     }
@@ -191,77 +390,6 @@ void colvar::coordnum::compute_tolerance_l2_max()
 }
 
 
-template <bool use_group1_com, bool use_group2_com, int flags>
-void inline colvar::coordnum::main_loop()
-{
-  size_t const group1_num_coords = use_group1_com ? 1 : group1->size();
-  size_t const group2_num_coords = use_group2_com ? 1 : group2->size();
-
-  cvm::atom_pos const group1_com = group1->center_of_mass();
-  cvm::atom_pos const group2_com = group2->center_of_mass();
-  cvm::rvector group1_com_grad, group2_com_grad;
-
-  bool *pairlist_elem = pairlist.get();
-
-  for (size_t i = 0; i < group1_num_coords; ++i) {
-
-    cvm::real const x1 = use_group1_com ? group1_com.x : group1->pos_x(i);
-    cvm::real const y1 = use_group1_com ? group1_com.y : group1->pos_y(i);
-    cvm::real const z1 = use_group1_com ? group1_com.z : group1->pos_z(i);
-
-    cvm::real &gx1 = use_group1_com ? group1_com_grad.x : group1->grad_x(i);
-    cvm::real &gy1 = use_group1_com ? group1_com_grad.y : group1->grad_y(i);
-    cvm::real &gz1 = use_group1_com ? group1_com_grad.z : group1->grad_z(i);
-
-    for (size_t j = 0; j < group2_num_coords; ++j) {
-
-      cvm::real const x2 = use_group2_com ? group2_com.x : group2->pos_x(j);
-      cvm::real const y2 = use_group2_com ? group2_com.y : group2->pos_y(j);
-      cvm::real const z2 = use_group2_com ? group2_com.z : group2->pos_z(j);
-
-      cvm::real &gx2 = use_group2_com ? group2_com_grad.x : group2->grad_x(j);
-      cvm::real &gy2 = use_group2_com ? group2_com_grad.y : group2->grad_y(j);
-      cvm::real &gz2 = use_group2_com ? group2_com_grad.z : group2->grad_z(j);
-
-      bool const within =
-          ((flags & ef_use_pairlist) && (*pairlist_elem || (flags & ef_rebuild_pairlist))) ||
-          !(flags & ef_use_pairlist);
-
-      cvm::real const partial = within ?
-        (b_use_internal_pbc ?
-         compute_pair_coordnum<flags | ef_use_internal_pbc>(inv_r0_vec, inv_r0sq_vec, en, ed,
-                                                            x1, y1, z1, x2, y2, z2,
-                                                            gx1, gy1, gz1, gx2, gy2, gz2,
-                                                            tolerance, tolerance_l2_max,
-                                                            cvmodule) :
-         compute_pair_coordnum<flags>(inv_r0_vec, inv_r0sq_vec, en, ed,
-                                      x1, y1, z1, x2, y2, z2,
-                                      gx1, gy1, gz1, gx2, gy2, gz2,
-                                      tolerance, tolerance_l2_max,
-                                      cvmodule) ) :
-        0.0;
-
-      if ((flags & ef_use_pairlist) && (flags & ef_rebuild_pairlist)) {
-        *pairlist_elem = partial > 0.0 ? true : false;
-      }
-
-      x.real_value += partial;
-
-      if (flags & ef_use_pairlist) {
-        pairlist_elem++;
-      }
-    }
-  }
-
-  if (use_group1_com) {
-    group1->set_weighted_gradient(group1_com_grad);
-  }
-  if (use_group2_com) {
-    group2->set_weighted_gradient(group2_com_grad);
-  }
-}
-
-
 template <bool use_group1_com, bool use_group2_com, int compute_flags>
 int colvar::coordnum::compute_coordnum()
 {
@@ -271,14 +399,14 @@ int colvar::coordnum::compute_coordnum()
   if (use_pairlist) {
     if (rebuild_pairlist) {
       constexpr int flags = compute_flags | ef_use_pairlist | ef_rebuild_pairlist;
-      main_loop<use_group1_com, use_group2_com, flags>();
+      main_loop<use_group1_com, use_group2_com, flags, 0, 0>();
     } else {
       constexpr int flags = compute_flags | ef_use_pairlist;
-      main_loop<use_group1_com, use_group2_com, flags>();
+      main_loop<use_group1_com, use_group2_com, flags, 0, 0>();
     }
   } else {
     constexpr int flags = compute_flags;
-    main_loop<use_group1_com, use_group2_com, flags>();
+    main_loop<use_group1_com, use_group2_com, flags, 0, 0>();
   }
 
   return COLVARS_OK;
@@ -288,42 +416,53 @@ int colvar::coordnum::compute_coordnum()
 void colvar::coordnum::calc_value()
 {
   x.real_value = 0.0;
-  if (is_enabled(f_cvc_gradient)) {
-
-    constexpr int flags = ef_gradients;
-
-    if (b_group1_center_only) {
-      if (b_group2_center_only) {
-        compute_coordnum<true, true, flags>();
-      } else {
-        compute_coordnum<true, false, flags>();
-      }
-    } else {
-      if (b_group2_center_only) {
-        compute_coordnum<false, true, flags>();
-      } else {
-        compute_coordnum<false, false, flags>();
-      }
-    }
-
-  } else {
-
-    constexpr int flags = ef_null;
-
-    if (b_group1_center_only) {
-      if (b_group2_center_only) {
-        compute_coordnum<true, true, flags>();
-      } else {
-        compute_coordnum<true, false, flags>();
-      }
-    } else {
-      if (b_group2_center_only) {
-        compute_coordnum<false, true, flags>();
-      } else {
-        compute_coordnum<false, false, flags>();
-      }
+  bool const use_pairlist = pairlist.get();
+  bool const rebuild_pairlist = use_pairlist && (cvmodule->step_relative() % pairlist_freq == 0);
+  const bool gradients = is_enabled(f_cvc_gradient);
+  const int options = \
+    int(gradients) + \
+    (int(b_use_internal_pbc) << 8) + \
+    (int(use_pairlist) << 9) + \
+    (int(rebuild_pairlist) << 10);
+#define CALL_KERNEL(G1C, G2C, N) do { \
+  constexpr auto const table_i = \
+    colvar::coordnum::static_function_table_impl::get_table_index<G1C, G2C>(); \
+  auto kernel = static_function_table->get_func<table_i, N>(en, ed); \
+  if (kernel) { /*printf("kernel = %p\n", (void*&)kernel);*/ ((*this).*kernel)(); } \
+  else {main_loop<G1C, G2C, N, 0, 0>();} \
+} while (0)
+#define CASE(N) \
+  case N: { \
+    if (b_group1_center_only) { \
+      if (b_group2_center_only) {CALL_KERNEL(true, true, N);}\
+      else {CALL_KERNEL(true, false, N);}  \
+    } else { \
+      if (b_group2_center_only) {CALL_KERNEL(false, true, N);}\
+      else {CALL_KERNEL(false, false, N);}  \
+    } break; }
+  switch (options) {
+    CASE(colvar::coordnum::ef_null);
+    CASE(colvar::coordnum::ef_gradients);
+    CASE(colvar::coordnum::ef_use_internal_pbc);
+    CASE(colvar::coordnum::ef_use_pairlist);
+    // CASE(colvar::coordnum::ef_rebuild_pairlist);
+    CASE(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_internal_pbc);
+    CASE(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_pairlist);
+    // CASE(colvar::coordnum::ef_gradients | colvar::coordnum::ef_rebuild_pairlist);
+    CASE(colvar::coordnum::ef_use_internal_pbc | colvar::coordnum::ef_use_pairlist);
+    // CASE(colvar::coordnum::ef_use_internal_pbc | colvar::coordnum::ef_rebuild_pairlist);
+    CASE(colvar::coordnum::ef_use_pairlist | colvar::coordnum::ef_rebuild_pairlist);
+    CASE(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_internal_pbc | colvar::coordnum::ef_use_pairlist);
+    // CASE(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_internal_pbc | colvar::coordnum::ef_rebuild_pairlist);
+    CASE(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_pairlist | colvar::coordnum::ef_rebuild_pairlist);
+    CASE(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_internal_pbc | colvar::coordnum::ef_use_pairlist | colvar::coordnum::ef_rebuild_pairlist);
+    default: {
+      cvmodule->error("Unknown kernel call in colvar::coordnum::calc_value\n",
+                      COLVARS_BUG_ERROR);
     }
   }
+#undef CASE
+#undef CALL_KERNEL
 }
 
 
@@ -482,67 +621,134 @@ void colvar::h_bond::calc_gradients()
                                          cvmodule);
 }
 
-
-colvar::selfcoordnum::selfcoordnum()
+template <int flags, int n, int m> inline void colvar::selfcoordnum::selfcoordnum_sequential_loop()
 {
-  set_function_type("selfCoordNum");
-}
+  constexpr bool use_pairlist = bool(flags & ef_use_pairlist);
+  constexpr bool rebuild_pairlist = bool(flags & ef_rebuild_pairlist);
+  if constexpr (rebuild_pairlist) {
+    static_assert(use_pairlist, "rebuild_pairlist requires use_pairlist.");
+  }
+  size_t const natoms = group1->size();
+  bool *pairlist_elem = pairlist.get();
 
+  /// Bravais lattice vectors
+  cvm::rvector unit_cell_x, unit_cell_y, unit_cell_z;
+  /// Reciprocal lattice vectors
+  cvm::rvector reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z;
 
-template <int flags> inline void colvar::selfcoordnum::selfcoordnum_sequential_loop()
-{
-  size_t const n = group1->size();
-  bool *pairlist_elem = pairlist.get();
+  if (flags & ef_use_internal_pbc) {
+    cvmodule->proxy->get_lattice(
+      unit_cell_x, unit_cell_y, unit_cell_z,
+      reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z);
+  }
+  colvarproxy_system::Boundaries_type boundaries_type = cvmodule->proxy->get_boundary_type();
+  if (flags & ef_use_internal_pbc) {
+    if (boundaries_type == colvarproxy_system::Boundaries_type::boundaries_unsupported) {
+      cvmodule->error("Error: unsupported boundary conditions.\n", COLVARS_INPUT_ERROR);
+    }
+  }
 
-  for (size_t i = 0; i < n - 1; i++) {
+  for (size_t i = 0; i < natoms - 1; i++) {
 
     cvm::real const x1 = group1->pos_x(i);
     cvm::real const y1 = group1->pos_y(i);
     cvm::real const z1 = group1->pos_z(i);
-    cvm::real &gx1 = group1->grad_x(i);
-    cvm::real &gy1 = group1->grad_y(i);
-    cvm::real &gz1 = group1->grad_z(i);
-
-    for (size_t j = i + 1; j < n; j++) {
+    cvm::rvector g1{0, 0, 0};
 
+    for (size_t j = i + 1; j < natoms; j++) {
       cvm::real const x2 = group1->pos_x(j);
       cvm::real const y2 = group1->pos_y(j);
       cvm::real const z2 = group1->pos_z(j);
       cvm::real &gx2 = group1->grad_x(j);
       cvm::real &gy2 = group1->grad_y(j);
       cvm::real &gz2 = group1->grad_z(j);
-
+      cvm::real partial = 0;
       bool const within =
         ((flags & ef_use_pairlist) && (*pairlist_elem || (flags & ef_rebuild_pairlist))) ||
         !(flags & ef_use_pairlist);
+      if (within) {
+        cvm::rvector diff{0, 0, 0};
+        if constexpr (flags & ef_use_internal_pbc) {
+          if (boundaries_type == colvarproxy_system::Boundaries_type::boundaries_non_periodic) {
+            diff = cvm::rvector{x2, y2, z2} - cvm::rvector{x1, y1, z1};
+          } else {
+            diff = colvarproxy_system::position_distance_kernel(
+              cvm::rvector{x1, y1, z1}, cvm::rvector{x2, y2, z2},
+              unit_cell_x, unit_cell_y, unit_cell_z,
+              reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z,
+              true, true, true);
+          }
+        } else {
+          diff = cvmodule->proxy->position_distance(cvm::rvector{x1, y1, z1}, cvm::rvector{x2, y2, z2});
+        }
+        partial = compute_pair_coordnum<flags, n, m>(
+          inv_r0_vec, inv_r0sq_vec, diff, en, ed,
+          g1.x, g1.y, g1.z, gx2, gy2, gz2,
+          tolerance, tolerance_l2_max);
+      }
 
-      cvm::real const partial = within ?
-        (b_use_internal_pbc ?
-         compute_pair_coordnum<flags | ef_use_internal_pbc>(inv_r0_vec, inv_r0sq_vec, en, ed,
-                                                            x1, y1, z1, x2, y2, z2,
-                                                            gx1, gy1, gz1, gx2, gy2, gz2,
-                                                            tolerance, tolerance_l2_max,
-                                                            cvmodule) :
-         compute_pair_coordnum<flags>(inv_r0_vec, inv_r0sq_vec, en, ed,
-                                      x1, y1, z1, x2, y2, z2,
-                                      gx1, gy1, gz1, gx2, gy2, gz2,
-                                      tolerance, tolerance_l2_max,
-                                      cvmodule) ) :
-        0.0;
-
-      if ((flags & ef_use_pairlist) && (flags & ef_rebuild_pairlist)) {
+      if constexpr ((flags & ef_use_pairlist) && (flags & ef_rebuild_pairlist)) {
         *pairlist_elem = partial > 0.0 ? true : false;
       }
 
       x.real_value += partial;
 
-      if (flags & ef_use_pairlist) {
+      if constexpr (flags & ef_use_pairlist) {
         pairlist_elem++;
       }
     }
+    if constexpr (flags & ef_gradients) {
+      group1->grad_x(i) += g1.x;
+      group1->grad_y(i) += g1.y;
+      group1->grad_z(i) += g1.z;
+    }
   }
 }
 
+class colvar::selfcoordnum::static_function_table_impl:
+  public static_function_table_impl_base<1, colvar::selfcoordnum> {
+public:
+  template <int n, int m>
+  void set_func() {
+    static_assert(n <= max_n, "n is larger than max_n!");
+    static_assert(m <= max_m, "m is larger than max_m!");
+#define SET_FUNC(FLAG) \
+    do { \
+      funcs_[0][select<FLAG>()][n-1][m-1] = &colvar::selfcoordnum::selfcoordnum_sequential_loop<FLAG, n, m>;\
+    } while (0)
+    SET_FUNC(colvar::coordnum::ef_null);
+    SET_FUNC(colvar::coordnum::ef_gradients);
+    SET_FUNC(colvar::coordnum::ef_use_internal_pbc);
+    SET_FUNC(colvar::coordnum::ef_use_pairlist);
+    // SET_FUNC(colvar::coordnum::ef_rebuild_pairlist);
+    SET_FUNC(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_internal_pbc);
+    SET_FUNC(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_pairlist);
+    // SET_FUNC(colvar::coordnum::ef_gradients | colvar::coordnum::ef_rebuild_pairlist);
+    SET_FUNC(colvar::coordnum::ef_use_internal_pbc | colvar::coordnum::ef_use_pairlist);
+    // SET_FUNC(colvar::coordnum::ef_use_internal_pbc | colvar::coordnum::ef_rebuild_pairlist);
+    SET_FUNC(colvar::coordnum::ef_use_pairlist | colvar::coordnum::ef_rebuild_pairlist);
+    SET_FUNC(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_internal_pbc | colvar::coordnum::ef_use_pairlist);
+    // SET_FUNC(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_internal_pbc | colvar::coordnum::ef_rebuild_pairlist);
+    SET_FUNC(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_pairlist | colvar::coordnum::ef_rebuild_pairlist);
+    SET_FUNC(colvar::coordnum::ef_gradients | colvar::coordnum::ef_use_internal_pbc | colvar::coordnum::ef_use_pairlist | colvar::coordnum::ef_rebuild_pairlist);
+#undef SET_FUNC
+  }
+};
+
+colvar::selfcoordnum::selfcoordnum()
+{
+  set_function_type("selfCoordNum");
+  static_function_table = std::make_unique<static_function_table_impl>();
+#ifdef COLVARS_COORDNUM_ENABLE_MORE_STATIC_TEMPLATE
+  // Enable some commonly used combinations of (n,m)
+  static_function_table->set_func<2, 4>();
+  static_function_table->set_func<4, 8>();
+  static_function_table->set_func<8, 16>();
+  static_function_table->set_func<10, 20>();
+#endif // COLVARS_COORDNUM_ENABLE_MORE_STATIC_TEMPLATE
+  static_function_table->set_func<6, 12>();
+}
+
 
 template<int compute_flags> int colvar::selfcoordnum::compute_selfcoordnum()
 {
@@ -551,15 +757,60 @@ template<int compute_flags> int colvar::selfcoordnum::compute_selfcoordnum()
 
   if (use_pairlist) {
     if (rebuild_pairlist) {
-      int constexpr flags = compute_flags | ef_use_pairlist | ef_rebuild_pairlist;
-      selfcoordnum_sequential_loop<flags>();
+      if (b_use_internal_pbc) {
+        int constexpr flags = compute_flags | ef_use_pairlist | ef_rebuild_pairlist | ef_use_internal_pbc;
+        auto kernel = static_function_table->get_func<0, flags>(en, ed);
+        if (kernel) {
+          ((*this).*kernel)();
+        } else {
+          selfcoordnum_sequential_loop<flags, 0, 0>();
+        }
+      } else {
+        int constexpr flags = compute_flags | ef_use_pairlist | ef_rebuild_pairlist;
+        auto kernel = static_function_table->get_func<0, flags>(en, ed);
+        if (kernel) {
+          ((*this).*kernel)();
+        } else {
+          selfcoordnum_sequential_loop<flags, 0, 0>();
+        }
+      }
     } else {
-      int constexpr flags = compute_flags | ef_use_pairlist;
-      selfcoordnum_sequential_loop<flags>();
+      if (b_use_internal_pbc) {
+        int constexpr flags = compute_flags | ef_use_pairlist | ef_use_internal_pbc;
+        auto kernel = static_function_table->get_func<0, flags>(en, ed);
+        if (kernel) {
+          ((*this).*kernel)();
+        } else {
+          selfcoordnum_sequential_loop<flags, 0, 0>();
+        }
+      } else {
+        int constexpr flags = compute_flags | ef_use_pairlist;
+        auto kernel = static_function_table->get_func<0, flags>(en, ed);
+        if (kernel) {
+          ((*this).*kernel)();
+        } else {
+          selfcoordnum_sequential_loop<flags, 0, 0>();
+        }
+      }
     }
   } else {
-    int constexpr flags = compute_flags | ef_null;
-    selfcoordnum_sequential_loop<flags>();
+    if (b_use_internal_pbc) {
+      int constexpr flags = compute_flags | ef_null | ef_use_internal_pbc;
+      auto kernel = static_function_table->get_func<0, flags>(en, ed);
+      if (kernel) {
+        ((*this).*kernel)();
+      } else {
+        selfcoordnum_sequential_loop<flags, 0, 0>();
+      }
+    } else {
+      int constexpr flags = compute_flags | ef_null;
+      auto kernel = static_function_table->get_func<0, flags>(en, ed);
+      if (kernel) {
+        ((*this).*kernel)();
+      } else {
+        selfcoordnum_sequential_loop<flags, 0, 0>();
+      }
+    }
   }
   return COLVARS_OK;
 }
diff --git a/src/colvarcomp_coordnums.h b/src/colvarcomp_coordnums.h
index 650fa5b2d..d1336a4dc 100644
--- a/src/colvarcomp_coordnums.h
+++ b/src/colvarcomp_coordnums.h
@@ -37,32 +37,42 @@ class colvar::coordnum : public colvar::cvc {
     ef_rebuild_pairlist = (1 << 10)
   };
 
+  /// Main kernel for the coordination number
+  template <int flags>
+  inline static cvm::real compute_pair_coordnum(cvm::rvector const &inv_r0_vec,
+                                         cvm::rvector const &inv_r0sq_vec, int en, int ed,
+                                         const cvm::real a1x, const cvm::real a1y, const cvm::real a1z,
+                                         const cvm::real a2x, const cvm::real a2y, const cvm::real a2z,
+                                         cvm::real &g1x, cvm::real &g1y, cvm::real &g1z,
+                                         cvm::real &g2x, cvm::real &g2y, cvm::real &g2z,
+                                         cvm::real pairlist_tol, cvm::real pairlist_tol_l2_max,
+                                         colvarmodule *cvmodule);
+
   /// Compute the switching function (1-l2**(en/2))/(1-l2**(ed/2)) for a given squared scaled distance
   /// @param[in] l2 Square norm of (Dx/r0x, Dy/r0y, Dz/r0z)
   /// @param[out] dFdl2 Derivative of the result with respect to l2
   /// @param en Numerator exponent
   /// @param ed Denominator exponent
   /// @param pairlist_tol Pairlist tolerance
-  template <int flags>
-  static cvm::real switching_function(cvm::real const &l2, cvm::real &dFdl2, int en, int ed,
+  template <int flags, int t_en, int t_ed>
+  inline static cvm::real switching_function(cvm::real const &l2, cvm::real &dFdl2,
+                                      int en, int ed,
                                       cvm::real pairlist_tol);
 
-  /// Main kernel for the coordination number
-  template <int flags>
-  static cvm::real compute_pair_coordnum(cvm::rvector const &inv_r0_vec,
-                                         cvm::rvector const &inv_r0sq_vec, int en, int ed,
-                                         const cvm::real a1x, const cvm::real a1y, const cvm::real a1z,
-                                         const cvm::real a2x, const cvm::real a2y, const cvm::real a2z,
+  template <int flags, int t_en, int t_ed>
+  inline static cvm::real compute_pair_coordnum(cvm::rvector const &inv_r0_vec,
+                                         cvm::rvector const &inv_r0sq_vec,
+                                         const cvm::rvector& diff,
+                                         int en, int ed,
                                          cvm::real &g1x, cvm::real &g1y, cvm::real &g1z,
                                          cvm::real &g2x, cvm::real &g2y, cvm::real &g2z,
-                                         cvm::real pairlist_tol, cvm::real pairlist_tol_l2_max,
-                                         colvarmodule *cvmodule);
+                                         cvm::real pairlist_tol, cvm::real pairlist_tol_l2_max);
 
   /// Workhorse function
   template <bool use_group1_com, bool use_group2_com, int flags> int compute_coordnum();
 
   /// Workhorse function
-  template <bool use_group1_com, bool use_group2_com, int flags> void main_loop();
+  template <bool use_group1_com, bool use_group2_com, int flags, int n, int m> void main_loop();
 
 protected:
   /// First atom group
@@ -114,12 +124,18 @@ class colvar::coordnum : public colvar::cvc {
   /// Pair list
   std::unique_ptr<bool []> pairlist;
 
+private:
+  class static_function_table_impl;
+  std::unique_ptr<static_function_table_impl> static_function_table;
 };
 
 
 /// \brief Colvar component: self-coordination number within a group
 /// (colvarvalue::type_scalar type, range [0:N*(N-1)/2])
 class colvar::selfcoordnum : public colvar::coordnum {
+private:
+  class static_function_table_impl;
+  std::unique_ptr<static_function_table_impl> static_function_table;
 public:
 
   selfcoordnum();
@@ -127,7 +143,7 @@ class colvar::selfcoordnum : public colvar::coordnum {
   virtual void calc_gradients();
 
   /// Workhorse function
-  template <int flags> void selfcoordnum_sequential_loop();
+  template <int flags, int n, int m> void selfcoordnum_sequential_loop();
 
   /// Main workhorse function
   template <int flags> int compute_selfcoordnum();
@@ -169,64 +185,85 @@ class colvar::h_bond : public colvar::cvc {
 };
 
 
-template <int flags>
-inline cvm::real colvar::coordnum::switching_function(cvm::real const &l2, cvm::real &dFdl2,
-                                               int en, int ed,
-                                               cvm::real pairlist_tol)
+template <int flags, int t_en, int t_ed>
+inline cvm::real colvar::coordnum::switching_function(
+  cvm::real const &l2, cvm::real &dFdl2, int en, int ed, cvm::real pairlist_tol)
 {
-  // Assume en and ed are even integers, and avoid sqrt in the following
-  int const en2 = en/2;
-  int const ed2 = ed/2;
-
-  cvm::real const xn = cvm::integer_power(l2, en2);
-  cvm::real const xd = cvm::integer_power(l2, ed2);
-  cvm::real const eps_l2 = 1.0e-7;
-  cvm::real const h = l2 - 1.0;
-  cvm::real const en2_r = (cvm::real) en2;
-  cvm::real const ed2_r = (cvm::real) ed2;
-  cvm::real func_no_pairlist;
-
-  if (std::abs(h) < eps_l2) {
-    // Order-2 Taylor expansion: c0 + c1*h + c2*h^2
-    cvm::real const c0 = en2_r / ed2_r;
-    cvm::real const c1 = (en2_r * (en2_r - ed2_r)) / (2.0 * ed2_r);
-    cvm::real const c2 = (en2_r * (en2_r - ed2_r) * (2.0 * en2_r - ed2_r - 3.0)) / (12.0 * ed2_r);
-    func_no_pairlist = c0 + h * (c1 + h * c2);
-  } else {
-    func_no_pairlist = (1.0 - xn) / (1.0 - xd);
-  }
-
-  cvm::real func, inv_one_pairlist_tol;
-  if (flags & ef_use_pairlist) {
-    inv_one_pairlist_tol = 1 / (1.0-pairlist_tol);
-    func = (func_no_pairlist - pairlist_tol) * inv_one_pairlist_tol;
+  constexpr bool ed_two_en = (t_ed == 2 * t_en);
+  if constexpr (ed_two_en && t_en != 0) {
+    static_assert(t_en % 2 == 0, "Unsupported instantiation of N (N % 2 != 0) in colvar::coordnum::switching_function.");
+    cvm::real func_no_pairlist, func, inv_one_pairlist_tol;
+    cvm::real xn = cvm::integer_power<t_en/2>(l2);
+    func_no_pairlist = 1.0 / (1.0 + xn);
+    if (flags & ef_use_pairlist) {
+      inv_one_pairlist_tol = 1 / (1.0-pairlist_tol);
+      func = (func_no_pairlist - pairlist_tol) * inv_one_pairlist_tol;
+    } else {
+      func = func_no_pairlist;
+    }
+    if (func < 0)
+      return 0;
+    if (flags & ef_gradients) {
+      if (flags & ef_use_pairlist) {
+        dFdl2 = -0.5 * (func_no_pairlist * func_no_pairlist) * t_en * xn / l2 * (inv_one_pairlist_tol);
+      } else {
+        dFdl2 = -0.5 * (func_no_pairlist * func_no_pairlist) * t_en * xn / l2;
+      }
+    }
+    return func;
   } else {
-    func = func_no_pairlist;
-  }
-
-  // If the value is too small and we are correcting for the tolerance, the result is negative
-  // and we need to exclude it rather than let it contribute to the sum or the gradients.
-  if (func < 0)
-    return 0;
+    // Assume en and ed are even integers, and avoid sqrt in the following
+    const int en2 = t_en != 0 ? t_en / 2 : en / 2;
+    const int ed2 = t_ed != 0 ? t_ed / 2 : ed / 2;
+    cvm::real const xn = t_en != 0 ? cvm::integer_power<t_en / 2>(l2) : cvm::integer_power(l2, en2);
+    cvm::real const xd = t_ed != 0 ? cvm::integer_power<t_ed / 2>(l2) : cvm::integer_power(l2, ed2);
+    cvm::real const eps_l2 = 1.0e-7;
+    cvm::real const h = l2 - 1.0;
+    cvm::real const en2_r = (cvm::real) en2;
+    cvm::real const ed2_r = (cvm::real) ed2;
+    cvm::real func_no_pairlist;
 
-  if (flags & ef_gradients) {
-    // Logarithmic derivative: 1st-order Taylor expansion around l2 = 1
-    cvm::real log_deriv;
     if (std::abs(h) < eps_l2) {
-      cvm::real const g0 = 0.5 * (en2_r - ed2_r);
-      cvm::real const g1 = ((en2_r - ed2_r) * (en2_r + ed2_r - 6.0)) / 12.0;
-      log_deriv = g0 + h * g1;
+      // Order-2 Taylor expansion: c0 + c1*h + c2*h^2
+      cvm::real const c0 = en2_r / ed2_r;
+      cvm::real const c1 = (en2_r * (en2_r - ed2_r)) / (2.0 * ed2_r);
+      cvm::real const c2 = (en2_r * (en2_r - ed2_r) * (2.0 * en2_r - ed2_r - 3.0)) / (12.0 * ed2_r);
+      func_no_pairlist = c0 + h * (c1 + h * c2);
     } else {
-      log_deriv = (ed2_r * xd / ((1.0 - xd) * l2)) - (en2_r * xn / ((1.0 - xn) * l2));
+      func_no_pairlist = (1.0 - xn) / (1.0 - xd);
     }
-    dFdl2 = (flags & ef_use_pairlist) ?
-      func_no_pairlist * inv_one_pairlist_tol * log_deriv :
-      func * log_deriv;
-  }
 
-  return func;
-}
+    cvm::real func, inv_one_pairlist_tol;
+    if (flags & ef_use_pairlist) {
+      inv_one_pairlist_tol = 1 / (1.0-pairlist_tol);
+      func = (func_no_pairlist - pairlist_tol) * inv_one_pairlist_tol;
+    } else {
+      func = func_no_pairlist;
+    }
+
+    // If the value is too small and we are correcting for the tolerance, the result is negative
+    // and we need to exclude it rather than let it contribute to the sum or the gradients.
+    if (func < 0)
+      return 0;
+
+    if (flags & ef_gradients) {
+      // Logarithmic derivative: 1st-order Taylor expansion around l2 = 1
+      cvm::real log_deriv;
+      if (std::abs(h) < eps_l2) {
+        cvm::real const g0 = 0.5 * (en2_r - ed2_r);
+        cvm::real const g1 = ((en2_r - ed2_r) * (en2_r + ed2_r - 6.0)) / 12.0;
+        log_deriv = g0 + h * g1;
+      } else {
+        log_deriv = (ed2_r * xd / ((1.0 - xd) * l2)) - (en2_r * xn / ((1.0 - xn) * l2));
+      }
+      dFdl2 = (flags & ef_use_pairlist) ?
+        func_no_pairlist * inv_one_pairlist_tol * log_deriv :
+        func * log_deriv;
+    }
 
+    return func;
+  }
+}
 
 template<int flags>
 inline cvm::real colvar::coordnum::compute_pair_coordnum(cvm::rvector const &inv_r0_vec,
@@ -255,6 +292,26 @@ inline cvm::real colvar::coordnum::compute_pair_coordnum(cvm::rvector const &inv
                                 ? cvmodule_in->proxy->position_distance_internal(pos1, pos2)
                                 : cvmodule_in->proxy->position_distance(pos1, pos2);
 
+  return compute_pair_coordnum<flags, 0, 0>(
+    inv_r0_vec, inv_r0sq_vec, diff, en, ed,
+    g1x, g1y, g1z, g2x, g2y, g2z,
+    pairlist_tol, pairlist_tol_l2_max);
+}
+
+template<int flags, int t_en, int t_ed>
+inline cvm::real colvar::coordnum::compute_pair_coordnum(cvm::rvector const &inv_r0_vec,
+                                                         cvm::rvector const &inv_r0sq_vec,
+                                                         const cvm::rvector& diff,
+                                                         int en, int ed,
+                                                         cvm::real& g1x,
+                                                         cvm::real& g1y,
+                                                         cvm::real& g1z,
+                                                         cvm::real& g2x,
+                                                         cvm::real& g2y,
+                                                         cvm::real& g2z,
+                                                         cvm::real pairlist_tol,
+                                                         cvm::real pairlist_tol_l2_max)
+{
   cvm::rvector const scal_diff(diff.x * inv_r0_vec.x,
                                diff.y * inv_r0_vec.y,
                                diff.z * inv_r0_vec.z);
@@ -267,7 +324,7 @@ inline cvm::real colvar::coordnum::compute_pair_coordnum(cvm::rvector const &inv
   }
 
   cvm::real dFdl2 = 0.0;
-  cvm::real F = switching_function<flags>(l2, dFdl2, en, ed, pairlist_tol);
+  const cvm::real F = switching_function<flags, t_en, t_ed>(l2, dFdl2, en, ed, pairlist_tol);
 
   if ((flags & ef_gradients) && (F > 0.0)) {
     cvm::rvector const dl2dx((2.0 * inv_r0sq_vec.x) * diff.x,
diff --git a/src/colvarmodule.h b/src/colvarmodule.h
index 16b0e0822..90bae5f1c 100644
--- a/src/colvarmodule.h
+++ b/src/colvarmodule.h
@@ -115,6 +115,20 @@ class colvarmodule {
     return (n > 0) ? yy : 1.0/yy;
   }
 
+  template <int n>
+  static inline real integer_power(real const &x)
+  {
+    // Original code: math_special.h in LAMMPS
+    double yy, ww;
+    if (x == 0.0) return 0.0;
+    int nn = (n > 0) ? n : -n;
+    ww = x;
+    for (yy = 1.0; nn != 0; nn >>= 1, ww *=ww) {
+      if (nn & 1) yy *= ww;
+    }
+    return (n > 0) ? yy : 1.0/yy;
+  }
+
   /// Reimplemented to work around MS compiler issues
   static inline real pow(real const &x, real const &y)
   {
diff --git a/src/colvarproxy_system.h b/src/colvarproxy_system.h
index 3746a5b66..1ebabc414 100644
--- a/src/colvarproxy_system.h
+++ b/src/colvarproxy_system.h
@@ -179,6 +179,37 @@ class colvarproxy_system {
     return false;
   }
 
+  void get_lattice(
+    cvm::rvector& unit_cell_x_out,
+    cvm::rvector& unit_cell_y_out,
+    cvm::rvector& unit_cell_z_out,
+    cvm::rvector& reciprocal_cell_x_out,
+    cvm::rvector& reciprocal_cell_y_out,
+    cvm::rvector& reciprocal_cell_z_out) const {
+    unit_cell_x_out = unit_cell_x;
+    unit_cell_y_out = unit_cell_y;
+    unit_cell_z_out = unit_cell_z;
+    reciprocal_cell_x_out = reciprocal_cell_x;
+    reciprocal_cell_y_out = reciprocal_cell_y;
+    reciprocal_cell_z_out = reciprocal_cell_z;
+  }
+
+  /// Type of boundary conditions defined for the current computation
+  ///
+  /// Orthogonal and triclinic cells are made available to objects.
+  /// For any other conditions (mixed periodicity, triclinic cells in LAMMPS)
+  /// minimum-image distances are computed by the host engine by default
+  enum Boundaries_type {
+    boundaries_non_periodic,
+    boundaries_pbc_ortho,
+    boundaries_pbc_triclinic,
+    boundaries_unsupported
+  };
+
+  Boundaries_type get_boundary_type() const {
+    return boundaries_type;
+  }
+
 protected:
 
   /// Next value of lambda to be sent to back-end
@@ -214,18 +245,6 @@ class colvarproxy_system {
   /// Use the PBC functions from the Colvars library (as opposed to MD engine)
   bool use_internal_pbc_ = false;
 
-  /// Type of boundary conditions defined for the current computation
-  ///
-  /// Orthogonal and triclinic cells are made available to objects.
-  /// For any other conditions (mixed periodicity, triclinic cells in LAMMPS)
-  /// minimum-image distances are computed by the host engine by default
-  enum Boundaries_type {
-    boundaries_non_periodic,
-    boundaries_pbc_ortho,
-    boundaries_pbc_triclinic,
-    boundaries_unsupported
-  };
-
   /// Type of boundary conditions
   Boundaries_type boundaries_type;