runtime switchable solid harmonics ordering

export/cmake/CMakeLists.txt.export

@@ -49,6 +49,7 @@ endif(ENABLE_FORTRAN)
 # these are known orderings, must match config.h
 set(LIBINT_SHGSHELL_ORDERING_STANDARD 1)
 set(LIBINT_SHGSHELL_ORDERING_GAUSSIAN 2)
+# August 2023: the following variable has an effect on `INT_SOLIDHARMINDEX(l, m)` but otherwise the choice can be deferred to runtime.
 set(LIBINT2_SHGAUSS_ORDERING "standard" CACHE STRING "Use one of the following known orderings for shells of solid harmonic Gaussians:
   standard -- standard ordering (-l, -l+1 ... l)
   gaussian -- the Gaussian ordering (0, 1, -1, 2, -2, ... l, -l)

include/libint2/initialize.h

@@ -34,6 +34,7 @@
 #include <libint2/util/deprecated.h>
 #include <libint2/util/singleton.h>
 #include <libint2/deriv_map.h>
+#include <libint2/shgshell_ordering.h>
 
 namespace libint2 {
 
@@ -56,6 +57,10 @@ namespace libint2 {
       static std::ostream* value = &std::clog;
       return value;
     }
+    inline std::atomic<SHGShellOrdering>& solid_harmonics_ordering_accessor() {
+      static std::atomic<SHGShellOrdering> value{libint2::SHGShellOrdering_Standard};
+      return value;
+    }
   } // namespace libint2::detail
 
   /// checks if the libint has been initialized.
@@ -74,6 +79,11 @@ namespace libint2 {
       (void) x;  // to suppress unused variable warning (not guaranteed to work) TODO revise when upgrade to C++17
       assert(x != nullptr);
       verbose_accessor() = verbose;
+
+      // initialize() functions that take `SHGShellOrdering sho` as an argument aren't provided because
+      // * (a) with casting, it's hard to disentangle from `bool verbose`
+      // * (b) a separate setter is needed anyways for cases like the Python module, where initialize(sho) lives in libint code
+      // * code in initializer would go here as `solid_harmonics_ordering_accessor() = sho;`
     }
   }
 
@@ -93,8 +103,18 @@ namespace libint2 {
       managed_singleton<__initializer>::delete_instance();
       verbose_accessor() = true;
       verbose_stream_accessor() = &std::clog;
+      solid_harmonics_ordering_accessor() = libint2::SHGShellOrdering_Standard;
     }
   }
