Merge pull request #273 from evaleev/feature/σpVσp

1-e σpVσp

.github/workflows/cmake.yml

@@ -16,12 +16,15 @@ jobs:
         include:
           - os: ubuntu-20.04
             cxx: g++-10
+            cc: gcc-10
           - os: macos-latest
             cxx: clang++
+            cc: clang
 
     name: "Repo • ${{ matrix.os }}: ${{ matrix.cxx }} ${{ matrix.build_type }}"
     runs-on: ${{ matrix.os }}
     env:
+      CC: ${{ matrix.cc }}
       CXX : ${{ matrix.cxx }}
       CCACHE_DIR : ${{github.workspace}}/build/.ccache
       CCACHE_COMPRESS : true
@@ -185,6 +188,7 @@ jobs:
               -DBUILD_SHARED_LIBS=OFF
               -DCMAKE_CXX_COMPILER=clang-cl
               -DCMAKE_C_COMPILER=clang-cl
+              -DCMAKE_CXX_FLAGS="/source-charset:utf-8 /execution-charset:utf-8 /DUNICODE /D_UNICODE"
             testargs: >
               -GNinja
               -DCMAKE_BUILD_TYPE=Release
@@ -313,4 +317,3 @@ jobs:
         cmake . -D CMAKE_PREFIX_PATH="${{github.workspace}}/installed" ${{ matrix.cfg.testargs }}
         cmake --build . --target libint2-python
         cmake --build . --target libint2-python-test
-

export/Makefile

