!> @brief Cartesian -> pure spherical-harmonic (c2s) transforms for shells. !> !> @details Provides the per-shell transform matrices B(l) that map the !> unit-normalized Cartesian components of a shell (OpenQP's !> canonical bf_names order; see constants::CART_X/Y/Z) onto the !> 2l+1 real solid harmonics in CCA/libint order (m = -l..+l). !> !> Convention is identical to the validated Python reference in !> pyoqp/oqp/library/symmetry.py (_solid_harmonic_coefficients): !> each column of C2S_x holds the Cartesian coefficients of one !> spherical component, orthonormal against the intra-shell metric !> S of unit-normalized Cartesian Gaussians (B S B^T = I). The !> matrices below were generated from that reference and are !> re-verified at runtime by c2s_selftest(). !> !> The transform is applied to integrals that are already in the !> unit-normalized Cartesian basis (e.g. 2e blocks AFTER the !> rotation/Rys/libint normalization in int2::shellquartet, where !> all backends agree). s and p shells are passed through unchanged !> (Cartesian == spherical up to the trivial 1:1 / 3:3 mapping). module cart2sph use precision, only: dp use constants, only: NUM_CART_BF, NUM_SPH_BF, BAS_MXANG implicit none private public :: c2s_ncomp public :: cart2sph_eri public :: cart2sph_mat public :: cart2sph_mat_unit public :: cart2sph_vec public :: c2s_expand_block public :: c2s_expansion_matrix public :: c2s_selftest ! l=2 (D): Cart(6) -> Sph(5); column i = Cartesian coeffs of spherical i (m=-l..+l) real(dp), parameter :: C2S_D(6,5) = reshape([ & 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 1.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, & 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 1.000000000000000e+00_dp, & -4.999999999999999e-01_dp, -4.999999999999999e-01_dp, 9.999999999999999e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, & 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 1.000000000000000e+00_dp, 0.000000000000000e+00_dp, & 8.660254037844386e-01_dp, -8.660254037844386e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp & ], shape=[6,5]) ! l=3 (F): Cart(10) -> Sph(7); column i = Cartesian coeffs of spherical i (m=-l..+l) real(dp), parameter :: C2S_F(10,7) = reshape([ & 0.000000000000000e+00_dp, -7.905694150420950e-01_dp, 0.000000000000000e+00_dp, 1.060660171779821e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, & 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 1.000000000000000e+00_dp, & 0.000000000000000e+00_dp, -6.123724356957946e-01_dp, 0.000000000000000e+00_dp, -2.738612787525830e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 1.095445115010332e+00_dp, 0.000000000000000e+00_dp, & 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 1.000000000000000e+00_dp, 0.000000000000000e+00_dp, -6.708203932499369e-01_dp, 0.000000000000000e+00_dp, -6.708203932499369e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, & -6.123724356957946e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, -2.738612787525830e-01_dp, 0.000000000000000e+00_dp, 1.095445115010332e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, & 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 8.660254037844385e-01_dp, 0.000000000000000e+00_dp, -8.660254037844385e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, & 7.905694150420950e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, -1.060660171779821e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp & ], shape=[10,7]) ! l=4 (G): Cart(15) -> Sph(9); column i = Cartesian coeffs of spherical i (m=-l..+l) real(dp), parameter :: C2S_G(15,9) = reshape([ & 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 1.118033988749895e+00_dp, 0.000000000000000e+00_dp, -1.118033988749895e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, & 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, -7.905694150420950e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 1.060660171779821e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, & 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, -4.225771273642583e-01_dp, 0.000000000000000e+00_dp, -4.225771273642583e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 1.133893419027681e+00_dp, & 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, -8.964214570007953e-01_dp, 0.000000000000000e+00_dp, 1.195228609334394e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, -4.008918628686365e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, & 3.749999999999999e-01_dp, 3.749999999999999e-01_dp, 9.999999999999998e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 2.195775164134199e-01_dp, -8.783100656536798e-01_dp, -8.783100656536798e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, & 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, -8.964214570007953e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 1.195228609334394e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, -4.008918628686365e-01_dp, 0.000000000000000e+00_dp, & -5.590169943749475e-01_dp, 5.590169943749475e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 9.819805060619657e-01_dp, -9.819805060619657e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, & 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 7.905694150420950e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, -1.060660171779821e+00_dp, 0.000000000000000e+00_dp, & 7.395099728874520e-01_dp, 7.395099728874520e-01_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, -1.299038105676658e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp, 0.000000000000000e+00_dp & ], shape=[15,9]) contains !> @brief AO component count for a shell, honoring its harmonic flag. elemental integer function c2s_ncomp(l, pure) result(n) integer, intent(in) :: l, pure if (pure == 1 .and. l >= 2) then n = NUM_SPH_BF(l) else n = NUM_CART_BF(l) end if end function c2s_ncomp !> @brief Return the c2s matrix B(l) (ncart x nsph) for a pure shell. !> @details Only l >= 2 carries a non-trivial transform. Caller guarantees !> l >= 2 (s/p never reach here because c2s_ncomp keeps them !> Cartesian). The result is a copy sized (NUM_CART_BF(l), 2l+1). subroutine c2s_get(l, b) integer, intent(in) :: l real(dp), allocatable, intent(out) :: b(:,:) real(dp), allocatable, save :: c2s_h(:,:), c2s_i(:,:) select case (l) case (2) b = C2S_D case (3) b = C2S_F case (4) b = C2S_G case (5) if (.not. allocated(c2s_h)) c2s_h = c2s_build(l) b = c2s_h case (6) if (.not. allocated(c2s_i)) c2s_i = c2s_build(l) b = c2s_i case default error stop 'cart2sph: requested angular momentum exceeds BAS_MXANG' end select end subroutine c2s_get !> @brief Build B(l) for higher shells from the closed-form real solid-harmonic !> expansion used to generate the checked d/f/g tables above. function c2s_build(l) result(b) use constants, only: CART_X, CART_Y, CART_Z integer, intent(in) :: l real(dp), allocatable :: b(:,:) integer :: nc, ns, col, m, am, t, u, k, k_start, ax, ay, az integer :: sign_pow, cidx, ic, jc real(dp), allocatable :: metric(:,:), row(:) real(dp) :: coeff, norm2 nc = NUM_CART_BF(l) ns = NUM_SPH_BF(l) allocate(b(nc, ns), metric(nc, nc), row(nc)) do ic = 1, nc do jc = 1, nc metric(ic,jc) = cart_overlap(CART_X(ic,l), CART_Y(ic,l), CART_Z(ic,l), & CART_X(jc,l), CART_Y(jc,l), CART_Z(jc,l)) end do end do col = 0 do m = -l, l col = col + 1 am = abs(m) row = 0.0_dp do t = 0, (l - am) / 2 do u = 0, t k_start = merge(0, 1, m >= 0) do k = k_start, am, 2 sign_pow = t + (k - k_start) / 2 coeff = merge(1.0_dp, -1.0_dp, mod(sign_pow, 2) == 0) & * (0.25_dp ** t) & * real(ibinom(l, t), dp) & * real(ibinom(l - t, am + t), dp) & * real(ibinom(t, u), dp) & * real(ibinom(am, k), dp) if (coeff == 0.0_dp) cycle ax = 2*t + am - 2*u - k ay = 2*u + k az = l - 2*t - am if (ax < 0 .or. ay < 0 .or. az < 0) cycle cidx = cart_index(l, ax, ay, az) row(cidx) = row(cidx) + coeff / cart_component_norm(ax, ay, az) end do end do end do norm2 = dot_product(row, matmul(metric, row)) b(:, col) = row / sqrt(norm2) end do end function c2s_build integer function cart_index(l, ax, ay, az) result(idx) use constants, only: CART_X, CART_Y, CART_Z integer, intent(in) :: l, ax, ay, az integer :: i do i = 1, NUM_CART_BF(l) if (CART_X(i,l) == ax .and. CART_Y(i,l) == ay .and. CART_Z(i,l) == az) then idx = i return end if end do error stop 'cart2sph: generated monomial is absent from CART_X/Y/Z' end function cart_index pure real(dp) function cart_component_norm(ax, ay, az) result(n) integer, intent(in) :: ax, ay, az n = 1.0_dp / sqrt(real(idfact(2*ax - 1) * idfact(2*ay - 1) * idfact(2*az - 1), dp)) end function cart_component_norm pure integer function ibinom(n, k) result(c) integer, intent(in) :: n, k integer :: i if (k < 0 .or. k > n) then c = 0 return end if c = 1 do i = 1, k c = c * (n - i + 1) / i end do end function ibinom !> @brief Contract one index of a 3-way-folded block: out(il,is,ir) = !> sum_ic B(ic,is) * a(il,ic,ir), with the index laid out as !> (left, n_cart, right) in column-major order. subroutine contract_index(a, left, ncart, right, b, nsph, out) integer, intent(in) :: left, ncart, right, nsph real(dp), intent(in) :: a(left, ncart, right) real(dp), intent(in) :: b(ncart, nsph) real(dp), intent(out) :: out(left, nsph, right) integer :: il, is, ic, ir real(dp) :: bval out = 0.0_dp do ir = 1, right do is = 1, nsph do ic = 1, ncart bval = b(ic, is) if (bval == 0.0_dp) cycle do il = 1, left out(il, is, ir) = out(il, is, ir) + bval * a(il, ic, ir) end do end do end do end do end subroutine contract_index !> @brief Transform a 2e shell-quartet block from unit-normalized Cartesian !> to pure spherical for any index whose shell is flagged harmonic. !> !> @param[inout] ints flat ERI buffer; on entry holds the Cartesian block !> with storage dims [nbf(1)..nbf(4)] (column-major, !> index 1 fastest); on exit the spherical block with !> dims [nbf_out(1)..nbf_out(4)]. !> @param[in] am angular momentum of the four shells, in storage order !> @param[in] pure per-shell harmonic flag (1=spherical), storage order !> @param[in] nbf Cartesian component counts, storage order !> @param[out] nbf_out spherical component counts, storage order !> !> Storage order means dimension k of `ints` corresponds to am(k)/pure(k). !> Callers in int2 pass these already in the flipped (stored) order. subroutine cart2sph_eri(ints, am, pure, nbf, nbf_out) real(dp), intent(inout) :: ints(:) integer, intent(in) :: am(4), pure(4), nbf(4) integer, intent(out) :: nbf_out(4) integer :: dims(4) ! running per-index sizes (Cartesian -> spherical) integer :: p, left, right, k, nc, ns real(dp), allocatable :: b(:,:), src(:), dst(:) nbf_out = nbf do p = 1, 4 if (pure(p) /= 1 .or. am(p) < 2) cycle nbf_out(p) = NUM_SPH_BF(am(p)) end do if (all(nbf_out == nbf)) return ! nothing pure -> leave Cartesian block dims = nbf allocate(src(product(nbf))) src(1:product(nbf)) = ints(1:product(nbf)) do p = 1, 4 if (pure(p) /= 1 .or. am(p) < 2) cycle nc = dims(p) ns = NUM_SPH_BF(am(p)) left = product(dims(1:p-1)) right = product(dims(p+1:4)) call c2s_get(am(p), b) allocate(dst(left * ns * right)) call contract_index(src, left, nc, right, b, ns, dst) call move_alloc(dst, src) dims(p) = ns deallocate(b) end do k = product(dims) ints(1:k) = src(1:k) deallocate(src) end subroutine cart2sph_eri !> @brief Transform a 1e shell-pair block (unit-normalized Cartesian) to !> pure spherical for any harmonic-flagged shell. !> !> @details The block is laid out with the "fast" shell varying quickest, !> i.e. blk(nn), nn over (slow outer, fast inner) -- the order !> consumed by update_triang_matrix/update_rectangular_matrix !> (fast = shj, slow = shi). On exit blk(1:n_out) holds the !> spherical block in the same fast/slow layout and n_fast_out/ !> n_slow_out give its extents. subroutine cart2sph_mat(blk, l_fast, pure_fast, l_slow, pure_slow, n_fast_out, n_slow_out, iandj, antisym) use constants, only: shells_pnrm2 real(dp), intent(inout) :: blk(:) integer, intent(in) :: l_fast, pure_fast, l_slow, pure_slow integer, intent(out), optional :: n_fast_out, n_slow_out logical, intent(in), optional :: iandj !< .true. for a same-shell block, !< stored as a lower triangle (fast<=slow) logical, intent(in), optional :: antisym !< .true. if the operator is !< antisymmetric (L, GIAO H10, SOC); !< affects the iandj unpacking only integer :: ncf, ncs, nsf, nss, k, ic, ir logical :: tri, anti real(dp), allocatable :: b(:,:), src(:), dst(:) ncf = NUM_CART_BF(l_fast) ncs = NUM_CART_BF(l_slow) nsf = c2s_ncomp(l_fast, pure_fast) nss = c2s_ncomp(l_slow, pure_slow) if (present(n_fast_out)) n_fast_out = nsf if (present(n_slow_out)) n_slow_out = nss if (nsf == ncf .and. nss == ncs) return ! nothing pure -> Cartesian block tri = .false. if (present(iandj)) tri = iandj anti = .false. if (present(antisym)) anti = antisym ! Same-shell blocks (update_triang_matrix iandj path) are stored as a lower ! triangle blk((i-1)i/2 + j), j<=i. Unpack to a full Cartesian block using ! the operator's parity, transform as a rectangle, then repack to the ! spherical triangle. if (tri) then call cart2sph_tri(blk, l_fast, pure_fast, l_slow, pure_slow, ncf, ncs, nsf, nss, anti) return end if allocate(src(ncf*ncs)) src(1:ncf*ncs) = blk(1:ncf*ncs) ! 1e blocks arrive in the pure-power Cartesian normalization (the ! shells_pnrm2 per-component factors are applied later by bas_norm_matrix ! for Cartesian shells). B is defined for unit-normalized Cartesians, so ! for each index we transform, first fold in shells_pnrm2 along that index; ! the resulting spherical components are unit-normalized (set_bfnorms then ! uses bfnrm = 1 for them). Non-pure indices are left pure-power untouched. ! Contract the fast index (storage layout (ncf, ncs), fast contiguous). if (pure_fast == 1 .and. l_fast >= 2) then do ir = 1, ncs do ic = 1, ncf src((ir-1)*ncf + ic) = src((ir-1)*ncf + ic) * shells_pnrm2(ic, l_fast) end do end do call c2s_get(l_fast, b) allocate(dst(nsf*ncs)) call contract_index(src, 1, ncf, ncs, b, nsf, dst) call move_alloc(dst, src) deallocate(b) end if ! Contract the slow index ((nsf, ncs, 1); the slow index is the outer one). if (pure_slow == 1 .and. l_slow >= 2) then do ir = 1, ncs do ic = 1, nsf src((ir-1)*nsf + ic) = src((ir-1)*nsf + ic) * shells_pnrm2(ir, l_slow) end do end do call c2s_get(l_slow, b) allocate(dst(nsf*nss)) call contract_index(src, nsf, ncs, 1, b, nss, dst) call move_alloc(dst, src) deallocate(b) end if k = nsf*nss blk(1:k) = src(1:k) deallocate(src) end subroutine cart2sph_mat !> @brief Transform a rectangular 1e shell-pair block that is already in the !> unit-normalized Cartesian convention. !> @details This is used by backends such as libecpint that return normalized !> Cartesian matrices directly. Unlike cart2sph_mat, this does not !> fold in shells_pnrm2 before applying the c2s coefficients. subroutine cart2sph_mat_unit(blk, l_fast, pure_fast, l_slow, pure_slow) real(dp), intent(inout) :: blk(:) integer, intent(in) :: l_fast, pure_fast, l_slow, pure_slow integer :: ncf, ncs, nsf, nss, k real(dp), allocatable :: b(:,:), src(:), dst(:) ncf = NUM_CART_BF(l_fast) ncs = NUM_CART_BF(l_slow) nsf = c2s_ncomp(l_fast, pure_fast) nss = c2s_ncomp(l_slow, pure_slow) if (nsf == ncf .and. nss == ncs) return allocate(src(ncf*ncs)) src(1:ncf*ncs) = blk(1:ncf*ncs) if (pure_fast == 1 .and. l_fast >= 2) then call c2s_get(l_fast, b) allocate(dst(nsf*ncs)) call contract_index(src, 1, ncf, ncs, b, nsf, dst) call move_alloc(dst, src) deallocate(b) end if if (pure_slow == 1 .and. l_slow >= 2) then call c2s_get(l_slow, b) allocate(dst(nsf*nss)) call contract_index(src, nsf, ncs, 1, b, nss, dst) call move_alloc(dst, src) deallocate(b) end if k = nsf*nss blk(1:k) = src(1:k) deallocate(src) end subroutine cart2sph_mat_unit !> @brief Per-shell density-expansion matrix B'(l) = B(l) * shells_pnrm2, !> shape (NUM_CART_BF(l), 2l+1). Maps a unit-spherical index back to !> the pure-power Cartesian index for contraction with derivative !> integrals: D_cart = B'_i D_sph B'_j^T (see c2s_expand_block). subroutine c2s_expansion_matrix(l, bp) use constants, only: shells_pnrm2 integer, intent(in) :: l real(dp), allocatable, intent(out) :: bp(:,:) integer :: nc, ns, c, s call c2s_get(l, bp) ! bp = B(l), shape (nc, ns) nc = NUM_CART_BF(l) ns = NUM_SPH_BF(l) do s = 1, ns do c = 1, nc bp(c, s) = bp(c, s) * shells_pnrm2(c, l) end do end do end subroutine c2s_expansion_matrix !> @brief Expand a spherical density block to the pure-power Cartesian !> ("effective") density used by the gradient/Hessian kernels: !> D_cart = B'_i D_sph B'_j^T. Pure shells (l>=2) use B'; otherwise !> the index passes through unchanged (Cartesian == spherical). !> @param[in] dsph (nsph_i, nsph_j) spherical density block (bfnrm-folded) !> @param[out] dcart (ncart_i, ncart_j) Cartesian-effective density block subroutine c2s_expand_block(dsph, dcart, l_i, pure_i, l_j, pure_j) real(dp), intent(in) :: dsph(:,:) real(dp), intent(out) :: dcart(:,:) integer, intent(in) :: l_i, pure_i, l_j, pure_j real(dp), allocatable :: bi(:,:), bj(:,:), tmp(:,:) integer :: nci, ncj, nsi, nsj logical :: pi, pj nci = NUM_CART_BF(l_i); nsi = c2s_ncomp(l_i, pure_i) ncj = NUM_CART_BF(l_j); nsj = c2s_ncomp(l_j, pure_j) pi = (pure_i == 1 .and. l_i >= 2) pj = (pure_j == 1 .and. l_j >= 2) if (.not. pi .and. .not. pj) then dcart(1:nci, 1:ncj) = dsph(1:nsi, 1:nsj) return end if ! Expand the i (row) index: tmp(nci, nsj) = B'_i (nci,nsi) . dsph (nsi,nsj) allocate(tmp(nci, nsj)) if (pi) then call c2s_expansion_matrix(l_i, bi) tmp = matmul(bi, dsph(1:nsi, 1:nsj)) else tmp = dsph(1:nci, 1:nsj) end if ! Expand the j (col) index: dcart(nci, ncj) = tmp (nci,nsj) . B'_j^T (nsj,ncj) if (pj) then call c2s_expansion_matrix(l_j, bj) dcart(1:nci, 1:ncj) = matmul(tmp, transpose(bj)) else dcart(1:nci, 1:ncj) = tmp(1:nci, 1:ncj) end if deallocate(tmp) end subroutine c2s_expand_block !> @brief Transform a 1-index AO vector (e.g. grid AO values or one !> derivative component) from pure-power Cartesian to pure spherical. !> @details sph(s) = sum_c B(c,s) * shells_pnrm2(c,l) * cart(c). The pnrm !> fold makes the spherical components unit-normalized (downstream !> bfnrm = 1 for them, matching set_bfnorms). For l < 2 this is a !> straight copy. cart has NUM_CART_BF(l) entries, sph has 2l+1. subroutine cart2sph_vec(cart, sph, l) use constants, only: shells_pnrm2 real(dp), intent(in) :: cart(:) real(dp), intent(out) :: sph(:) integer, intent(in) :: l real(dp), allocatable :: b(:,:) integer :: nc, ns, c, s real(dp) :: acc nc = NUM_CART_BF(l) if (l < 2) then sph(1:nc) = cart(1:nc) return end if ns = NUM_SPH_BF(l) call c2s_get(l, b) do s = 1, ns acc = 0.0_dp do c = 1, nc acc = acc + b(c, s) * shells_pnrm2(c, l) * cart(c) end do sph(s) = acc end do deallocate(b) end subroutine cart2sph_vec !> @brief Same-shell (iandj) variant: the block is a packed lower triangle !> blk((i-1)i/2 + j), j<=i. Unpack to a full Cartesian block with the !> operator's parity (full(j>i) = +/- blk, zero diagonal when !> antisymmetric), transform as a rectangle, repack to the spherical !> triangle. subroutine cart2sph_tri(blk, l_fast, pure_fast, l_slow, pure_slow, ncf, ncs, nsf, nss, antisym) real(dp), intent(inout) :: blk(:) integer, intent(in) :: l_fast, pure_fast, l_slow, pure_slow, ncf, ncs, nsf, nss logical, intent(in) :: antisym real(dp), allocatable :: full(:) integer :: i, j, nn real(dp) :: mirror mirror = 1.0_dp if (antisym) mirror = -1.0_dp allocate(full(ncf*ncs)) do i = 1, ncs ! slow index (shi) do j = 1, ncf ! fast index (shj) if (j < i) then full((i-1)*ncf + j) = blk(i*(i-1)/2 + j) else if (j > i) then full((i-1)*ncf + j) = mirror*blk(j*(j-1)/2 + i) ! mirrored counterpart else if (antisym) then full((i-1)*ncf + j) = 0.0_dp ! antisymmetric diagonal is exact zero else full((i-1)*ncf + j) = blk(i*(i-1)/2 + j) end if end do end do call cart2sph_mat(full, l_fast, pure_fast, l_slow, pure_slow) do i = 1, nss do j = 1, i blk(i*(i-1)/2 + j) = full((i-1)*nsf + j) end do end do deallocate(full) end subroutine cart2sph_tri !> @brief Self-test: rebuild the intra-shell metric S of unit-normalized !> Cartesian Gaussians from the canonical exponents and verify !> B(l) S B(l)^T = I for every supported pure shell. Returns the !> worst deviation. subroutine c2s_selftest(max_err) use constants, only: CART_X, CART_Y, CART_Z real(dp), intent(out) :: max_err integer :: l, nc, ns, i, j, is, js real(dp), allocatable :: b(:,:), s(:,:), g(:,:) real(dp) :: err max_err = 0.0_dp do l = 2, BAS_MXANG nc = NUM_CART_BF(l) ns = NUM_SPH_BF(l) call c2s_get(l, b) allocate(s(nc, nc)) do i = 1, nc do j = 1, nc s(i,j) = cart_overlap(CART_X(i,l), CART_Y(i,l), CART_Z(i,l), & CART_X(j,l), CART_Y(j,l), CART_Z(j,l)) end do end do ! g = B^T S B (ns x ns), should be identity allocate(g(ns, ns)) g = matmul(matmul(transpose(b), s), b) do is = 1, ns do js = 1, ns err = abs(g(is,js) - merge(1.0_dp, 0.0_dp, is == js)) if (err > max_err) max_err = err end do end do deallocate(b, s, g) end do end subroutine c2s_selftest !> @brief Overlap of two unit-normalized Cartesian Gaussians of the same !> shell (same exponent), i.e. the intra-shell metric element. !> For a normalized x^a y^b z^c: <i|j> = prod_k (a_k+b_k-1)!! / !> sqrt((2a_k-1)!! (2b_k-1)!!) over k in {x,y,z}; zero if any sum odd. pure real(dp) function cart_overlap(ax, ay, az, bx, by, bz) result(s) integer, intent(in) :: ax, ay, az, bx, by, bz if (mod(ax+bx,2) /= 0 .or. mod(ay+by,2) /= 0 .or. mod(az+bz,2) /= 0) then s = 0.0_dp return end if s = ratio(ax, bx) * ratio(ay, by) * ratio(az, bz) end function cart_overlap !> @brief (a+b-1)!! / sqrt((2a-1)!! (2b-1)!!) for one Cartesian axis. pure real(dp) function ratio(a, b) result(r) integer, intent(in) :: a, b r = real(idfact(a+b-1), dp) / sqrt(real(idfact(2*a-1), dp) * real(idfact(2*b-1), dp)) end function ratio !> @brief Integer double factorial n!! with (-1)!! = 0!! = 1. pure integer function idfact(n) result(r) integer, intent(in) :: n integer :: k r = 1 k = n do while (k > 1) r = r * k k = k - 2 end do end function idfact end module cart2sph