Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Test of Clebsch-Gordon Coefficients
10 : !> \par History
11 : !> none
12 : !> \author JGH (28.02.2002)
13 : ! **************************************************************************************************
14 : MODULE cg_test
15 :
16 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
17 : USE kinds, ONLY: dp
18 : USE lebedev, ONLY: deallocate_lebedev_grids,&
19 : get_number_of_lebedev_grid,&
20 : init_lebedev_grids,&
21 : lebedev_grid
22 : USE machine, ONLY: m_walltime
23 : USE mathconstants, ONLY: pi
24 : USE spherical_harmonics, ONLY: clebsch_gordon,&
25 : clebsch_gordon_deallocate,&
26 : clebsch_gordon_init,&
27 : y_lm
28 : #include "../base/base_uses.f90"
29 :
30 : IMPLICIT NONE
31 :
32 : PRIVATE
33 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cg_test'
34 : PUBLIC :: clebsch_gordon_test
35 :
36 : CONTAINS
37 :
38 : ! **************************************************************************************************
39 : !> \brief ...
40 : ! **************************************************************************************************
41 2 : SUBROUTINE clebsch_gordon_test()
42 :
43 : INTEGER, PARAMETER :: l = 7
44 :
45 2 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: a1, a2, a3
46 : INTEGER :: il, iw, l1, l2, ll, lp, m1, m2, mm, mp, &
47 : na
48 : REAL(KIND=dp) :: ca, cga(10), cn, rga(10, 21), tend, &
49 : tstart
50 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: b1, b2, b3, wa
51 :
52 2 : iw = cp_logger_get_default_io_unit()
53 :
54 2 : IF (iw > 0) THEN
55 :
56 1 : WRITE (iw, '(/,A,/)') " Test of Clebsch-Gordon Coefficients"
57 1 : WRITE (iw, '(T40,A,T77,I4)') " Maximum l value tested:", l
58 :
59 1 : na = 500
60 1 : CALL init_lebedev_grids
61 1 : ll = get_number_of_lebedev_grid(n=na)
62 1 : na = lebedev_grid(ll)%n
63 3 : ALLOCATE (wa(na))
64 5 : ALLOCATE (a1(na), a2(na), a3(na))
65 4 : ALLOCATE (b1(na), b2(na), b3(na))
66 :
67 591 : wa(1:na) = 4.0_dp*pi*lebedev_grid(ll)%w(1:na)
68 :
69 1 : tstart = m_walltime()
70 1 : CALL clebsch_gordon_init(l)
71 1 : tend = m_walltime()
72 1 : tend = tend - tstart
73 1 : WRITE (iw, '(T30,A,T71,F10.3)') " Time for Clebsch-Gordon Table [s] ", tend
74 : lp = (l**4 + 6*l**3 + 15*l**2 + 18*l + 8)/8
75 1 : lp = 2*lp*(l + 1)
76 1 : WRITE (iw, '(T30,A,T71,I10)') " Size of Clebsch-Gordon Table ", lp
77 1 : WRITE (iw, '(/,A)') " Start Test for Complex Spherical Harmonics "
78 :
79 9 : DO l1 = 0, l
80 72 : DO m1 = -l1, l1
81 64 : CALL y_lm(lebedev_grid(ll)%r, a1, l1, m1)
82 584 : DO l2 = 0, l
83 4672 : DO m2 = -l2, l2
84 4096 : CALL y_lm(lebedev_grid(ll)%r, a2, l2, m2)
85 4096 : CALL clebsch_gordon(l1, m1, l2, m2, cga)
86 27408 : DO lp = MOD(l1 + l2, 2), l1 + l2, 2
87 22800 : mp = m1 + m2
88 22800 : IF (lp < ABS(mp)) CYCLE
89 15240 : CALL y_lm(lebedev_grid(ll)%r, a3, lp, mp)
90 9006840 : cn = REAL(SUM(a1*a2*CONJG(a3)*wa), KIND=dp)
91 15240 : il = lp/2 + 1
92 15240 : ca = cga(il)
93 19336 : IF (ABS(ca - cn) > 1.e-10_dp) THEN
94 0 : WRITE (*, '(A,3I5,A,F20.12)') " l ", l1, l2, lp, " A ", ca
95 0 : WRITE (*, '(A,3I5,A,F20.12)') " m ", m1, m2, mp, " N ", cn
96 0 : WRITE (*, *)
97 : END IF
98 : END DO
99 : END DO
100 : END DO
101 : END DO
102 9 : WRITE (iw, '(A,i2,A)') " Test for l = ", l1, " done"
103 : END DO
104 :
105 1 : WRITE (iw, '(/,A)') " Start Test for Real Spherical Harmonics "
106 9 : DO l1 = 0, l
107 72 : DO m1 = -l1, l1
108 64 : CALL y_lm(lebedev_grid(ll)%r, b1, l1, m1)
109 584 : DO l2 = 0, l
110 4672 : DO m2 = -l2, l2
111 4096 : CALL y_lm(lebedev_grid(ll)%r, b2, l2, m2)
112 4096 : CALL clebsch_gordon(l1, m1, l2, m2, rga)
113 4096 : mp = m1 + m2
114 4096 : mm = m1 - m2
115 4096 : IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN
116 2016 : mp = -ABS(mp)
117 2016 : mm = -ABS(mm)
118 : ELSE
119 2080 : mp = ABS(mp)
120 2080 : mm = ABS(mm)
121 : END IF
122 27408 : DO lp = MOD(l1 + l2, 2), l1 + l2, 2
123 22800 : IF (ABS(mp) <= lp) THEN
124 15240 : CALL y_lm(lebedev_grid(ll)%r, b3, lp, mp)
125 9006840 : cn = SUM(b1*b2*b3*wa)
126 15240 : il = lp/2 + 1
127 15240 : ca = rga(il, 1)
128 15240 : IF (ABS(ca - cn) > 1.e-10_dp) THEN
129 0 : WRITE (*, '(A,3I5,A,F20.12)') " l ", l1, l2, lp, " A ", ca
130 0 : WRITE (*, '(A,3I5,A,F20.12)') " m ", m1, m2, mp, " N ", cn
131 0 : WRITE (*, *)
132 : END IF
133 : END IF
134 26896 : IF (mp /= mm .AND. ABS(mm) <= lp) THEN
135 11832 : CALL y_lm(lebedev_grid(ll)%r, b3, lp, mm)
136 6992712 : cn = SUM(b1*b2*b3*wa)
137 11832 : il = lp/2 + 1
138 11832 : ca = rga(il, 2)
139 11832 : IF (ABS(ca - cn) > 1.e-10_dp) THEN
140 0 : WRITE (*, '(A,3I5,A,F20.12)') " l ", l1, l2, lp, " A ", ca
141 0 : WRITE (*, '(A,3I5,A,F20.12)') " m ", m1, m2, mm, " N ", cn
142 0 : WRITE (*, *)
143 : END IF
144 : END IF
145 : END DO
146 : END DO
147 : END DO
148 : END DO
149 9 : WRITE (iw, '(A,i2,A)') " Test for l = ", l1, " done"
150 : END DO
151 :
152 1 : DEALLOCATE (wa)
153 1 : DEALLOCATE (a1, a2, a3)
154 1 : DEALLOCATE (b1, b2, b3)
155 :
156 1 : CALL deallocate_lebedev_grids()
157 1 : CALL clebsch_gordon_deallocate()
158 :
159 : END IF
160 :
161 2 : END SUBROUTINE clebsch_gordon_test
162 :
163 : END MODULE cg_test
164 :
|
