added gaunt and 3j coefficients

CMakeLists.txt

@@ -101,6 +101,11 @@ if (${LBASIS_COMPONENT})
    add_subdirectory(GX-LocalizedBasis)
 endif()
 
+option (NAO_COMPONENT "Enable the localized basis component" ON) 
+if (${NAO_COMPONENT})
+   add_subdirectory(GX-NAOIntegrals)
+endif()
+
 # Elements to compile the libraries inside submodules
 option(COMPILE_SUBMODULES "Compile GreenX component contained inside submodules" ON)
 if (${COMPILE_SUBMODULES})

GX-NAOIntegrals/CMakeLists.txt

@@ -0,0 +1,75 @@
+# ***************************************************************************************************
+#  Copyright (C) 2020-2024 GreenX library
+#  This file is distributed under the terms of the APACHE2 License.
+#
+# ***************************************************************************************************
+
+# -----------------------------------------------
+# CMakeLists for NAO Integrals Library
+# -----------------------------------------------
+add_library(LibGXNAO "")
+
+set_target_properties(LibGXNAO
+  PROPERTIES
+  VERSION 0.0.1
+  SOVERSION 0.0.1
+  LIBRARY_OUTPUT_NAME gx_nao
+  ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_Fortran_LIB_DIRECTORY}
+  LIBRARY_OUTPUT_DIRECTORY ${CMAKE_Fortran_LIB_DIRECTORY})
+
+target_include_directories(LibGXNAO PUBLIC src/)
+
+target_sources(LibGXNAO PRIVATE src/wigner.f90 src/gaunt.f90)
+target_link_libraries(LibGXNAO GXCommon)
+
+# -----------------------------------------------
+# Library Installation
+# -----------------------------------------------
+# Install library
+# Destination relative to ${CMAKE_INSTALL_PREFIX}, defined in top-level CMake
+install(TARGETS LibGXNAO ARCHIVE DESTINATION lib LIBRARY DESTINATION lib)
+
+# Install modules
+# Destination relative to ${CMAKE_INSTALL_PREFIX}, defined in top-level CMake
+install(DIRECTORY ${CMAKE_Fortran_MODULE_DIRECTORY}  DESTINATION include)
+
+# -----------------------------------------------
+# Application Testing Set-Up
+# -----------------------------------------------
+# Include cmake custom function
+include(../cmake/testFunctions.cmake)
+
+# Set name of test sub-directory in the build directory
+set(TEST_TARGET_DIR "nao_integrals")
+
+target_include_directories(LibGXNAO PUBLIC test/)
+
+# -----------------------------------------------
+# Application Tests
+# -----------------------------------------------
+
+# Add a test target
+add_executable(test_nao)
+
+# Set binary name
+set_target_properties(test_nao
+        PROPERTIES
+        RUNTIME_OUTPUT_NAME test_nao.exe)
+
+# Define source that comprise the binary
+target_sources(test_nao
+        PRIVATE
+        test/test_nao.f90
+        )
+
+# Libraries that the binary links to
+target_link_libraries(test_nao
+        PUBLIC
+	    LibGXNAO
+        )
+
+# Build location of the binary
+# TODO(Alex) Consider move this to test/, to sit with the python drivers
+set_target_properties(test_nao
+        PROPERTIES
+        RUNTIME_OUTPUT_DIRECTORY ${CMAKE_Fortran_BIN_DIRECTORY})

GX-NAOIntegrals/src/gaunt.f90

