diff --git a/CMakeLists.txt b/CMakeLists.txt index f78341d28..244259c26 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -109,64 +109,67 @@ option_with_default(CMAKE_BUILD_TYPE "Build type (Release or Debug)" Release) ### compiler-only option_with_print(LIBINT2_BUILD_LIBRARY_AS_SUBPROJECT - "[EXPERT] Build generated library as a subproject: if FALSE will configure and build separately" OFF) + "[EXPERT] Build generated library as a subproject: if FALSE will configure and build separately" OFF) ### library-only option_with_print(LIBINT2_REQUIRE_CXX_API - "C++11 Libint API: define library targets + test (requires Eigen3, Boost is optional but strongly recommended)" ON) + "C++11 Libint API: define library targets + test (requires Eigen3, Boost is optional but strongly recommended)" ON) option_with_print(LIBINT2_REQUIRE_CXX_API_COMPILED - "Build C++11 Compiled (not just header-only) targets (requires Eigen3, Boost strongly recommended)" ON) + "Build C++11 Compiled (not just header-only) targets (requires Eigen3, Boost strongly recommended)" ON) option_with_print(LIBINT2_ENABLE_FORTRAN - "Build Fortran03+ Libint interface (requires Fortran)" OFF) + "Build Fortran03+ Libint interface (requires Fortran)" OFF) option_with_print(LIBINT2_ENABLE_PYTHON - "Build Python bindings (requires Python and Pybind11 and Eigen3)" OFF) + "Build Python bindings (requires Python and Pybind11 and Eigen3)" OFF) option_with_print(LIBINT2_PREFIX_PYTHON_INSTALL - "For LIBINT2_ENABLE_PYTHON=ON, whether to install the Python module in the Linux manner to CMAKE_INSTALL_PREFIX or to not install it. See target libint2-python-wheel for alternate installation in the Python manner to Python_EXECUTABLE's site-packages." OFF) + "For LIBINT2_ENABLE_PYTHON=ON, whether to install the Python module in the Linux manner to CMAKE_INSTALL_PREFIX or to not install it. See target libint2-python-wheel for alternate installation in the Python manner to Python_EXECUTABLE's site-packages." OFF) option_with_print(BUILD_SHARED_LIBS - "Build Libint library as shared, not static" OFF) + "Build Libint library as shared, not static" OFF) option_with_print(LIBINT2_BUILD_SHARED_AND_STATIC_LIBS - "Build both shared and static Libint libraries in one shot. Uses -fPIC." OFF) + "Build both shared and static Libint libraries in one shot. Uses -fPIC." OFF) option_with_print(LIBINT2_ENABLE_MPFR - "Use GNU MPFR library for high-precision testing (EXPERTS ONLY). Consumed at library build-time." OFF) + "Use GNU MPFR library for high-precision testing (EXPERTS ONLY). Consumed at library build-time." OFF) # <<< Which Integrals Classes, Which Derivative Levels >>> option_with_default(LIBINT2_ENABLE_ONEBODY - "Compile with support for up to N-th derivatives of 1-body integrals (-1 for OFF)" 0) + "Compile with support for up to N-th derivatives of 1-body integrals (-1 for OFF)" 0) option_with_default(LIBINT2_ENABLE_ERI - "Compile with support for up to N-th derivatives of 4-center electron repulsion integrals (-1 for OFF)" 0) + "Compile with support for up to N-th derivatives of 4-center electron repulsion integrals (-1 for OFF)" 0) option_with_default(LIBINT2_ENABLE_ERI3 - "Compile with support for up to N-th derivatives of 3-center electron repulsion integrals (-1 for OFF)" -1) + "Compile with support for up to N-th derivatives of 3-center electron repulsion integrals (-1 for OFF)" -1) option_with_default(LIBINT2_ENABLE_ERI2 - "Compile with support for up to N-th derivatives of 2-center electron repulsion integrals (-1 for OFF)" -1) + "Compile with support for up to N-th derivatives of 2-center electron repulsion integrals (-1 for OFF)" -1) +option_with_default(LIBINT2_ENABLE_RKB_ERI + "Compile with support for up to N-th derivatives of relativistic restricted kinetic + balance (RKB) 4-center electron repulsion integrals (-1 for OFF)" 0) option_with_default(LIBINT2_ENABLE_G12 - "Compile with support for N-th derivatives of MP2-F12 energies with Gaussian factors (-1 for OFF)" -1) + "Compile with support for N-th derivatives of MP2-F12 energies with Gaussian factors (-1 for OFF)" -1) option_with_default(LIBINT2_ENABLE_G12DKH - "Compile with support for N-th derivatives of DKH-MP2-F12 energies with Gaussian factors (-1 for OFF)" -1) + "Compile with support for N-th derivatives of DKH-MP2-F12 energies with Gaussian factors (-1 for OFF)" -1) option_with_print(LIBINT2_DISABLE_ONEBODY_PROPERTY_DERIVS - "Disable geometric derivatives of 1-body property integrals (all but overlap, kinetic, elecpot). + "Disable geometric derivatives of 1-body property integrals (all but overlap, kinetic, elecpot). These derivatives are disabled by default to save compile time. (enable with OFF) Note that the libtool build won't enable this- if forcibly enabled, build_libint balks." ON) option_with_print(LIBINT2_ENABLE_T1G12 - "Enable [Ti,G12] integrals when G12 integrals are enabled. Irrelevant when `LIBINT2_ENABLE_G12=OFF`. (disable with OFF)" ON) + "Enable [Ti,G12] integrals when G12 integrals are enabled. Irrelevant when `LIBINT2_ENABLE_G12=OFF`. (disable with OFF)" ON) # <<< Ordering Conventions >>> option_with_default(LIBINT2_SHGAUSS_ORDERING - "Ordering for shells of solid harmonic Gaussians: + "Ordering for shells of solid harmonic Gaussians: standard -- standard ordering (-l, -l+1 ... l) gaussian -- the Gaussian ordering (0, 1, -1, 2, -2, ... l, -l) See https://github.com/evaleev/libint/blob/master/INSTALL.md#solid-harmonic-ordering-scope-and-history ." standard) option_with_default(LIBINT2_CARTGAUSS_ORDERING - "Orderings for shells of cartesian Gaussians: + "Orderings for shells of cartesian Gaussians: standard -- standard ordering (xxx, xxy, xxz, xyy, xyz, xzz, yyy, ...) intv3 -- intv3 ordering (yyy, yyz, yzz, zzz, xyy, xyz, xzz, xxy, xxz, xxx) gamess -- GAMESS ordering (xxx, yyy, zzz, xxy, xxz, yyx, yyz, zzx, zzy, xyz) orca -- ORCA ordering (hydrid between GAMESS and standard) bagel -- axis-permuted version of intv3 (xxx, xxy, xyy, yyy, xxz, xyz, yyz, xzz, yzz, zzz)" standard) option_with_default(LIBINT2_SHELL_SET - "Support computation of shell sets sets subject to these restrictions: + "Support computation of shell sets sets subject to these restrictions: standard -- standard ordering: for (ab|cd): l(a) >= l(b), @@ -195,99 +198,107 @@ option_with_default(LIBINT2_SHELL_SET # `export CMAKE_BUILD_PARALLEL_LEVEL=N`. option_with_default(LIBINT2_MAX_AM - "Support Gaussians of angular momentum up to N. + "Support Gaussians of angular momentum up to N. Can specify values for each derivative level as a semicolon-separated string. If ERI3 ints are enabled, this option also controls the AM of the paired centers." 4) option_with_default(LIBINT2_OPT_AM - "Optimize maximally for up to angular momentum N (N <= max-am). + "Optimize maximally for up to angular momentum N (N <= max-am). Can specify values for each derivative level as a semicolon-separated string. (default: (libint_max_am/2)+1)" -1) option_with_default(LIBINT2_MULTIPOLE_MAX_ORDER - "Maximum order of spherical multipole integrals. There is no maximum" 4) + "Maximum order of spherical multipole integrals. There is no maximum" 4) option_with_default(LIBINT2_ONEBODY_MAX_AM - "Support 1-body ints for Gaussians of angular momentum up to N. + "Support 1-body ints for Gaussians of angular momentum up to N. Can specify values for each derivative level as a semicolon-separated string. (default: max_am)" -1) option_with_default(LIBINT2_ONEBODY_OPT_AM - "Optimize 1-body ints maximally for up to angular momentum N (N <= max-am). + "Optimize 1-body ints maximally for up to angular momentum N (N <= max-am). Can specify values for each derivative level as a semicolon-separated string (default: (max_am/2)+1)" -1) option_with_default(LIBINT2_ERI_MAX_AM - "Support 4-center ERIs for Gaussians of angular momentum up to N. + "Support 4-center ERIs for Gaussians of angular momentum up to N. Can specify values for each derivative level as a semicolon-separated string. (default: max_am)" -1) option_with_default(LIBINT2_ERI_OPT_AM - "Optimize 4-center ERIs maximally for up to angular momentum N (N <= max-am). + "Optimize 4-center ERIs maximally for up to angular momentum N (N <= max-am). Can specify values for each derivative level as a semicolon-separated string (default: (max_am/2)+1)" -1) +option_with_default(LIBINT2_RKB_ERI_MAX_AM + "Support relativistic restricted kinetic balance (RKB) 4-center ERIs for Gaussians of angular momentum up to N. + Can specify values for each derivative level as a semicolon-separated string. (default: max_am)" -1) +option_with_default(LIBINT2_RKB_ERI_OPT_AM + "Optimize relativistic restricted kinetic balance (RKB) 4-center ERIs maximally for up to angular momentum N (N <= max-am). + Can specify values for each derivative level as a semicolon-separated string (default: (max_am/2)+1)" -1) + + option_with_default(LIBINT2_ERI3_MAX_AM - "Support 3-center ERIs for Gaussians of angular momentum up to N. + "Support 3-center ERIs for Gaussians of angular momentum up to N. Can specify values for each derivative level as a semicolon-separated string. (default: max_am) This option controls only the single fitting center. The paired centers use LIBINT2_MAX_AM." -1) option_with_default(LIBINT2_ERI3_OPT_AM - "Optimize 3-center ERIs maximally for up to angular momentum N (N <= max-am). + "Optimize 3-center ERIs maximally for up to angular momentum N (N <= max-am). Can specify values for each derivative level as a semicolon-separated string. (default: (max_am/2)+1)" -1) option_with_print(LIBINT2_ERI3_PURE_SH - "Assume the 'unpaired' center of 3-center ERIs will be transformed to pure solid harmonics" OFF) + "Assume the 'unpaired' center of 3-center ERIs will be transformed to pure solid harmonics" OFF) option_with_default(LIBINT2_ERI2_MAX_AM - "Support 2-center ERIs for Gaussians of angular momentum up to N. + "Support 2-center ERIs for Gaussians of angular momentum up to N. Can specify values for each derivative level as a semicolon-separated string. (default: max_am)" -1) option_with_default(LIBINT2_ERI2_OPT_AM - "Optimize 2-center ERIs maximally for up to angular momentum N (N <= max-am). + "Optimize 2-center ERIs maximally for up to angular momentum N (N <= max-am). Can specify values for each derivative level as a semicolon-separated string. (default: (max_am/2)+1)" -1) option_with_print(LIBINT2_ERI2_PURE_SH - "Assume the 2-center ERIs will be transformed to pure solid harmonics" OFF) + "Assume the 2-center ERIs will be transformed to pure solid harmonics" OFF) option_with_default(LIBINT2_G12_MAX_AM - "Support integrals for G12 methods of angular momentum up to N. (default: max_am)" -1) + "Support integrals for G12 methods of angular momentum up to N. (default: max_am)" -1) option_with_default(LIBINT2_G12_OPT_AM - "Optimize G12 integrals for up to angular momentum N (N <= max-am). (default: (max_am/2)+1)" -1) + "Optimize G12 integrals for up to angular momentum N (N <= max-am). (default: (max_am/2)+1)" -1) option_with_default(LIBINT2_G12DKH_MAX_AM - "Support integrals for relativistic G12 methods of angular momentum up to N. (default: max_am)" -1) + "Support integrals for relativistic G12 methods of angular momentum up to N. (default: max_am)" -1) option_with_default(LIBINT2_G12DKH_OPT_AM - "Optimize G12DKH integrals for up to angular momentum N (N <= max-am). (default: (max_am/2)+1)" -1) + "Optimize G12DKH integrals for up to angular momentum N (N <= max-am). (default: (max_am/2)+1)" -1) # <<< Miscellaneous >>> option_with_print(LIBINT2_CONTRACTED_INTS - "Turn on support for contracted integrals." ON) + "Turn on support for contracted integrals." ON) option_with_default(LIBINT2_ERI_STRATEGY - "(EXPERT) Compute ERIs using the following strategy. (0 for OS, 1 for HGP, 2 for HL)" 1) + "(EXPERT) Compute ERIs using the following strategy. (0 for OS, 1 for HGP, 2 for HL)" 1) option_with_print(LIBINT2_USE_COMPOSITE_EVALUATORS - "Libint will use composite evaluators (i.e. every evaluator will compute one integral type only)" ON) + "Libint will use composite evaluators (i.e. every evaluator will compute one integral type only)" ON) option_with_print(LIBINT2_SINGLE_EVALTYPE - "Generate single evaluator type (i.e. all tasks use the same evaluator). OFF is NYI" ON) + "Generate single evaluator type (i.e. all tasks use the same evaluator). OFF is NYI" ON) option_with_default(LIBINT2_ENABLE_UNROLLING - "Unroll shell sets into integrals (will unroll shell sets larger than N) (0 for never, N for N, 1000000000 for always)" 100) + "Unroll shell sets into integrals (will unroll shell sets larger than N) (0 for never, N for N, 1000000000 for always)" 100) option_with_default(LIBINT2_ALIGN_SIZE - "(EXPERT) if posix_memalign is available, this will specify alignment of Libint data, in units of + "(EXPERT) if posix_memalign is available, this will specify alignment of Libint data, in units of sizeof(LIBINT2_REALTYPE). Default is to use built-in heuristics: system-determined for vectorization off (default) or veclen * sizeof(LIBINT2_REALTYPE) for vectorization on." 0) mark_as_advanced(LIBINT2_ALIGN_SIZE) option_with_default(LIBINT2_REALTYPE - "Specifies the floating-point data type used by the library. Consumed at library build-time." double) + "Specifies the floating-point data type used by the library. Consumed at library build-time." double) option_with_print(LIBINT2_USER_DEFINED_REAL_INCLUDES - "Additional #includes necessary to use the real type." OFF) + "Additional #includes necessary to use the real type." OFF) include(int_userreal) option_with_print(LIBINT2_GENERATE_FMA - "Generate FMA (fused multiply-add) instructions (to benefit must have FMA-capable hardware and compiler)" OFF) + "Generate FMA (fused multiply-add) instructions (to benefit must have FMA-capable hardware and compiler)" OFF) option_with_print(LIBINT2_ENABLE_GENERIC_CODE "Use manually-written generic code" OFF) option_with_default(LIBINT2_API_PREFIX "Prepend this string to every name in the library API (except for the types)." "") option_with_print(LIBINT2_VECTOR_LENGTH - "Compute integrals in vectors of length N." OFF) + "Compute integrals in vectors of length N." OFF) option_with_default(LIBINT2_VECTOR_METHOD - "Specifies how to vectorize integrals. Irrelevant when `LIBINT2_VECTOR_LENGTH=OFF. Allowed values are 'block' (default), and 'line'." block) + "Specifies how to vectorize integrals. Irrelevant when `LIBINT2_VECTOR_LENGTH=OFF. Allowed values are 'block' (default), and 'line'." block) option_with_print(LIBINT2_ACCUM_INTS - "Accumulate integrals to the buffer, rather than copy (OFF for copy, ON for accum)." OFF) + "Accumulate integrals to the buffer, rather than copy (OFF for copy, ON for accum)." OFF) option_with_print(LIBINT2_FLOP_COUNT - "Support (approximate) FLOP counting by the library. (Generated code will require C++11!)" OFF) + "Support (approximate) FLOP counting by the library. (Generated code will require C++11!)" OFF) option_with_print(LIBINT2_PROFILE - "Turn on profiling instrumentation of the library. (Generated code will require C++11!)" OFF) + "Turn on profiling instrumentation of the library. (Generated code will require C++11!)" OFF) option_with_print(LIBINT2_ENABLE_MPFR - "Use GNU MPFR library for high-precision testing (EXPERTS ONLY). Consumed at library build-time." OFF) + "Use GNU MPFR library for high-precision testing (EXPERTS ONLY). Consumed at library build-time." OFF) option_with_default(LIBINT2_EXPORT_COMPRESSOR - "Export tarball with compression gzip or bzip2" gzip) + "Export tarball with compression gzip or bzip2" gzip) # next one defined by `include(CTest)` message(STATUS "Showing option BUILD_TESTING: ${BUILD_TESTING}") @@ -304,13 +315,13 @@ include(int_am) check_function_exists(posix_memalign HAVE_POSIX_MEMALIGN) if (NOT HAVE_POSIX_MEMALIGN) message(FATAL_ERROR "did not find posix_memalign ... this SHOULD NOT happen. Cannot proceed.") -endif() +endif () check_include_file_cxx(stdint.h HAVE_STDINT_H) # limits.h? if (cxx_std_11 IN_LIST CMAKE_CXX_COMPILE_FEATURES) set(LIBINT_HAS_CXX11 1) -endif() +endif () booleanize01(LIBINT2_ERI3_PURE_SH) booleanize01(LIBINT2_ERI2_PURE_SH) @@ -331,9 +342,9 @@ if (LIBINT2_EXPORT_COMPRESSOR STREQUAL "gzip") elseif (LIBINT2_EXPORT_COMPRESSOR STREQUAL "bzip2") set(LIBINT_EXPORT_COMPRESSOR_CMD "jcf") set(LIBINT_EXPORT_COMPRESSOR_EXT "tbz2") -else() +else () message(FATAL_ERROR "No valid compressor; invoke CMake with -DLIBINT2_EXPORT_COMPRESSOR=gzip|bzip2") -endif() +endif () ################################## Dependencies ################################# @@ -344,9 +355,9 @@ if (LIBINT2_ENABLE_MPFR) # mpfr detected in CMakeLists.txt.export at appropriate time for library, but prechecking here find_package(Multiprecision MODULE REQUIRED COMPONENTS gmpxx mpfr) set(LIBINT_HAS_MPFR 1) -else() +else () find_package(Multiprecision MODULE REQUIRED COMPONENTS gmpxx) -endif() +endif () get_property(_loc TARGET Multiprecision::gmp PROPERTY LOCATION) message(VERBOSE "${Cyan}Found GMP${ColourReset}: ${_loc}") @@ -355,12 +366,12 @@ message(VERBOSE "${Cyan}Found GMPXX${ColourReset}: ${_loc}") if (TARGET Multiprecision::mpfr) get_property(_loc TARGET Multiprecision::mpfr PROPERTY LOCATION) message(VERBOSE "${Cyan}Found MPFR${ColourReset}: ${_loc} (found version ${MPFR_VERSION})") -endif() +endif () find_package(Boost 1.57 REQUIRED) if (TARGET Boost::headers) set(LIBINT_HAS_SYSTEM_BOOST_PREPROCESSOR_VARIADICS 1) -endif() +endif () # deferring find_package(Eigen3) to library (CMakeLists.txt.export) @@ -370,9 +381,9 @@ endif() set(EXPORT_STAGE_DIR ${PROJECT_BINARY_DIR}/libint-${LIBINT_EXT_VERSION}) configure_file( - cmake/modules/int_computed.cmake.in - cmake/modules/int_computed.cmake - @ONLY) + cmake/modules/int_computed.cmake.in + cmake/modules/int_computed.cmake + @ONLY) # CMake data transmitted to C++ via config.h for generator/compiler (_EXPORT_MODE=0). # Same info is positioned for the library export, but _EXPORT_MODE=1 turns on @@ -380,28 +391,28 @@ configure_file( # library build time. set(_EXPORT_MODE 0) # convert user-facing LIBINT2_ variables to LIBINT_ internal variables -foreach(_var API_PREFIX;ERI3_PURE_SH;ERI2_PURE_SH;DISABLE_ONEBODY_PROPERTY_DERIVS;ENABLE_UNROLLING;ENABLE_GENERIC_CODE;VECTOR_LENGTH;VECTOR_METHOD;ALIGN_SIZE;USER_DEFINED_REAL;USER_DEFINED_REAL_INCLUDES;GENERATE_FMA;ACCUM_INTS;FLOP_COUNT;PROFILE;CONTRACTED_INTS;SINGLE_EVALTYPE;USE_COMPOSITE_EVALUATORS;ERI_STRATEGY;MULTIPOLE_MAX_ORDER) +foreach (_var API_PREFIX;ERI3_PURE_SH;ERI2_PURE_SH;DISABLE_ONEBODY_PROPERTY_DERIVS;ENABLE_UNROLLING;ENABLE_GENERIC_CODE;VECTOR_LENGTH;VECTOR_METHOD;ALIGN_SIZE;USER_DEFINED_REAL;USER_DEFINED_REAL_INCLUDES;GENERATE_FMA;ACCUM_INTS;FLOP_COUNT;PROFILE;CONTRACTED_INTS;SINGLE_EVALTYPE;USE_COMPOSITE_EVALUATORS;ERI_STRATEGY;MULTIPOLE_MAX_ORDER) if (DEFINED LIBINT2_${_var}) if (DEFINED LIBINT_${_var}) message(FATAL_ERROR "renaming user-facing LIBINT2_${_var} variable but internal variable LIBINT_${_var} already exists") else () set(LIBINT_${_var} ${LIBINT2_${_var}}) - endif() - endif() + endif () + endif () endforeach () configure_file( - include/libint2/config.h.cmake.in - include/libint2/config.h - @ONLY) + include/libint2/config.h.cmake.in + include/libint2/config.h + @ONLY) set(_EXPORT_MODE 1) configure_file( - include/libint2/config.h.cmake.in - ${EXPORT_STAGE_DIR}/include/libint2/config.h - @ONLY) + include/libint2/config.h.cmake.in + ${EXPORT_STAGE_DIR}/include/libint2/config.h + @ONLY) configure_file( - include/libint2/config2.h.cmake.in - ${EXPORT_STAGE_DIR}/include/libint2/config2.h.cmake.in - COPYONLY) + include/libint2/config2.h.cmake.in + ${EXPORT_STAGE_DIR}/include/libint2/config2.h.cmake.in + COPYONLY) add_subdirectory(src) diff --git a/cmake/modules/int_am.cmake b/cmake/modules/int_am.cmake index c1d61cd55..cc86b7aa7 100644 --- a/cmake/modules/int_am.cmake +++ b/cmake/modules/int_am.cmake @@ -357,6 +357,7 @@ endmacro() process_integrals_class(ONEBODY) process_integrals_class(ERI) +process_integrals_class(RKB_ERI) process_integrals_class(ERI3) process_integrals_class(ERI2) # unlike above, these classes (1) don't do AM_LIST and (2) require value in config.h if enabled @@ -396,7 +397,7 @@ list(REVERSE _amlist) list(APPEND Libint2_ERI_COMPONENTS "${_amlist}") message(VERBOSE "setting components ${_amlist}") -foreach(_cls ONEBODY;ERI;ERI3;ERI2;G12;G12DKH) +foreach(_cls ONEBODY;ERI;RKB_ERI;ERI3;ERI2;G12;G12DKH) if((_cls STREQUAL G12) OR (_cls STREQUAL G12DKH)) add_feature_info( "integral class ${_cls}" @@ -433,6 +434,8 @@ foreach(_cls ONEBODY;ERI;ERI3;ERI2;G12;G12DKH) list(APPEND _amlist "onebody_${_am${_l}}${_am${_l}}_d${_d}") elseif (_cls STREQUAL "G12") list(APPEND _amlist "g12_${_am${_l}}${_am${_l}}${_am${_l}}${_am${_l}}_d${_d}") + elseif (_cls STREQUAL "RKB_ERI") + list(APPEND _amlist "rkb_eri_${_am${_l}}${_am${_l}}${_am${_l}}${_am${_l}}_d${_d}") endif() endforeach() if (_cls STREQUAL "ERI3") diff --git a/export/tests/unit/test-2body.cc b/export/tests/unit/test-2body.cc index fd602a910..2f0383fc1 100644 --- a/export/tests/unit/test-2body.cc +++ b/export/tests/unit/test-2body.cc @@ -1,5 +1,5 @@ /* - * Copyright (C) 2004-2024 Edward F. Valeev + * Copyright (C) 2004-2026 Edward F. Valeev * * This file is part of Libint library. * @@ -344,6 +344,262 @@ TEST_CASE("eri geometric derivatives", "[engine][2-body]") { } } +TEST_CASE("RKB Coulomb integrals", "[engine][2-body]") { + std::vector obs{// pseudorandom s + Shell{{1.0}, {{0, false, {1.0}}}, {{0.0, 0.0, 0.0}}}, + // pseudorandom p + Shell{{2.0}, {{1, false, {1.0}}}, {{1.0, 1.0, 1.0}}}}; + + const auto max_nprim = libint2::max_nprim(obs); + const auto max_l = libint2::max_l(obs); + typedef std::array der_idx; + + // e.g. d_xx maps the derivative index of derivative w.r.t x + // coord of ket1 and x coord of ket2 in Chemist notation. + // deriv indices for (LL|SS) + der_idx d_xx = {0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0}; + der_idx d_yy = {0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0}; + der_idx d_zz = {0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1}; + der_idx d_yz = {0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1}; + der_idx d_zy = {0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0}; + der_idx d_zx = {0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0}; + der_idx d_xz = {0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1}; + der_idx d_xy = {0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0}; + der_idx d_yx = {0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0}; + + // deriv indices for (SS|SS) + // 0th component + der_idx xxxx = {1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0}; + der_idx yyxx = {0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0}; + der_idx zzxx = {0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0}; + der_idx yxyx = {0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0}; + der_idx xyyx = {1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0}; + der_idx yxxy = {0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0}; + der_idx xyxy = {1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0}; + der_idx xxyy = {1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0}; + der_idx yyyy = {0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0}; + der_idx zzyy = {0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0}; + der_idx xxzz = {1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1}; + der_idx yyzz = {0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1}; + der_idx zzzz = {0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1}; + + // x-component + der_idx zxzx = {0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0}; + der_idx xzzx = {1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0}; + der_idx zyzy = {0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0}; + der_idx yzzy = {0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0}; + der_idx zxxz = {0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1}; + der_idx xzxz = {1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1}; + der_idx zyyz = {0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1}; + der_idx yzyz = {0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1}; + + // y-component + der_idx zyzx = {0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0}; + der_idx yzzx = {0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0}; + der_idx zxzy = {0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0}; + der_idx xzzy = {1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0}; + der_idx zyxz = {0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1}; + der_idx yzxz = {0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1}; + der_idx zxyz = {0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1}; + der_idx xzyz = {1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1}; + + // z-component + der_idx yxxx = {0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0}; + der_idx xyxx = {1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0}; + der_idx xxyx = {1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0}; + der_idx yyyx = {0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0}; + der_idx zzyx = {0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0}; + der_idx xxxy = {1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0}; + der_idx yyxy = {0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0}; + der_idx zzxy = {0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0}; + der_idx yxyy = {0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0}; + der_idx xyyy = {1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0}; + der_idx yxzz = {0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1}; + der_idx xyzz = {1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1}; + + SECTION("Coulombσpσp and σpσpCoulombσpσp") { + Engine engine_llss, engine_ssss; + try { + engine_llss = Engine(Operator::coulomb_opop, max_nprim, max_l, 0); + engine_ssss = Engine(Operator::opop_coulomb_opop, max_nprim, max_l, 0); + // TODO: need another unit test for derivatives of RKB ERIs + } catch ( + Engine::lmax_exceeded &) { // skip the test if lmax exceeded or libint2 + // not configured with RKB support + return; + } + + const auto nshell = obs.size(); + for (int s0 = 0; s0 != nshell; ++s0) { + for (int s1 = 0; s1 != nshell; ++s1) { + for (int s2 = 0; s2 != nshell; ++s2) { + for (int s3 = 0; s3 != nshell; ++s3) { + const auto &results_llss = + engine_llss.compute(obs[s0], obs[s1], obs[s2], obs[s3]); + const auto &results_ssss = + engine_ssss.compute(obs[s0], obs[s1], obs[s2], obs[s3]); + assert(results_llss.size() == + 4); // we get 4 buffers for each quaternion component + + LIBINT2_REF_REALTYPE Aref[3]; + for (int i = 0; i < 3; ++i) Aref[i] = obs[s0].O[i]; + LIBINT2_REF_REALTYPE Bref[3]; + for (int i = 0; i < 3; ++i) Bref[i] = obs[s1].O[i]; + LIBINT2_REF_REALTYPE Cref[3]; + for (int i = 0; i < 3; ++i) Cref[i] = obs[s2].O[i]; + LIBINT2_REF_REALTYPE Dref[3]; + for (int i = 0; i < 3; ++i) Dref[i] = obs[s3].O[i]; + + int ijkl = 0; + + int l0, m0, n0; + FOR_CART(l0, m0, n0, obs[s0].contr[0].l) + + int l1, m1, n1; + FOR_CART(l1, m1, n1, obs[s1].contr[0].l) + + int l2, m2, n2; + FOR_CART(l2, m2, n2, obs[s2].contr[0].l) + + int l3, m3, n3; + FOR_CART(l3, m3, n3, obs[s3].contr[0].l) + + std::array ref_coulomb_opop{0.0, 0.0, 0.0, + 0.0}; + std::array ref_opop_coulomb_opop{0.0, 0.0, + 0.0, 0.0}; + uint p0123 = 0; + for (uint p0 = 0; p0 < obs[s0].nprim(); p0++) { + for (uint p1 = 0; p1 < obs[s1].nprim(); p1++) { + for (uint p2 = 0; p2 < obs[s2].nprim(); p2++) { + for (uint p3 = 0; p3 < obs[s3].nprim(); p3++, p0123++) { + const LIBINT2_REF_REALTYPE alpha0 = obs[s0].alpha[p0]; + const LIBINT2_REF_REALTYPE alpha1 = obs[s1].alpha[p1]; + const LIBINT2_REF_REALTYPE alpha2 = obs[s2].alpha[p2]; + const LIBINT2_REF_REALTYPE alpha3 = obs[s3].alpha[p3]; + + const LIBINT2_REF_REALTYPE c0 = obs[s0].contr[0].coeff[p0]; + const LIBINT2_REF_REALTYPE c1 = obs[s1].contr[0].coeff[p1]; + const LIBINT2_REF_REALTYPE c2 = obs[s2].contr[0].coeff[p2]; + const LIBINT2_REF_REALTYPE c3 = obs[s3].contr[0].coeff[p3]; + const LIBINT2_REF_REALTYPE c0123 = c0 * c1 * c2 * c3; + + auto eri_drrrr = [&](der_idx d_rrrr) { + return eri(d_rrrr.data(), l0, m0, n0, alpha0, Aref, l1, + m1, n1, alpha1, Bref, l2, m2, n2, alpha2, Cref, + l3, m3, n3, alpha3, Dref, 0); + }; + + // (LL|SS) + ref_coulomb_opop[0] += + c0123 * + (eri_drrrr(d_xx) + eri_drrrr(d_yy) + eri_drrrr(d_zz)); + ref_coulomb_opop[1] += + c0123 * (eri_drrrr(d_yz) - eri_drrrr(d_zy)); + ref_coulomb_opop[2] += + c0123 * (eri_drrrr(d_zx) - eri_drrrr(d_xz)); + ref_coulomb_opop[3] += + c0123 * (eri_drrrr(d_xy) - eri_drrrr(d_yx)); + + // (SS|SS) + ref_opop_coulomb_opop[0] += + c0123 * + (eri_drrrr(xxxx) + eri_drrrr(yyxx) + eri_drrrr(zzxx) - + eri_drrrr(yxyx) + eri_drrrr(xyyx) + eri_drrrr(yxxy) - + eri_drrrr(xyxy) + eri_drrrr(xxyy) + eri_drrrr(yyyy) + + eri_drrrr(zzyy) + eri_drrrr(xxzz) + eri_drrrr(yyzz) + + eri_drrrr(zzzz)); + ref_opop_coulomb_opop[1] += + c0123 * + (eri_drrrr(zxzx) - eri_drrrr(xzzx) - eri_drrrr(zyzy) + + eri_drrrr(yzzy) - eri_drrrr(zxxz) + eri_drrrr(xzxz) + + eri_drrrr(zyyz) - eri_drrrr(yzyz)); + ref_opop_coulomb_opop[2] += + c0123 * + (-eri_drrrr(zyzx) + eri_drrrr(yzzx) - eri_drrrr(zxzy) + + eri_drrrr(xzzy) + eri_drrrr(zyxz) - eri_drrrr(yzxz) + + eri_drrrr(zxyz) - eri_drrrr(xzyz)); + ref_opop_coulomb_opop[3] += + c0123 * + (-eri_drrrr(yxxx) + eri_drrrr(xyxx) - eri_drrrr(xxyx) - + eri_drrrr(yyyx) - eri_drrrr(zzyx) + eri_drrrr(xxxy) + + eri_drrrr(yyxy) + eri_drrrr(zzxy) - eri_drrrr(yxyy) + + eri_drrrr(xyyy) - eri_drrrr(yxzz) + eri_drrrr(xyzz)); + } + } + } + } + + const double ABSOLUTE_DEVIATION_THRESHOLD = 5.0E-14; + const double RELATIVE_DEVIATION_THRESHOLD = + 1.0E-9; // For more detail on choice of these thresholds, see + // the comments in the TEST_CASE "eri geometric + // derivatives" + + std::array abs_errs_llss; + std::array rel_abs_errs_llss; + + std::array abs_errs_ssss; + std::array rel_abs_errs_ssss; + + for (auto comp = 0; comp < 4; ++comp) { + abs_errs_llss[comp] = + abs(ref_coulomb_opop[comp] - results_llss[comp][ijkl]); + rel_abs_errs_llss[comp] = + abs(abs_errs_llss[comp] / ref_coulomb_opop[comp]); + + abs_errs_ssss[comp] = + abs(ref_opop_coulomb_opop[comp] - results_ssss[comp][ijkl]); + rel_abs_errs_ssss[comp] = + abs(abs_errs_ssss[comp] / ref_opop_coulomb_opop[comp]); + + bool llss_not_ok = + rel_abs_errs_llss[comp] > RELATIVE_DEVIATION_THRESHOLD && + abs_errs_llss[comp] > ABSOLUTE_DEVIATION_THRESHOLD; + + bool ssss_not_ok = + rel_abs_errs_ssss[comp] > RELATIVE_DEVIATION_THRESHOLD && + abs_errs_ssss[comp] > ABSOLUTE_DEVIATION_THRESHOLD; + + // no 3^n prefactor here since the intrinsic deriv order is 2 + if (llss_not_ok) { + std::cout << "(l0 l1| l2 l3) = " + << "(" << s0 << " " << s1 << " | " << s2 << " " << s3 + << ") " + << "Elem " << ijkl << " comp= " << comp + << " : ref = " << ref_coulomb_opop[comp] + << " libint = " << results_llss[comp][ijkl] + << " relabs_error = " << rel_abs_errs_llss[comp] + << " abs_error = " << abs_errs_llss[comp] + << std::endl; + } + if (ssss_not_ok) { + std::cout << "(l0 l1| l2 l3) = " + << "(" << s0 << " " << s1 << " | " << s2 << " " << s3 + << ") " + << "Elem " << ijkl << " comp= " << comp + << " : ref = " << ref_opop_coulomb_opop[comp] + << " libint = " << results_ssss[comp][ijkl] + << " relabs_error = " << rel_abs_errs_ssss[comp] + << " abs_error = " << abs_errs_ssss[comp] + << std::endl; + } + REQUIRE(!llss_not_ok); + REQUIRE(!ssss_not_ok); + } + + ++ijkl; + END_FOR_CART + END_FOR_CART + END_FOR_CART + END_FOR_CART + } + } + } + } + } +} + TEST_CASE("Erfx_Coulomb integrals", "[engine][2-body]") { // pseudorandom s shells std::vector obs{ @@ -374,12 +630,12 @@ TEST_CASE("Erfx_Coulomb integrals", "[engine][2-body]") { REQUIRE(results[0] != nullptr); switch (k) { /* VALIDATION WOLFRAM CODE: -(* Integral of Coulomb kernel damped by (\[Lambda] Erf[\[Omega] r] + \ -\[Sigma] Erfc[\[Omega] r]), over unit-normalized s functions, \ -see Eq 52 in DOI 10.1039/b605188j *) -F0[T_] := If[T == 0, 1, Sqrt[\[Pi]/T]*Erf[Sqrt[T]]/2]; -sN[a_] := ((2 a)/\[Pi])^(3/4); -VVeeErfx[\[Alpha]1_, A1_List, \[Alpha]2_, A2_List, \[Beta]1_, + (* Integral of Coulomb kernel damped by (\[Lambda] Erf[\[Omega] r] + \ + \[Sigma] Erfc[\[Omega] r]), over unit-normalized s functions, \ + see Eq 52 in DOI 10.1039/b605188j *) + F0[T_] := If[T == 0, 1, Sqrt[\[Pi]/T]*Erf[Sqrt[T]]/2]; + sN[a_] := ((2 a)/\[Pi])^(3/4); + VVeeErfx[\[Alpha]1_, A1_List, \[Alpha]2_, A2_List, \[Beta]1_, B1_List, \[Beta]2_, B2_List, \[Omega]_, \[Lambda]_, \[Sigma]_] := Module[{\[Gamma]1, \[Gamma]2, P1, P2, K1, K2, T, result, \[Rho]}, \[Gamma]1 = \[Alpha]1 + \[Beta]1; @@ -397,13 +653,13 @@ VVeeErfx[\[Alpha]1_, A1_List, \[Alpha]2_, A2_List, \[Beta]1_, T]) sN[\[Alpha]1] sN[\[Alpha]2] sN[\[Beta]1] sN[\[Beta]2]; Return[result]; ]; -Print[CForm[ + Print[CForm[ N[VVeeErfx[1, {0, 0, 0}, 3, {2, 2, 2}, 2, {1, 1, 1}, 4, {3, 3, 3}, 1.1, 1, 0], 20]]] -Print[CForm[ + Print[CForm[ N[VVeeErfx[1, {0, 0, 0}, 3, {2, 2, 2}, 2, {1, 1, 1}, 4, {3, 3, 3}, 1.1, 0, 1], 20]]] -Print[CForm[ + Print[CForm[ N[VVeeErfx[1, {0, 0, 0}, 3, {2, 2, 2}, 2, {1, 1, 1}, 4, {3, 3, 3}, 1.1, 2, 3], 20]]] */ diff --git a/include/libint2.h b/include/libint2.h index 59d17ef65..25a231be7 100644 --- a/include/libint2.h +++ b/include/libint2.h @@ -22,18 +22,19 @@ #define _libint2_header_ #define LIBINT_T_SS_EREP_SS(mValue) \ - _aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_##mValue + _aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_##mValue #define LIBINT_T_SS_Km1G12_SS(mValue) \ - _aB_s___0__s___1___r12_minus_1_g12_s___0__s___1___Ab__up_##mValue + _aB_s____0__s____1___r12_minus_1_g12_s____0__s____1___Ab__up_##mValue #define LIBINT_T_SS_K0G12_SS_0 \ - _aB_s___0__s___1___r12_0_g12_s___0__s___1___Ab__up_0 + _aB_s____0__s____1___r12_0_g12_s____0__s____1___Ab__up_0 #define LIBINT_T_SS_K2G12_SS_0 \ - _aB_s___0__s___1___r12_2_g12_s___0__s___1___Ab__up_0 + _aB_s____0__s____1___r12_2_g12_s____0__s____1___Ab__up_0 #define LIBINT_T_SS_K4G12_SS_0 \ - _aB_s___0__s___1___r12_4_g12_s___0__s___1___Ab__up_0 -#define LIBINT_T_S_OVERLAP_S _aB_s___0___Overlap_s___0___Ab__up_ -#define LIBINT_T_S_KINETIC_S _aB_s___0___Kinetic_s___0___Ab__up_ -#define LIBINT_T_S_ELECPOT_S(mValue) _aB_s___0___ElecPot_s___0___Ab__up_##mValue + _aB_s____0__s____1___r12_4_g12_s____0__s____1___Ab__up_0 +#define LIBINT_T_S_OVERLAP_S _aB_s____0___Overlap_s____0___Ab__up_ +#define LIBINT_T_S_KINETIC_S _aB_s____0___Kinetic_s____0___Ab__up_ +#define LIBINT_T_S_ELECPOT_S(mValue) \ + _aB_s____0___ElecPot_s____0___Ab__up_##mValue #include #include diff --git a/include/libint2/config.h.cmake.in b/include/libint2/config.h.cmake.in index 640e68099..3eb32e2f0 100644 --- a/include/libint2/config.h.cmake.in +++ b/include/libint2/config.h.cmake.in @@ -71,6 +71,13 @@ #undef LIBINT_INCLUDE_ERI #endif +/* Support ERI derivatives up to this order */ +#define LIBINT_INCLUDE_RKB_ERI @LIBINT_INCLUDE_RKB_ERI@ +#if @LIBINT_INCLUDE_RKB_ERI@ == -1 +#undef LIBINT_INCLUDE_RKB_ERI +#endif + + /* Support 3-center ERI derivatives up to this order */ #define LIBINT_INCLUDE_ERI3 @LIBINT_INCLUDE_ERI3@ #if @LIBINT_INCLUDE_ERI3@ == -1 @@ -122,6 +129,18 @@ /* Max optimized AM for ERI and its derivatives */ #cmakedefine LIBINT_ERI_OPT_AM_LIST "@LIBINT_ERI_OPT_AM_LIST@" +/* Max AM for RKB_ERI (same for all derivatives; if not defined see LIBINT_RKB_ERI_MAX_AM_LIST) */ +#cmakedefine LIBINT_RKB_ERI_MAX_AM @LIBINT_RKB_ERI_MAX_AM@ + +/* Max AM for RKB_ERI and its derivatives */ +#cmakedefine LIBINT_RKB_ERI_MAX_AM_LIST "@LIBINT_RKB_ERI_MAX_AM_LIST@" + +/* Max optimized AM for ERI (same for all derivatives; if not defined see LIBINT_RKB_ERI_OPT_AM_LIST) */ +#cmakedefine LIBINT_RKB_ERI_OPT_AM @LIBINT_RKB_ERI_OPT_AM@ + +/* Max optimized AM for ERI and its derivatives */ +#cmakedefine LIBINT_RKB_ERI_OPT_AM_LIST "@LIBINT_RKB_ERI_OPT_AM_LIST@" + /* Max AM for 3-center ERI (same for all derivatives; if not defined see LIBINT_ERI3_MAX_AM_LIST) */ #cmakedefine LIBINT_ERI3_MAX_AM @LIBINT_ERI3_MAX_AM@ diff --git a/include/libint2/cxxapi.h b/include/libint2/cxxapi.h index 0aceb509b..22686f958 100644 --- a/include/libint2/cxxapi.h +++ b/include/libint2/cxxapi.h @@ -35,9 +35,9 @@ #if !defined(LIBINT_INCLUDE_ONEBODY) || \ !(defined(LIBINT_INCLUDE_ERI) || defined(LIBINT_INCLUDE_ERI3) || \ - defined(LIBINT_INCLUDE_ERI2)) + defined(LIBINT_INCLUDE_ERI2) || defined(LIBINT_INCLUDE_RKB_ERI)) #error \ - "C++ API is only supported if both 1-body and some (eri, eri3, eri2) 2-body integrals are enabled" + "C++ API is only supported if both 1-body and some (eri, eri3, eri2, rkb_eri) 2-body integrals are enabled" #endif #include diff --git a/include/libint2/engine.h b/include/libint2/engine.h index 36db41821..94077da63 100644 --- a/include/libint2/engine.h +++ b/include/libint2/engine.h @@ -153,6 +153,13 @@ enum class Operator { coulomb, /// alias for Operator::coulomb r12_m1 = coulomb, + /// (2-body) \f$ r_{12}^{-1} (σ.p_{k1})(σ.p_{k2})\f$ where k1 & k2 are + /// centers of ket1 and ket2, respectively + coulomb_opop, + /// (2-body) \f$ (σ.p_{b1})(σ.p_{b2}) r_{12}^{-1} (σ.p_{k1})(σ.p_{k2})\f$ + /// where b1 & b2 are centers of bra1 and bra2 and k1 & k2 are centers of + /// ket1 and ket2, respectively + opop_coulomb_opop, /// contracted Gaussian geminal cgtg, /// contracted Gaussian geminal times Coulomb @@ -246,6 +253,7 @@ struct operator_traits typedef const libint2::FmEval_Reference core_eval_type; #endif }; + template <> struct operator_traits : public operator_traits { @@ -346,6 +354,20 @@ struct operator_traits typedef const libint2::FmEval_Reference core_eval_type; #endif }; + +template <> +struct operator_traits + : public operator_traits { + static constexpr auto nopers = 4; + static constexpr auto intrinsic_deriv_order = 2; +}; +template <> +struct operator_traits + : public operator_traits { + static constexpr auto nopers = 4; + static constexpr auto intrinsic_deriv_order = 4; +}; + namespace detail { template struct cgtg_operator_traits : public detail::default_operator_traits { @@ -839,16 +861,16 @@ class Engine { const Shell& ket2, const ShellPair* spbra, const ShellPair* spket); // clang-format off - /** this specifies target precision for computing the integrals, i.e. - * the target absolute (i.e., not relative) error of the integrals. - * It is used to screen out primitive integrals. For some screening - * methods precision can be almost guaranteed (due to finite precision - * of the precomputed interpolation tables used to evaluate the core integrals - * it is not in general possible to guarantee precision rigorously). - * - * @param[in] prec the target precision - * @sa ScreeningMethod - */ + /** this specifies target precision for computing the integrals, i.e. + * the target absolute (i.e., not relative) error of the integrals. + * It is used to screen out primitive integrals. For some screening + * methods precision can be almost guaranteed (due to finite precision + * of the precomputed interpolation tables used to evaluate the core integrals + * it is not in general possible to guarantee precision rigorously). + * + * @param[in] prec the target precision + * @sa ScreeningMethod + */ // clang-format on Engine& set_precision(scalar_type prec) { if (prec <= 0.) { diff --git a/include/libint2/engine.impl.h b/include/libint2/engine.impl.h index a97fd914a..a39244e7e 100644 --- a/include/libint2/engine.impl.h +++ b/include/libint2/engine.impl.h @@ -70,30 +70,32 @@ typename std::remove_all_extents::type* to_ptr1(T (&a)[N]) { /// These MUST appear in the same order as in Operator. /// You must also update BOOST_PP_NBODY_OPERATOR_LAST_ONEBODY_INDEX when you add /// one-body ints -#define BOOST_PP_NBODY_OPERATOR_LIST \ - (overlap, /* overlap */ \ - (kinetic, /* kinetic */ \ - (elecpot, /* nuclear */ \ - (elecpot, /* erf_nuclear */ \ - (elecpot, /* erfc_nuclear */ \ - (elecpot, /* erfx_nuclear */ \ - (1emultipole, /* emultipole1 */ \ - (2emultipole, /* emultipole2 */ \ - (3emultipole, /* emultipole3 */ \ - (sphemultipole, /* sphemultipole */ \ - (opVop, /* opVop */ \ - (eri, /* delta */ \ - (eri, /* coulomb */ \ - (eri, /* cgtg */ \ - (eri, /* cgtg_x_coulomb */ \ - (eri, /* delcgtg2 */ \ - (eri, /* r12 */ \ - (eri, /* erf_coulomb */ \ - (eri, /* erfc_coulomb */ \ - (eri, /* erfx_coulomb */ \ - (eri, /* stg */ \ - (eri, /* yukawa */ \ - BOOST_PP_NIL)))))))))))))))))))))) +#define BOOST_PP_NBODY_OPERATOR_LIST \ + (overlap, /* overlap */ \ + (kinetic, /* kinetic */ \ + (elecpot, /* nuclear */ \ + (elecpot, /* erf_nuclear */ \ + (elecpot, /* erfc_nuclear */ \ + (elecpot, /* erfx_nuclear */ \ + (1emultipole, /* emultipole1 */ \ + (2emultipole, /* emultipole2 */ \ + (3emultipole, /* emultipole3 */ \ + (sphemultipole, /* sphemultipole */ \ + (opVop, /* opVop */ \ + (eri, /* delta */ \ + (eri, /* coulomb */ \ + (coulomb_opop, /* coulomb_opop */ \ + (opop_coulomb_opop, /* coulomb_opop */ \ + (eri, /* cgtg */ \ + (eri, /* cgtg_x_coulomb */ \ + (eri, /* delcgtg2 */ \ + (eri, /* r12 */ \ + (eri, /* erf_coulomb */ \ + (eri, /* erfc_coulomb */ \ + (eri, /* erfx_coulomb */ \ + (eri, /* stg */ \ + (eri, /* yukawa */ \ + BOOST_PP_NIL)))))))))))))))))))))))) #define BOOST_PP_NBODY_OPERATOR_INDEX_TUPLE \ BOOST_PP_MAKE_TUPLE(BOOST_PP_LIST_SIZE(BOOST_PP_NBODY_OPERATOR_LIST)) @@ -628,7 +630,8 @@ __libint2_engine_inline void Engine::_initialize() { if (lmax_ >= lmax) { \ throw Engine::lmax_exceeded(BOOST_PP_STRINGIZE(BOOST_PP_NBODYENGINE_MCR3_TASK(product)), lmax, lmax_); \ } \ - if (stack_size_ > 0) LIBINT2_PREFIXED_NAME(libint2_cleanup_default)(&primdata_[0]); \ + if (stack_size_ > 0) \ + LIBINT2_PREFIXED_NAME(libint2_cleanup_default)(&primdata_[0]); \ stack_size_ = LIBINT2_PREFIXED_NAME(BOOST_PP_CAT( \ libint2_need_memory_, BOOST_PP_NBODYENGINE_MCR3_TASK(product)))( \ lmax_); \ @@ -663,23 +666,23 @@ __libint2_engine_inline void Engine::initialize(size_t max_nprim) { // validate braket #ifndef LIBINT_INCLUDE_ONEBODY assert(braket_ != BraKet::x_x && - "this braket type not supported by the library; give --enable-1body " - "to configure"); + "this braket type not supported by the library; configure with " + "-DLIBINT2_ENABLE_ONEBODY >= 0"); #endif #ifndef LIBINT_INCLUDE_ERI assert(braket_ != BraKet::xx_xx && - "this braket type not supported by the library; give --enable-eri to " - "configure"); + "this braket type not supported by the library; configure with " + "-DLIBINT2_ENABLE_ERI >= 0"); #endif #ifndef LIBINT_INCLUDE_ERI3 assert((braket_ != BraKet::xs_xx && braket_ != BraKet::xx_xs) && - "this braket type not supported by the library; give --enable-eri3 to " - "configure"); + "this braket type not supported by the library; configure with " + "-DLIBINT2_ENABLE_ERI3 >= 0"); #endif #ifndef LIBINT_INCLUDE_ERI2 assert(braket_ != BraKet::xs_xs && - "this braket type not supported by the library; give --enable-eri2 to " - "configure"); + "this braket type not supported by the library; configure with " + "-DLIBINT2_ENABLE_ERI2 >= 0"); #endif // make sure it's no default initialized @@ -701,6 +704,7 @@ __libint2_engine_inline void Engine::initialize(size_t max_nprim) { // target indices. const auto permutable_targets = deriv_order_ > 0 && + (braket_ == BraKet::xx_xx || braket_ == BraKet::xs_xx || braket_ == BraKet::xx_xs); if (permutable_targets) @@ -1212,19 +1216,31 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( #if LIBINT2_SHELLQUARTET_SET == \ LIBINT2_SHELLQUARTET_SET_STANDARD // standard angular momentum ordering - const auto swap_tbra = (tbra1.contr[0].l < tbra2.contr[0].l); - const auto swap_tket = (tket1.contr[0].l < tket2.contr[0].l); - const auto swap_braket = - ((braket_ == BraKet::xx_xx) && (tbra1.contr[0].l + tbra2.contr[0].l > - tket1.contr[0].l + tket2.contr[0].l)) || - braket_ == BraKet::xx_xs; + const auto swap_braket = ((braket_ == BraKet::xx_xx) && + (tbra1.contr[0].l + tbra2.contr[0].l > + tket1.contr[0].l + tket2.contr[0].l) && + (oper_ != Operator::coulomb_opop)) || + braket_ == BraKet::xx_xs; + bool swap_tbra, swap_tket; + if (oper_ == Operator::opop_coulomb_opop) { + bool swap_p1p2 = swap_braket ? (tbra1.contr[0].l < tbra2.contr[0].l) + : (tket1.contr[0].l < tket2.contr[0].l); + swap_tbra = swap_tket = swap_p1p2; + } else { + swap_tbra = (tbra1.contr[0].l < tbra2.contr[0].l); + swap_tket = (tket1.contr[0].l < tket2.contr[0].l); + } + + // N.B. cannot swap bra and ket for coulomb_opop since the ket is mutated by + // this operator #else // orca angular momentum ordering const auto swap_tbra = (tbra1.contr[0].l > tbra2.contr[0].l); const auto swap_tket = (tket1.contr[0].l > tket2.contr[0].l); - const auto swap_braket = - ((braket_ == BraKet::xx_xx) && (tbra1.contr[0].l + tbra2.contr[0].l < - tket1.contr[0].l + tket2.contr[0].l)) || - braket_ == BraKet::xx_xs; + const auto swap_braket = ((braket_ == BraKet::xx_xx) && + (tbra1.contr[0].l + tbra2.contr[0].l < + tket1.contr[0].l + tket2.contr[0].l) && + (oper_ != Operator::coulomb_opop)) || + braket_ == BraKet::xx_xs; assert(false && "feature not implemented"); abort(); #endif @@ -1421,7 +1437,7 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( const scalar_type rho = gammap * gammaq * oogammapq; const scalar_type T = PQ2 * rho; auto* gm_ptr = &(primdata.LIBINT_T_SS_EREP_SS(0)[0]); - const auto mmax = l + deriv_order_; + const auto mmax = l + deriv_order_ + intrinsic_deriv_order(); if (!skip_core_ints) { switch (oper_) { @@ -1432,6 +1448,20 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( .first(); core_eval_ptr->eval(gm_ptr, T, mmax); } break; + case Operator::coulomb_opop: { + const auto& core_eval_ptr = + any_cast&>(core_eval_pack_) + .first(); + core_eval_ptr->eval(gm_ptr, T, mmax); + } break; + case Operator::opop_coulomb_opop: { + const auto& core_eval_ptr = + any_cast&>(core_eval_pack_) + .first(); + core_eval_ptr->eval(gm_ptr, T, mmax); + } break; case Operator::cgtg_x_coulomb: { const auto& core_eval_ptr = any_cast 0 || lmax_bra > 0) { + if (deriv_order_ + intrinsic_deriv_order() > 0 || lmax_bra > 0) { #if LIBINT2_DEFINED(eri, WP_x) primdata.WP_x[0] = Wx - P[0]; #endif @@ -1662,7 +1692,7 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( primdata.WP_z[0] = Wz - P[2]; #endif } - if (deriv_order_ > 0 || lmax_ket > 0) { + if (deriv_order_ + intrinsic_deriv_order() > 0 || lmax_ket > 0) { #if LIBINT2_DEFINED(eri, WQ_x) primdata.WQ_x[0] = Wx - Q[0]; #endif @@ -1746,7 +1776,7 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( #endif // prefactors for derivative ERI relations - if (deriv_order_ > 0) { + if (deriv_order_ + intrinsic_deriv_order() > 0) { #if LIBINT2_DEFINED(eri, alpha1_rho_over_zeta2) primdata.alpha1_rho_over_zeta2[0] = alpha0 * (oogammap * gammaq_o_gammapgammaq); @@ -1829,7 +1859,8 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( } // compute directly (ss|ss) - const auto compute_directly = lmax == 0 && deriv_order_ == 0; + const auto compute_directly = + lmax == 0 && deriv_order_ == 0 & intrinsic_deriv_order() == 0; if (compute_directly) { #ifdef LIBINT2_ENGINE_TIMERS @@ -1907,8 +1938,10 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( "the angular momentum limit is exceeded"); assert(ket2.contr[0].l <= ket_lmax && "the angular momentum limit is exceeded"); + buildfnidx = (bra1.contr[0].l * ket_lmax + ket1.contr[0].l) * ket_lmax + ket2.contr[0].l; + #ifdef LIBINT_ERI3_PURE_SH if (bra1.contr[0].l > 1) assert(bra1.contr[0].pure && @@ -2101,9 +2134,19 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( const auto tgt_row_idx = !swap_tbra ? r1 * nr2 + r2 : r2 * nr1 + r1; Map tgt_blk_mat(tgt_ptr + tgt_row_idx * ncol, nc1_tgt, nc2_tgt); - if (swap_tket) - tgt_blk_mat = src_blk_mat.transpose(); - else + if (swap_tket) { + Shell::real_t oper_cart_component_phase = 1.0; + if (oper_ == Operator::opop_coulomb_opop && s == 3) + oper_cart_component_phase = + -1.0; // z quaternion components flip sign on + // swapping ket + if (oper_ == Operator::coulomb_opop && s > 0) + oper_cart_component_phase = + -1.0; // x,y,z quaternion components flip sign on + // swapping ket for coulomb_opop + tgt_blk_mat = + oper_cart_component_phase * src_blk_mat.transpose(); + } else tgt_blk_mat = src_blk_mat; } } // end of loop diff --git a/src/bin/libint/build_libint.cc b/src/bin/libint/build_libint.cc index 8b5ece9ea..c2a9de6cd 100644 --- a/src/bin/libint/build_libint.cc +++ b/src/bin/libint/build_libint.cc @@ -73,18 +73,37 @@ enum ShellSetType { template struct ShellQuartetSetPredicate { // return true if this set of angular momenta is included - static bool value(int la, int lb, int lc, int ld); + static bool value(int la, int lb, int lc, int ld, bool p1p2_swappable = true); }; + +/** + * standard ordering for angular momenta la, lb, lc, ld + * @param p1p2_swappable whether operator allows swaps of particle 1 and 2 + * functions (e.g., not allowed for Coulombσpσp but allowed + * for Coulomb (TwoPRep)). + * @param bra_ket_coswappable whether need to swap within both bra and ket. + * Not individually swapping of either ket of bra allowed + * ( e.g., for σpσpCoulombσpσp) + */ template <> struct ShellQuartetSetPredicate { - static bool value(int la, int lb, int lc, int ld) { - return la >= lb && lc >= ld && la + lb <= lc + ld; + static bool value(int la, int lb, int lc, int ld, bool p1p2_swappable = true, + bool bra_ket_coswappable = false) { + if (bra_ket_coswappable) + return (la + lb <= lc + ld) && lc >= ld; + else + return la >= lb && lc >= ld && (!p1p2_swappable || la + lb <= lc + ld); } }; template <> struct ShellQuartetSetPredicate { - static bool value(int la, int lb, int lc, int ld) { - return la <= lb && lc <= ld && (la < lc || (la == lc && lb <= ld)); + static bool value(int la, int lb, int lc, int ld, bool p1p2_swappable = true, + bool bra_ket_coswappable = false) { + if (bra_ket_coswappable) + return (la < lc || (la == lc && lb <= ld)); + else + return la <= lb && lc <= ld && + (!p1p2_swappable || (la < lc || (la == lc && lb <= ld))); } }; template @@ -193,12 +212,16 @@ static void config_to_api(const std::shared_ptr& cparams, #ifdef LIBINT_INCLUDE_ERI #define USE_GENERIC_ERI_BUILD 1 #if !USE_GENERIC_ERI_BUILD +template static void build_TwoPRep_2b_2k( - std::ostream& os, const std::shared_ptr& cparams, + std::ostream& os, std::string label, + const std::shared_ptr& cparams, std::shared_ptr& iface); #else +template static void build_TwoPRep_2b_2k( - std::ostream& os, const std::shared_ptr& cparams, + std::ostream& os, std::string label, + const std::shared_ptr& cparams, std::shared_ptr& iface, unsigned int deriv_level); #endif #endif @@ -274,6 +297,17 @@ template <> σpVσp_Descr make_descr<σpVσp_Descr>(int p, int, int) { return σpVσp_Descr(p); } + +template <> +Coulombσpσp_Descr make_descr(int p, int, int) { + return Coulombσpσp_Descr(p); +} + +template <> +σpσpCoulombσpσp_Descr make_descr<σpσpCoulombσpσp_Descr>(int p, int, int) { + return σpσpCoulombσpσp_Descr(p); +} + } // namespace template @@ -567,6 +601,24 @@ void try_main(int argc, char* argv[]) { taskmgr.add(task_label("eri", d)); } #endif + +#ifdef LIBINT_INCLUDE_RKB_ERI +#define BOOST_PP_RKB_ERI_TASK_TUPLE (coulomb_opop, opop_coulomb_opop) +#define BOOST_PP_RKB_ERI_TASK_OPER_TUPLE (CoulombσpσpOper, σpσpCoulombσpσpOper) +#define BOOST_PP_RKB_ERI_TASK_LIST \ + BOOST_PP_TUPLE_TO_LIST(BOOST_PP_RKB_ERI_TASK_TUPLE) +#define BOOST_PP_RKB_ERI_TASK_OPER_LIST \ + BOOST_PP_TUPLE_TO_LIST(BOOST_PP_RKB_ERI_TASK_OPER_TUPLE) + + for (unsigned int d = 0; d <= LIBINT_INCLUDE_RKB_ERI; ++d) { +#define BOOST_PP_RKB_ERI_MCR1(r, data, elem) \ + taskmgr.add(task_label(BOOST_PP_STRINGIZE(elem), d)); + + BOOST_PP_LIST_FOR_EACH(BOOST_PP_RKB_ERI_MCR1, _, BOOST_PP_RKB_ERI_TASK_LIST) +#undef BOOST_PP_RKB_ERI_MCR1 + } +#endif + #ifdef LIBINT_INCLUDE_ERI3 for (unsigned int d = 0; d <= LIBINT_INCLUDE_ERI3; ++d) { taskmgr.add(task_label("3eri", d)); @@ -668,6 +720,46 @@ void try_main(int argc, char* argv[]) { cparams->num_bf(task_label("eri", d), 4); } #endif + +#ifdef LIBINT_INCLUDE_RKB_ERI + for (unsigned int d = 0; d <= LIBINT_INCLUDE_RKB_ERI; ++d) { +#if defined(LIBINT_RKB_ERI_MAX_AM_LIST) +#define BOOST_PP_RKB_ERI_MCR2(r, data, elem) \ + cparams->max_am( \ + task_label(BOOST_PP_STRINGIZE(elem), d), \ + token(LIBINT_RKB_ERI_MAX_AM_LIST, ',', d)); + BOOST_PP_LIST_FOR_EACH(BOOST_PP_RKB_ERI_MCR2, _, BOOST_PP_RKB_ERI_TASK_LIST) +#undef BOOST_PP_RKB_ERI_MCR2 +#elif defined(LIBINT_RKB_ERI_MAX_AM) +#define BOOST_PP_RKB_ERI_MCR3(r, data, elem) \ + cparams->max_am(task_label(BOOST_PP_STRINGIZE(elem), d), \ + LIBINT_RKB_ERI_MAX_AM); + BOOST_PP_LIST_FOR_EACH(BOOST_PP_RKB_ERI_MCR3, _, BOOST_PP_RKB_ERI_TASK_LIST) +#undef BOOST_PP_RKB_ERI_MCR3 +#endif +#if defined(LIBINT_RKB_ERI_OPT_AM_LIST) +#define BOOST_PP_RKB_ERI_MCR4(r, data, elem) \ + cparams->max_am_opt( \ + task_label(BOOST_PP_STRINGIZE(elem), d), \ + token(LIBINT_RKB_ERI_OPT_AM_LIST, ',', d)); + BOOST_PP_LIST_FOR_EACH(BOOST_PP_RKB_ERI_MCR4, _, BOOST_PP_RKB_ERI_TASK_LIST) +#undef BOOST_PP_RKB_ERI_MCR4 +#elif defined(LIBINT_RKB_ERI_OPT_AM) +#define BOOST_PP_RKB_ERI_MCR5(r, data, elem) \ + cparams->max_am_opt(task_label(BOOST_PP_STRINGIZE(elem), d), \ + LIBINT_RKB_ERI_OPT_AM); + BOOST_PP_LIST_FOR_EACH(BOOST_PP_RKB_ERI_MCR5, _, BOOST_PP_RKB_ERI_TASK_LIST) +#undef BOOST_PP_RKB_ERI_MCR5 +#endif + } + for (unsigned int d = 0; d <= LIBINT_INCLUDE_RKB_ERI; ++d) { +#define BOOST_PP_RKB_ERI_MCR6(r, data, elem) \ + cparams->num_bf(task_label(BOOST_PP_STRINGIZE(elem), d), 4); + BOOST_PP_LIST_FOR_EACH(BOOST_PP_RKB_ERI_MCR6, _, BOOST_PP_RKB_ERI_TASK_LIST) +#undef BOOST_PP_RKB_ERI_MCR6 + } +#endif // LIBINT_INCLUDE_RKB_ERI + #ifdef LIBINT_INCLUDE_ERI3 for (unsigned int d = 0; d <= LIBINT_INCLUDE_ERI3; ++d) { #if defined(LIBINT_ERI3_MAX_AM_LIST) @@ -852,6 +944,9 @@ void try_main(int argc, char* argv[]) { #ifdef LIBINT_INCLUDE_ERI max_deriv = std::max(LIBINT_INCLUDE_ERI, max_deriv); #endif +#ifdef LIBINT_INCLUDE_RKB_ERI + max_deriv = std::max(LIBINT_INCLUDE_RKB_ERI, max_deriv); +#endif #ifdef LIBINT_INCLUDE_ERI3 max_deriv = std::max(LIBINT_INCLUDE_ERI3, max_deriv); #endif @@ -879,13 +974,25 @@ void try_main(int argc, char* argv[]) { #endif #ifdef LIBINT_INCLUDE_ERI #if !USE_GENERIC_ERI_BUILD - build_TwoPRep_2b_2k(os, cparams, iface); + build_TwoPRep_2b_2k(os, "eri", cparams, iface); #else for (unsigned int d = 0; d <= LIBINT_INCLUDE_ERI; ++d) { - build_TwoPRep_2b_2k(os, cparams, iface, d); + build_TwoPRep_2b_2k(os, "eri", cparams, iface, d); } #endif #endif + +#ifdef LIBINT_INCLUDE_RKB_ERI + for (unsigned int d = 0; d <= LIBINT_INCLUDE_RKB_ERI; ++d) { +#define BOOST_PP_RKB_ERI_MCR7(r, data, i, elem) \ + build_TwoPRep_2b_2k( \ + os, BOOST_PP_STRINGIZE(elem), cparams, iface, d); + BOOST_PP_LIST_FOR_EACH_I(BOOST_PP_RKB_ERI_MCR7, _, + BOOST_PP_RKB_ERI_TASK_LIST) +#undef BOOST_PP_RKB_ERI_MCR7 + } +#endif + #ifdef LIBINT_INCLUDE_ERI3 for (unsigned int d = 0; d <= LIBINT_INCLUDE_ERI3; ++d) { build_TwoPRep_1b_2k(os, cparams, iface, d); @@ -985,15 +1092,26 @@ void print_config(std::ostream& os) { #ifdef LIBINT_INCLUDE_G12DKH os << "Will support G12DKH" << endl; #endif +#ifdef LIBINT_INCLUDE_RKB_ERI + os << "Will support restricted kinetically balance (RKB) 4-center ERIs " + << std::endl; + if (LIBINT_INCLUDE_RKB_ERI > 0) + os << "(deriv order = " << LIBINT_INCLUDE_RKB_ERI << ")"; + os << endl; +#endif } #ifdef LIBINT_INCLUDE_ERI -void build_TwoPRep_2b_2k(std::ostream& os, +template +void build_TwoPRep_2b_2k(std::ostream& os, std::string label, const std::shared_ptr& cparams, std::shared_ptr& iface, unsigned int deriv_level) { - const std::string task = task_label("eri", deriv_level); - typedef TwoPRep_11_11_sq TwoPRep_sh_11_11; + typedef GenIntegralSet_11_11 TwoBody_sh_11_11; + typedef typename OperType::Descriptor OperDescrType; + + const std::string task = task_label(label, deriv_level); + vector shells; unsigned int lmax = cparams->max_am(task); for (unsigned int l = 0; l <= lmax; l++) { @@ -1005,6 +1123,7 @@ void build_TwoPRep_2b_2k(std::ostream& os, taskmgr.current(task); iface->to_params(iface->macro_define(std::string("MAX_AM_") + task, lmax)); + const auto nullaux = typename TwoBody_sh_11_11::AuxIndexType(0u); // // Construct graphs for each desired target integral and // 1) generate source code for the found traversal path @@ -1019,12 +1138,18 @@ void build_TwoPRep_2b_2k(std::ostream& os, std::shared_ptr context(new CppCodeContext(cparams)); std::shared_ptr memman(new WorstFitMemoryManager()); + bool p1_p2_swappable = !std::is_same::value; + bool bra_ket_coswappable = std::is_same::value; + + // Note: la, lb, lc, ld generate code for chemist notation (ab|O|cd), where O + // is a two-body operator. for (unsigned int la = 0; la <= lmax; la++) { for (unsigned int lb = 0; lb <= lmax; lb++) { for (unsigned int lc = 0; lc <= lmax; lc++) { for (unsigned int ld = 0; ld <= lmax; ld++) { if (!ShellQuartetSetPredicate( - LIBINT_SHELL_SET)>::value(la, lb, lc, ld)) + LIBINT_SHELL_SET)>::value(la, lb, lc, ld, p1_p2_swappable, + bra_ket_coswappable)) continue; // std::shared_ptr tactic(new ParticleDirectionTactic(la+lb > @@ -1036,13 +1161,30 @@ void build_TwoPRep_2b_2k(std::ostream& os, const int lim = 1; if (!(la == lim && lb == lim && lc == lim && ld == lim)) continue; #endif + // this will hold all target shell sets + std::vector> targets; + + ///////////////////////////////// + // loop over operator components + ///////////////////////////////// + std::vector descrs(1); + if (std::is_same::value || + std::is_same::value) { + // reset descriptors array + descrs.resize(0); + // iterate over quaternion components + for (int p = 0; p != 4; ++p) { + descrs.emplace_back(make_descr(p)); + } + } // unroll only if max_am <= cparams->max_am_opt(task) using std::max; const unsigned int max_am = max(max(la, lb), max(lc, ld)); const bool need_to_optimize = (max_am <= cparams->max_am_opt(task)); + const auto nopers = descrs.size(); const bool need_to_unroll = - l_to_cgshellsize(la) * l_to_cgshellsize(lb) * + nopers * l_to_cgshellsize(la) * l_to_cgshellsize(lb) * l_to_cgshellsize(lc) * l_to_cgshellsize(ld) <= cparams->unroll_threshold(); const unsigned int unroll_threshold = @@ -1067,7 +1209,7 @@ void build_TwoPRep_2b_2k(std::ostream& os, //////////// // NB translational invariance is now handled by CR_DerivGauss CartesianDerivIterator<4> diter(deriv_level); - std::vector> targets; + bool last_deriv = false; do { CGShell a(la); @@ -1084,18 +1226,22 @@ void build_TwoPRep_2b_2k(std::ostream& os, } } - std::shared_ptr abcd = - TwoPRep_sh_11_11::Instance(a, b, c, d, mType(0u)); - targets.push_back(abcd); + // operator component loop + for (unsigned int op = 0; op != descrs.size(); ++op) { + OperType oper(descrs[op]); + + std::shared_ptr abcd = + TwoBody_sh_11_11::Instance(a, b, c, d, nullaux, oper); + targets.push_back(abcd); + } + last_deriv = diter.last(); if (!last_deriv) diter.next(); } while (!last_deriv); // append all derivatives as targets to the graph - for (std::vector>::const_iterator - t = targets.begin(); - t != targets.end(); ++t) { + for (auto it = targets.begin(); it != targets.end(); ++it) { std::shared_ptr t_ptr = - std::dynamic_pointer_cast(*t); + std::dynamic_pointer_cast(*it); dg_xxxx->append_target(t_ptr); } @@ -1107,33 +1253,45 @@ void build_TwoPRep_2b_2k(std::ostream& os, CGShell b(lb); CGShell c(lc); CGShell d(ld); - std::shared_ptr abcd = - TwoPRep_sh_11_11::Instance(a, b, c, d, mType(0u)); - abcd_label = abcd->label(); + + if constexpr (std::is_same::value) { + OperType oper; + oper = OperType(descrs[0]); + std::shared_ptr abcd = + TwoBody_sh_11_11::Instance(a, b, c, d, nullaux, oper); + abcd_label = abcd->label(); + } else { + std::ostringstream oss; + oss << "_" << a.label() << "_" << b.label(); + oss << "_" << label; + oss << "_" << c.label() << "_" << d.label(); + abcd_label = oss.str(); + } } // + derivative level (if deriv_level > 0) - std::string label; + std::string eval_label; { - label = ""; + eval_label = ""; if (deriv_level != 0) { std::ostringstream oss; oss << "deriv" << deriv_level; - label += oss.str(); + eval_label += oss.str(); } - label += abcd_label; + eval_label += abcd_label; } - std::cout << "working on " << label << " ... "; + std::cout << "working on " << eval_label << " ... "; std::cout.flush(); std::string prefix(cparams->source_directory()); std::deque decl_filenames; std::deque def_filenames; - // this will generate code for these targets, and potentially generate - // code for its prerequisites + // this will generate code for these targets, and potentially + // generate code for its prerequisites GenerateCode(dg_xxxx, context, cparams, strat, tactic, memman, - decl_filenames, def_filenames, prefix, label, false); + decl_filenames, def_filenames, prefix, eval_label, + false); // update max stack size and # of targets const std::shared_ptr& tparams = @@ -1147,16 +1305,14 @@ void build_TwoPRep_2b_2k(std::ostream& os, ostringstream oss; oss << context->label_to_function_name("libint2_build_" + task) << "[" << la << "][" << lb << "][" << lc << "][" << ld - << "] = " << context->label_to_function_name(label) + << "] = " << context->label_to_function_name(eval_label) << context->end_of_stat() << endl; iface->to_static_init(oss.str()); // need to declare this function internally - for (std::deque::const_iterator i = - decl_filenames.begin(); - i != decl_filenames.end(); ++i) { + for (auto& decl_filename : decl_filenames) { oss.str(""); - oss << "#include <" << *i << ">" << endl; + oss << "#include <" << decl_filename << ">" << endl; iface->to_int_iface(oss.str()); } @@ -1363,10 +1519,9 @@ void build_TwoPRep_1b_2k(std::ostream& os, iface->to_static_init(oss.str()); // need to declare this function internally - for (std::deque::const_iterator i = decl_filenames.begin(); - i != decl_filenames.end(); ++i) { + for (auto& decl_filename : decl_filenames) { oss.str(""); - oss << "#include <" << *i << ">" << endl; + oss << "#include <" << decl_filename << ">" << endl; iface->to_int_iface(oss.str()); } @@ -1375,7 +1530,7 @@ void build_TwoPRep_1b_2k(std::ostream& os, #endif dg_xxx->reset(); memman->reset(); - + std::cout << "done" << std::endl; } // end of d loop } // end of c loop } // end of bra loop @@ -1539,8 +1694,8 @@ void build_TwoPRep_1b_1k(std::ostream& os, std::deque decl_filenames; std::deque def_filenames; - // this will generate code for this targets, and potentially generate code - // for its prerequisites + // this will generate code for this targets, and potentially generate + // code for its prerequisites GenerateCode(dg_xxx, context, cparams, strat, tactic, memman, decl_filenames, def_filenames, prefix, label, false); @@ -1549,7 +1704,8 @@ void build_TwoPRep_1b_1k(std::ostream& os, taskmgr.current().params(); tparams->max_stack_size(max_am, memman->max_memory_used()); tparams->max_ntarget(targets.size()); - // os << " Max memory used = " << memman->max_memory_used() << std::endl; + // os << " Max memory used = " << memman->max_memory_used() << + // std::endl; // set pointer to the top-level evaluator function ostringstream oss; @@ -1560,10 +1716,9 @@ void build_TwoPRep_1b_1k(std::ostream& os, iface->to_static_init(oss.str()); // need to declare this function internally - for (std::deque::const_iterator i = decl_filenames.begin(); - i != decl_filenames.end(); ++i) { + for (auto& decl_filename : decl_filenames) { oss.str(""); - oss << "#include <" << *i << ">" << endl; + oss << "#include <" << decl_filename << ">" << endl; iface->to_int_iface(oss.str()); } @@ -1572,7 +1727,7 @@ void build_TwoPRep_1b_1k(std::ostream& os, #endif dg_xxx->reset(); memman->reset(); - + std::cout << "done" << std::endl; } // end of ket loop } // end of bra loop } @@ -1741,8 +1896,8 @@ void build_R12kG12_2b_2k(std::ostream& os, std::deque decl_filenames; std::deque def_filenames; - // this will generate code for this targets, and potentially generate - // code for its prerequisites + // this will generate code for this targets, and potentially + // generate code for its prerequisites GenerateCode(dg_xxxx, context, cparams, strat, tactic, memman, decl_filenames, def_filenames, prefix, label, false); @@ -2094,11 +2249,11 @@ void build_G12DKH_2b_2k(std::ostream& os, oss << "#include <" << decl_filename << ">" << endl; iface->to_int_iface(oss.str()); - // For the most expensive (i.e. presumably complete) graph extract all - // precomputed quantities -- these will be members of the evaluator - // structure also extract all RRs -- need to keep track of these to - // figure out which external symbols appearing in RR code belong to - // this task also + // For the most expensive (i.e. presumably complete) graph extract + // all precomputed quantities -- these will be members of the + // evaluator structure also extract all RRs -- need to keep track of + // these to figure out which external symbols appearing in RR code + // belong to this task also if (la == lmax && lb == lmax && lc == lmax && ld == lmax) extract_symbols(dg_xxxx); @@ -2136,6 +2291,12 @@ void config_to_api(const std::shared_ptr& cparams, iface->to_params(iface->macro_define("DERIV_ERI_ORDER", LIBINT_INCLUDE_ERI)); max_deriv_order = std::max(max_deriv_order, LIBINT_INCLUDE_ERI); #endif +#ifdef LIBINT_INCLUDE_RKB_ERI + iface->to_params(iface->macro_define("SUPPORT_RKB_ERI", 1)); + iface->to_params( + iface->macro_define("DERIV_RKB_ERI_ORDER", LIBINT_INCLUDE_RKB_ERI)); + max_deriv_order = std::max(max_deriv_order, LIBINT_INCLUDE_RKB_ERI); +#endif #ifdef LIBINT_INCLUDE_ERI3 iface->to_params(iface->macro_define("SUPPORT_ERI3", 1)); iface->to_params( @@ -2165,7 +2326,8 @@ void config_to_api(const std::shared_ptr& cparams, // generated tasks declare all tasks in a range of valid tasks as defined or // not LibraryTaskManager& taskmgr = LibraryTaskManager::Instance(); - // the range is defined by max # of centers, max deriv order, and operator set + // the range is defined by max # of centers, max deriv order, and operator + // set const size_t max_ncenter = 4; for (unsigned int ncenter = 0; ncenter <= max_ncenter; ++ncenter) { std::stringstream oss; @@ -2191,8 +2353,9 @@ void config_to_api(const std::shared_ptr& cparams, { // 2-body ints -#define BOOST_PP_TWOBODY_TASKOPER_TUPLE \ - ("eri", "r12kg12", "r12_0_g12", "r12_2_g12", "g12_T1_g12", "g12dkh") +#define BOOST_PP_TWOBODY_TASKOPER_TUPLE \ + ("eri", "coulomb_opop", "opop_coulomb_opop", "r12kg12", "r12_0_g12", \ + "r12_2_g12", "g12_T1_g12", "g12dkh") #define BOOST_PP_TWOBODY_TASKOPER_LIST \ BOOST_PP_TUPLE_TO_LIST(BOOST_PP_TWOBODY_TASKOPER_TUPLE) diff --git "a/src/bin/libint/comp_11_Coulomb\317\203p\317\203p_11.h" "b/src/bin/libint/comp_11_Coulomb\317\203p\317\203p_11.h" new file mode 100644 index 000000000..315135994 --- /dev/null +++ "b/src/bin/libint/comp_11_Coulomb\317\203p\317\203p_11.h" @@ -0,0 +1,164 @@ +/* + * Copyright (C) 2004-2026 Edward F. Valeev + * + * This file is part of Libint compiler. + * + * Libint compiler is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * Libint compiler is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with Libint compiler. If not, see . + * + */ + +#ifndef LIBINT_COMP_11_COULOMBΣPΣP_11_H +#define LIBINT_COMP_11_COULOMBΣPΣP_11_H + +#include +#include +#include + +namespace libint2 { + +/** + * this computes integral of + * \f$ \frac{1}{r_{ij}} \sigma \cdot \hat{p}_1 \sigma \cdot \hat{p}_2 \f$ over + * CGShell/CGF by rewriting it as a linear combination of integrals over + * derivatives of \frac{1}{r_{ij}} + * @tparam F basis function type. valid choices are CGShell or CGF + */ +template +class CR_11_Coulombσpσp_11 + : public GenericRecurrenceRelation< + CR_11_Coulombσpσp_11, F, + GenIntegralSet_11_11> { + public: + typedef CR_11_Coulombσpσp_11 ThisType; + typedef F BasisFunctionType; + typedef CoulombσpσpOper OperType; + typedef GenIntegralSet_11_11 TargetType; + typedef GenericRecurrenceRelation + ParentType; + friend class GenericRecurrenceRelation; + static const unsigned int max_nchildren = 100; // TODO figure out + + using ParentType::Instance; + + static bool directional() { return false; } + + private: + using ParentType::is_simple; + using ParentType::target_; + using ParentType::RecurrenceRelation::expr_; + using ParentType::RecurrenceRelation::nflops_; + + /// Constructor is private, used by ParentType::Instance that maintains + /// registry of these objects + CR_11_Coulombσpσp_11(const std::shared_ptr &, unsigned int = 0); + + static std::string descr() { return "CR"; } +}; + +template +CR_11_Coulombσpσp_11::CR_11_Coulombσpσp_11( + const std::shared_ptr &Tint, unsigned int) + : ParentType(Tint, 0) { + assert(Tint->num_func_bra(/* particle */ 0) == 1); + assert(Tint->num_func_bra(/* particle */ 1) == 1); + assert(Tint->num_func_ket(/* particle */ 0) == 1); + assert(Tint->num_func_ket(/* particle */ 1) == 1); + + F a(Tint->bra(0, 0)); + F b(Tint->ket(0, 0)); + F c(Tint->bra(1, 0)); + F d(Tint->ket(1, 0)); + + const auto &oper = Tint->oper(); + + if (a.contracted() || b.contracted() || c.contracted() || d.contracted()) + return; + + using namespace libint2::algebra; + using namespace libint2::prefactor; + using libint2::algebra::operator*; + + const mType zero_m(0u); + + ChildFactory> + factory(this); + + constexpr auto x = 0; + constexpr auto y = 1; + constexpr auto z = 2; + + F c_x{c}; + c_x.deriv().inc(x); // d(c)/dx = c_x + F c_y{c}; + c_y.deriv().inc(y); // d(c)/dy = c_y + F c_z{c}; + c_z.deriv().inc(z); // d(c)/dz = c_z + + F d_x{d}; + d_x.deriv().inc(x); // d(d)/dx = d_x + F d_y{d}; + d_y.deriv().inc(y); // d(d)/dy = d_y + F d_z{d}; + d_z.deriv().inc(z); // d(d)/dz = d_z + + // Component wise generation for quaternion ( a b | 1/r12 | (σ.p) c (σ.p) d ) + switch (oper->descr().quaternion_index()) { + case 0: { + // zeroth component = (a b | c_x d_x) + (a b | c_y d_y) + (a b | c_z d_z) + auto a_b_cx_dx = factory.make_child(a, b, c_x, d_x, zero_m); + auto a_b_cy_dy = factory.make_child(a, b, c_y, d_y, zero_m); + auto a_b_cz_dz = factory.make_child(a, b, c_z, d_z, zero_m); + if (is_simple()) { + expr_ = a_b_cx_dx + a_b_cy_dy + a_b_cz_dz; + nflops_ += 2; + } + } break; + case 1: { + // x component = (a b | c_y d_z) - (a b | c_z d_y) + auto a_b_cy_dz = factory.make_child(a, b, c_y, d_z, zero_m); + auto a_b_cz_dy = factory.make_child(a, b, c_z, d_y, zero_m); + if (is_simple()) { + expr_ = a_b_cy_dz - a_b_cz_dy; + nflops_ += 1; + } + } break; + case 2: { + // y component = (a b | c_z d_x) - (a b | c_x d_z) + auto a_b_cz_dx = factory.make_child(a, b, c_z, d_x, zero_m); + auto a_b_cx_dz = factory.make_child(a, b, c_x, d_z, zero_m); + if (is_simple()) { + expr_ = a_b_cz_dx - a_b_cx_dz; + nflops_ += 1; + } + } break; + case 3: { + // z component = (a b | c_x d_y) - (a b | c_y d_x) + auto a_b_cx_dy = factory.make_child(a, b, c_x, d_y, zero_m); + auto a_b_cy_dx = factory.make_child(a, b, c_y, d_x, zero_m); + if (is_simple()) { + expr_ = a_b_cx_dy - a_b_cy_dx; + nflops_ += 1; + } + } break; + default: + throw std::runtime_error( + "CR_11_Coulombσpσp_11: invalid quaternionic index"); + } + +} // CR_11_Coulombσpσp_11::CR_11_Coulombσpσp_11 +}; // namespace libint2 + +#endif // LIBINT_COMP_11_COULOMBΣPΣP_11_H diff --git "a/src/bin/libint/comp_11_\317\203p\317\203pCoulomb\317\203p\317\203p_11.h" "b/src/bin/libint/comp_11_\317\203p\317\203pCoulomb\317\203p\317\203p_11.h" new file mode 100644 index 000000000..7bf0a4b9a --- /dev/null +++ "b/src/bin/libint/comp_11_\317\203p\317\203pCoulomb\317\203p\317\203p_11.h" @@ -0,0 +1,196 @@ +/* + * Copyright (C) 2004-2026 Edward F. Valeev + * + * This file is part of Libint compiler. + * + * Libint compiler is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * Libint compiler is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with Libint compiler. If not, see . + * + */ + +#ifndef LIBINT_COMP_11_ΣPΣPCOULOMBΣPΣP_11_H +#define LIBINT_COMP_11_ΣPΣPCOULOMBΣPΣP_11_H + +#include +#include +#include + +namespace libint2 { + +/** + * this computes integral of + * \sigma \cdot \hat{p}_1 \sigma \cdot \hat{p}_2 \f$ \frac{1}{r_{ij}} \sigma + * \cdot \hat{p}_3 \sigma \cdot \hat{p}_4 \f$ over CGShell/CGF by rewriting it + * as a linear combination of integrals over derivatives of \frac{1}{r_{ij}} + * @tparam F basis function type. valid choices are CGShell or CGF + */ +template +class CR_11_σpσpCoulombσpσp_11 + : public GenericRecurrenceRelation< + CR_11_σpσpCoulombσpσp_11, F, + GenIntegralSet_11_11> { + public: + typedef CR_11_σpσpCoulombσpσp_11 ThisType; + typedef F BasisFunctionType; + typedef σpσpCoulombσpσpOper OperType; + typedef GenIntegralSet_11_11 TargetType; + typedef GenericRecurrenceRelation + ParentType; + friend class GenericRecurrenceRelation; + static const unsigned int max_nchildren = 100; // TODO figure out + + using ParentType::Instance; + + static bool directional() { return false; } + + private: + using ParentType::is_simple; + using ParentType::target_; + using ParentType::RecurrenceRelation::expr_; + using ParentType::RecurrenceRelation::nflops_; + + /// Constructor is private, used by ParentType::Instance that maintains + /// registry of these objects + CR_11_σpσpCoulombσpσp_11(const std::shared_ptr &, + unsigned int = 0); + + static std::string descr() { return "CR"; } +}; + +template +CR_11_σpσpCoulombσpσp_11::CR_11_σpσpCoulombσpσp_11( + const std::shared_ptr &Tint, unsigned int) + : ParentType(Tint, 0) { + assert(Tint->num_func_bra(/* particle */ 0) == 1); + assert(Tint->num_func_bra(/* particle */ 1) == 1); + assert(Tint->num_func_ket(/* particle */ 0) == 1); + assert(Tint->num_func_ket(/* particle */ 1) == 1); + + F a(Tint->bra(0, 0)); + F b(Tint->ket(0, 0)); + F c(Tint->bra(1, 0)); + F d(Tint->ket(1, 0)); + + const auto &oper = Tint->oper(); + + if (a.contracted() || b.contracted() || c.contracted() || d.contracted()) + return; + + using namespace libint2::algebra; + using namespace libint2::prefactor; + using libint2::algebra::operator*; + + const mType zero_m(0u); + + ChildFactory> + factory(this); + + constexpr auto x = 0; + constexpr auto y = 1; + constexpr auto z = 2; + + auto mc = [&](const int r1, const int r2, const int r3, const int r4) { + F a_r1{a}; + a_r1.deriv().inc(r1); + F b_r2{b}; + b_r2.deriv().inc(r2); + F c_r3{c}; + c_r3.deriv().inc(r3); + F d_r4{d}; + d_r4.deriv().inc(r4); + return factory.make_child(a_r1, b_r2, c_r3, d_r4, zero_m); + }; + + // Component wise generation for quaternion : + // ( (σ.p) a (σ.p) b | 1/r12 | (σ.p) c (σ.p) d ) + switch (oper->descr().quaternion_index()) { + case 0: { + auto xxxx = mc(x, x, x, x); + auto yyxx = mc(y, y, x, x); + auto zzxx = mc(z, z, x, x); + auto yxyx = mc(y, x, y, x); + auto xyyx = mc(x, y, y, x); + auto yxxy = mc(y, x, x, y); + auto xyxy = mc(x, y, x, y); + auto xxyy = mc(x, x, y, y); + auto yyyy = mc(y, y, y, y); + auto zzyy = mc(z, z, y, y); + auto xxzz = mc(x, x, z, z); + auto yyzz = mc(y, y, z, z); + auto zzzz = mc(z, z, z, z); + if (is_simple()) { + expr_ = xxxx + yyxx + zzxx - yxyx + xyyx + yxxy - xyxy + xxyy + yyyy + + zzyy + xxzz + yyzz + zzzz; + nflops_ += 12; + } + } break; + case 1: { + auto zxzx = mc(z, x, z, x); + auto xzzx = mc(x, z, z, x); + auto zyzy = mc(z, y, z, y); + auto yzzy = mc(y, z, z, y); + auto zxxz = mc(z, x, x, z); + auto xzxz = mc(x, z, x, z); + auto zyyz = mc(z, y, y, z); + auto yzyz = mc(y, z, y, z); + if (is_simple()) { + expr_ = zxzx - xzzx - zyzy + yzzy - zxxz + xzxz + zyyz - yzyz; + nflops_ += 7; + } + } break; + case 2: { + auto zyzx = mc(z, y, z, x); + auto yzzx = mc(y, z, z, x); + auto zxzy = mc(z, x, z, y); + auto xzzy = mc(x, z, z, y); + auto zyxz = mc(z, y, x, z); + auto yzxz = mc(y, z, x, z); + auto zxyz = mc(z, x, y, z); + auto xzyz = mc(x, z, y, z); + if (is_simple()) { + // swapped order of first two terms compiler does not like negative sign + // in front of first term + expr_ = yzzx - zyzx - zxzy + xzzy + zyxz - yzxz + zxyz - xzyz; + nflops_ += 7; + } + } break; + case 3: { + auto yxxx = mc(y, x, x, x); + auto xyxx = mc(x, y, x, x); + auto xxyx = mc(x, x, y, x); + auto yyyx = mc(y, y, y, x); + auto zzyx = mc(z, z, y, x); + auto xxxy = mc(x, x, x, y); + auto yyxy = mc(y, y, x, y); + auto zzxy = mc(z, z, x, y); + auto yxyy = mc(y, x, y, y); + auto xyyy = mc(x, y, y, y); + auto yxzz = mc(y, x, z, z); + auto xyzz = mc(x, y, z, z); + if (is_simple()) { + expr_ = xyxx - yxxx - xxyx - yyyx - zzyx + xxxy + yyxy + zzxy - yxyy + + xyyy - yxzz + xyzz; + nflops_ += 11; + } + } break; + default: + throw std::runtime_error( + "CR_11_σpσpCoulombσpσp_11: invalid quaternionic index"); + } + +} // CR_11_σpσpCoulombσpσp_11::CR_11_σpσpCoulombσpσp_11 +}; // namespace libint2 + +#endif // LIBINT_COMP_11_ΣPΣPCOULOMBΣPΣP_11_H diff --git "a/src/bin/libint/comp_1_\317\203pV\317\203p_1.h" "b/src/bin/libint/comp_1_\317\203pV\317\203p_1.h" index cb131ebb8..9fbdef361 100644 --- "a/src/bin/libint/comp_1_\317\203pV\317\203p_1.h" +++ "b/src/bin/libint/comp_1_\317\203pV\317\203p_1.h" @@ -107,7 +107,7 @@ CR_1_σpVσp_1::CR_1_σpVσp_1(const std::shared_ptr &Tint, // (a|W0|b) = (d a/dAx | V | d b/dBx) + (d a/dAy | V | d b/dBy) + (d a/dAz | V // | d b/dBz) - switch (oper->descr().pauli_index()) { + switch (oper->descr().quaternion_index()) { case 0: { auto Dx_a_V_Dx_b = factory.make_child(Dx_a, Dx_b, zero_m); auto Dy_a_V_Dy_b = factory.make_child(Dy_a, Dy_b, zero_m); @@ -146,7 +146,7 @@ CR_1_σpVσp_1::CR_1_σpVσp_1(const std::shared_ptr &Tint, } } break; default: - throw std::runtime_error("CR_1_σpVσp_1: invalid Pauli index"); + throw std::runtime_error("CR_1_σpVσp_1: invalid quaternionic index"); } } // CR_1_σpVσp_1::CR_1_σpVσp_1 diff --git a/src/bin/libint/gauss.cc b/src/bin/libint/gauss.cc index 3899283ab..aa60b1de7 100644 --- a/src/bin/libint/gauss.cc +++ b/src/bin/libint/gauss.cc @@ -115,8 +115,15 @@ std::string CGF::label() const { unsigned int am = qn_[0] + qn_[1] + qn_[2]; std::string deriv_label; if (!deriv_.zero()) deriv_label = deriv_.label(); - const std::string am_string = am_to_symbol(am, contracted()); + std::string am_string = am_to_symbol(am, contracted()); std::ostringstream oss; + + // Some OSs can have case-insensitive filesystem e.g., MacOS. So here we add + // additional identifier for primitive function labels + if (!this->contracted()) { + am_string += "_"; + } + oss << (pure_sh_ && am > 0 ? "W" : "") << am_string << deriv_label << "_"; if (am == 0) return oss.str(); @@ -223,8 +230,16 @@ CGShell::~CGShell() {} std::string CGShell::label() const { if (is_unit()) return "unit"; - std::string result = std::string(pure_sh_ && qn_[0] > 0 ? "W" : "") + - am_to_symbol(qn_[0], contracted()); + std::string am_symbol = am_to_symbol(qn_[0], contracted()); + + // Some OSs can have case-insensitive filesystem e.g., MacOS. So here we add + // additional identifier for primitive shell labels + if (!this->contracted()) { + am_symbol += "_"; + } + + std::string result = + std::string(pure_sh_ && qn_[0] > 0 ? "W" : "") + am_symbol; if (!deriv_.zero()) result += deriv_.label(); return result; } diff --git a/src/bin/libint/master_ints_list.h b/src/bin/libint/master_ints_list.h index 1aa8c3e64..37bfa29a7 100644 --- a/src/bin/libint/master_ints_list.h +++ b/src/bin/libint/master_ints_list.h @@ -106,6 +106,13 @@ typedef GenIntegralSet_1_1, CartesianMultipoleOper<1u>, ////////////////////////// typedef GenIntegralSet_11_11 TwoPRep_11_11_sq; typedef GenIntegralSet_11_11 TwoPRep_11_11_int; +typedef GenIntegralSet_11_11 + Coulombσpσp_11_11_sq; +typedef GenIntegralSet_11_11 Coulombσpσp_11_11_int; +typedef GenIntegralSet_11_11 + σpσpCoulombσpσp_11_11_sq; +typedef GenIntegralSet_11_11 + σpσpCoulombσpσp_11_11_int; typedef GenIntegralSet_11_11 R12kG12_11_11_sq; typedef GenIntegralSet_11_11 R12kG12_11_11_int; typedef GenIntegralSet_11_11 @@ -144,10 +151,11 @@ typedef boost::mpl::list< CMultipole_1_1_int_y, CMultipole_1_1_int_z, SMultipole_1_1_sh, SMultipole_1_1_int, #endif - TwoPRep_11_11_sq, TwoPRep_11_11_int, R12kG12_11_11_sq, R12kG12_11_11_int, - R12kR12lG12_11_11_sq, R12kR12lG12_11_11_int, TiG12_11_11_sq, - TiG12_11_11_int, G12TiG12_11_11_sq, G12TiG12_11_11_int, - DivG12prime_xTx_11_11_sq, DivG12prime_xTx_11_11_int, + TwoPRep_11_11_sq, TwoPRep_11_11_int, Coulombσpσp_11_11_sq, + Coulombσpσp_11_11_int, σpσpCoulombσpσp_11_11_sq, σpσpCoulombσpσp_11_11_int, + R12kG12_11_11_sq, R12kG12_11_11_int, R12kR12lG12_11_11_sq, + R12kR12lG12_11_11_int, TiG12_11_11_sq, TiG12_11_11_int, G12TiG12_11_11_sq, + G12TiG12_11_11_int, DivG12prime_xTx_11_11_sq, DivG12prime_xTx_11_11_int, DummySymmIntegral_11_11_sq, DummySymmIntegral_11_11_int> MasterIntegralTypeList; diff --git a/src/bin/libint/master_rrs_list.h b/src/bin/libint/master_rrs_list.h index f3ec4e2d2..d55cfa301 100644 --- a/src/bin/libint/master_rrs_list.h +++ b/src/bin/libint/master_rrs_list.h @@ -21,10 +21,12 @@ #ifndef _libint2_src_bin_libint_masterrrslist_h_ #define _libint2_src_bin_libint_masterrrslist_h_ +#include #include #include #include #include +#include #include #include #include @@ -266,6 +268,11 @@ typedef CR_DerivGauss Deriv_d_11_TwoPRep_11_int; +typedef CR_11_Coulombσpσp_11 CR_11_Coulombσpσp_11_sh; +typedef CR_11_Coulombσpσp_11 CR_11_Coulombσpσp_11_int; + +typedef CR_11_σpσpCoulombσpσp_11 CR_11_σpσpCoulombσpσp_11_sh; +typedef CR_11_σpσpCoulombσpσp_11 CR_11_σpσpCoulombσpσp_11_int; }; // namespace libint2 #endif // header guard diff --git a/src/bin/libint/oper.h b/src/bin/libint/oper.h index 091df3ac0..cecbda72b 100644 --- a/src/bin/libint/oper.h +++ b/src/bin/libint/oper.h @@ -289,23 +289,23 @@ BOOST_PP_LIST_FOR_EACH(BOOST_PP_DECLARE_HERMITIAN_ONEBODY_DESCRIPTOR, struct σpVσp_Descr : public Contractable<σpVσp_Descr> { typedef MultiplicativeODep1Body_Props Properties; - σpVσp_Descr() : pauli_index_(0) {} - σpVσp_Descr(int pauli_index) : pauli_index_(pauli_index) { - assert(pauli_index <= 3); + σpVσp_Descr() : quaternion_index_(0) {} + σpVσp_Descr(int quaternion_index) : quaternion_index_(quaternion_index) { + assert(quaternion_index <= 3); } static const unsigned int max_key = 4; - unsigned int key() const { return pauli_index(); } + unsigned int key() const { return quaternion_index(); } std::string description() const { std::string descr("opVop["); - if (pauli_index() == 0) + if (quaternion_index() == 0) descr += "0"; - else if (pauli_index() == 1) - descr += "Z"; - else if (pauli_index() == 2) + else if (quaternion_index() == 1) descr += "X"; - else if (pauli_index() == 3) + else if (quaternion_index() == 2) descr += "Y"; + else if (quaternion_index() == 3) + descr += "Z"; else abort(); return descr + "]"; @@ -314,10 +314,10 @@ struct σpVσp_Descr : public Contractable<σpVσp_Descr> { int psymm(int i, int j) const { abort(); } int hermitian(int i) const { return +1; } - int pauli_index() const { return pauli_index_; } + int quaternion_index() const { return quaternion_index_; } private: - const int pauli_index_ = -1; + const int quaternion_index_ = -1; }; typedef GenOper<σpVσp_Descr> σpVσpOper; @@ -399,6 +399,80 @@ struct TwoPRep_Descr : public Contractable { }; typedef GenOper TwoPRep; +/** Coulombσpσp is the two-body repulsion operator. + */ +struct Coulombσpσp_Descr : public Contractable { + typedef MultiplicativeSymm2Body_Props Properties; + + Coulombσpσp_Descr() : quaternion_index_(0) {} + Coulombσpσp_Descr(int quaternion_index) + : quaternion_index_(quaternion_index) { + assert(quaternion_index <= 3); + } + + static const unsigned int max_key = 4; + unsigned int key() const { return quaternion_index(); } + std::string description() const { + std::string descr("coulomb_opop["); + if (quaternion_index() == 0) + descr += "0"; + else if (quaternion_index() == 1) + descr += "X"; + else if (quaternion_index() == 2) + descr += "Y"; + else if (quaternion_index() == 3) + descr += "Z"; + else + abort(); + return descr + "]"; + } + std::string label() const { return description(); } + int psymm(int i, int j) const { abort(); } + int hermitian(int i) const { return +1; } + + int quaternion_index() const { return quaternion_index_; } + + private: + const int quaternion_index_ = -1; +}; +typedef GenOper CoulombσpσpOper; + +struct σpσpCoulombσpσp_Descr : public Contractable<σpσpCoulombσpσp_Descr> { + typedef MultiplicativeSymm2Body_Props Properties; + + σpσpCoulombσpσp_Descr() : quaternion_index_(0) {} + σpσpCoulombσpσp_Descr(int quaternion_index) + : quaternion_index_(quaternion_index) { + assert(quaternion_index <= 3); + } + + static const unsigned int max_key = 4; + unsigned int key() const { return quaternion_index(); } + std::string description() const { + std::string descr("opop_coulomb_opop["); + if (quaternion_index() == 0) + descr += "0"; + else if (quaternion_index() == 1) + descr += "X"; + else if (quaternion_index() == 2) + descr += "Y"; + else if (quaternion_index() == 3) + descr += "Z"; + else + abort(); + return descr + "]"; + } + std::string label() const { return description(); } + int psymm(int i, int j) const { abort(); } + int hermitian(int i) const { return +1; } + + int quaternion_index() const { return quaternion_index_; } + + private: + const int quaternion_index_ = -1; +}; +typedef GenOper<σpσpCoulombσpσp_Descr> σpσpCoulombσpσpOper; + /** GTG_1d is the two-body 1-dimensional Gaussian geminal */ struct GTG_1d_Descr : public Contractable { diff --git a/src/bin/libint/strategy.cc b/src/bin/libint/strategy.cc index bbfc3fe21..58fb8d2bd 100644 --- a/src/bin/libint/strategy.cc +++ b/src/bin/libint/strategy.cc @@ -115,6 +115,23 @@ struct MasterStrategy { }; #endif +template <> +struct MasterStrategy { + typedef boost::mpl::list value; +}; +template <> +struct MasterStrategy { + typedef boost::mpl::list value; +}; +template <> +struct MasterStrategy<σpσpCoulombσpσp_11_11_sq> { + typedef boost::mpl::list value; +}; +template <> +struct MasterStrategy<σpσpCoulombσpσp_11_11_int> { + typedef boost::mpl::list value; +}; + #if LIBINT_SHELLQUARTET_STRATEGY == LIBINT_SHELLQUARTET_STRATEGY_A0C0 template <> struct MasterStrategy {