@@ -102,7 +102,5 @@ exportdir::
 	-$(INSTALL) $(INSTALLLIBOPT) $(SRCDIR)/cmake/hftest.cmake $(TOPDIR)/$(EXPORTDIR)/cmake/hftest.cmake
 	$(INSTALL) $(INSTALLDIROPT) $(TOPDIR)/$(EXPORTDIR)/cmake/modules
 	-$(INSTALL) $(INSTALLLIBOPT) $(SRCDIR)/cmake/modules/*.cmake $(TOPDIR)/$(EXPORTDIR)/cmake/modules
-	echo "set(LIBINT2_LIBRARY_CXX_SRC" > $(TOPDIR)/$(EXPORTDIR)/srclist.cmake
 	ls $(TOPDIR)/$(EXPORTDIR)/src >> $(TOPDIR)/$(EXPORTDIR)/srclist.cmake
-	echo ")" >> $(TOPDIR)/$(EXPORTDIR)/srclist.cmake
 	cp -Rap $(SRCTOPDIR)/python $(TOPDIR)/$(EXPORTDIR)/

export/cmake/CMakeLists.txt.export

@@ -188,10 +188,24 @@ set(LIBINT2_INSTALL_CMAKEDIR "lib/cmake/libint2"
 
 # Libint library =======================================================================================================
 
-include(srclist.cmake)
+file(STRINGS "${CMAKE_CURRENT_SOURCE_DIR}/srclist.cmake" LIBINT2_LIBRARY_CXX_SRC ENCODING UTF-8)
 set(LIB_CXX_SRC )
 foreach(FN IN LISTS LIBINT2_LIBRARY_CXX_SRC)
-    list(APPEND LIB_CXX_SRC "src/${FN}")
+    set(srcfile "${CMAKE_CURRENT_SOURCE_DIR}/src/${FN}")
+    if (MSVC)
+         string(REGEX REPLACE "[a-zA-Z0-9_.]" "" NONASCII_FN "${FN}")
+         if (NONASCII_FN)
+             message(WARNING "Non-ASCII characters in filename ${FN}")
+             string(MAKE_C_IDENTIFIER "${FN}" C_FN)
+             if (FN STREQUAL C_FN)
+                 message(FATAL_ERROR "Non-ASCII characters in filename ${FN} but C identifier is the same")
+             endif()
+             file(CREATE_LINK "${CMAKE_CURRENT_SOURCE_DIR}/src/${FN}" "${CMAKE_CURRENT_SOURCE_DIR}/src/${C_FN}" SYMBOLIC)
+             message(WARNING "Created symbolic link ${CMAKE_CURRENT_SOURCE_DIR}/src/${C_FN} -> ${CMAKE_CURRENT_SOURCE_DIR}/src/${FN}")
+             set(srcfile "${CMAKE_CURRENT_SOURCE_DIR}/src/${C_FN}")
+         endif()
+    endif(MSVC)
+    list(APPEND LIB_CXX_SRC "${srcfile}")
 endforeach()
 # Create object files to use for static and shared libraries
 add_library(libint2_obj OBJECT ${LIB_CXX_SRC} "src/configuration.cc")
@@ -203,7 +217,7 @@ target_compile_features(libint2_obj PUBLIC "cxx_std_11")
 set_target_properties(
   libint2_obj
   PROPERTIES
-    UNITY_BUILD TRUE
+    UNITY_BUILD FALSE
   )
 if (TARGET MPFR::GMPXX)
   target_link_libraries(libint2_obj PUBLIC MPFR::GMPXX)

include/libint2/engine.h

@@ -128,6 +128,8 @@ enum class Operator {
   //! \f$ \mathcal{N}^+_{0,0} , \mathcal{N}^-_{1,1}, \mathcal{N}^+_{1,0}, \mathcal{N}^+_{1,1}, \mathcal{N}^-_{2,2}, \mathcal{N}^-_{2,1}, \mathcal{N}^+_{2,0}, \mathcal{N}^+_{2,1}, \mathcal{N}^+_{2,2}. \dots \f$ .
   //! Previous to cdbb9f3 released in v2.8.0, Standard -or- Gaussian ordering could be be specified at generator/compiler configure time.
   sphemultipole,
+  /// The four components of σp . V . σp, where V is the nuclear potential.
+  opVop,
   /// \f$ \delta(\vec{r}_1 - \vec{r}_2) \f$
   delta,
   /// (2-body) Coulomb operator = \f$ r_{12}^{-1} \f$
@@ -160,7 +162,7 @@ enum class Operator {
   invalid = -1,
   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!keep this updated!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   first_1body_oper = overlap,
-  last_1body_oper = sphemultipole,
+  last_1body_oper = opVop,
   first_2body_oper = delta,
   last_2body_oper = stg_x_coulomb,
   first_oper = first_1body_oper,
@@ -185,6 +187,9 @@ struct default_operator_traits {
   } oper_params_type;
   static oper_params_type default_params() { return oper_params_type{}; }
   static constexpr auto nopers = 1u;
+  // N.B.: Below field means we *should* template specialize operator_traits for
+  // Operator::kinetic, but L2 doesn't use that anywhere.
+  static constexpr auto intrinsic_deriv_order = 0u;
   struct _core_eval_type {
     template <typename... params>
     static std::shared_ptr<const _core_eval_type> instance(params...) {
@@ -216,6 +221,12 @@ struct operator_traits<Operator::nuclear>
   typedef const libint2::FmEval_Reference<scalar_type> core_eval_type;
 #endif
 };
+template <>
+struct operator_traits<Operator::opVop>
+    : public operator_traits<Operator::nuclear> {
+  static constexpr auto nopers = 4;
+  static constexpr auto intrinsic_deriv_order = 2;
+};
 
 template <>
 struct operator_traits<Operator::erf_nuclear>
@@ -993,6 +1004,7 @@ class Engine {
   //-------
   unsigned int nparams() const;
   unsigned int nopers() const;
+  unsigned int intrinsic_deriv_order() const;
   /// if Params == operator_traits<oper>::oper_params_type, will return
   /// \c any(params)
   /// else will set return \c any initialized with default value for

include/libint2/engine.impl.h

@@ -80,14 +80,15 @@ typename std::remove_all_extents<T>::type* to_ptr1(T (&a)[N]) {
         (2emultipole,                \
          (3emultipole,               \
            (sphemultipole,           \
-          (eri, (eri, (eri, (eri, (eri, (eri, (eri, (eri, (eri, (eri, BOOST_PP_NIL)))))))))))))))))))
+            (opVop,                  \
+          (eri, (eri, (eri, (eri, (eri, (eri, (eri, (eri, (eri, (eri, BOOST_PP_NIL))))))))))))))))))))
 
 #define BOOST_PP_NBODY_OPERATOR_INDEX_TUPLE \
   BOOST_PP_MAKE_TUPLE(BOOST_PP_LIST_SIZE(BOOST_PP_NBODY_OPERATOR_LIST))
 #define BOOST_PP_NBODY_OPERATOR_INDEX_LIST \
   BOOST_PP_TUPLE_TO_LIST(BOOST_PP_NBODY_OPERATOR_INDEX_TUPLE)
 #define BOOST_PP_NBODY_OPERATOR_LAST_ONEBODY_INDEX \
-  8  // sphemultipole, the 9th member of BOOST_PP_NBODY_OPERATOR_LIST, is the last
+  9  // opVop, the 10th member of BOOST_PP_NBODY_OPERATOR_LIST, is the last
      // 1-body operator
 
 // make list of braket indices for n-body ints
@@ -179,7 +180,7 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute1(
 
   const auto oper_is_nuclear =
       (oper_ == Operator::nuclear || oper_ == Operator::erf_nuclear ||
-       oper_ == Operator::erfc_nuclear);
+       oper_ == Operator::erfc_nuclear || oper_ == Operator::opVop);
 
   const auto l1 = s1.contr[0].l;
   const auto l2 = s2.contr[0].l;
@@ -190,7 +191,7 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute1(
   // user likely forgot to call set_params
   if (oper_is_nuclear && nparams() == 0)
     throw std::logic_error(
-        "Engine<*nuclear>, but no charges found; forgot to call "
+        "Engine requires charges, but no charges found; forgot to call "
         "set_params()?");
 
   const auto n1 = s1.size();
@@ -247,7 +248,8 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute1(
   // element of stack
   const auto compute_directly =
       lmax == 0 && deriv_order_ == 0 &&
-      (oper_ == Operator::overlap || oper_is_nuclear);
+      (oper_ == Operator::overlap || oper_is_nuclear) &&
+      oper_ != Operator::opVop;
   if (compute_directly) {
     primdata_[0].stack[0] = 0;
     targets_[0] = primdata_[0].stack;
@@ -748,6 +750,7 @@ Engine::compute2_ptrs() const {
 __libint2_engine_inline unsigned int Engine::nparams() const {
   switch (oper_) {
     case Operator::nuclear:
+    case Operator::opVop:
       return any_cast<const operator_traits<Operator::nuclear>::oper_params_type&>(params_)
           .size();
     case Operator::erf_nuclear:
@@ -772,6 +775,19 @@ __libint2_engine_inline unsigned int Engine::nopers() const {
   assert(false && "missing case in switch");  // unreachable
   abort();
 }
+__libint2_engine_inline unsigned int Engine::intrinsic_deriv_order() const {
+  switch (static_cast<int>(oper_)) {
+#define BOOST_PP_NBODYENGINE_MCR9(r, data, i, elem) \
+  case i:                                           \
+    return operator_traits<static_cast<Operator>(i)>::intrinsic_deriv_order;
+    BOOST_PP_LIST_FOR_EACH_I(BOOST_PP_NBODYENGINE_MCR9, _,
+                             BOOST_PP_NBODY_OPERATOR_LIST)
+    default:
+      break;
+  }
+  assert(false && "missing case in switch");  // unreachable
+  abort();
+}
 
 template <>
 __libint2_engine_inline any Engine::enforce_params_type<any>(
@@ -825,20 +841,26 @@ __libint2_engine_inline any Engine::enforce_params_type(
   return result;
 }
 
+/// @param[in] oper the Operator for which to return core-evaluator objects
+/// @return a core-integral evaluator and its scratch for @c oper
+/// @throw @c oper not registered in BOOST_PP_NBODY_OPERATOR_LIST
+/// @throw @c oper has no core_eval_type specialization in operator_traits.
 __libint2_engine_inline any Engine::make_core_eval_pack(Operator oper) const {
   any result;
   switch (static_cast<int>(oper)) {
-#define BOOST_PP_NBODYENGINE_MCR6(r, data, i, elem)                          \
-  case i:                                                                    \
-    result = libint2::detail::make_compressed_pair(                          \
-        operator_traits<static_cast<Operator>(i)>::core_eval_type::instance( \
-            braket_rank() * lmax_ + deriv_order_,                            \
-            std::numeric_limits<scalar_type>::epsilon()),                    \
-        libint2::detail::CoreEvalScratch<                                    \
-            operator_traits<static_cast<Operator>(i)>::core_eval_type>(      \
-            braket_rank() * lmax_ + deriv_order_));                          \
-    assert(any_cast<detail::core_eval_pack_type<static_cast<Operator>(i)>>(  \
-               &result) != nullptr);                                         \
+#define BOOST_PP_NBODYENGINE_MCR6(r, data, i, elem)                             \
+  case i:                                                                       \
+    result = libint2::detail::make_compressed_pair(                             \
+        operator_traits<static_cast<Operator>(i)>::core_eval_type::instance(    \
+            braket_rank() * lmax_ + deriv_order_ +                              \
+            operator_traits<static_cast<Operator>(i)>::intrinsic_deriv_order,   \
+            std::numeric_limits<scalar_type>::epsilon()),                       \
+        libint2::detail::CoreEvalScratch<                                       \
+            operator_traits<static_cast<Operator>(i)>::core_eval_type>(         \
+            braket_rank() * lmax_ + deriv_order_ +                              \
+            operator_traits<static_cast<Operator>(i)>::intrinsic_deriv_order)); \
+    assert(any_cast<detail::core_eval_pack_type<static_cast<Operator>(i)>>(     \
+               &result) != nullptr);                                            \
     break;
 
     BOOST_PP_LIST_FOR_EACH_I(BOOST_PP_NBODYENGINE_MCR6, _,
@@ -907,7 +929,7 @@ __libint2_engine_inline void Engine::compute_primdata(Libint_t& primdata, const
 
   const auto oper_is_nuclear =
       (oper_ == Operator::nuclear || oper_ == Operator::erf_nuclear ||
-       oper_ == Operator::erfc_nuclear);
+       oper_ == Operator::erfc_nuclear || oper_ == Operator::opVop);
 
   // need to use HRR? see strategy.cc
   const auto l1 = s1.contr[0].l;
@@ -1010,7 +1032,7 @@ __libint2_engine_inline void Engine::compute_primdata(Libint_t& primdata, const
   primdata._0_Overlap_0_y[0] = ovlp_ss_y;
   primdata._0_Overlap_0_z[0] = ovlp_ss_z;
 
-  if (oper_ == Operator::kinetic || (deriv_order_ > 0)) {
+  if (oper_ == Operator::kinetic || (deriv_order_ > 0) || oper_ == Operator::opVop) {
 #if LIBINT2_DEFINED(eri, two_alpha0_bra)
     primdata.two_alpha0_bra[0] = 2.0 * alpha1;
 #endif
@@ -1021,7 +1043,7 @@ __libint2_engine_inline void Engine::compute_primdata(Libint_t& primdata, const
 
   if (oper_is_nuclear) {
 
-    const auto& params = (oper_ == Operator::nuclear) ?
+    const auto& params = (oper_ == Operator::nuclear || oper_ == Operator::opVop) ?
         any_cast<const operator_traits<Operator::nuclear>::oper_params_type&>(params_) :
         std::get<1>(any_cast<const operator_traits<Operator::erfc_nuclear>::oper_params_type&>(params_));
 
@@ -1039,7 +1061,7 @@ __libint2_engine_inline void Engine::compute_primdata(Libint_t& primdata, const
 
 #if LIBINT2_DEFINED(eri, rho12_over_alpha1) || \
     LIBINT2_DEFINED(eri, rho12_over_alpha2)
-    if (deriv_order_ > 0) {
+    if (deriv_order_ + intrinsic_deriv_order() > 0) {
 #if LIBINT2_DEFINED(eri, rho12_over_alpha1)
       primdata.rho12_over_alpha1[0] = rhop_over_alpha1;
 #endif
@@ -1054,9 +1076,9 @@ __libint2_engine_inline void Engine::compute_primdata(Libint_t& primdata, const
                      primdata.PC_y[0] * primdata.PC_y[0] +
                      primdata.PC_z[0] * primdata.PC_z[0];
     const scalar_type U = gammap * PC2;
-    const auto mmax = s1.contr[0].l + s2.contr[0].l + deriv_order_;
+    const auto mmax = s1.contr[0].l + s2.contr[0].l + deriv_order_ + intrinsic_deriv_order();
     auto* fm_ptr = &(primdata.LIBINT_T_S_ELECPOT_S(0)[0]);
-    if (oper_ == Operator::nuclear) {
+    if (oper_ == Operator::nuclear || oper_ == Operator::opVop) {
       const auto& fm_engine_ptr =
           any_cast<const detail::core_eval_pack_type<Operator::nuclear>&>(core_eval_pack_)
           .first();

src/bin/libint/build_libint.cc

@@ -244,7 +244,7 @@ namespace {
     typedef EmptySet type;
   };
 
-  template <typename OperDescrType> OperDescrType make_descr(int, int, int = 0) {
+  template <typename OperDescrType> OperDescrType make_descr(int, int = 0, int = 0) {
     return OperDescrType();
   }
   template <> CartesianMultipole_Descr<3u> make_descr<CartesianMultipole_Descr<3u>>(int x, int y, int z) {
@@ -259,6 +259,9 @@ namespace {
     SphericalMultipole_Descr result(l,m);
     return result;
   }
+  template <> σpVσp_Descr make_descr<σpVσp_Descr>(int p, int, int) {
+    return σpVσp_Descr(p);
+  }
 }
 
 template <typename _OperType>
@@ -371,6 +374,14 @@ build_onebody_1b_1k(std::ostream& os, std::string label, const SafePtr<Compilati
               END_FOR_SOLIDHARM
             }
           }
+          if (std::is_same<_OperType,σpVσpOper>::value) {
+            // reset descriptors array
+            descrs.resize(0);
+            // iterate over pauli components
+            for (int p = 0; p != 4; ++p) {
+              descrs.emplace_back(make_descr<OperDescrType>(p));
+            }
+          }
 
           // derivative index is the outermost (slowest running)
           // operator component is second slowest
@@ -384,7 +395,7 @@ build_onebody_1b_1k(std::ostream& os, std::string label, const SafePtr<Compilati
           do {
             BFType a(la);
             BFType b(lb);
-            
+
             for(unsigned int c=0; c!=2; ++c) {
               const unsigned int ndir = std::is_same<BFType,CGShell>::value ? 3 : 1;
               for(unsigned int xyz=0; xyz<ndir; ++xyz) {
@@ -400,7 +411,7 @@ build_onebody_1b_1k(std::ostream& os, std::string label, const SafePtr<Compilati
               SafePtr<Onebody_sh_1_1> target = Onebody_sh_1_1::Instance(a,b,nullaux,oper);
               targets.push_back(target);
             } // loop over operator components
-            
+
             last_deriv = diter.last();
             if (!last_deriv) diter.next();
           } while (!last_deriv); // loop over derivatives
@@ -506,15 +517,17 @@ void try_main (int argc, char* argv[])
                                      1emultipole,           \
                                      2emultipole,           \
                                      3emultipole,           \
-                                     sphemultipole    \
+                                     sphemultipole,         \
+                                     opVop                  \
                                     )
 #define BOOST_PP_ONEBODY_TASK_OPER_TUPLE (OverlapOper,                    \
                                           KineticOper,                    \
                                           ElecPotOper,                    \
                                           CartesianMultipoleOper<3u>,     \
                                           CartesianMultipoleOper<3u>,     \
                                           CartesianMultipoleOper<3u>,     \
-                                          SphericalMultipoleOper          \
+                                          SphericalMultipoleOper,         \
+                                          σpVσpOper                       \
                                          )
 #define BOOST_PP_ONEBODY_TASK_LIST BOOST_PP_TUPLE_TO_LIST( BOOST_PP_ONEBODY_TASK_TUPLE )
 #define BOOST_PP_ONEBODY_TASK_OPER_LIST BOOST_PP_TUPLE_TO_LIST( BOOST_PP_ONEBODY_TASK_OPER_TUPLE )