Line data Source code
1 : MODULE m_rotMMPmat
2 :
3 : USE m_constants
4 : USE m_types_sym
5 :
6 : IMPLICIT NONE
7 :
8 : PRIVATE
9 : PUBLIC :: rotMMPmat
10 :
11 : INTERFACE rotMMPmat
12 : PROCEDURE :: rotMMPmat_dwgn, rotMMPmat_angle, rotMMPmat_sym_op
13 : PROCEDURE :: rotMMPmat_angle_one_spin, rotMMPmat_dwgn_one_spin, rotMMPmat_sym_op_one_spin
14 : END INTERFACE
15 :
16 : CONTAINS
17 :
18 5051502 : PURE FUNCTION rotMMPmat_dwgn(mmpmat,dwgn,dwgnp,su, spin_rotation, real_space_rotation, inverse) Result(mmpmatOut)
19 :
20 : COMPLEX, INTENT(IN) :: mmpmat(-lmaxU_const:,-lmaxU_const:,:)
21 : COMPLEX, OPTIONAL, INTENT(IN) :: dwgn(-lmaxU_const:,-lmaxU_const:)
22 : COMPLEX, OPTIONAL, INTENT(IN) :: dwgnp(-lmaxU_const:,-lmaxU_const:)
23 : COMPLEX, OPTIONAL, INTENT(IN) :: su(:,:)
24 : LOGICAL, OPTIONAL, INTENT(IN) :: spin_rotation
25 : LOGICAL, OPTIONAL, INTENT(IN) :: real_space_rotation
26 : LOGICAL, OPTIONAL, INTENT(IN) :: inverse
27 :
28 : COMPLEX, ALLOCATABLE :: mmpmatOut(:,:,:)
29 :
30 : COMPLEX :: d(2,2)
31 : INTEGER :: ispin,m,mp
32 : LOGICAL :: spin_rotation_arg, real_space_rotation_arg, inverseArg
33 :
34 5051502 : spin_rotation_arg = .FALSE.
35 5051502 : IF(PRESENT(spin_rotation)) spin_rotation_arg = spin_rotation
36 :
37 5051502 : real_space_rotation_arg = .TRUE.
38 5051502 : IF(PRESENT(real_space_rotation)) real_space_rotation_arg = real_space_rotation
39 :
40 5051502 : inverseArg = .FALSE.
41 5051502 : IF(PRESENT(inverse)) inverseArg = inverse
42 :
43 25257510 : IF(.NOT.ALLOCATED(mmpmatOut)) ALLOCATE(mmpmatOut,mold=mmpmat)
44 303539802 : mmpmatOut = mmpmat
45 :
46 :
47 5051502 : IF(spin_rotation_arg.AND.PRESENT(su)) THEN
48 0 : DO m = -lmaxU_const, lmaxU_const
49 0 : DO mp = -lmaxU_const, lmaxU_const
50 :
51 0 : d = cmplx_0
52 0 : d(1,1) = mmpmatOut(m,mp,1)
53 0 : d(2,2) = mmpmatOut(m,mp,2)
54 0 : IF(SIZE(mmpmat,3)>=3) THEN
55 0 : d(2,1) = mmpmatOut(m,mp,3)
56 0 : IF(SIZE(mmpmat,3)==3) THEN
57 0 : d(1,2) = conjg(mmpmatOut(mp,m,3))
58 : ELSE
59 0 : d(1,2) = mmpmatOut(m,mp,4)
60 : ENDIF
61 : ENDIF
62 :
63 0 : IF(inverseArg) THEN
64 0 : d = matmul(conjg(transpose(su)),d)
65 0 : d = matmul(d,su)
66 : ELSE
67 0 : d = matmul(su,d)
68 0 : d = matmul(d,conjg(transpose(su)))
69 : ENDIF
70 :
71 0 : mmpmatOut(m,mp,1) = d(1,1)
72 0 : mmpmatOut(m,mp,2) = d(2,2)
73 0 : IF(SIZE(mmpmat,3)>=3) THEN
74 0 : mmpmatOut(m,mp,3) = d(2,1)
75 0 : IF(SIZE(mmpmat,3)==4) THEN
76 0 : mmpmatOut(mp,m,4) = d(1,2)
77 : ENDIF
78 : ENDIF
79 :
80 : ENDDO
81 : ENDDO
82 : ENDIF
83 :
84 5051502 : IF(real_space_rotation_arg.AND.PRESENT(dwgn)) THEN
85 10199516 : DO ispin = 1, SIZE(mmpmat,3)
86 10199516 : IF(inverseArg) THEN
87 0 : mmpmatOut(:,:,ispin) = matmul(dwgn,mmpmatOut(:,:,ispin))
88 0 : IF(PRESENT(dwgnp)) THEN
89 0 : mmpmatOut(:,:,ispin) = matmul(mmpmatOut(:,:,ispin),conjg(transpose(dwgnp)))
90 : ELSE
91 0 : mmpmatOut(:,:,ispin) = matmul(mmpmatOut(:,:,ispin),conjg(transpose(dwgn)))
92 : ENDIF
93 : ELSE
94 7948533616 : mmpmatOut(:,:,ispin) = matmul(conjg(transpose(dwgn)),mmpmatOut(:,:,ispin))
95 5148014 : IF(PRESENT(dwgnp)) THEN
96 7948496560 : mmpmatOut(:,:,ispin) = matmul(mmpmatOut(:,:,ispin),dwgnp)
97 : ELSE
98 37056 : mmpmatOut(:,:,ispin) = matmul(mmpmatOut(:,:,ispin),dwgn)
99 : ENDIF
100 : ENDIF
101 : ENDDO
102 : ENDIF
103 :
104 5051502 : END FUNCTION rotMMPmat_dwgn
105 :
106 3958 : PURE FUNCTION rotMMPmat_angle(mmpmat,alpha,beta,gamma,l,lp,spin_rotation,real_space_rotation,inverse) Result(mmpmatOut)
107 : use m_types_nococonv
108 : COMPLEX, INTENT(IN) :: mmpmat(-lmaxU_const:,-lmaxU_const:,:)
109 : REAL, INTENT(IN) :: alpha,beta,gamma !Euler angles
110 : INTEGER, INTENT(IN) :: l
111 : INTEGER, OPTIONAL, INTENT(IN) :: lp
112 : LOGICAL, OPTIONAL, INTENT(IN) :: spin_rotation
113 : LOGICAL, OPTIONAL, INTENT(IN) :: real_space_rotation
114 : LOGICAL, OPTIONAL, INTENT(IN) :: inverse
115 :
116 : COMPLEX, ALLOCATABLE :: mmpmatOut(:,:,:)
117 : COMPLEX :: su(2,2), eia
118 : REAL :: co_bh,si_bh, alphaArg, betaArg, gammaArg
119 : COMPLEX :: d(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const)
120 : COMPLEX :: dp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const)
121 : LOGICAL :: spin_rotation_arg, inverseArg
122 :
123 15832 : TYPE(t_nococonv):: nococonv
124 :
125 19790 : IF(.NOT.ALLOCATED(mmpmatOut)) ALLOCATE(mmpmatOut,mold=mmpmat)
126 236714 : mmpmatOut = mmpmat
127 :
128 3958 : IF(ABS(alpha)<1e-10.AND.ABS(beta)<1e-10.AND.ABS(gamma)<1e-10) RETURN
129 :
130 3108 : inverseArg = .FALSE.
131 3108 : IF(PRESENT(inverse)) inverseArg = inverse
132 :
133 3084 : IF(inverseArg) THEN
134 3084 : alphaArg = -gamma; betaArg = -beta; gammaArg = -alpha
135 : ELSE
136 24 : alphaArg = alpha; betaArg = beta; gammaArg = gamma
137 : ENDIF
138 :
139 3108 : d = d_wigner_mat(alphaArg,betaArg,gammaArg,l)
140 3108 : IF(PRESENT(lp)) THEN
141 3084 : dp = d_wigner_mat(alphaArg,betaArg,gammaArg,lp)
142 : ENDIF
143 :
144 3108 : co_bh = cos(betaArg*0.5)
145 3108 : si_bh = sin(betaArg*0.5)
146 3108 : eia = exp( ImagUnit * alphaArg/2.0 )
147 3108 : su(1,1) = eia*co_bh
148 3108 : su(2,1) = -conjg(eia)*si_bh
149 3108 : su(1,2) = eia*si_bh
150 3108 : su(2,2) = conjg(eia)*co_bh
151 3108 : su=nococonv%chi(alphaArg,betaArg)
152 :
153 3108 : IF(PRESENT(lp)) THEN
154 : mmpmatOut = rotMMPmat_dwgn(mmpmat,d,dwgnp=dp,su=su,spin_rotation=spin_rotation,&
155 181956 : real_space_rotation=real_space_rotation)
156 : ELSE
157 : mmpmatOut = rotMMPmat_dwgn(mmpmat,d,su=su,spin_rotation=spin_rotation,&
158 1416 : real_space_rotation=real_space_rotation)
159 : ENDIF
160 :
161 : END FUNCTION rotMMPmat_angle
162 :
163 5363864 : PURE FUNCTION rotMMPmat_sym_op(mmpmat,sym,iop,l,lp,inverse,reciprocal,spin_rotation,real_space_rotation) Result(mmpmatOut)
164 :
165 : COMPLEX, INTENT(IN) :: mmpmat(-lmaxU_const:,-lmaxU_const:,:)
166 : TYPE(t_sym), INTENT(IN) :: sym
167 : INTEGER, INTENT(IN) :: iop
168 : INTEGER, INTENT(IN) :: l
169 : INTEGER, OPTIONAL, INTENT(IN) :: lp
170 : LOGICAL, OPTIONAL, INTENT(IN) :: inverse, reciprocal, spin_rotation
171 : LOGICAL, OPTIONAL, INTENT(IN) :: real_space_rotation
172 :
173 : COMPLEX, ALLOCATABLE :: mmpmatOut(:,:,:)
174 : COMPLEX :: dwgn(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const)
175 : COMPLEX :: dwgnp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const)
176 : INTEGER :: iopArg, lpArg
177 : LOGICAL :: reciprocalArg, inverseArg
178 :
179 26819320 : IF(.NOT.ALLOCATED(mmpmatOut)) ALLOCATE(mmpmatOut,mold=mmpmat)
180 321969160 : mmpmatOut = mmpmat
181 :
182 5363864 : IF(iop==1.or.l==0) RETURN
183 :
184 5048394 : reciprocalArg = .FALSE.
185 5048394 : IF(PRESENT(reciprocal)) reciprocalArg = reciprocal
186 :
187 5048394 : inverseArg = .FALSE.
188 5048394 : IF(PRESENT(inverse)) inverseArg = inverse
189 :
190 5048394 : lpArg = l
191 5048394 : IF(PRESENT(lp)) lpArg = lp
192 :
193 5048394 : iopArg = iop
194 5048394 : IF(reciprocalArg) THEN
195 0 : IF(iop <= sym%nop) THEN
196 0 : iopArg = sym%invtab(iop)
197 : ELSE
198 0 : iopArg = sym%invtab(iop-sym%nop)
199 : ENDIF
200 5048394 : ELSE IF(inverseArg) THEN
201 0 : iopArg = sym%invtab(iop)
202 : ENDIF
203 :
204 287758458 : dwgn = sym%d_wgn(:,:,l,iopArg)
205 287758458 : dwgnp = sym%d_wgn(:,:,lpArg,iopArg)
206 :
207 5048394 : IF(reciprocalArg) THEN
208 0 : IF(iop <= sym%nop) THEN
209 0 : dwgn = transpose(dwgn)
210 0 : dwgnp = transpose(dwgnp)
211 : ELSE
212 0 : dwgn = -transpose(dwgn)
213 0 : dwgnp = -transpose(dwgnp)
214 : ENDIF
215 : ENDIF
216 :
217 : mmpmatOut = rotMMPmat_dwgn(mmpmat,dwgn=dwgn,dwgnp=dwgnp,spin_rotation=spin_rotation,&
218 303356430 : real_space_rotation=real_space_rotation)
219 :
220 : END FUNCTION rotMMPmat_sym_op
221 :
222 0 : PURE FUNCTION rotMMPmat_sym_op_one_spin(mmpmat,sym,iop,l,lp,inverse,reciprocal) Result(mmpmatOut)
223 :
224 : COMPLEX, INTENT(IN) :: mmpmat(-lmaxU_const:,-lmaxU_const:)
225 : TYPE(t_sym), INTENT(IN) :: sym
226 : INTEGER, INTENT(IN) :: iop
227 : INTEGER, INTENT(IN) :: l
228 : INTEGER, OPTIONAL, INTENT(IN) :: lp
229 : LOGICAL, OPTIONAL, INTENT(IN) :: inverse
230 : LOGICAL, OPTIONAL, INTENT(IN) :: reciprocal
231 :
232 0 : COMPLEX, ALLOCATABLE :: mmpmatOut(:,:), mmpmatOutsplit(:,:,:)
233 : COMPLEX :: mmpmatsplit(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,1)
234 :
235 0 : ALLOCATE(mmpmatOut,mold=mmpMat)
236 0 : mmpmatsplit(:,:,1) = mmpmat
237 0 : mmpmatOutsplit = rotMMPmat_sym_op(mmpmatsplit,sym,iop,l,lp=lp,inverse=inverse,reciprocal=reciprocal)
238 0 : mmpmatOut = mmpmatOutsplit(:,:,1)
239 :
240 0 : END FUNCTION rotMMPmat_sym_op_one_spin
241 :
242 3798 : PURE FUNCTION rotMMPmat_angle_one_spin(mmpmat,alpha,beta,gamma,l,lp,inverse) Result(mmpmatOut)
243 :
244 : COMPLEX, INTENT(IN) :: mmpmat(-lmaxU_const:,-lmaxU_const:)
245 : REAL, INTENT(IN) :: alpha,beta,gamma !Euler angles
246 : INTEGER, INTENT(IN) :: l
247 : INTEGER, OPTIONAL, INTENT(IN) :: lp
248 : LOGICAL, OPTIONAL, INTENT(IN) :: inverse
249 :
250 3798 : COMPLEX, ALLOCATABLE :: mmpmatOut(:,:), mmpmatOutsplit(:,:,:)
251 : COMPLEX :: mmpmatsplit(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,1)
252 :
253 15192 : ALLOCATE(mmpmatOut,mold=mmpMat)
254 216486 : mmpmatsplit(:,:,1) = mmpmat
255 220284 : mmpmatOutsplit = rotMMPmat_angle(mmpmatsplit,alpha,beta,gamma,l,lp=lp, inverse=inverse)
256 220284 : mmpmatOut = mmpmatOutsplit(:,:,1)
257 :
258 3798 : END FUNCTION rotMMPmat_angle_one_spin
259 :
260 0 : PURE FUNCTION rotMMPmat_dwgn_one_spin(mmpmat,dwgn,dwgnp,inverse) Result(mmpmatOut)
261 :
262 : COMPLEX, INTENT(IN) :: mmpmat(-lmaxU_const:,-lmaxU_const:)
263 : COMPLEX, INTENT(IN) :: dwgn(-lmaxU_const:,-lmaxU_const:)
264 : COMPLEX, OPTIONAL, INTENT(IN) :: dwgnp(-lmaxU_const:,-lmaxU_const:)
265 : LOGICAL, OPTIONAL, INTENT(IN) :: inverse
266 :
267 0 : COMPLEX, ALLOCATABLE :: mmpmatOut(:,:), mmpmatOutsplit(:,:,:)
268 : COMPLEX :: mmpmatsplit(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,1)
269 :
270 0 : ALLOCATE(mmpmatOut,mold=mmpMat)
271 0 : mmpmatsplit(:,:,1) = mmpmat
272 0 : mmpmatOutsplit = rotMMPmat_dwgn(mmpmatsplit,dwgn=dwgn,dwgnp=dwgnp,inverse=inverse)
273 0 : mmpmatOut = mmpmatOutsplit(:,:,1)
274 :
275 0 : END FUNCTION rotMMPmat_dwgn_one_spin
276 :
277 6192 : PURE FUNCTION d_wigner_mat(alpha, beta, gamma, l) RESULT(d)
278 :
279 : REAL, INTENT(IN) :: alpha,beta,gamma !Euler angles
280 : INTEGER, INTENT(IN) :: l
281 :
282 : COMPLEX :: d(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const)
283 :
284 : INTEGER :: m,mp,x_lo,x_up,x,e_c,e_s
285 : REAL :: fac_l_m,fac_l_mp,fac_lmpx,fac_lmx,fac_x,fac_xmpm
286 : REAL :: co_bh,si_bh,zaehler,nenner,cp,sp
287 : COMPLEX :: phase_g,phase_a,bas
288 :
289 :
290 6192 : co_bh = cos(beta*0.5)
291 6192 : si_bh = -sin(beta*0.5)
292 352944 : d = cmplx_0
293 :
294 46788 : DO m = -l,l
295 121788 : fac_l_m = fac(l+m) * fac(l-m)
296 40596 : phase_g = exp( - ImagUnit * gamma * m )
297 :
298 317220 : DO mp = -l,l
299 811296 : fac_l_mp = fac(l+mp) * fac(l-mp)
300 :
301 270432 : zaehler = sqrt( real(fac_l_m * fac_l_mp) )
302 270432 : phase_a = exp( - ImagUnit * alpha * mp )
303 270432 : x_lo = max(0, m-mp)
304 270432 : x_up = min(l-mp, l+m)
305 :
306 270432 : bas = zaehler * phase_a * phase_g
307 763830 : DO x = x_lo,x_up
308 905604 : fac_lmpx = fac(l-mp-x)
309 905604 : fac_lmx = fac(l+m-x)
310 452802 : fac_x = fac(x)
311 905604 : fac_xmpm = fac(x+mp-m)
312 452802 : nenner = fac_lmpx * fac_lmx * fac_x * fac_xmpm
313 452802 : e_c = 2*l + m - mp - 2*x
314 452802 : e_s = 2*x + mp - m
315 452802 : IF (e_c.EQ.0) THEN
316 : cp = 1.0
317 : ELSE
318 412206 : cp = co_bh ** e_c
319 : ENDIF
320 452802 : IF (e_s.EQ.0) THEN
321 : sp = 1.0
322 : ELSE
323 412206 : sp = si_bh ** e_s
324 : ENDIF
325 723234 : d(m,mp) = d(m,mp) + bas * (-1)**x * cp * sp / nenner
326 : ENDDO
327 :
328 : ENDDO ! loop over mp
329 : ENDDO ! loop over m
330 46788 : DO m = -l,l
331 317220 : DO mp = -l,l
332 311028 : d( m,mp ) = d( m,mp ) * (-1)**(m-mp)
333 : ENDDO
334 : ENDDO
335 :
336 6192 : END FUNCTION d_wigner_mat
337 :
338 2433264 : ELEMENTAL REAL FUNCTION fac(n)
339 :
340 : INTEGER, INTENT (IN) :: n
341 : INTEGER :: i
342 :
343 1669434 : fac = 0
344 1669434 : IF (n.LT.0) RETURN
345 2433264 : fac = 1
346 2230962 : IF (n.EQ.0) RETURN
347 4404360 : DO i = 2,n
348 4404360 : fac = fac * i
349 : ENDDO
350 :
351 : END FUNCTION fac
352 :
353 :
354 : END MODULE m_rotMMPmat
|
