Line data Source code
1 : MODULE m_ylm
2 : use iso_c_binding
3 : IMPLICIT NONE
4 : PRIVATE
5 : PUBLIC ylm4, ylm4_batched
6 : !************************************************************
7 : ! generate the spherical harmonics for the vector v
8 : ! using a stable upward recursion in l. (see notes
9 : ! by m. weinert.)
10 : ! m.weinert january 1982
11 : ! modified by R. Podloucky (added in ynorm); July 1989
12 : ! cleaned up mw 1995
13 : !
14 : ! modified to make use of f90 constructs. note that
15 : ! the normalization is an internal subroutine and hence
16 : ! can only be called from here. also, no need to dimension
17 : ! arrays for ynorm, done dynamically. mw 1999
18 : !
19 : ! GPU version added
20 : ! U.Alekseeva Oktober 2018
21 : !************************************************************
22 : CONTAINS
23 21359881 : SUBROUTINE ylm4(lmax, v, ylm)
24 : ! USE m_juDFT
25 : INTEGER, VALUE, INTENT(IN) :: lmax
26 : REAL, INTENT(IN), target :: v(3)
27 : COMPLEX, INTENT(OUT), target :: ylm((lmax + 1)**2)
28 :
29 : real, pointer :: v_2d(:,:)
30 21359881 : complex, pointer :: ylm_2d(:,:)
31 :
32 64079643 : call c_f_pointer(c_loc(v), v_2d, [3,1])
33 64079643 : call c_f_pointer(c_loc(ylm), ylm_2d, [(lmax + 1)**2, 1])
34 21359881 : call ylm4_batched(lmax, v_2d, ylm_2d)
35 21359881 : END SUBROUTINE ylm4
36 :
37 21409490 : subroutine ylm4_batched(lmax, v, ylm)
38 : implicit none
39 : INTEGER, VALUE, INTENT(IN) :: lmax
40 : REAL, INTENT(IN) :: v(:,:)
41 : COMPLEX, INTENT(OUT):: ylm(:,:)
42 :
43 : INTEGER :: l, lm0, m, i
44 : REAL :: fac, x, y, z, xy, r, rxy, cth, sth, cph, sph, cph2
45 21409490 : REAL :: c(0:max(lmax,1)), s(0:max(lmax,1))
46 21409490 : REAL :: p(0:lmax, 0:lmax)
47 : COMPLEX :: ylms
48 21409490 : REAL :: ynorm_dev((lmax + 1)**2)
49 : REAL :: a, cd
50 : real, parameter :: fpi = 4.0*3.1415926535897932, small = 1.0e-12
51 :
52 : !---> calculate norm
53 168436220 : DO l = 0, lmax
54 147026730 : lm0 = l*(l + 1) + 1
55 147026730 : cd = 1.0
56 147026730 : a = sqrt((2*l + 1)/fpi)
57 147026730 : ynorm_dev(lm0) = a
58 756959171 : DO m = 1, l
59 588522951 : cd = cd/((l + 1 - m)*(l + m))
60 588522951 : ynorm_dev(lm0 + m) = a*sqrt(cd)
61 735549681 : ynorm_dev(lm0 - m) = ((-1.0)**m)*ynorm_dev(lm0 + m)
62 : ENDDO
63 : ENDDO
64 :
65 :
66 : !---> calculate sin and cos of theta and phi
67 : !$OMP parallel do default(none) &
68 : !$OMP private(i, x, y, z, xy, r, rxy, cth, sth, fac, m, l, p, c ,s, ylms, lm0, cph, sph, cph2)&
69 21409490 : !$OMP shared(ylm, ynorm_dev, lmax, v)
70 : do i = 1,size(v,2)
71 : x = v(1,i)
72 : y = v(2,i)
73 : z = v(3,i)
74 : xy = x*x + y*y
75 : r = sqrt(xy + z*z)
76 : rxy = sqrt(xy)
77 :
78 : IF (r .GT. small) THEN
79 : cth = z/r
80 : sth = rxy/r
81 : ELSE
82 : sth = 0.0
83 : cth = 1.0
84 : ENDIF
85 : IF (rxy .GT. small) THEN
86 : cph = x/rxy
87 : sph = y/rxy
88 : ELSE
89 : cph = 1.0
90 : sph = 0.0
91 : ENDIF
92 :
93 : !---> generate associated legendre functions for m.ge.0
94 : fac = 1.0
95 : !---> loop over m values
96 : DO m = 0, lmax - 1
97 : fac = -(m + m - 1)*fac
98 : p(m, m) = fac
99 : p(m + 1, m) = (m + m + 1)*cth*fac
100 : !---> recurse upward in l
101 : DO l = m + 2, lmax
102 : p(l, m) = ((l + l - 1)*cth*p(l - 1, m) - (l + m - 1)*p(l - 2, m))/(l - m)
103 : ENDDO
104 : fac = fac*sth
105 : ENDDO
106 : p(lmax, lmax) = -(lmax + lmax - 1)*fac
107 :
108 : !---> determine sin and cos of phi
109 : s(0) = 0.0
110 : s(1) = sph
111 : c(0) = 1.0
112 : c(1) = cph
113 : cph2 = cph + cph
114 : DO m = 2, lmax
115 : s(m) = cph2*s(m - 1) - s(m - 2)
116 : c(m) = cph2*c(m - 1) - c(m - 2)
117 : ENDDO
118 :
119 : !---> multiply in the normalization factors
120 : DO l = 0, lmax
121 : ylm(l*(l + 1) + 1, i) = ynorm_dev(l*(l + 1) + 1)*cmplx(p(l, 0), 0.0)
122 : ENDDO
123 : DO m = 1, lmax
124 : DO l = m, lmax
125 : lm0 = l*(l + 1) + 1
126 : ylms = p(l, m)*cmplx(c(m), s(m))
127 : ylm(lm0 + m, i) = ynorm_dev(lm0 + m)*ylms
128 : ylm(lm0 - m, i) = conjg(ylms)*ynorm_dev(lm0 - m)
129 : ENDDO
130 : ENDDO
131 : enddo
132 : !$OMP end parallel do
133 21409490 : end subroutine ylm4_batched
134 : END MODULE m_ylm
|