+  /// Setter for the SHGShellOrdering
+  inline void set_solid_harmonics_ordering(SHGShellOrdering sho) {
+    detail::solid_harmonics_ordering_accessor() = sho;
+  }
+  /// Accessor for the SHGShellOrdering
+  /// @return the val for Operator::nucleus
+  inline SHGShellOrdering solid_harmonics_ordering() {
+    return detail::solid_harmonics_ordering_accessor();
+  }
   /// Accessor for the disgnostics stream
   /// @return the stream to which the diagnostics will be written if verbose() returns true
   inline std::ostream& verbose_stream() {

include/libint2/shgshell_ordering.h

@@ -39,54 +39,41 @@ enum SHGShellOrdering {
 // Macros that define orderings
 //
 
-#if LIBINT_SHGSHELL_ORDERING == LIBINT_SHGSHELL_ORDERING_STANDARD
-
 /* Computes an index to a Cartesian function within a shell given
  * l = total angular momentum
  * m = real solid harmonic index (|m| = the absolute value of the projection of
  * the angular momentum on the z axis) m runs from -l to l
  */
 namespace libint2 {
-inline int INT_SOLIDHARMINDEX(int l, int m) { return m + l; }
+inline int INT_SOLIDHARMINDEX_STANDARD(int l, int m) { return m + l; }
 }
-LIBINT_DEPRECATED("please use libint2::INT_SOLIDHARMINDEX instead")
-inline int INT_SOLIDHARMINDEX(int l, int m) { return libint2::INT_SOLIDHARMINDEX(l, m); }
 
 /* This sets up the above loop over cartesian exponents as follows
  * int m;
- * FOR_SOLIDHARM(l,m)
+ * FOR_SOLIDHARM_STANDARD(l,m)
  * END_FOR_SOLIDHARM
  */
-#define FOR_SOLIDHARM(l, m) for ((m) = -(l); (m) <= (l); ++(m)) {
+#define FOR_SOLIDHARM_STANDARD(l, m) for ((m) = -(l); (m) <= (l); ++(m)) {
 #define END_FOR_SOLIDHARM }
 
-#endif  // Standard ordering
-
-#if LIBINT_SHGSHELL_ORDERING == LIBINT_SHGSHELL_ORDERING_GAUSSIAN
-
 /* Computes an index to a Cartesian function within a shell given
  * l = total angular momentum
  * m = real solid harmonic index (|m| = the absolute value of the projection of
  * the angular momentum on the z axis) m runs as 0, +1, -1, +2, -2 ... +l, -l
  */
 namespace libint2 {
-inline int INT_SOLIDHARMINDEX(int l, int m) {
+inline int INT_SOLIDHARMINDEX_GAUSSIAN(int l, int m) {
   return 2 * std::abs(m) + (m > 0 ? -1 : 0);
 }
 }
-LIBINT_DEPRECATED("please use libint2::INT_SOLIDHARMINDEX instead")
-inline int INT_SOLIDHARMINDEX(int l, int m) { return libint2::INT_SOLIDHARMINDEX(l, m); }
 
 /* This sets up the above loop over cartesian exponents as follows
  * int m;
- * FOR_SOLIDHARM(l,m)
+ * FOR_SOLIDHARM_GAUSSIAN(l,m)
  * END_FOR_SOLIDHARM
  */
-#define FOR_SOLIDHARM(l, m) \
+#define FOR_SOLIDHARM_GAUSSIAN(l, m) \
   for ((m) = 0; (m) != (l) + 1; (m) = ((m) > 0 ? -(m) : 1 - (m))) {
-#define END_FOR_SOLIDHARM }
-
-#endif  // Gaussian ordering
 
 /// these always-available macros encode orderings assumed by Molden
 
@@ -95,11 +82,38 @@ inline int INT_SOLIDHARMINDEX_MOLDEN(int l, int m) {
   return 2 * std::abs(m) + (m > 0 ? -1 : 0);
 }
 }
-LIBINT_DEPRECATED("please use libint2::INT_SOLIDHARMINDEX_MOLDEN instead")
-inline int INT_SOLIDHARMINDEX_MOLDEN(int l, int m) { return libint2::INT_SOLIDHARMINDEX_MOLDEN(l, m); }
 
 #define FOR_SOLIDHARM_MOLDEN(l, m) \
   for ((m) = 0; (m) != (l) + 1; (m) = ((m) > 0 ? -(m) : 1 - (m))) {
 #define END_FOR_SOLIDHARM_MOLDEN }
 
+namespace libint2 {
+  LIBINT_DEPRECATED("please use libint2::INT_SOLIDHARMINDEX(libint2::solid_harmonics_ordering(), l, m) instead. Current function returns the standard or gaussian index based on build configuration, not on libint2::set_solid_harmonics_ordering() argument.")
+  inline int INT_SOLIDHARMINDEX(int l, int m) {
+#if LIBINT_SHGSHELL_ORDERING == LIBINT_SHGSHELL_ORDERING_STANDARD
+    return libint2::INT_SOLIDHARMINDEX_STANDARD(l, m);
+#elif LIBINT_SHGSHELL_ORDERING == LIBINT_SHGSHELL_ORDERING_GAUSSIAN
+    return libint2::INT_SOLIDHARMINDEX_GAUSSIAN(l, m);
+#else
+#  error "unknown value of macro LIBINT_SHGSHELL_ORDERING"
+#endif
+  }
+
+  inline int INT_SOLIDHARMINDEX(int sho, int l, int m) {
+    if (sho == LIBINT_SHGSHELL_ORDERING_STANDARD)
+      return libint2::INT_SOLIDHARMINDEX_STANDARD(l, m);
+    else
+      return libint2::INT_SOLIDHARMINDEX_GAUSSIAN(l, m);
+  }
+}
+
+// continue providing a generic FOR_SOLIDHARM?
+//#if LIBINT_SHGSHELL_ORDERING == LIBINT_SHGSHELL_ORDERING_STANDARD
+//#define FOR_SOLIDHARM FOR_SOLIDHARM_STANDARD
+//#elif LIBINT_SHGSHELL_ORDERING == LIBINT_SHGSHELL_ORDERING_GAUSSIAN
+//#define FOR_SOLIDHARM FOR_SOLIDHARM_GAUSSIAN
+//#else
+//#  error "unknown value of macro LIBINT_SHGSHELL_ORDERING"
+//#endif
+
 #endif  // _libint2_src_bin_libint_shgshellordering_h_

include/libint2/solidharmonics.h

@@ -171,13 +171,19 @@ namespace libint2 {
           const unsigned short ncart = (l_ + 1) * (l_ + 2) / 2;
           std::vector<Real> full_coeff(npure * ncart);
 
-#if LIBINT_SHGSHELL_ORDERING == LIBINT_SHGSHELL_ORDERING_STANDARD
-          for(signed char pure_idx=0, m=-l_; pure_idx!=npure; ++pure_idx, ++m) {
-#elif LIBINT_SHGSHELL_ORDERING == LIBINT_SHGSHELL_ORDERING_GAUSSIAN
-          for(signed char pure_idx=0, m=0; pure_idx!=npure; ++pure_idx, m=(m>0?-m:1-m)) {
-#else
-#  error "unknown value of macro LIBINT_SHGSHELL_ORDERING"
-#endif
+          std::vector<int> shg_indices;
+          if (libint2::solid_harmonics_ordering() == libint2::SHGShellOrdering_Standard) {
+            for(signed char pure_idx=0, m=-l_; pure_idx!=npure; ++pure_idx, ++m)
+              shg_indices.push_back(m);
+          } else if (libint2::solid_harmonics_ordering() == libint2::SHGShellOrdering_Gaussian) {
+            for(signed char pure_idx=0, m=0; pure_idx!=npure; ++pure_idx, m=(m>0?-m:1-m))
+              shg_indices.push_back(m);
+          } else {
+            throw std::invalid_argument(std::string("libint2::solid_harmonics_ordering() value not recognized."));
+          }
+
+          for(signed char pure_idx=0; pure_idx!=npure; ++pure_idx) {
+            int m = shg_indices[pure_idx];
             signed char cart_idx = 0;
             signed char lx, ly, lz;
             FOR_CART(lx, ly, lz, l_)

python/src/libint2/libint2.cc

@@ -137,6 +137,12 @@ std::vector<double> coeffs_normalized(const Shell &s) {
 
 PYBIND11_MODULE(libint2, m) {
 
+  py::enum_<SHGShellOrdering>(m, "SHGShellOrdering")
+    .value("SHGShellOrdering_Standard", libint2::SHGShellOrdering_Standard)
+    .value("SHGShellOrdering_Gaussian", libint2::SHGShellOrdering_Gaussian)
+    .value("SHGShellOrdering_MOLDEN", libint2::SHGShellOrdering_Gaussian)
+    ;
+
   libint2::initialize();
 
   m.attr("MAX_AM") = LIBINT2_MAX_AM;

src/bin/libint/build_libint.cc

@@ -366,7 +366,13 @@ build_onebody_1b_1k(std::ostream& os, std::string label, const SafePtr<Compilati
             for(int l=0; l<=MULTIPOLE_MAX_ORDER; ++l) {
               // we iterate over them same way as over solid harmonic Gaussian shells
               int m;
-              FOR_SOLIDHARM(l,m)
+#if LIBINT_SHGSHELL_ORDERING == LIBINT_SHGSHELL_ORDERING_STANDARD
+              FOR_SOLIDHARM_STANDARD(l,m)
+#elif LIBINT_SHGSHELL_ORDERING == LIBINT_SHGSHELL_ORDERING_GAUSSIAN
+              FOR_SOLIDHARM_GAUSSIAN(l,m)
+#else
+#  error "unknown value of macro LIBINT_SHGSHELL_ORDERING"
+#endif
                 descrs.push_back(make_descr<OperDescrType>(l,m));
               END_FOR_SOLIDHARM
             }