@@ -0,0 +1,142 @@
+!> @brief calculation of R-Gaunt coefficients
+!!
+!! see Talman, Int. J. of Quant. Chem., 111, 2221–2227 (2011)
+!! and  Homeier, H. H. H.; Steinborn, E. O. JMol Struct (THEOCHEM) 1996, 368, 31.
+module gaunt
+
+    use wigner, only: threej
+    use kinds, only: dp
+    use constants, only: pi
+
+    implicit none
+
+    private
+    public :: r_gaunt
+
+    contains
+
+    !> @prief get the R-Gaunt coefficient
+    !!
+    !! @param[in] l1 -- angular quantum number 
+    !! @param[in] l2 -- angular quantum number 
+    !! @param[in] l3 -- angular quantum number 
+    !! @param[in] m1 -- magnetic quantum number 
+    !! @param[in] m2 -- magnetic quantum number 
+    !! @param[in] m3 -- magnetic quantum number 
+    !! @return        R-Gaunt coefficient
+    real(kind=dp) function r_gaunt(l1, l2, l3, m1, m2, m3) result(coeff)
+        integer, intent(in) :: l1, l2, l3, m1, m2, m3 
+
+        ! internal variables
+        integer :: sign
+        integer :: l1_aux, l2_aux, l3_aux, m1_aux, m2_aux, m3_aux 
+        integer :: m1_abs, m2_abs, m3_abs 
+
+        l1_aux = l1
+        l2_aux = l2
+        l3_aux = l3
+        m1_aux = m1
+        m2_aux = m2
+        m3_aux = m3
+        m1_abs = abs(m1)
+        m2_abs = abs(m2)
+        m3_abs = abs(m3)
+
+        ! create permutation to ensure m1_abs < m2_abs < m3_abs 
+        call order_for_m_abs(l1_aux, l2_aux, l3_aux, m1_aux, m2_aux, m3_aux, m1_abs, m2_abs, m3_abs)
+
+        if ((m1_abs.eq.0) .and. (m2_aux.ne.m3_aux)) then 
+            coeff = 0.0_dp
+        else if ((m1_abs.eq.0) .and. (m2_aux.eq.m3_aux)) then
+            coeff = 4.0_dp * pi * solid_harm_prefactor(l1_aux)   &
+                    * solid_harm_prefactor(l2_aux)          &
+                    * solid_harm_prefactor(l3_aux)          &
+                    * alpha_function(l1_aux, l2_aux, l3_aux)        &
+                    * wigner_3j(l1_aux, l2_aux, l3_aux, 0, m2_aux, -m3_aux)
+        else 
+            if (m1_aux+m2_aux+m3_aux .gt. 0) then 
+                sign = 1
+            else 
+                sign = -1
+            end if 
+            coeff = sign * 4.0_dp * pi * 1.0_dp/sqrt(2.0_dp) &
+                    * solid_harm_prefactor(l1_aux)          &
+                    * solid_harm_prefactor(l2_aux)          &
+                    * solid_harm_prefactor(l3_aux)          &
+                    * alpha_function(l1_aux, l2_aux, l3_aux)        &
+                    * wigner_3j(l1_aux, l2_aux, l3_aux, -m1_abs, -m2_abs, m3_abs)
+            ! consistency check 
+            if (m3_abs .ne. (m1_abs+m2_abs)) then 
+                print *, "inconsistency in R-Gaunt routine"
+                stop
+            end if 
+        end if 
+
+    end function r_gaunt
+
+
+
+    !> @brief reorder quantum numbers so that m1_abs < m2_abs < m3_abs
+    !!
+    !! @param[inout] l -- angular quantum number
+    !! @param[inout] m -- magnetic quantum number
+    !! @param[inout] m_abs -- abs(m)
+    subroutine order_for_m_abs(l1, l2, l3, m1, m2, m3, m1_abs, m2_abs, m3_abs)
+        integer, intent(inout) :: l1, l2, l3, m1, m2, m3, m1_abs, m2_abs, m3_abs
+
+        ! internal variables
+        integer :: min_idx, mid_idx, max_idx
+        integer, dimension(3):: l_tmp, m_tmp, m_abs_tmp
+
+        if ((m1_abs .le. m2_abs) .and. (m2_abs .le. m3_abs)) return 
+
+        l_tmp = (/l1, l2, l3/)
+        m_tmp = (/m1, m2, m3/)
+        m_abs_tmp = (/m1_abs, m2_abs, m3_abs/)
+
+        max_idx = maxval(m_abs_tmp)
+        min_idx = minval(m_abs_tmp)
+        mid_idx = 6 - (max_idx + min_idx)
+
+        l1 = l_tmp(min_idx)
+        l2 = l_tmp(mid_idx)
+        l3 = l_tmp(max_idx)
+
+        m1 = m_tmp(min_idx)
+        m2 = m_tmp(mid_idx)
+        m3 = m_tmp(max_idx)
+
+        m1_abs = m_abs_tmp(min_idx)
+        m2_abs = m_abs_tmp(mid_idx)
+        m3_abs = m_abs_tmp(max_idx)
+    end subroutine order_for_m_abs
+
+    !> @brief returns the prefactor of a solid harmonic function Y_lm
+    !!
+    !! @param[in] l -- angular quantum number
+    !! @result         prefator
+    real(kind=dp) function solid_harm_prefactor(l) result(prefac)
+        integer, intent(in) :: l 
+        prefac = sqrt((2*l + 1)/(4*pi))
+    end function solid_harm_prefactor
+
+    !> @brief get the wigner-3j symbol where all magnetic quantum n. are zero
+    !!
+    !! @param[in] l -- angular quantum number
+    !! @return         3j symbol
+    real(kind=dp) function alpha_function(l1, l2, l3) result(alpha)
+        integer, intent(in) :: l1, l2, l3 
+        alpha = threej(2*l1, 2*l2, 2*l3, 0, 0, 0)
+    end function alpha_function
+
+    !> @brief interface to the wigner-3j module where the input is defined as 2*l and 2*m
+    !!
+    !! @param[inout] l -- angular quantum number
+    !! @param[inout] m -- magnetic quantum number
+    !! @return         3j symbol
+    real(kind=dp) function wigner_3j(l1, l2, l3, m1, m2, m3) result(symbol)
+        integer :: l1, l2, l3, m1, m2, m3
+        symbol = threej(2*l1, 2*l2, 2*l3, 2*m1, 2*m2, 2*m3)
+    end function wigner_3j
+
+end module gaunt

GX-NAOIntegrals/src/logsbt.f90

@@ -0,0 +1,8 @@
+module logsbt
+
+    implicit none 
+
+    contains
+
+
+end module logsbt