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 Test of the Gaunt coefficients
10 : !> \note
11 : !> Checks that the Gaunt coefficients returned by libwignernj agree with a
12 : !> Lebedev quadrature of the corresponding triple product of the spherical
13 : !> harmonics of spherical_harmonics::y_lm. This guards the phase and
14 : !> normalisation conventions on both sides; the values themselves are
15 : !> exact by construction in libwignernj.
16 : !> \par History
17 : !> JGH 1-Mar-2002 : Clebsch-Gordon coefficients
18 : !> 09.2026 rewritten against libwignernj [Susi Lehtola]
19 : !> \author JGH 1-Mar-2002
20 : ! **************************************************************************************************
21 : MODULE cg_test
22 :
23 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
24 : USE kinds, ONLY: dp
25 : USE lebedev, ONLY: deallocate_lebedev_grids,&
26 : get_number_of_lebedev_grid,&
27 : init_lebedev_grids,&
28 : lebedev_grid
29 : USE mathconstants, ONLY: pi
30 : USE spherical_harmonics, ONLY: y_lm
31 : USE wignernj_interface, ONLY: complex_gaunt,&
32 : real_gaunt
33 : #include "../base/base_uses.f90"
34 :
35 : IMPLICIT NONE
36 : PRIVATE
37 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cg_test'
38 : PUBLIC :: gaunt_coefficient_test
39 :
40 : CONTAINS
41 :
42 : ! **************************************************************************************************
43 : !> \brief Compare the Gaunt coefficients against a Lebedev quadrature of the
44 : !> triple product of spherical harmonics
45 : ! **************************************************************************************************
46 2 : SUBROUTINE gaunt_coefficient_test()
47 :
48 : INTEGER, PARAMETER :: lmax = 7
49 : REAL(KIND=dp), PARAMETER :: eps = 1.0E-10_dp
50 :
51 2 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: yc
52 : INTEGER :: i1, i2, i3, iw, l, ll, m, n, na, nfun, &
53 : nfun1
54 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: lof, mof
55 : REAL(KIND=dp) :: dmax_complex, dmax_real, ref
56 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wa
57 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: yr
58 :
59 2 : iw = cp_logger_get_default_io_unit()
60 :
61 2 : IF (iw <= 0) RETURN
62 :
63 1 : WRITE (iw, '(/,A,/)') " Test of the Gaunt coefficients"
64 1 : WRITE (iw, '(T40,A,T77,I4)') " Maximum l value tested:", lmax
65 :
66 1 : na = 500
67 1 : CALL init_lebedev_grids()
68 1 : ll = get_number_of_lebedev_grid(n=na)
69 1 : na = lebedev_grid(ll)%n
70 3 : ALLOCATE (wa(na))
71 591 : wa(1:na) = 4.0_dp*pi*lebedev_grid(ll)%w(1:na)
72 :
73 : ! The two harmonics that are multiplied run up to lmax, their product
74 : ! is expanded in harmonics up to 2*lmax
75 1 : nfun1 = (lmax + 1)**2
76 1 : nfun = (2*lmax + 1)**2
77 1 : ALLOCATE (lof(nfun), mof(nfun))
78 5 : ALLOCATE (yr(na, nfun), yc(na, nfun))
79 1 : n = 0
80 16 : DO l = 0, 2*lmax
81 241 : DO m = -l, l
82 225 : n = n + 1
83 225 : lof(n) = l
84 225 : mof(n) = m
85 225 : CALL y_lm(lebedev_grid(ll)%r, yr(:, n), l, m)
86 240 : CALL y_lm(lebedev_grid(ll)%r, yc(:, n), l, m)
87 : END DO
88 : END DO
89 1 : CPASSERT(n == nfun)
90 :
91 1 : dmax_real = 0.0_dp
92 1 : dmax_complex = 0.0_dp
93 65 : DO i1 = 1, nfun1
94 4160 : DO i2 = 1, nfun1
95 925760 : DO i3 = 1, nfun
96 : ! Real spherical harmonics
97 921600 : ref = real_gaunt(lof(i1), mof(i1), lof(i2), mof(i2), lof(i3), mof(i3))
98 : dmax_real = MAX(dmax_real, &
99 544665600 : ABS(ref - SUM(wa*yr(:, i1)*yr(:, i2)*yr(:, i3))))
100 : ! Complex spherical harmonics. Unlike the real ones, the complex
101 : ! harmonics of y_lm are defined without the Condon-Shortley phase
102 : ! (see cvy_lm: P_l^|m| with no (-1)^m), so they are related to the
103 : ! Condon-Shortley harmonics that libwignernj uses by
104 : ! Y_cp2k(l,m) = (-1)^max(m,0) Y_cs(l,m).
105 : ! The Gaunt coefficient vanishes unless m1 + m2 + m3 = 0, in which
106 : ! case the integral is real.
107 921600 : ref = complex_gaunt(lof(i1), mof(i1), lof(i2), mof(i2), lof(i3), mof(i3))
108 921600 : ref = ref*(-1.0_dp)**(MAX(mof(i1), 0) + MAX(mof(i2), 0) + MAX(mof(i3), 0))
109 : dmax_complex = MAX(dmax_complex, &
110 : ABS(CMPLX(ref, 0.0_dp, KIND=dp) - &
111 544669696 : SUM(wa*yc(:, i1)*yc(:, i2)*yc(:, i3))))
112 : END DO
113 : END DO
114 65 : WRITE (iw, '(A,I2,A)') " Test for l = ", lof(i1), " done"
115 : END DO
116 :
117 1 : WRITE (iw, '(/,T30,A,T71,ES10.3)') " Real spherical harmonics, max deviation ", dmax_real
118 1 : WRITE (iw, '(T30,A,T71,ES10.3)') " Complex spherical harmonics, max deviation ", dmax_complex
119 :
120 1 : DEALLOCATE (wa, lof, mof, yr, yc)
121 1 : CALL deallocate_lebedev_grids()
122 :
123 1 : IF (MAX(dmax_real, dmax_complex) > eps) THEN
124 0 : CPABORT("The Gaunt coefficients disagree with the numerical integration")
125 : END IF
126 :
127 : END SUBROUTINE gaunt_coefficient_test
128 :
129 : END MODULE cg_test
|
