(git:374b731)
Loading...
Searching...
No Matches
spherical_harmonics.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculate spherical harmonics
10!> \note
11!> Spherical Harmonics
12!> Numerical Stability up to L=15
13!> Accuracy > 1.E-12 up to L=15 tested
14!> Definition is consistent with orbital_transformation_matrices
15!> Clebsch-Gordon Coefficients
16!> Tested up to l=7 (i.e. L=14)
17!> \todo
18!> 1) Check if this definition is consistent with the
19!> Slater-Koster module
20!> \par History
21!> JGH 28-Feb-2002 : Change of sign convention (-1^m)
22!> JGH 1-Mar-2002 : Clebsch-Gordon Coefficients
23!> - Clebsch-Gordon coefficients and Wigner 3-j symbols added as generic routines using the
24!> standard normalization (19.09.2022, MK)
25!> \author JGH 6-Oct-2000, MK
26! **************************************************************************************************
27MODULE spherical_harmonics
28
29 USE kinds, ONLY: dp
30 USE mathconstants, ONLY: fac,&
31 maxfac,&
32 pi
33#include "../base/base_uses.f90"
34
35 IMPLICIT NONE
36
37 PRIVATE
38
39 PUBLIC :: y_lm, dy_lm
40 PUBLIC :: clebsch_gordon_deallocate, clebsch_gordon_init, &
41 clebsch_gordon
42 PUBLIC :: legendre, dlegendre
43 PUBLIC :: clebsch_gordon_coefficient, wigner_3j
44
45 INTERFACE y_lm
46 MODULE PROCEDURE rvy_lm, rry_lm, cvy_lm, ccy_lm
47 END INTERFACE
48
49 INTERFACE dy_lm
50 MODULE PROCEDURE dry_lm, dcy_lm
51 END INTERFACE
52
53 INTERFACE clebsch_gordon
54 MODULE PROCEDURE clebsch_gordon_real, clebsch_gordon_complex
55 END INTERFACE
56
57 INTERFACE clebsch_gordon_getm
58 MODULE PROCEDURE getm
59 END INTERFACE
60
61 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'spherical_harmonics'
62
63 REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: cg_table
64 INTEGER :: lmax = -1
65 REAL(KIND=dp) :: osq2, sq2
66
67CONTAINS
68
69! Clebsch-Gordon Coefficients
70
71! **************************************************************************************************
72!> \brief ...
73!> \param l1 ...
74!> \param m1 ...
75!> \param l2 ...
76!> \param m2 ...
77!> \param clm ...
78! **************************************************************************************************
79 SUBROUTINE clebsch_gordon_complex(l1, m1, l2, m2, clm)
80 INTEGER, INTENT(IN) :: l1, m1, l2, m2
81 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: clm
82
83 INTEGER :: icase, ind, l, lm, lp, n
84
85 l = l1 + l2
86 IF (l > lmax) CALL clebsch_gordon_init(l)
87 n = l/2 + 1
88 IF (n > SIZE(clm)) cpabort("Array too small")
89
90 IF ((m1 >= 0 .AND. m2 >= 0) .OR. (m1 < 0 .AND. m2 < 0)) THEN
91 icase = 1
92 ELSE
93 icase = 2
94 END IF
95 ind = order(l1, m1, l2, m2)
96
97 DO lp = mod(l, 2), l, 2
98 lm = lp/2 + 1
99 clm(lm) = cg_table(ind, lm, icase)
100 END DO
101
102 END SUBROUTINE clebsch_gordon_complex
103
104! **************************************************************************************************
105!> \brief ...
106!> \param l1 ...
107!> \param m1 ...
108!> \param l2 ...
109!> \param m2 ...
110!> \param rlm ...
111! **************************************************************************************************
112 SUBROUTINE clebsch_gordon_real(l1, m1, l2, m2, rlm)
113 INTEGER, INTENT(IN) :: l1, m1, l2, m2
114 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: rlm
115
116 INTEGER :: icase1, icase2, ind, l, lm, lp, mm(2), n
117 REAL(KIND=dp) :: xsi
118
119 l = l1 + l2
120 IF (l > lmax) CALL clebsch_gordon_init(l)
121 n = l/2 + 1
122 IF (n > SIZE(rlm, 1)) cpabort("Array too small")
123
124 ind = order(l1, m1, l2, m2)
125 mm = getm(m1, m2)
126 IF ((m1 >= 0 .AND. m2 >= 0) .OR. (m1 < 0 .AND. m2 < 0)) THEN
127 icase1 = 1
128 icase2 = 2
129 ELSE
130 icase1 = 2
131 icase2 = 1
132 END IF
133
134 DO lp = mod(l, 2), l, 2
135 lm = lp/2 + 1
136 xsi = get_factor(m1, m2, mm(1))
137 rlm(lm, 1) = xsi*cg_table(ind, lm, icase1)
138 xsi = get_factor(m1, m2, mm(2))
139 rlm(lm, 2) = xsi*cg_table(ind, lm, icase2)
140 END DO
141
142 END SUBROUTINE clebsch_gordon_real
143
144! **************************************************************************************************
145!> \brief ...
146!> \param m1 ...
147!> \param m2 ...
148!> \return ...
149! **************************************************************************************************
150 FUNCTION getm(m1, m2) RESULT(m)
151 INTEGER, INTENT(IN) :: m1, m2
152 INTEGER, DIMENSION(2) :: m
153
154 INTEGER :: mm, mp
155
156 mp = m1 + m2
157 mm = m1 - m2
158 IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN
159 mp = -abs(mp)
160 mm = -abs(mm)
161 ELSE
162 mp = abs(mp)
163 mm = abs(mm)
164 END IF
165 m(1) = mp
166 m(2) = mm
167 END FUNCTION getm
168
169! **************************************************************************************************
170!> \brief ...
171!> \param m1 ...
172!> \param m2 ...
173!> \param m ...
174!> \return ...
175! **************************************************************************************************
176 FUNCTION get_factor(m1, m2, m) RESULT(f)
177 INTEGER, INTENT(IN) :: m1, m2, m
178 REAL(KIND=dp) :: f
179
180 INTEGER :: mx, my
181
182 f = 1.0_dp
183 IF (abs(m1) >= abs(m2)) THEN
184 mx = m1
185 my = m2
186 ELSE
187 mx = m2
188 my = m1
189 END IF
190 IF (mx*my == 0) THEN
191 f = 1.0_dp
192 ELSE IF (m == 0) THEN
193 IF (abs(mx) /= abs(my)) WRITE (*, '(A,3I6)') " 1) Illegal Case ", m1, m2, m
194 IF (mx*my > 0) THEN
195 f = 1.0_dp
196 ELSE
197 f = 0.0_dp
198 END IF
199 ELSE IF (abs(mx) + abs(my) == m) THEN
200 f = osq2
201 IF (mx < 0) f = -osq2
202 ELSE IF (abs(mx) + abs(my) == -m) THEN
203 f = osq2
204 ELSE IF (abs(mx) - abs(my) == -m) THEN
205 IF (mx*my > 0) WRITE (*, '(A,3I6)') " 2) Illegal Case ", m1, m2, m
206 IF (mx > 0) f = -osq2
207 IF (mx < 0) f = osq2
208 ELSE IF (abs(mx) - abs(my) == m) THEN
209 IF (mx*my < 0) WRITE (*, '(A,3I6)') " 3) Illegal Case ", m1, m2, m
210 f = osq2
211 ELSE
212 WRITE (*, '(A,3I6)') " 4) Illegal Case ", m1, m2, m
213 END IF
214 END FUNCTION get_factor
215
216! **************************************************************************************************
217!> \brief ...
218! **************************************************************************************************
219 SUBROUTINE clebsch_gordon_deallocate()
220 CHARACTER(len=*), PARAMETER :: routinen = 'clebsch_gordon_deallocate'
221
222 INTEGER :: handle
223
224 CALL timeset(routinen, handle)
225
226 IF (ALLOCATED(cg_table)) THEN
227 DEALLOCATE (cg_table)
228 END IF
229 CALL timestop(handle)
230
231 END SUBROUTINE clebsch_gordon_deallocate
232
233! **************************************************************************************************
234!> \brief ...
235!> \param l ...
236! **************************************************************************************************
237 SUBROUTINE clebsch_gordon_init(l)
238 INTEGER, INTENT(IN) :: l
239
240 CHARACTER(len=*), PARAMETER :: routinen = 'clebsch_gordon_init'
241
242 INTEGER :: handle, i1, i2, ix, iy, l1, l2, lp, m, &
243 m1, m2, ml, mp, n
244
245 CALL timeset(routinen, handle)
246
247 sq2 = sqrt(2.0_dp)
248 osq2 = 1.0_dp/sq2
249
250 IF (l < 0) cpabort("l < 0")
251 IF (ALLOCATED(cg_table)) THEN
252 DEALLOCATE (cg_table)
253 END IF
254 ! maximum size of table
255 n = (l**4 + 6*l**3 + 15*l**2 + 18*l + 8)/8
256 m = l + 1
257 ALLOCATE (cg_table(n, m, 2))
258 lmax = l
259
260 DO l1 = 0, lmax
261 DO m1 = 0, l1
262 iy = (l1*(l1 + 1))/2 + m1 + 1
263 DO l2 = l1, lmax
264 ml = 0
265 IF (l1 == l2) ml = m1
266 DO m2 = ml, l2
267 ix = (l2*(l2 + 1))/2 + m2 + 1
268 i1 = (ix*(ix - 1))/2 + iy
269 DO lp = mod(l1 + l2, 2), l1 + l2, 2
270 i2 = lp/2 + 1
271 mp = m2 + m1
272 cg_table(i1, i2, 1) = cgc(l1, m1, l2, m2, lp, mp)
273 mp = abs(m2 - m1)
274 IF (m2 >= m1) THEN
275 cg_table(i1, i2, 2) = cgc(l1, m1, lp, mp, l2, m2)
276 ELSE
277 cg_table(i1, i2, 2) = cgc(l2, m2, lp, mp, l1, m1)
278 END IF
279 END DO
280 END DO
281 END DO
282 END DO
283 END DO
284
285 CALL timestop(handle)
286
287 END SUBROUTINE clebsch_gordon_init
288
289! **************************************************************************************************
290!> \brief ...
291!> \param l1 ...
292!> \param m1 ...
293!> \param l2 ...
294!> \param m2 ...
295!> \param lp ...
296!> \param mp ...
297!> \return ...
298! **************************************************************************************************
299 FUNCTION cgc(l1, m1, l2, m2, lp, mp)
300 INTEGER, INTENT(IN) :: l1, m1, l2, m2, lp, mp
301 REAL(kind=dp) :: cgc
302
303 INTEGER :: la, lb, ll, ma, mb, s, t, tmax, tmin, &
304 z1, z2, z3, z4, z5
305 REAL(kind=dp) :: f1, f2, pref
306
307! m1 >= 0; m2 >= 0; mp >= 0
308
309 IF (m1 < 0 .OR. m2 < 0 .OR. mp < 0) THEN
310 WRITE (*, *) l1, l2, lp
311 WRITE (*, *) m1, m2, mp
312 cpabort("Illegal input values")
313 END IF
314 IF (l2 < l1) THEN
315 la = l2
316 ma = m2
317 lb = l1
318 mb = m1
319 ELSE
320 la = l1
321 ma = m1
322 lb = l2
323 mb = m2
324 END IF
325
326 IF (mod(la + lb + lp, 2) == 0 .AND. la + lb >= lp .AND. lp >= lb - la &
327 .AND. lb - mb >= 0) THEN
328 ll = (2*lp + 1)*(2*la + 1)*(2*lb + 1)
329 pref = 1.0_dp/sqrt(4.0_dp*pi)*0.5_dp*sqrt(real(ll, dp)* &
330 (sfac(lp - mp)/sfac(lp + mp))* &
331 (sfac(la - ma)/sfac(la + ma))*(sfac(lb - mb)/sfac(lb + mb)))
332 s = (la + lb + lp)/2
333 tmin = max(0, -lb + la - mp)
334 tmax = min(lb + la - mp, lp - mp, la - ma)
335 f1 = real(2*(-1)**(s - lb - ma), kind=dp)*(sfac(lb + mb)/sfac(lb - mb))* &
336 sfac(la + ma)/(sfac(s - lp)*sfac(s - lb))*sfac(2*s - 2*la)/sfac(s - la)* &
337 (sfac(s)/sfac(2*s + 1))
338 f2 = 0.0_dp
339 DO t = tmin, tmax
340 z1 = lp + mp + t
341 z2 = la + lb - mp - t
342 z3 = lp - mp - t
343 z4 = lb - la + mp + t
344 z5 = la - ma - t
345 f2 = f2 + (-1)**t*(sfac(z1)/(sfac(t)*sfac(z3)))*(sfac(z2)/(sfac(z4)*sfac(z5)))
346 END DO
347 cgc = pref*f1*f2
348 ELSE
349 cgc = 0.0_dp
350 END IF
351
352 END FUNCTION cgc
353
354! **************************************************************************************************
355!> \brief Calculate factorial even for integer values larger than the tabulated (pre-computed)
356!> values of up to 30!
357!> \param n Integer number for which the factorial has to be returned
358!> \return Factorial n!
359! **************************************************************************************************
360 FUNCTION sfac(n)
361 INTEGER, INTENT(IN) :: n
362 REAL(kind=dp) :: sfac
363
364 INTEGER :: i
365
366 IF (n > 170) THEN
367 cpabort("Factorials greater than 170! cannot be computed with double-precision")
368 ELSE IF (n > maxfac) THEN
369 sfac = fac(maxfac)
370 DO i = maxfac + 1, n
371 sfac = real(i, kind=dp)*sfac
372 END DO
373 ELSE IF (n >= 0) THEN
374 sfac = fac(n)
375 ELSE
376 sfac = 1.0_dp
377 END IF
378
379 END FUNCTION sfac
380
381! **************************************************************************************************
382!> \brief ...
383!> \param l1 ...
384!> \param m1 ...
385!> \param l2 ...
386!> \param m2 ...
387!> \return ...
388! **************************************************************************************************
389 FUNCTION order(l1, m1, l2, m2) RESULT(ind)
390 INTEGER, INTENT(IN) :: l1, m1, l2, m2
391 INTEGER :: ind
392
393 INTEGER :: i1, i2, ix, iy
394
395 i1 = (l1*(l1 + 1))/2 + abs(m1) + 1
396 i2 = (l2*(l2 + 1))/2 + abs(m2) + 1
397 ix = max(i1, i2)
398 iy = min(i1, i2)
399 ind = (ix*(ix - 1))/2 + iy
400 END FUNCTION order
401
402! Calculation of Spherical Harmonics
403
404! **************************************************************************************************
405!> \brief ...
406!> \param r ...
407!> \param y ...
408!> \param l ...
409!> \param m ...
410! **************************************************************************************************
411 SUBROUTINE rvy_lm(r, y, l, m)
412!
413! Real Spherical Harmonics
414! _ _
415! | [(2l+1)(l-|m|)!] |1/2 m cos(m p) m>=0
416! Y_lm ( t, p ) = |---------------------| P_l(cos(t))
417! |[2Pi(1+d_m0)(l+|m|)!]| sin(|m| p) m<0
418!
419! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1
420!
421 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r
422 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: y
423 INTEGER, INTENT(IN) :: l, m
424
425 INTEGER :: i
426 REAL(KIND=dp) :: cp, lmm, lpm, pf, plm, rxy, sp, t, z
427
428 SELECT CASE (l)
429 CASE (:-1)
430 cpabort("Negative l value")
431 CASE (0)
432 pf = sqrt(1.0_dp/(4.0_dp*pi))
433 IF (m /= 0) cpabort("l = 0 and m value out of bounds")
434 y(:) = pf
435 CASE (1)
436 SELECT CASE (m)
437 CASE DEFAULT
438 cpabort("l = 1 and m value out of bounds")
439 CASE (1)
440 pf = sqrt(3.0_dp/(4.0_dp*pi))
441 y(:) = pf*r(1, :)
442 CASE (0)
443 pf = sqrt(3.0_dp/(4.0_dp*pi))
444 y(:) = pf*r(3, :)
445 CASE (-1)
446 pf = sqrt(3.0_dp/(4.0_dp*pi))
447 y(:) = pf*r(2, :)
448 END SELECT
449 CASE (2)
450 SELECT CASE (m)
451 CASE DEFAULT
452 cpabort("l = 2 and m value out of bounds")
453 CASE (2)
454 pf = sqrt(15.0_dp/(16.0_dp*pi))
455 y(:) = pf*(r(1, :)*r(1, :) - r(2, :)*r(2, :))
456 CASE (1)
457 pf = sqrt(15.0_dp/(4.0_dp*pi))
458 y(:) = pf*r(3, :)*r(1, :)
459 CASE (0)
460 pf = sqrt(5.0_dp/(16.0_dp*pi))
461 y(:) = pf*(3.0_dp*r(3, :)*r(3, :) - 1.0_dp)
462 CASE (-1)
463 pf = sqrt(15.0_dp/(4.0_dp*pi))
464 y(:) = pf*r(3, :)*r(2, :)
465 CASE (-2)
466 pf = sqrt(15.0_dp/(16.0_dp*pi))
467 y(:) = pf*2.0_dp*r(1, :)*r(2, :)
468 END SELECT
469 CASE (3)
470 SELECT CASE (m)
471 CASE DEFAULT
472 cpabort("l = 3 and m value out of bounds")
473 CASE (3)
474 pf = sqrt(35.0_dp/(32.0_dp*pi))
475 y(:) = pf*r(1, :)*(r(1, :)**2 - 3.0_dp*r(2, :)**2)
476 CASE (2)
477 pf = sqrt(105.0_dp/(16.0_dp*pi))
478 y(:) = pf*r(3, :)*(r(1, :)**2 - r(2, :)**2)
479 CASE (1)
480 pf = sqrt(21.0_dp/(32.0_dp*pi))
481 y(:) = pf*r(1, :)*(5.0_dp*r(3, :)*r(3, :) - 1.0_dp)
482 CASE (0)
483 pf = sqrt(7.0_dp/(16.0_dp*pi))
484 y(:) = pf*r(3, :)*(5.0_dp*r(3, :)*r(3, :) - 3.0_dp)
485 CASE (-1)
486 pf = sqrt(21.0_dp/(32.0_dp*pi))
487 y(:) = pf*r(2, :)*(5.0_dp*r(3, :)*r(3, :) - 1.0_dp)
488 CASE (-2)
489 pf = sqrt(105.0_dp/(16.0_dp*pi))
490 y(:) = pf*2.0_dp*r(1, :)*r(2, :)*r(3, :)
491 CASE (-3)
492 pf = sqrt(35.0_dp/(32.0_dp*pi))
493 y(:) = pf*r(2, :)*(3.0_dp*r(1, :)**2 - r(2, :)**2)
494 END SELECT
495 CASE DEFAULT
496 IF (m < -l .OR. m > l) cpabort("m value out of bounds")
497 lpm = fac(l + abs(m))
498 lmm = fac(l - abs(m))
499 IF (m == 0) THEN
500 t = 4.0_dp*pi
501 ELSE
502 t = 2.0_dp*pi
503 END IF
504 IF (abs(lpm) < epsilon(1.0_dp)) THEN
505 pf = real(2*l + 1, kind=dp)/t
506 ELSE
507 pf = (real(2*l + 1, kind=dp)*lmm)/(t*lpm)
508 END IF
509 pf = sqrt(pf)
510 DO i = 1, SIZE(r, 2)
511 z = r(3, i)
512! plm = legendre ( z, l, m )
513 plm = legendre(z, l, abs(m))
514 IF (m == 0) THEN
515 y(i) = pf*plm
516 ELSE
517 rxy = sqrt(r(1, i)**2 + r(2, i)**2)
518 IF (rxy < epsilon(1.0_dp)) THEN
519 y(i) = 0.0_dp
520 ELSE
521 cp = r(1, i)/rxy
522 sp = r(2, i)/rxy
523 IF (m > 0) THEN
524 y(i) = pf*plm*cosn(m, cp, sp)
525 ELSE
526 y(i) = pf*plm*sinn(-m, cp, sp)
527 END IF
528 END IF
529 END IF
530 END DO
531 END SELECT
532
533 END SUBROUTINE rvy_lm
534
535! **************************************************************************************************
536!> \brief ...
537!> \param r ...
538!> \param y ...
539!> \param l ...
540!> \param m ...
541! **************************************************************************************************
542 SUBROUTINE rry_lm(r, y, l, m)
543!
544! Real Spherical Harmonics
545! _ _
546! | [(2l+1)(l-|m|)!] |1/2 m cos(m p) m>=0
547! Y_lm ( t, p ) = |---------------------| P_l(cos(t))
548! |[2Pi(1+d_m0)(l+|m|)!]| sin(|m| p) m<0
549!
550! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1
551!
552 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r
553 REAL(KIND=dp), INTENT(OUT) :: y
554 INTEGER, INTENT(IN) :: l, m
555
556 REAL(KIND=dp) :: cp, lmm, lpm, pf, plm, rxy, sp, t, z
557
558 SELECT CASE (l)
559 CASE (:-1)
560 cpabort("Negative l value")
561 CASE (0)
562 pf = sqrt(1.0_dp/(4.0_dp*pi))
563 IF (m /= 0) cpabort("l = 0 and m value out of bounds")
564 y = pf
565 CASE (1)
566 SELECT CASE (m)
567 CASE DEFAULT
568 cpabort("l = 1 and m value out of bounds")
569 CASE (1)
570 pf = sqrt(3.0_dp/(4.0_dp*pi))
571 y = pf*r(1)
572 CASE (0)
573 pf = sqrt(3.0_dp/(4.0_dp*pi))
574 y = pf*r(3)
575 CASE (-1)
576 pf = sqrt(3.0_dp/(4.0_dp*pi))
577 y = pf*r(2)
578 END SELECT
579 CASE (2)
580 SELECT CASE (m)
581 CASE DEFAULT
582 cpabort("l = 2 and m value out of bounds")
583 CASE (2)
584 pf = sqrt(15.0_dp/(16.0_dp*pi))
585 y = pf*(r(1)*r(1) - r(2)*r(2))
586 CASE (1)
587 pf = sqrt(15.0_dp/(4.0_dp*pi))
588 y = pf*r(3)*r(1)
589 CASE (0)
590 pf = sqrt(5.0_dp/(16.0_dp*pi))
591 y = pf*(3.0_dp*r(3)*r(3) - 1.0_dp)
592 CASE (-1)
593 pf = sqrt(15.0_dp/(4.0_dp*pi))
594 y = pf*r(3)*r(2)
595 CASE (-2)
596 pf = sqrt(15.0_dp/(16.0_dp*pi))
597 y = pf*2.0_dp*r(1)*r(2)
598 END SELECT
599 CASE (3)
600 SELECT CASE (m)
601 CASE DEFAULT
602 cpabort("l = 3 and m value out of bounds")
603 CASE (3)
604 pf = sqrt(35.0_dp/(32.0_dp*pi))
605 y = pf*r(1)*(r(1)**2 - 3.0_dp*r(2)**2)
606 CASE (2)
607 pf = sqrt(105.0_dp/(16.0_dp*pi))
608 y = pf*r(3)*(r(1)**2 - r(2)**2)
609 CASE (1)
610 pf = sqrt(21.0_dp/(32.0_dp*pi))
611 y = pf*r(1)*(5.0_dp*r(3)*r(3) - 1.0_dp)
612 CASE (0)
613 pf = sqrt(7.0_dp/(16.0_dp*pi))
614 y = pf*r(3)*(5.0_dp*r(3)*r(3) - 3.0_dp)
615 CASE (-1)
616 pf = sqrt(21.0_dp/(32.0_dp*pi))
617 y = pf*r(2)*(5.0_dp*r(3)*r(3) - 1.0_dp)
618 CASE (-2)
619 pf = sqrt(105.0_dp/(16.0_dp*pi))
620 y = pf*2.0_dp*r(1)*r(2)*r(3)
621 CASE (-3)
622 pf = sqrt(35.0_dp/(32.0_dp*pi))
623 y = pf*r(2)*(3.0_dp*r(1)**2 - r(2)**2)
624 END SELECT
625 CASE DEFAULT
626 IF (m < -l .OR. m > l) cpabort("m value out of bounds")
627 lpm = fac(l + abs(m))
628 lmm = fac(l - abs(m))
629 IF (m == 0) THEN
630 t = 4.0_dp*pi
631 ELSE
632 t = 2.0_dp*pi
633 END IF
634 IF (abs(lpm) < epsilon(1.0_dp)) THEN
635 pf = real(2*l + 1, kind=dp)/t
636 ELSE
637 pf = (real(2*l + 1, kind=dp)*lmm)/(t*lpm)
638 END IF
639 pf = sqrt(pf)
640 z = r(3)
641 plm = legendre(z, l, m)
642 IF (m == 0) THEN
643 y = pf*plm
644 ELSE
645 rxy = sqrt(r(1)**2 + r(2)**2)
646 IF (rxy < epsilon(1.0_dp)) THEN
647 y = 0.0_dp
648 ELSE
649 cp = r(1)/rxy
650 sp = r(2)/rxy
651 IF (m > 0) THEN
652 y = pf*plm*cosn(m, cp, sp)
653 ELSE
654 y = pf*plm*sinn(-m, cp, sp)
655 END IF
656 END IF
657 END IF
658 END SELECT
659
660 END SUBROUTINE rry_lm
661
662! **************************************************************************************************
663!> \brief ...
664!> \param r ...
665!> \param y ...
666!> \param l ...
667!> \param m ...
668! **************************************************************************************************
669 SUBROUTINE cvy_lm(r, y, l, m)
670!
671! Complex Spherical Harmonics
672! _ _
673! | [(2l+1)(l-|m|)!] |1/2 m
674! Y_lm ( t, p ) = |------------------| P_l(cos(t)) Exp[ i m p ]
675! | [4Pi(l+|m|)!]| |
676!
677! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1
678!
679 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r
680 COMPLEX(KIND=dp), DIMENSION(:), INTENT(OUT) :: y
681 INTEGER, INTENT(IN) :: l, m
682
683 INTEGER :: i, n
684 REAL(KIND=dp) :: cp, lmm, lpm, pf, plm, rxy, sp, t, ym, &
685 yp, z
686
687 n = SIZE(r, 2)
688 SELECT CASE (l)
689 CASE (:-1)
690 cpabort("Negative l value")
691 CASE (0)
692 pf = sqrt(1.0_dp/(4.0_dp*pi))
693 IF (m /= 0) cpabort("l = 0 and m value out of bounds")
694 y(:) = pf
695 CASE (1)
696 SELECT CASE (m)
697 CASE DEFAULT
698 cpabort("l = 1 and m value out of bounds")
699 CASE (1)
700 pf = sqrt(3.0_dp/(8.0_dp*pi))
701 DO i = 1, n
702 yp = r(1, i)
703 ym = r(2, i)
704 y(i) = pf*cmplx(yp, ym, kind=dp)
705 END DO
706 CASE (0)
707 pf = sqrt(3.0_dp/(4.0_dp*pi))
708 y(:) = pf*r(3, :)
709 CASE (-1)
710 pf = sqrt(3.0_dp/(8.0_dp*pi))
711 DO i = 1, n
712 yp = r(1, i)
713 ym = r(2, i)
714 y(i) = pf*cmplx(yp, -ym, kind=dp)
715 END DO
716 END SELECT
717 CASE (2)
718 SELECT CASE (m)
719 CASE DEFAULT
720 cpabort("l = 2 and m value out of bounds")
721 CASE (2)
722 pf = sqrt(15.0_dp/(32.0_dp*pi))
723 DO i = 1, n
724 yp = (r(1, i)*r(1, i) - r(2, i)*r(2, i))
725 ym = 2.0_dp*r(1, i)*r(2, i)
726 y(i) = pf*cmplx(yp, ym, kind=dp)
727 END DO
728 CASE (1)
729 pf = sqrt(15.0_dp/(8.0_dp*pi))
730 DO i = 1, n
731 yp = r(3, i)*r(1, i)
732 ym = r(3, i)*r(2, i)
733 y(i) = pf*cmplx(yp, ym, kind=dp)
734 END DO
735 CASE (0)
736 pf = sqrt(5.0_dp/(16.0_dp*pi))
737 y(:) = pf*(3.0_dp*r(3, :)*r(3, :) - 1.0_dp)
738 CASE (-1)
739 pf = sqrt(15.0_dp/(8.0_dp*pi))
740 DO i = 1, n
741 yp = r(3, i)*r(1, i)
742 ym = r(3, i)*r(2, i)
743 y(i) = pf*cmplx(yp, -ym, kind=dp)
744 END DO
745 CASE (-2)
746 pf = sqrt(15.0_dp/(32.0_dp*pi))
747 DO i = 1, n
748 yp = (r(1, i)*r(1, i) - r(2, i)*r(2, i))
749 ym = 2.0_dp*r(1, i)*r(2, i)
750 y(i) = pf*cmplx(yp, -ym, kind=dp)
751 END DO
752 END SELECT
753 CASE (3)
754 SELECT CASE (m)
755 CASE DEFAULT
756 cpabort("l = 3 and m value out of bounds")
757 CASE (3)
758 pf = sqrt(35.0_dp/(64.0_dp*pi))
759 DO i = 1, n
760 yp = r(1, i)*(r(1, i)**2 - 3.0_dp*r(2, i)**2)
761 ym = r(2, i)*(3.0_dp*r(1, i)**2 - r(2, i)**2)
762 y(i) = pf*cmplx(yp, ym, kind=dp)
763 END DO
764 CASE (2)
765 pf = sqrt(105.0_dp/(32.0_dp*pi))
766 DO i = 1, n
767 yp = r(3, i)*(r(1, i)**2 - r(2, i)**2)
768 ym = 2.0_dp*r(1, i)*r(2, i)*r(3, i)
769 y(i) = pf*cmplx(yp, ym, kind=dp)
770 END DO
771 CASE (1)
772 pf = sqrt(21.0_dp/(64.0_dp*pi))
773 DO i = 1, n
774 yp = r(1, i)*(5.0_dp*r(3, i)*r(3, i) - 1.0_dp)
775 ym = r(2, i)*(5.0_dp*r(3, i)*r(3, i) - 1.0_dp)
776 y(i) = pf*cmplx(yp, ym, kind=dp)
777 END DO
778 CASE (0)
779 pf = sqrt(7.0_dp/(16.0_dp*pi))
780 y(:) = pf*r(3, :)*(5.0_dp*r(3, :)*r(3, :) - 3.0_dp)
781 CASE (-1)
782 pf = sqrt(21.0_dp/(64.0_dp*pi))
783 DO i = 1, n
784 yp = r(1, i)*(5.0_dp*r(3, i)*r(3, i) - 1.0_dp)
785 ym = r(2, i)*(5.0_dp*r(3, i)*r(3, i) - 1.0_dp)
786 y(i) = pf*cmplx(yp, -ym, kind=dp)
787 END DO
788 CASE (-2)
789 pf = sqrt(105.0_dp/(32.0_dp*pi))
790 DO i = 1, n
791 yp = r(3, i)*(r(1, i)**2 - r(2, i)**2)
792 ym = 2.0_dp*r(1, i)*r(2, i)*r(3, i)
793 y(i) = pf*cmplx(yp, -ym, kind=dp)
794 END DO
795 CASE (-3)
796 pf = sqrt(35.0_dp/(64.0_dp*pi))
797 DO i = 1, n
798 yp = r(1, i)*(r(1, i)**2 - 3.0_dp*r(2, i)**2)
799 ym = r(2, i)*(3.0_dp*r(1, i)**2 - r(2, i)**2)
800 y(i) = pf*cmplx(yp, -ym, kind=dp)
801 END DO
802 END SELECT
803 CASE DEFAULT
804 IF (m < -l .OR. m > l) cpabort("m value out of bounds")
805 lpm = fac(l + abs(m))
806 lmm = fac(l - abs(m))
807 t = 4.0_dp*pi
808 IF (abs(lpm) < epsilon(1.0_dp)) THEN
809 pf = real(2*l + 1, kind=dp)/t
810 ELSE
811 pf = (real(2*l + 1, kind=dp)*lmm)/(t*lpm)
812 END IF
813 pf = sqrt(pf)
814 DO i = 1, n
815 z = r(3, i)
816 plm = legendre(z, l, m)
817 IF (m == 0) THEN
818 y(i) = pf*plm
819 ELSE
820 rxy = sqrt(r(1, i)**2 + r(2, i)**2)
821 IF (rxy < epsilon(1.0_dp)) THEN
822 y(i) = 0.0_dp
823 ELSE
824 cp = r(1, i)/rxy
825 sp = r(2, i)/rxy
826 IF (m > 0) THEN
827 yp = cosn(m, cp, sp)
828 ym = sinn(m, cp, sp)
829 y(i) = pf*plm*cmplx(yp, ym, kind=dp)
830 ELSE
831 yp = cosn(-m, cp, sp)
832 ym = sinn(-m, cp, sp)
833 y(i) = pf*plm*cmplx(yp, -ym, kind=dp)
834 END IF
835 END IF
836 END IF
837 END DO
838 END SELECT
839
840 END SUBROUTINE cvy_lm
841
842! **************************************************************************************************
843!> \brief ...
844!> \param r ...
845!> \param y ...
846!> \param l ...
847!> \param m ...
848! **************************************************************************************************
849 SUBROUTINE ccy_lm(r, y, l, m)
850!
851! Complex Spherical Harmonics
852! _ _
853! | [(2l+1)(l-|m|)!] |1/2 m
854! Y_lm ( t, p ) = |------------------| P_l(cos(t)) Exp[ i m p ]
855! | [4Pi(l+|m|)!]| |
856!
857! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1
858!
859 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r
860 COMPLEX(KIND=dp), INTENT(OUT) :: y
861 INTEGER, INTENT(IN) :: l, m
862
863 REAL(KIND=dp) :: cp, lmm, lpm, pf, plm, rxy, sp, t, ym, &
864 yp, z
865
866 SELECT CASE (l)
867 CASE (:-1)
868 cpabort("Negative l value")
869 CASE (0)
870 pf = sqrt(1.0_dp/(4.0_dp*pi))
871 IF (m /= 0) cpabort("l = 0 and m value out of bounds")
872 y = pf
873 CASE (1)
874 SELECT CASE (m)
875 CASE DEFAULT
876 cpabort("l = 1 and m value out of bounds")
877 CASE (1)
878 pf = sqrt(3.0_dp/(8.0_dp*pi))
879 yp = r(1)
880 ym = r(2)
881 y = pf*cmplx(yp, ym, kind=dp)
882 CASE (0)
883 pf = sqrt(3.0_dp/(4.0_dp*pi))
884 y = pf*r(3)
885 CASE (-1)
886 pf = sqrt(3.0_dp/(8.0_dp*pi))
887 yp = r(1)
888 ym = r(2)
889 y = pf*cmplx(yp, -ym, kind=dp)
890 END SELECT
891 CASE (2)
892 SELECT CASE (m)
893 CASE DEFAULT
894 cpabort("l = 2 and m value out of bounds")
895 CASE (2)
896 pf = sqrt(15.0_dp/(32.0_dp*pi))
897 yp = (r(1)*r(1) - r(2)*r(2))
898 ym = 2.0_dp*r(1)*r(2)
899 y = pf*cmplx(yp, ym, kind=dp)
900 CASE (1)
901 pf = sqrt(15.0_dp/(8.0_dp*pi))
902 yp = r(3)*r(1)
903 ym = r(3)*r(2)
904 y = pf*cmplx(yp, ym, kind=dp)
905 CASE (0)
906 pf = sqrt(5.0_dp/(16.0_dp*pi))
907 y = pf*(3.0_dp*r(3)*r(3) - 1.0_dp)
908 CASE (-1)
909 pf = sqrt(15.0_dp/(8.0_dp*pi))
910 yp = r(3)*r(1)
911 ym = r(3)*r(2)
912 y = pf*cmplx(yp, -ym, kind=dp)
913 CASE (-2)
914 pf = sqrt(15.0_dp/(32.0_dp*pi))
915 yp = (r(1)*r(1) - r(2)*r(2))
916 ym = 2.0_dp*r(1)*r(2)
917 y = pf*cmplx(yp, -ym, kind=dp)
918 END SELECT
919 CASE (3)
920 SELECT CASE (m)
921 CASE DEFAULT
922 cpabort("l = 3 and m value out of bounds")
923 CASE (3)
924 pf = sqrt(35.0_dp/(64.0_dp*pi))
925 yp = r(1)*(r(1)**2 - 3.0_dp*r(2)**2)
926 ym = r(2)*(3.0_dp*r(1)**2 - r(2)**2)
927 y = pf*cmplx(yp, ym, kind=dp)
928 CASE (2)
929 pf = sqrt(105.0_dp/(32.0_dp*pi))
930 yp = r(3)*(r(1)**2 - r(2)**2)
931 ym = 2.0_dp*r(1)*r(2)*r(3)
932 y = pf*cmplx(yp, ym, kind=dp)
933 CASE (1)
934 pf = sqrt(21.0_dp/(64.0_dp*pi))
935 yp = r(1)*(5.0_dp*r(3)*r(3) - 1.0_dp)
936 ym = r(2)*(5.0_dp*r(3)*r(3) - 1.0_dp)
937 y = pf*cmplx(yp, ym, kind=dp)
938 CASE (0)
939 pf = sqrt(7.0_dp/(16.0_dp*pi))
940 y = pf*r(3)*(5.0_dp*r(3)*r(3) - 3.0_dp)
941 CASE (-1)
942 pf = sqrt(21.0_dp/(64.0_dp*pi))
943 yp = r(1)*(5.0_dp*r(3)*r(3) - 1.0_dp)
944 ym = r(2)*(5.0_dp*r(3)*r(3) - 1.0_dp)
945 y = pf*cmplx(yp, -ym, kind=dp)
946 CASE (-2)
947 pf = sqrt(105.0_dp/(32.0_dp*pi))
948 yp = r(3)*(r(1)**2 - r(2)**2)
949 ym = 2.0_dp*r(1)*r(2)*r(3)
950 y = pf*cmplx(yp, -ym, kind=dp)
951 CASE (-3)
952 pf = sqrt(35.0_dp/(64.0_dp*pi))
953 yp = r(1)*(r(1)**2 - 3.0_dp*r(2)**2)
954 ym = r(2)*(3.0_dp*r(1)**2 - r(2)**2)
955 y = pf*cmplx(yp, -ym, kind=dp)
956 END SELECT
957 CASE DEFAULT
958 IF (m < -l .OR. m > l) cpabort("m value out of bounds")
959 lpm = fac(l + abs(m))
960 lmm = fac(l - abs(m))
961 t = 4.0_dp*pi
962 IF (abs(lpm) < epsilon(1.0_dp)) THEN
963 pf = real(2*l + 1, kind=dp)/t
964 ELSE
965 pf = (real(2*l + 1, kind=dp)*lmm)/(t*lpm)
966 END IF
967 pf = sqrt(pf)
968 z = r(3)
969 plm = legendre(z, l, m)
970 IF (m == 0) THEN
971 y = pf*plm
972 ELSE
973 rxy = sqrt(r(1)**2 + r(2)**2)
974 IF (rxy < epsilon(1.0_dp)) THEN
975 y = 0.0_dp
976 ELSE
977 cp = r(1)/rxy
978 sp = r(2)/rxy
979 IF (m > 0) THEN
980 yp = cosn(m, cp, sp)
981 ym = sinn(m, cp, sp)
982 y = pf*plm*cmplx(yp, ym, kind=dp)
983 ELSE
984 yp = cosn(-m, cp, sp)
985 ym = sinn(-m, cp, sp)
986 y = pf*plm*cmplx(yp, -ym, kind=dp)
987 END IF
988 END IF
989 END IF
990 END SELECT
991
992 END SUBROUTINE ccy_lm
993
994! Calculation of derivatives of Spherical Harmonics
995
996! **************************************************************************************************
997!> \brief ...
998!> \param c ...
999!> \param dy ...
1000!> \param l ...
1001!> \param m ...
1002! **************************************************************************************************
1003 SUBROUTINE dry_lm(c, dy, l, m)
1004!
1005! Real Spherical Harmonics
1006! _ _
1007! | [(2l+1)(l-|m|)!] |1/2 m cos(m p) m>=0
1008! Y_lm ( t, p ) = |---------------------| P_l(cos(t))
1009! |[2Pi(1+d_m0)(l+|m|)!]| sin(|m| p) m<0
1010!
1011! Input: c == (t,p)
1012! Output: dy == (dy/dt, dy/dp)
1013!
1014! x == sin(t)*cos(p)
1015! y == sin(t)*sin(p)
1016! z == cos(t)
1017!
1018
1019 REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: c
1020 REAL(KIND=dp), DIMENSION(2), INTENT(OUT) :: dy
1021 INTEGER, INTENT(IN) :: l, m
1022
1023 REAL(KIND=dp) :: cp, ct, dplm, lmm, lpm, p, pf, rxy, sp, &
1024 st, t, tt, y, z
1025 REAL(KIND=dp), DIMENSION(3) :: r
1026
1027 t = c(1)
1028 ct = cos(t)
1029 st = sin(t)
1030 p = c(2)
1031 cp = cos(p)
1032 sp = sin(p)
1033 r(1) = st*cp
1034 r(2) = st*sp
1035 r(3) = ct
1036
1037! dY/dp
1038 IF (m == 0) THEN
1039 dy(2) = 0.0_dp
1040 ELSE
1041 CALL y_lm(r, y, l, -m)
1042 dy(2) = -real(m, kind=dp)*y
1043 END IF
1044
1045! dY/dt
1046 SELECT CASE (l)
1047 CASE (:-1)
1048 cpabort("Negative l value")
1049 CASE (0)
1050 IF (m /= 0) cpabort("l = 0 and m value out of bounds")
1051 dy(1) = 0.0_dp
1052 CASE (1)
1053 SELECT CASE (m)
1054 CASE DEFAULT
1055 cpabort("l = 1 and m value out of bounds")
1056 CASE (1)
1057 pf = sqrt(3.0_dp/(4.0_dp*pi))
1058 dy(1) = pf*ct*cp
1059 CASE (0)
1060 pf = sqrt(3.0_dp/(4.0_dp*pi))
1061 dy(1) = -pf*st
1062 CASE (-1)
1063 pf = sqrt(3.0_dp/(4.0_dp*pi))
1064 dy(1) = pf*ct*sp
1065 END SELECT
1066 CASE (2)
1067 SELECT CASE (m)
1068 CASE DEFAULT
1069 cpabort("l = 2 and m value out of bounds")
1070 CASE (2)
1071 pf = sqrt(15.0_dp/(16.0_dp*pi))
1072 dy(1) = pf*2.0_dp*st*ct*cos(2._dp*p)
1073 CASE (1)
1074 pf = sqrt(15.0_dp/(4.0_dp*pi))
1075 dy(1) = pf*cp*(ct*ct - st*st)
1076 CASE (0)
1077 pf = sqrt(5.0_dp/(16.0_dp*pi))
1078 dy(1) = -pf*6.0_dp*ct*st
1079 CASE (-1)
1080 pf = sqrt(15.0_dp/(4.0_dp*pi))
1081 dy(1) = pf*sp*(ct*ct - st*st)
1082 CASE (-2)
1083 pf = sqrt(15.0_dp/(16.0_dp*pi))
1084 dy(1) = pf*2.0_dp*st*ct*sin(2._dp*p)
1085 END SELECT
1086 CASE (3)
1087 SELECT CASE (m)
1088 CASE DEFAULT
1089 cpabort("l = 3 and m value out of bounds")
1090 CASE (3)
1091 pf = sqrt(35.0_dp/(32.0_dp*pi))
1092 dy(1) = pf*3.0_dp*cos(3._dp*p)*ct*st*st
1093 CASE (2)
1094 pf = sqrt(105.0_dp/(16.0_dp*pi))
1095 dy(1) = pf*2.0_dp*cos(2._dp*p)*ct*st
1096 CASE (1)
1097 pf = sqrt(21.0_dp/(32.0_dp*pi))
1098 dy(1) = pf*cp*(ct*(5.0_dp*ct - 1.0_dp) - 5.0_dp*st*st)
1099 CASE (0)
1100 pf = sqrt(7.0_dp/(16.0_dp*pi))
1101 dy(1) = pf*r(3)*(3.0_dp - 15.0_dp*ct*ct)*st
1102 CASE (-1)
1103 pf = sqrt(21.0_dp/(32.0_dp*pi))
1104 dy(1) = pf*sp*(ct*(5.0_dp*ct - 1.0_dp) - 5.0_dp*st*st)
1105 CASE (-2)
1106 pf = sqrt(105.0_dp/(16.0_dp*pi))
1107 dy(1) = pf*2.0_dp*sin(2._dp*p)*ct*st
1108 CASE (-3)
1109 pf = sqrt(35.0_dp/(32.0_dp*pi))
1110 dy(1) = pf*3.0_dp*sin(3._dp*p)*ct*st*st
1111 END SELECT
1112 CASE DEFAULT
1113 IF (m < -l .OR. m > l) cpabort("m value out of bounds")
1114 lpm = fac(l + abs(m))
1115 lmm = fac(l - abs(m))
1116 IF (m == 0) THEN
1117 tt = 4.0_dp*pi
1118 ELSE
1119 tt = 2.0_dp*pi
1120 END IF
1121 IF (abs(lpm) < epsilon(1.0_dp)) THEN
1122 pf = real(2*l + 1, kind=dp)/tt
1123 ELSE
1124 pf = (real(2*l + 1, kind=dp)*lmm)/(tt*lpm)
1125 END IF
1126 pf = sqrt(pf)
1127 z = ct
1128 dplm = dlegendre(z, l, m)
1129 IF (m == 0) THEN
1130 y = pf*dplm
1131 ELSE
1132 rxy = sqrt(r(1)**2 + r(2)**2)
1133 IF (rxy < epsilon(1.0_dp)) THEN
1134 y = 0.0_dp
1135 ELSE
1136 IF (m > 0) THEN
1137 y = pf*dplm*cosn(m, cp, sp)
1138 ELSE
1139 y = pf*dplm*sinn(-m, cp, sp)
1140 END IF
1141 END IF
1142 END IF
1143 END SELECT
1144
1145 END SUBROUTINE dry_lm
1146
1147! **************************************************************************************************
1148!> \brief ...
1149!> \param c ...
1150!> \param dy ...
1151!> \param l ...
1152!> \param m ...
1153! **************************************************************************************************
1154 SUBROUTINE dcy_lm(c, dy, l, m)
1155!
1156! Complex Spherical Harmonics
1157! _ _
1158! | [(2l+1)(l-|m|)!] |1/2 m
1159! Y_lm ( t, p ) = |------------------| P_l(cos(t)) Exp[ i m p ]
1160! | [4Pi(l+|m|)!]| |
1161!
1162! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1
1163!
1164! Input: c == (t,p)
1165! Output: dy == (dy/dt, dy/dp)
1166!
1167! x == sin(t)*cos(p)
1168! y == sin(t)*sin(p)
1169! z == cos(t)
1170!
1171
1172 REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: c
1173 COMPLEX(KIND=dp), DIMENSION(2), INTENT(OUT) :: dy
1174 INTEGER, INTENT(IN) :: l, m
1175
1176 REAL(KIND=dp), DIMENSION(2) :: dd
1177
1178 cpabort("Not implemented")
1179 CALL dry_lm(c, dd, l, m)
1180 dy(1) = cmplx(dd(1), 0.0_dp, kind=dp)
1181 dy(2) = cmplx(dd(2), 0.0_dp, kind=dp)
1182
1183 END SUBROUTINE dcy_lm
1184
1185! **************************************************************************************************
1186!> \brief ...
1187!> \param x ...
1188!> \param l ...
1189!> \param m ...
1190!> \return ...
1191! **************************************************************************************************
1192 FUNCTION legendre(x, l, m) RESULT(plm)
1193
1194 REAL(kind=dp), INTENT(IN) :: x
1195 INTEGER, INTENT(IN) :: l, m
1196 REAL(kind=dp) :: plm
1197
1198 INTEGER :: il, im, mm
1199 REAL(kind=dp) :: fact, pll, pmm, pmmp1, somx2
1200
1201 IF (abs(x) > 1.0_dp) cpabort("x value > 1")
1202 SELECT CASE (l)
1203 CASE (:-1)
1204 cpabort("Negative l value")
1205 CASE (0)
1206 plm = 1.0_dp
1207 CASE (1)
1208 SELECT CASE (abs(m))
1209 CASE DEFAULT
1210 cpabort("l = 1 and m value out of bounds")
1211 CASE (1)
1212 plm = sqrt(1.0_dp - x*x)
1213 CASE (0)
1214 plm = x
1215 END SELECT
1216 CASE (2)
1217 SELECT CASE (abs(m))
1218 CASE DEFAULT
1219 cpabort("l = 2 and m value out of bounds")
1220 CASE (2)
1221 plm = 3.0_dp*(1.0_dp - x*x)
1222 CASE (1)
1223 plm = 3.0_dp*x*sqrt(1.0_dp - x*x)
1224 CASE (0)
1225 plm = 1.5_dp*x*x - 0.5_dp
1226 END SELECT
1227 CASE (3)
1228 SELECT CASE (abs(m))
1229 CASE DEFAULT
1230 cpabort("l = 3 and m value out of bounds")
1231 CASE (3)
1232 plm = 15.0_dp*(1.0_dp - x*x)**1.5_dp
1233 CASE (2)
1234 plm = 15.0_dp*x*(1.0_dp - x*x)
1235 CASE (1)
1236 plm = (7.5_dp*x*x - 1.5_dp)*sqrt(1.0_dp - x*x)
1237 CASE (0)
1238 plm = 2.5_dp*x**3 - 1.5_dp*x
1239 END SELECT
1240 CASE (4)
1241 SELECT CASE (abs(m))
1242 CASE DEFAULT
1243 cpabort("l = 4 and m value out of bounds")
1244 CASE (4)
1245 plm = 105.0_dp*(1.0_dp - x*x)**2
1246 CASE (3)
1247 plm = 105.0_dp*x*(1.0_dp - x*x)**1.5_dp
1248 CASE (2)
1249 plm = (52.5_dp*x*x - 7.5_dp)*(1.0_dp - x*x)
1250 CASE (1)
1251 plm = (17.5_dp*x**3 - 7.5_dp*x)*sqrt(1.0_dp - x*x)
1252 CASE (0)
1253 plm = 4.375_dp*x**4 - 3.75_dp*x**2 + 0.375_dp
1254 END SELECT
1255 CASE (5)
1256 SELECT CASE (abs(m))
1257 CASE DEFAULT
1258 cpabort("l = 5 and m value out of bounds")
1259 CASE (5)
1260 plm = 945.0_dp*(1.0_dp - x*x)**2.5_dp
1261 CASE (4)
1262 plm = 945.0_dp*x*(1.0_dp - x*x)**2
1263 CASE (3)
1264 plm = -(-472.5_dp*x*x + 52.5_dp)*(1.0_dp - x*x)**1.5_dp
1265 CASE (2)
1266 plm = (157.5_dp*x**3 - 52.5_dp*x)*(1.0_dp - x*x)
1267 CASE (1)
1268 plm = -(-39.375_dp*x**4 + 26.25_dp*x*x - &
1269 1.875_dp)*sqrt(1.0_dp - x*x)
1270 CASE (0)
1271 plm = 7.875_dp*x**5 - 8.75_dp*x**3 + 1.875_dp*x
1272 END SELECT
1273 CASE (6)
1274 SELECT CASE (abs(m))
1275 CASE DEFAULT
1276 cpabort("l = 6 and m value out of bounds")
1277 CASE (6)
1278 plm = 10395.0_dp*(1.0_dp - x*x)**3
1279 CASE (5)
1280 plm = 10395.0_dp*x*(1.0_dp - x*x)**2.5_dp
1281 CASE (4)
1282 plm = (5197.5_dp*x*x - 472.5_dp)*(1.0_dp - x*x)**2
1283 CASE (3)
1284 plm = -(-1732.5_dp*x**3 + 472.5_dp*x)* &
1285 (1.0_dp - x*x)**1.5_dp
1286 CASE (2)
1287 plm = (433.125_dp*x**4 - 236.25_dp*x**2 + &
1288 13.125_dp)*(1.0_dp - x*x)
1289 CASE (1)
1290 plm = -(-86.625_dp*x**5 + 78.75_dp*x**3 - &
1291 13.125_dp*x)*sqrt(1.0_dp - x*x)
1292 CASE (0)
1293 plm = 14.4375_dp*x**6 - 19.6875_dp*x**4 + &
1294 6.5625_dp*x**2 - 0.3125_dp
1295 END SELECT
1296 CASE DEFAULT
1297 mm = abs(m)
1298 IF (mm > l) cpabort("m out of bounds")
1299! use recurence from numerical recipies
1300 pmm = 1.0_dp
1301 IF (mm > 0) THEN
1302 somx2 = sqrt((1.0_dp - x)*(1.0_dp + x))
1303 fact = 1.0_dp
1304 DO im = 1, mm
1305 pmm = pmm*fact*somx2
1306 fact = fact + 2.0_dp
1307 END DO
1308 END IF
1309 IF (l == mm) THEN
1310 plm = pmm
1311 ELSE
1312 pmmp1 = x*real(2*mm + 1, kind=dp)*pmm
1313 IF (l == mm + 1) THEN
1314 plm = pmmp1
1315 ELSE
1316 DO il = mm + 2, l
1317 pll = (x*real(2*il - 1, kind=dp)*pmmp1 - &
1318 REAL(il + mm - 1, kind=dp)*pmm)/real(il - mm, kind=dp)
1319 pmm = pmmp1
1320 pmmp1 = pll
1321 END DO
1322 plm = pll
1323 END IF
1324 END IF
1325 END SELECT
1326
1327 END FUNCTION legendre
1328
1329! **************************************************************************************************
1330!> \brief ...
1331!> \param x ...
1332!> \param l ...
1333!> \param m ...
1334!> \return ...
1335! **************************************************************************************************
1336 FUNCTION dlegendre(x, l, m) RESULT(dplm)
1337 REAL(kind=dp), INTENT(IN) :: x
1338 INTEGER, INTENT(IN) :: l, m
1339 REAL(kind=dp) :: dplm
1340
1341 INTEGER :: mm
1342
1343 IF (abs(x) > 1.0_dp) cpabort("x value > 1")
1344 SELECT CASE (l)
1345 CASE (0)
1346 dplm = 0.0_dp
1347 CASE (1)
1348 SELECT CASE (abs(m))
1349 CASE DEFAULT
1350 cpabort("l = 1 and m value out of bounds")
1351 CASE (1)
1352 dplm = -x/sqrt(1.0_dp - x*x)
1353 CASE (0)
1354 dplm = 1.0_dp
1355 END SELECT
1356 CASE (2)
1357 SELECT CASE (abs(m))
1358 CASE DEFAULT
1359 cpabort("l = 2 and m value out of bounds")
1360 CASE (2)
1361 dplm = -6.0_dp*x
1362 CASE (1)
1363 dplm = 3.0_dp*sqrt(1.0_dp - x*x) - 3.0_dp*x*x/sqrt(1.0_dp - x*x)
1364 CASE (0)
1365 dplm = 3.0_dp*x
1366 END SELECT
1367 CASE (3)
1368 SELECT CASE (abs(m))
1369 CASE DEFAULT
1370 cpabort("l = 3 and m value out of bounds")
1371 CASE (3)
1372 dplm = -45.0_dp*sqrt(1.0_dp - x*x)*x
1373 CASE (2)
1374 dplm = 15.0_dp*(1.0_dp - x*x) - 30.0_dp*x*x
1375 CASE (1)
1376 dplm = 15.0_dp*x*sqrt(1.0_dp - x*x) - (x*(7.5_dp*x*x - 1.5_dp))/sqrt(1.0_dp - x*x)
1377 CASE (0)
1378 dplm = 7.5_dp*x*x - 1.5_dp
1379 END SELECT
1380 CASE (4)
1381 SELECT CASE (abs(m))
1382 CASE DEFAULT
1383 cpabort("l = 4 and m value out of bounds")
1384 CASE (4)
1385 dplm = -420*x*(1 - x*x)
1386 CASE (3)
1387 dplm = 105.0_dp*((1.0_dp - x*x)**1.5_dp - 3.0_dp*x*x*(1.0_dp - x*x)**0.5_dp)
1388 CASE (2)
1389 dplm = 105.0_dp*x*(1.0_dp - x*x) - 2.0_dp*x*(52.5_dp*x*x - 7.5_dp)
1390 CASE (1)
1391 IF (abs(x) - 1.0_dp < epsilon(1.0_dp)) THEN
1392 dplm = 0.0_dp
1393 ELSE
1394 dplm = (17.5_dp*3.0_dp*x**2 - 7.5_dp)*sqrt(1.0_dp - x*x) - &
1395 x*(17.5_dp*x**3 - 7.5_dp*x)/sqrt(1.0_dp - x*x)
1396 END IF
1397 CASE (0)
1398 dplm = 4.375_dp*4.0_dp*x**3 - 3.75_dp*2.0_dp*x
1399 END SELECT
1400 CASE (5)
1401 SELECT CASE (abs(m))
1402 CASE DEFAULT
1403 cpabort("l = 5 and m value out of bounds")
1404 CASE (5)
1405 dplm = -945.0_dp*5.0_dp*x*(1.0_dp - x*x)**1.5_dp
1406 CASE (4)
1407 dplm = 945.0_dp*((1.0_dp - x*x)**2 - 2.0_dp*x*x*(1.0_dp - x*x))
1408 CASE (3)
1409 dplm = 945.0_dp*x*(1.0_dp - x*x)**1.5_dp - &
1410 3.0_dp*x*(472.5_dp*x*x - 52.5_dp)*(1.0_dp - x*x)**0.5_dp
1411 CASE (2)
1412 dplm = (3.0_dp*157.5_dp*x**2 - 52.5_dp)*(1.0_dp - x*x) - &
1413 (157.5_dp*x**3 - 52.5_dp*x)*(-2.0_dp*x)
1414 CASE (1)
1415 IF (abs(x) - 1.0_dp < epsilon(1.0_dp)) THEN
1416 dplm = 0.0_dp
1417 ELSE
1418 dplm = -(-39.375_dp*4.0_dp*x*x*x + 2.0_dp*26.25_dp*x)*sqrt(1.0_dp - x*x) + &
1419 x*(-39.375_dp*x**4 + 26.25_dp*x*x - 1.875_dp)/sqrt(1.0_dp - x*x)
1420 END IF
1421 CASE (0)
1422 dplm = 5.0_dp*7.875_dp*x**4 - 3.0_dp*8.75_dp*x**2 + 1.875_dp
1423 END SELECT
1424 CASE (6)
1425 SELECT CASE (abs(m))
1426 CASE DEFAULT
1427 cpabort("l = 6 and m value out of bounds")
1428 CASE (6)
1429 dplm = -10395.0_dp*6.0_dp*x*(1.0_dp - x*x)**2
1430 CASE (5)
1431 dplm = 10395.0_dp*((1.0_dp - x*x)**2.5_dp - 5.0_dp*x*x*(1.0_dp - x*x)**1.5_dp)
1432 CASE (4)
1433 dplm = 2.0_dp*5197.5_dp*x*(1.0_dp - x*x)**2 - &
1434 4.0_dp*x*(5197.5_dp*x*x - 472.5_dp)*(1.0_dp - x*x)
1435 CASE (3)
1436 dplm = -(-3.0_dp*1732.5_dp*x*x + 472.5_dp)*(1.0_dp - x*x)**1.5_dp + &
1437 (-1732.5_dp*x**3 + 472.5_dp*x)*3.0_dp*x*(1.0_dp - x*x)**0.5_dp
1438 CASE (2)
1439 dplm = (433.125_dp*4.0_dp*x**3 - 2.0_dp*236.25_dp*x)*(1.0_dp - x*x) - &
1440 2.0_dp*x*(433.125_dp*x**4 - 236.25_dp*x**2 + 13.125_dp)
1441 CASE (1)
1442 IF (abs(x) - 1.0_dp < epsilon(1.0_dp)) THEN
1443 dplm = 0.0_dp
1444 ELSE
1445 dplm = -(-5.0_dp*86.625_dp*x**4 + 3.0_dp*78.75_dp**2 - 13.125_dp)*sqrt(1.0_dp - x*x) + &
1446 x*(-86.625_dp*x**5 + 78.75_dp*x**3 - 13.125_dp*x)/sqrt(1.0_dp - x*x)
1447 END IF
1448 CASE (0)
1449 dplm = 14.4375_dp*6.0_dp*x**5 - 19.6875_dp*4.0_dp*x**3 + &
1450 6.5625_dp*2.0_dp*x
1451 END SELECT
1452 CASE DEFAULT
1453 mm = abs(m)
1454 IF (mm > l) cpabort("m out of bounds")
1455
1456 !From Wikipedia: dPlm(x) = -(l+1)x/(x^2-1)*Plm(x) + (l-m+1)/(x^2-1)Pl+1m(x)
1457 IF (abs(x) - 1.0_dp < epsilon(1.0_dp)) THEN
1458 dplm = 0.0_dp
1459 ELSE
1460 dplm = -real(l + 1, dp)*x/(x**2 - 1.0_dp)*legendre(x, l, m) &
1461 + real(l - m + 1, dp)/(x**2 - 1.0_dp)*legendre(x, l + 1, m)
1462 END IF
1463 END SELECT
1464
1465 END FUNCTION dlegendre
1466
1467! **************************************************************************************************
1468!> \brief ...
1469!> \param x ...
1470!> \param l ...
1471!> \return ...
1472! **************************************************************************************************
1473 FUNCTION dpof1(x, l)
1474
1475 REAL(kind=dp), INTENT(IN) :: x
1476 INTEGER, INTENT(IN) :: l
1477 REAL(kind=dp) :: dpof1
1478
1479 IF (abs(x) - 1.0_dp > epsilon(1.0_dp)) THEN
1480 cpabort("|x| is not 1")
1481 END IF
1482 IF (x > 0.0_dp) THEN
1483 SELECT CASE (l)
1484 CASE (:-1)
1485 cpabort("Negative l value")
1486 CASE (0)
1487 dpof1 = 0.0_dp
1488 CASE (1)
1489 dpof1 = 1.0_dp
1490 CASE (2)
1491 dpof1 = 3.0_dp
1492 CASE (3)
1493 dpof1 = 6.0_dp
1494 CASE (4)
1495 dpof1 = 10.0_dp
1496 CASE (5)
1497 dpof1 = 15.0_dp
1498 CASE (6)
1499 dpof1 = 21.0_dp
1500 CASE (7:)
1501 cpabort("Not implemented")
1502 END SELECT
1503 ELSE
1504 SELECT CASE (l)
1505 CASE (:-1)
1506 cpabort("Negative l value")
1507 CASE (0)
1508 dpof1 = 0.0_dp
1509 CASE (1)
1510 dpof1 = 1.0_dp
1511 CASE (2)
1512 dpof1 = -3.0_dp
1513 CASE (3)
1514 dpof1 = 6.0_dp
1515 CASE (4)
1516 dpof1 = -10.0_dp
1517 CASE (5)
1518 dpof1 = 15.0_dp
1519 CASE (6)
1520 dpof1 = -21.0_dp
1521 CASE (7:)
1522 cpabort("Not implemented")
1523 END SELECT
1524 END IF
1525
1526 END FUNCTION dpof1
1527
1528! **************************************************************************************************
1529!> \brief ...
1530!> \param n ...
1531!> \param k ...
1532!> \return ...
1533! **************************************************************************************************
1534 FUNCTION choose(n, k)
1535
1536 INTEGER, INTENT(IN) :: n, k
1537 REAL(kind=dp) :: choose
1538
1539 IF (n >= k) THEN
1540 choose = real(nint(fac(n)/(fac(k)*fac(n - k))), kind=dp)
1541 ELSE
1542 choose = 0.0_dp
1543 END IF
1544
1545 END FUNCTION choose
1546
1547! **************************************************************************************************
1548!> \brief ...
1549!> \param n ...
1550!> \param c ...
1551!> \param s ...
1552!> \return ...
1553! **************************************************************************************************
1554 FUNCTION cosn(n, c, s)
1555
1556 INTEGER, INTENT(IN) :: n
1557 REAL(kind=dp), INTENT(IN) :: c, s
1558 REAL(kind=dp) :: cosn
1559
1560 INTEGER :: i, j
1561
1562 cosn = 0.0_dp
1563 IF (abs(c) < epsilon(1.0_dp) .OR. n == 0) THEN
1564 IF (mod(n, 2) == 0) THEN
1565 cosn = (-1.0_dp)**int(n/2)
1566 ELSE
1567 cosn = 0.0_dp
1568 END IF
1569 ELSE
1570 DO i = n, 0, -2
1571 IF (i == 0) THEN
1572 j = n
1573 IF (j == 0) THEN
1574 cosn = cosn + choose(n, j)
1575 ELSE IF (mod(j/2, 2) == 0) THEN
1576 cosn = cosn + choose(n, j)*s**j
1577 ELSE
1578 cosn = cosn - choose(n, j)*s**j
1579 END IF
1580 ELSE
1581 j = n - i
1582 IF (j == 0) THEN
1583 cosn = cosn + choose(n, j)*c**i
1584 ELSE IF (mod(j/2, 2) == 0) THEN
1585 cosn = cosn + choose(n, j)*c**i*s**j
1586 ELSE
1587 cosn = cosn - choose(n, j)*c**i*s**j
1588 END IF
1589 END IF
1590 END DO
1591 END IF
1592
1593 END FUNCTION cosn
1594
1595! **************************************************************************************************
1596!> \brief ...
1597!> \param n ...
1598!> \param c ...
1599!> \param s ...
1600!> \return ...
1601! **************************************************************************************************
1602 FUNCTION dcosn_dcp(n, c, s)
1603
1604 INTEGER, INTENT(IN) :: n
1605 REAL(kind=dp), INTENT(IN) :: c, s
1606 REAL(kind=dp) :: dcosn_dcp
1607
1608 INTEGER :: i, j
1609
1610 dcosn_dcp = 0.0_dp
1611
1612 IF (s < epsilon(1.0_dp)) THEN
1613 dcosn_dcp = 0.0_dp
1614 ELSE
1615 DO i = n, 0, -2
1616 IF (i == 0) THEN
1617 dcosn_dcp = dcosn_dcp
1618 ELSE
1619 j = n - i
1620 IF (mod(j/2, 2) == 0) THEN
1621 dcosn_dcp = dcosn_dcp + choose(n, j)*real(i, dp)*c**(i - 1)*s**j
1622 ELSE
1623 dcosn_dcp = dcosn_dcp - choose(n, j)*real(i, dp)*c**(i - 1)*s**j
1624 END IF
1625 END IF
1626 END DO
1627 END IF
1628
1629 END FUNCTION dcosn_dcp
1630
1631! **************************************************************************************************
1632!> \brief ...
1633!> \param n ...
1634!> \param c ...
1635!> \param s ...
1636!> \return ...
1637! **************************************************************************************************
1638 FUNCTION dcosn_dsp(n, c, s)
1639
1640 INTEGER, INTENT(IN) :: n
1641 REAL(kind=dp), INTENT(IN) :: c, s
1642 REAL(kind=dp) :: dcosn_dsp
1643
1644 INTEGER :: i, j
1645
1646 dcosn_dsp = 0.0_dp
1647 IF (c < epsilon(1.0_dp) .OR. s < epsilon(1.0_dp)) THEN
1648 dcosn_dsp = 0.0_dp
1649 ELSE
1650 DO i = n, 0, -2
1651 j = n - i
1652 IF (j == 0) THEN
1653 dcosn_dsp = dcosn_dsp
1654 ELSE
1655 IF (mod(j/2, 2) == 0) THEN
1656 dcosn_dsp = dcosn_dsp + choose(n, j)*real(j, dp)*s**(j - 1)*c**i
1657 ELSE
1658 dcosn_dsp = dcosn_dsp - choose(n, j)*real(j, dp)*s**(j - 1)*c**i
1659 END IF
1660 END IF
1661 END DO
1662 END IF
1663
1664 END FUNCTION dcosn_dsp
1665
1666! **************************************************************************************************
1667!> \brief ...
1668!> \param n ...
1669!> \param c ...
1670!> \param s ...
1671!> \return ...
1672! **************************************************************************************************
1673 FUNCTION sinn(n, c, s)
1674
1675 INTEGER, INTENT(IN) :: n
1676 REAL(kind=dp), INTENT(IN) :: c, s
1677 REAL(kind=dp) :: sinn
1678
1679 INTEGER :: i, j
1680
1681 sinn = 0.0_dp
1682
1683 IF (abs(s) < epsilon(1.0_dp) .OR. n == 0) THEN
1684 sinn = 0.0_dp
1685 ELSE IF (abs(c) < epsilon(1.0_dp)) THEN
1686 IF (mod(n, 2) == 0) THEN
1687 sinn = 0.0_dp
1688 ELSE
1689 sinn = s*(-1.0_dp)**int((n - 1)/2)
1690 END IF
1691 ELSE
1692 DO i = n - 1, 0, -2
1693 IF (i == 0) THEN
1694 j = n
1695 IF (mod(j/2, 2) == 0) THEN
1696 sinn = sinn + choose(n, j)*s**j
1697 ELSE
1698 sinn = sinn - choose(n, j)*s**j
1699 END IF
1700 ELSE
1701 j = n - i
1702 IF (mod(j/2, 2) == 0) THEN
1703 sinn = sinn + choose(n, j)*c**i*s**j
1704 ELSE
1705 sinn = sinn - choose(n, j)*c**i*s**j
1706 END IF
1707 END IF
1708 END DO
1709 END IF
1710
1711 END FUNCTION sinn
1712
1713! **************************************************************************************************
1714!> \brief ...
1715!> \param n ...
1716!> \param c ...
1717!> \param s ...
1718!> \return ...
1719! **************************************************************************************************
1720 FUNCTION dsinn_dcp(n, c, s)
1721
1722 INTEGER, INTENT(IN) :: n
1723 REAL(kind=dp), INTENT(IN) :: c, s
1724 REAL(kind=dp) :: dsinn_dcp
1725
1726 INTEGER :: i, j
1727
1728 dsinn_dcp = 0.0_dp
1729
1730 IF (c < epsilon(1.0_dp) .OR. s < epsilon(1.0_dp)) THEN
1731 dsinn_dcp = 0.0_dp
1732 ELSE
1733 DO i = n - 1, 0, -2
1734 IF (i == 0) THEN
1735 dsinn_dcp = dsinn_dcp
1736 ELSE
1737 j = n - i
1738 IF (mod(j/2, 2) == 0) THEN
1739 dsinn_dcp = dsinn_dcp + choose(n, j)*real(i, dp)*c**(i - 1)*s**j
1740 ELSE
1741 dsinn_dcp = dsinn_dcp - choose(n, j)*real(i, dp)*c**(i - 1)*s**j
1742 END IF
1743 END IF
1744 END DO
1745 END IF
1746
1747 END FUNCTION dsinn_dcp
1748
1749! **************************************************************************************************
1750!> \brief ...
1751!> \param n ...
1752!> \param c ...
1753!> \param s ...
1754!> \return ...
1755! **************************************************************************************************
1756 FUNCTION dsinn_dsp(n, c, s)
1757
1758 INTEGER, INTENT(IN) :: n
1759 REAL(kind=dp), INTENT(IN) :: c, s
1760 REAL(kind=dp) :: dsinn_dsp
1761
1762 INTEGER :: i, j
1763
1764 dsinn_dsp = 0.0_dp
1765
1766 IF (c < epsilon(1.0_dp) .OR. s < epsilon(1.0_dp)) THEN
1767 dsinn_dsp = 0.0_dp
1768 ELSE
1769 DO i = n - 1, 0, -2
1770 j = n - i
1771 IF (j == 0) THEN
1772 dsinn_dsp = dsinn_dsp
1773 ELSE
1774 IF (mod(j/2, 2) == 0) THEN
1775 dsinn_dsp = dsinn_dsp + choose(n, j)*c**i*real(j, dp)*s**(j - 1)
1776 ELSE
1777 dsinn_dsp = dsinn_dsp - choose(n, j)*c**i*real(j, dp)*s**(j - 1)
1778 END IF
1779 END IF
1780 END DO
1781 END IF
1782
1783 END FUNCTION dsinn_dsp
1784
1785! **************************************************************************************************
1786!> \brief Compute the Clebsch-Gordon coefficient C = < j1 m1 j2 m2 | J M > = < j1 j2; m1 m2 | J M >
1787!> \param j1 Angular momentum quantum number of the first state | j1 m1 >
1788!> \param m1 Magnetic quantum number of the first first state | j1 m1 >
1789!> \param j2 Angular momentum quantum number of the second state | j2 m2 >
1790!> \param m2 Magnetic quantum number of the second state | j2 m2 >
1791!> \param J Angular momentum quantum number of the coupled state | J M >
1792!> \param M Magnetic quantum number of the coupled state | J M >
1793!> \param CG_coeff Clebsch-Gordon coefficient C^{JM}_{j1 m1 j2 m2}
1794!> \author Matthias Krack (16.09.2022, based on a program by D. G. Simpson)
1795!> \note Generic routine allowing also for fractional arguments. It should return CG coefficients
1796!> consistent with the standard definition and normalization, e.g. of Wolfram Mathematica.
1797!> The input parameters have to be integer or half-integer.
1798! **************************************************************************************************
1799 SUBROUTINE clebsch_gordon_coefficient(j1, m1, j2, m2, J, M, CG_coeff)
1800
1801 REAL(kind=dp), INTENT(IN) :: j1, m1, j2, m2, j, m
1802 REAL(kind=dp), INTENT(OUT) :: cg_coeff
1803
1804 INTEGER :: k, kmax
1805 REAL(kind=dp) :: sumk, t
1806
1807 ! Check validity of the input parameters
1808 IF (j1 < 0.0_dp) &
1809 cpabort("The angular momentum quantum number j1 has to be nonnegative")
1810 IF (.NOT. (is_integer(j1) .OR. is_integer(2.0_dp*j1))) &
1811 cpabort("The angular momentum quantum number j1 has to be integer or half-integer")
1812 IF (j2 < 0.0_dp) &
1813 cpabort("The angular momentum quantum number j2 has to be nonnegative")
1814 IF (.NOT. (is_integer(j2) .OR. is_integer(2.0_dp*j2))) &
1815 cpabort("The angular momentum quantum number j2 has to be integer or half-integer")
1816 IF (j < 0.0_dp) &
1817 cpabort("The angular momentum quantum number J has to be nonnegative")
1818 IF (.NOT. (is_integer(j) .OR. is_integer(2.0_dp*j))) &
1819 cpabort("The angular momentum quantum number J has to be integer or half-integer")
1820 IF ((abs(m1) - j1) > epsilon(m1)) &
1821 cpabort("The angular momentum quantum number m1 has to satisfy -j1 <= m1 <= j1")
1822 IF (.NOT. (is_integer(m1) .OR. is_integer(2.0_dp*m1))) &
1823 cpabort("The angular momentum quantum number m1 has to be integer or half-integer")
1824 IF ((abs(m2) - j2) > epsilon(m2)) &
1825 cpabort("The angular momentum quantum number m2 has to satisfy -j2 <= m1 <= j2")
1826 IF (.NOT. (is_integer(m2) .OR. is_integer(2.0_dp*m2))) &
1827 cpabort("The angular momentum quantum number m2 has to be integer or half-integer")
1828 IF ((abs(m) - j) > epsilon(m)) &
1829 cpabort("The angular momentum quantum number M has to satisfy -J <= M <= J")
1830 IF (.NOT. (is_integer(m) .OR. is_integer(2.0_dp*m))) &
1831 cpabort("The angular momentum quantum number M has to be integer or half-integer")
1832
1833 IF (is_integer(j1 + j2 + j) .AND. &
1834 is_integer(j1 + m1) .AND. &
1835 is_integer(j2 + m2) .AND. &
1836 is_integer(j + m) .AND. &
1837 is_integer(j - j1 - m2) .AND. &
1838 is_integer(j - j2 + m1)) THEN
1839 IF ((j < abs(j1 - j2)) .OR. &
1840 (j > j1 + j2) .OR. &
1841 (abs(m1) > j1) .OR. &
1842 (abs(m2) > j2) .OR. &
1843 (abs(m) > j)) THEN
1844 cg_coeff = 0.0_dp
1845 ELSE
1846 ! Compute Clebsch-Gordan coefficient
1847 sumk = 0.0_dp
1848 kmax = nint(max(j1 + j2 - j, j1 - j + m2, j2 - j - m1, j1 - m1, j2 + m2))
1849 DO k = 0, kmax
1850 IF (j1 + j2 - j - k < 0.0_dp) cycle
1851 IF (j - j1 - m2 + k < 0.0_dp) cycle
1852 IF (j - j2 + m1 + k < 0.0_dp) cycle
1853 IF (j1 - m1 - k < 0.0_dp) cycle
1854 IF (j2 + m2 - k < 0.0_dp) cycle
1855 t = sfac(nint(j1 + j2 - j - k))*sfac(nint(j - j1 - m2 + k))* &
1856 sfac(nint(j - j2 + m1 + k))*sfac(nint(j1 - m1 - k))* &
1857 sfac(nint(j2 + m2 - k))*sfac(k)
1858 IF (mod(k, 2) == 1) t = -t
1859 sumk = sumk + 1.0_dp/t
1860 END DO
1861 cg_coeff = sqrt((2.0_dp*j + 1)/sfac(nint(j1 + j2 + j + 1.0_dp)))* &
1862 sqrt(sfac(nint(j1 + j2 - j))*sfac(nint(j2 + j - j1))*sfac(nint(j + j1 - j2)))* &
1863 sqrt(sfac(nint(j1 + m1))*sfac(nint(j1 - m1))* &
1864 sfac(nint(j2 + m2))*sfac(nint(j2 - m2))* &
1865 sfac(nint(j + m))*sfac(nint(j - m)))*sumk
1866 END IF
1867 ELSE
1868 cpabort("Invalid set of input parameters < j1 m1 j2 m2 | J M > specified")
1869 END IF
1870
1871 END SUBROUTINE clebsch_gordon_coefficient
1872
1873! **************************************************************************************************
1874!> \brief Compute the Wigner 3-j symbol
1875!> / j1 j2 j3 \
1876!> \ m1 m2 m3 /
1877!> using the Clebsch-Gordon coefficients
1878!> \param j1 Angular momentum quantum number of the first state | j1 m1 >
1879!> \param m1 Magnetic quantum number of the first first state | j1 m1 >
1880!> \param j2 Angular momentum quantum number of the second state | j2 m2 >
1881!> \param m2 Magnetic quantum number of the second state | j2 m2 >
1882!> \param j3 Angular momentum quantum number of the third state | j3 m3 >
1883!> \param m3 Magnetic quantum number of the third state | j3 m3 >
1884!> \param W_3j Wigner 3-j symbol
1885!> \author Matthias Krack (16.09.2022, MK)
1886! **************************************************************************************************
1887 SUBROUTINE wigner_3j(j1, m1, j2, m2, j3, m3, W_3j)
1888
1889 REAL(kind=dp), INTENT(IN) :: j1, m1, j2, m2, j3, m3
1890 REAL(kind=dp), INTENT(OUT) :: w_3j
1891
1892 REAL(kind=dp) :: cg_coeff
1893
1894 IF ((abs(m1 + m2 + m3) < epsilon(m1)) .AND. &
1895 (abs(j1 - j2) <= j3) .AND. (j3 <= abs(j1 + j2)) .AND. &
1896 is_integer(j1 + j2 + j3)) THEN
1897 IF (is_integer(j1 - j2 - m3)) THEN
1898 CALL clebsch_gordon_coefficient(j1, m1, j2, m2, j3, -m3, cg_coeff)
1899 w_3j = (-1.0_dp)**int(j1 - j2 - m3)*cg_coeff/sqrt(2.0_dp*j3 + 1.0_dp)
1900 ELSE
1901 cpabort("j1 - j2 - m3 results in an invalid non-integer exponent")
1902 END IF
1903 ELSE
1904 w_3j = 0.0_dp
1905 END IF
1906
1907 END SUBROUTINE wigner_3j
1908
1909! **************************************************************************************************
1910!> \brief Check if the input value is an integer number within double-precision accuracy
1911!> \param x Input value to be checked
1912!> \return True if the input value is integer otherwise false
1913!> \author Matthias Krack (16.09.2022, MK)
1914! **************************************************************************************************
1915 FUNCTION is_integer(x)
1916
1917 REAL(kind=dp), INTENT(IN) :: x
1918 LOGICAL :: is_integer
1919
1920 IF ((abs(x) - aint(abs(x), kind=dp)) > epsilon(x)) THEN
1921 is_integer = .false.
1922 ELSE
1923 is_integer = .true.
1924 END IF
1925
1926 END FUNCTION is_integer
1927
1928END MODULE spherical_harmonics
spherical_harmonics::clebsch_gordon
Definition spherical_harmonics.F:53
spherical_harmonics::dy_lm
Definition spherical_harmonics.F:49
spherical_harmonics::y_lm
Definition spherical_harmonics.F:45
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::sp
integer, parameter, public sp
Definition kinds.F:33
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::maxfac
integer, parameter, public maxfac
Definition mathconstants.F:36
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
spherical_harmonics
Calculate spherical harmonics.
Definition spherical_harmonics.F:27
spherical_harmonics::legendre
real(kind=dp) function, public legendre(x, l, m)
...
Definition spherical_harmonics.F:1193
spherical_harmonics::clebsch_gordon_coefficient
subroutine, public clebsch_gordon_coefficient(j1, m1, j2, m2, j, m, cg_coeff)
Compute the Clebsch-Gordon coefficient C = < j1 m1 j2 m2 | J M > = < j1 j2; m1 m2 | J M >
Definition spherical_harmonics.F:1800
spherical_harmonics::wigner_3j
subroutine, public wigner_3j(j1, m1, j2, m2, j3, m3, w_3j)
Compute the Wigner 3-j symbol / j1 j2 j3 \ \ m1 m2 m3 / using the Clebsch-Gordon coefficients.
Definition spherical_harmonics.F:1888
spherical_harmonics::clebsch_gordon_init
subroutine, public clebsch_gordon_init(l)
...
Definition spherical_harmonics.F:238
spherical_harmonics::clebsch_gordon_deallocate
subroutine, public clebsch_gordon_deallocate()
...
Definition spherical_harmonics.F:220
spherical_harmonics::dlegendre
real(kind=dp) function, public dlegendre(x, l, m)
...
Definition spherical_harmonics.F:1337