Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Expansion of the product of two real spherical harmonics in real
10 : !> spherical harmonics
11 : !> \note
12 : !> The expansion coefficients are the Gaunt coefficients of the real
13 : !> spherical harmonics; they are obtained from libwignernj, which evaluates
14 : !> them in exact integer arithmetic. Historically CP2K called these
15 : !> "Clebsch-Gordon coefficients", which they are not, and built them from a
16 : !> table of factorial ratios that lost about one digit per two units of l.
17 : !> \par History
18 : !> 09.2026 created out of the five copies of this loop that lived in
19 : !> qs_rho_atom_methods, qs_tddfpt2_soc, qs_cneo_utils,
20 : !> xas_tdp_atom and generic_shg_integrals_init [Susi Lehtola]
21 : !> \author Susi Lehtola
22 : ! **************************************************************************************************
23 : MODULE gaunt_coefficients
24 :
25 : USE kinds, ONLY: dp
26 : USE memory_utilities, ONLY: reallocate
27 : USE orbital_pointers, ONLY: indso,&
28 : nsoset
29 : USE wignernj_interface, ONLY: real_gaunt
30 : #include "../base/base_uses.f90"
31 :
32 : IMPLICIT NONE
33 :
34 : PRIVATE
35 :
36 : PUBLIC :: real_gaunt_coefficients
37 :
38 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gaunt_coefficients'
39 :
40 : ! Coefficients below this threshold are not entered in the non-zero list
41 : REAL(KIND=dp), PARAMETER, PRIVATE :: cg_eps = 1.0E-8_dp
42 :
43 : CONTAINS
44 :
45 : ! **************************************************************************************************
46 : !> \brief Calculate the expansion of the product of two real spherical harmonics
47 : !> in real spherical harmonics,
48 : !> S(l1,m1) S(l2,m2) = sum_lm my_cg(iso1,iso2,iso) S(l,m)
49 : !> \param my_cg the expansion coefficients, allocated here; the pointer is
50 : !> INTENT(OUT), so the caller must not rely on its previous contents
51 : !> \param maxl1 maximal angular momentum of the first harmonic
52 : !> \param maxl2 maximal angular momentum of the second harmonic
53 : !> \param llmax maximal angular momentum retained in the expansion
54 : !> \param cg_none0_list list of the harmonics with a non-zero coefficient,
55 : !> allocated here if present
56 : !> \param ncg_none0 number of harmonics with a non-zero coefficient,
57 : !> allocated here if present
58 : ! **************************************************************************************************
59 1852 : SUBROUTINE real_gaunt_coefficients(my_cg, maxl1, maxl2, llmax, cg_none0_list, ncg_none0)
60 :
61 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT), &
62 : POINTER :: my_cg
63 : INTEGER, INTENT(IN) :: maxl1, maxl2, llmax
64 : INTEGER, DIMENSION(:, :, :), INTENT(OUT), &
65 : OPTIONAL, POINTER :: cg_none0_list
66 : INTEGER, DIMENSION(:, :), INTENT(OUT), OPTIONAL, &
67 : POINTER :: ncg_none0
68 :
69 : CHARACTER(len=*), PARAMETER :: routineN = 'real_gaunt_coefficients'
70 :
71 : INTEGER :: handle, ilist, im, iso, iso1, iso2, l1, &
72 : l2, lp, m1, m2, mp(4), nlist, &
73 : nlist_max, nsfunc, nsfunc1, nsfunc2
74 : LOGICAL :: do_list
75 :
76 1852 : CALL timeset(routineN, handle)
77 :
78 1852 : CPASSERT(maxl1 >= 0)
79 1852 : CPASSERT(maxl2 >= 0)
80 1852 : CPASSERT(llmax >= 0)
81 :
82 1852 : nsfunc1 = nsoset(maxl1)
83 1852 : nsfunc2 = nsoset(maxl2)
84 1852 : nsfunc = nsoset(llmax)
85 :
86 9260 : ALLOCATE (my_cg(nsfunc1, nsfunc2, nsfunc))
87 6609510 : my_cg = 0.0_dp
88 :
89 15904 : DO iso1 = 1, nsfunc1
90 14052 : l1 = indso(1, iso1)
91 14052 : m1 = indso(2, iso1)
92 161980 : DO iso2 = 1, nsfunc2
93 146076 : l2 = indso(1, iso2)
94 146076 : m2 = indso(2, iso2)
95 : ! The product of two real harmonics only contains the azimuthal
96 : ! frequencies |m1+m2| and |m1-m2|; the cosine part goes to m >= 0
97 : ! and the sine part to m < 0, so of the four candidates below at
98 : ! most two have a non-vanishing coefficient. Rather than working
99 : ! out which one it is, we evaluate all four and let the selection
100 : ! rules zero out the rest.
101 146076 : mp(1) = ABS(m1 + m2)
102 146076 : mp(2) = -mp(1)
103 146076 : mp(3) = ABS(m1 - m2)
104 146076 : mp(4) = -mp(3)
105 538290 : DO lp = MOD(l1 + l2, 2), MIN(l1 + l2, llmax), 2
106 2036886 : DO im = 1, 4
107 1512648 : IF (ABS(mp(im)) > lp) CYCLE
108 1056312 : iso = nsoset(lp - 1) + lp + 1 + mp(im)
109 1890810 : my_cg(iso1, iso2, iso) = real_gaunt(l1, m1, l2, m2, lp, mp(im))
110 : END DO
111 : END DO
112 : END DO
113 : END DO
114 :
115 1852 : do_list = PRESENT(cg_none0_list) .OR. PRESENT(ncg_none0)
116 1852 : IF (do_list) THEN
117 72 : CPASSERT(PRESENT(cg_none0_list))
118 72 : CPASSERT(PRESENT(ncg_none0))
119 72 : nlist_max = 6
120 288 : ALLOCATE (ncg_none0(nsfunc1, nsfunc2))
121 360 : ALLOCATE (cg_none0_list(nsfunc1, nsfunc2, nlist_max))
122 15110 : ncg_none0 = 0
123 90732 : cg_none0_list = 0
124 682 : DO iso1 = 1, nsfunc1
125 14252 : DO iso2 = 1, nsfunc2
126 13570 : nlist = 0
127 790740 : DO ilist = 1, nsfunc
128 790740 : IF (ABS(my_cg(iso1, iso2, ilist)) > cg_eps) THEN
129 37734 : nlist = nlist + 1
130 37734 : IF (nlist > nlist_max) THEN
131 8 : CALL reallocate(cg_none0_list, 1, nsfunc1, 1, nsfunc2, 1, nlist)
132 8 : nlist_max = nlist
133 : END IF
134 37734 : cg_none0_list(iso1, iso2, nlist) = ilist
135 : END IF
136 : END DO
137 14180 : ncg_none0(iso1, iso2) = nlist
138 : END DO
139 : END DO
140 : END IF
141 :
142 1852 : CALL timestop(handle)
143 :
144 1852 : END SUBROUTINE real_gaunt_coefficients
145 :
146 : END MODULE gaunt_coefficients
|
