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 Interface to libwignernj for the angular momentum algebra
10 : !> \note
11 : !> libwignernj carries out all intermediate arithmetic exactly, in a prime
12 : !> factorization representation, and rounds only once at the end. The
13 : !> returned coefficients are therefore correct to the last bit, whereas the
14 : !> factorial-ratio expressions that CP2K used before lost roughly one digit
15 : !> per two units of l.
16 : !>
17 : !> The library takes all angular momenta as twice their value so that
18 : !> half-integers stay exact; the wrappers here take the plain integer l and
19 : !> m of a spherical harmonic and do the doubling.
20 : !>
21 : !> The real spherical harmonics S(l,m) follow the Condon-Shortley
22 : !> convention
23 : !> S(l, 0) = Y(l,0)
24 : !> S(l, m>0) = ( Y(l,-m) + (-1)^m Y(l,m) )/sqrt(2)
25 : !> S(l, m<0) = ( Y(l,m) - (-1)^|m| Y(l,-m) )*i/sqrt(2)
26 : !> which is the same convention as spherical_harmonics::y_lm, so no
27 : !> additional phase is needed on either side.
28 : !>
29 : !> The interfaces below bind directly to the C API, so no compiler-specific
30 : !> Fortran module from libwignernj is required.
31 : !> \par History
32 : !> 09.2026 created [Susi Lehtola]
33 : !> \author Susi Lehtola
34 : ! **************************************************************************************************
35 : MODULE wignernj_interface
36 :
37 : USE ISO_C_BINDING, ONLY: C_DOUBLE,&
38 : C_INT
39 : USE kinds, ONLY: dp
40 : #include "../base/base_uses.f90"
41 :
42 : IMPLICIT NONE
43 :
44 : PRIVATE
45 :
46 : PUBLIC :: real_gaunt, complex_gaunt
47 :
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'wignernj_interface'
49 :
50 : INTERFACE
51 : ! **************************************************************************************************
52 : !> \brief Gaunt coefficient of three complex spherical harmonics
53 : !> \param tl1 twice the angular momentum of the first harmonic
54 : !> \param tm1 twice the magnetic quantum number of the first harmonic
55 : !> \param tl2 twice the angular momentum of the second harmonic
56 : !> \param tm2 twice the magnetic quantum number of the second harmonic
57 : !> \param tl3 twice the angular momentum of the third harmonic
58 : !> \param tm3 twice the magnetic quantum number of the third harmonic
59 : !> \return integral of Y(l1,m1)*Y(l2,m2)*Y(l3,m3) over the unit sphere
60 : ! **************************************************************************************************
61 : FUNCTION gaunt_c(tl1, tm1, tl2, tm2, tl3, tm3) BIND(C, name="gaunt")
62 : IMPORT :: C_DOUBLE, C_INT
63 : INTEGER(KIND=C_INT), VALUE :: tl1, tm1, tl2, tm2, tl3, tm3
64 : REAL(KIND=C_DOUBLE) :: gaunt_c
65 : END FUNCTION gaunt_c
66 :
67 : ! **************************************************************************************************
68 : !> \brief Gaunt coefficient of three real spherical harmonics
69 : !> \param tl1 twice the angular momentum of the first harmonic
70 : !> \param tm1 twice the magnetic quantum number of the first harmonic
71 : !> \param tl2 twice the angular momentum of the second harmonic
72 : !> \param tm2 twice the magnetic quantum number of the second harmonic
73 : !> \param tl3 twice the angular momentum of the third harmonic
74 : !> \param tm3 twice the magnetic quantum number of the third harmonic
75 : !> \return integral of S(l1,m1)*S(l2,m2)*S(l3,m3) over the unit sphere
76 : ! **************************************************************************************************
77 : FUNCTION gaunt_real_c(tl1, tm1, tl2, tm2, tl3, tm3) BIND(C, name="gaunt_real")
78 : IMPORT :: C_DOUBLE, C_INT
79 : INTEGER(KIND=C_INT), VALUE :: tl1, tm1, tl2, tm2, tl3, tm3
80 : REAL(KIND=C_DOUBLE) :: gaunt_real_c
81 : END FUNCTION gaunt_real_c
82 : END INTERFACE
83 :
84 : CONTAINS
85 :
86 : ! **************************************************************************************************
87 : !> \brief Gaunt coefficient of three real spherical harmonics,
88 : !> int S(l1,m1) S(l2,m2) S(l3,m3) dOmega
89 : !> \param l1 angular momentum of the first harmonic
90 : !> \param m1 magnetic quantum number of the first harmonic
91 : !> \param l2 angular momentum of the second harmonic
92 : !> \param m2 magnetic quantum number of the second harmonic
93 : !> \param l3 angular momentum of the third harmonic
94 : !> \param m3 magnetic quantum number of the third harmonic
95 : !> \return the Gaunt coefficient, zero if the selection rules are violated
96 : ! **************************************************************************************************
97 1977912 : FUNCTION real_gaunt(l1, m1, l2, m2, l3, m3) RESULT(g)
98 : INTEGER, INTENT(IN) :: l1, m1, l2, m2, l3, m3
99 : REAL(KIND=dp) :: g
100 :
101 1977912 : g = REAL(gaunt_real_c(2*l1, 2*m1, 2*l2, 2*m2, 2*l3, 2*m3), KIND=dp)
102 :
103 1977912 : END FUNCTION real_gaunt
104 :
105 : ! **************************************************************************************************
106 : !> \brief Gaunt coefficient of three complex spherical harmonics,
107 : !> int Y(l1,m1) Y(l2,m2) Y(l3,m3) dOmega
108 : !> \param l1 angular momentum of the first harmonic
109 : !> \param m1 magnetic quantum number of the first harmonic
110 : !> \param l2 angular momentum of the second harmonic
111 : !> \param m2 magnetic quantum number of the second harmonic
112 : !> \param l3 angular momentum of the third harmonic
113 : !> \param m3 magnetic quantum number of the third harmonic
114 : !> \return the Gaunt coefficient, zero if the selection rules are violated
115 : ! **************************************************************************************************
116 921600 : FUNCTION complex_gaunt(l1, m1, l2, m2, l3, m3) RESULT(g)
117 : INTEGER, INTENT(IN) :: l1, m1, l2, m2, l3, m3
118 : REAL(KIND=dp) :: g
119 :
120 921600 : g = REAL(gaunt_c(2*l1, 2*m1, 2*l2, 2*m2, 2*l3, 2*m3), KIND=dp)
121 :
122 921600 : END FUNCTION complex_gaunt
123 :
124 : END MODULE wignernj_interface
|
