Line data Source code
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 Calculation of the spherical harmonics and the corresponding orbital
10 : !> transformation matrices.
11 : !> \par Literature
12 : !> H. B. Schlegel, M. J. Frisch, Int. J. Quantum Chem. 54, 83 (1995)
13 : !> \par History
14 : !> - restructured and cleaned (20.05.2004,MK)
15 : !> \author Matthias Krack (08.06.2000)
16 : ! **************************************************************************************************
17 : MODULE orbital_transformation_matrices
18 :
19 : USE kinds, ONLY: dp
20 : USE mathconstants, ONLY: dfac,&
21 : fac,&
22 : pi
23 : USE mathlib, ONLY: binomial
24 : USE orbital_pointers, ONLY: co,&
25 : nco,&
26 : ncoset,&
27 : nso,&
28 : nsoset
29 : USE orbital_symbols, ONLY: cgf_symbol,&
30 : sgf_symbol
31 :
32 : !$ USE OMP_LIB, ONLY: omp_get_level
33 :
34 : #include "../base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 :
38 : PRIVATE
39 :
40 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'orbital_transformation_matrices'
41 :
42 : TYPE orbtramat_type
43 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: c2s => NULL(), slm => NULL(), &
44 : slm_inv => NULL(), s2c => NULL()
45 : END TYPE orbtramat_type
46 :
47 : TYPE(orbtramat_type), DIMENSION(:), POINTER :: orbtramat => NULL()
48 :
49 : TYPE orbrotmat_type
50 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: mat => NULL()
51 : END TYPE orbrotmat_type
52 :
53 : INTEGER, SAVE :: current_maxl = -1
54 :
55 : PUBLIC :: orbtramat
56 : PUBLIC :: orbrotmat_type, calculate_rotmat, release_rotmat
57 : PUBLIC :: deallocate_spherical_harmonics, init_spherical_harmonics
58 :
59 : CONTAINS
60 :
61 : ! **************************************************************************************************
62 : !> \brief Allocate and initialize the orbital transformation matrices for
63 : !> the transformation and for the backtransformation of orbitals
64 : !> from the Cartesian representation to the spherical representation.
65 : !> \param maxl ...
66 : !> \date 20.05.2004
67 : !> \author MK
68 : !> \version 1.0
69 : ! **************************************************************************************************
70 6202 : SUBROUTINE create_spherical_harmonics(maxl)
71 :
72 : INTEGER, INTENT(IN) :: maxl
73 :
74 : INTEGER :: expo, i, ic, ic1, ic2, is, j, k, l, lx, &
75 : lx1, lx2, ly, ly1, ly2, lz, lz1, lz2, &
76 : m, ma, nc, ns
77 : REAL(KIND=dp) :: s, s1, s2
78 :
79 6202 : !$ IF (omp_get_level() > 0) &
80 0 : !$ CPABORT("create_spherical_harmonics is not thread-safe")
81 :
82 6202 : IF (current_maxl > -1) THEN
83 : CALL cp_abort(__LOCATION__, &
84 : "Spherical harmonics are already allocated. "// &
85 0 : "Use the init routine for an update")
86 : END IF
87 :
88 6202 : IF (maxl < 0) THEN
89 : CALL cp_abort(__LOCATION__, &
90 : "A negative maximum angular momentum quantum "// &
91 0 : "number is invalid")
92 : END IF
93 :
94 51614 : ALLOCATE (orbtramat(0:maxl))
95 6202 : nc = ncoset(maxl)
96 6202 : ns = nsoset(maxl)
97 :
98 39210 : DO l = 0, maxl
99 33008 : nc = nco(l)
100 33008 : ns = nso(l)
101 132032 : ALLOCATE (orbtramat(l)%c2s(ns, nc))
102 2799550 : orbtramat(l)%c2s = 0.0_dp
103 99024 : ALLOCATE (orbtramat(l)%s2c(ns, nc))
104 2799550 : orbtramat(l)%s2c = 0.0_dp
105 99024 : ALLOCATE (orbtramat(l)%slm(ns, nc))
106 2799550 : orbtramat(l)%slm = 0.0_dp
107 99024 : ALLOCATE (orbtramat(l)%slm_inv(ns, nc))
108 2805752 : orbtramat(l)%slm_inv = 0.0_dp
109 : END DO
110 :
111 : ! Loop over all angular momentum quantum numbers l
112 :
113 39210 : DO l = 0, maxl
114 :
115 : ! Build the orbital transformation matrix for the transformation
116 : ! from Cartesian to spherical orbitals (c2s, formula 15)
117 :
118 142325 : DO lx = 0, l
119 429099 : DO ly = 0, l - lx
120 286774 : lz = l - lx - ly
121 286774 : ic = co(lx, ly, lz)
122 2875859 : DO m = -l, l
123 2479768 : is = l + m + 1
124 2479768 : ma = ABS(m)
125 2479768 : j = lx + ly - ma
126 2766542 : IF ((j >= 0) .AND. (MODULO(j, 2) == 0)) THEN
127 1017772 : j = j/2
128 1017772 : s1 = 0.0_dp
129 3006704 : DO i = 0, (l - ma)/2
130 1988932 : s2 = 0.0_dp
131 5782228 : DO k = 0, j
132 3111045 : IF (((m < 0) .AND. (MODULO(ABS(ma - lx), 2) == 1)) .OR. &
133 3793296 : ((m > 0) .AND. (MODULO(ABS(ma - lx), 2) == 0))) THEN
134 1455522 : expo = (ma - lx + 2*k)/2
135 1455522 : s = (-1.0_dp)**expo*SQRT(2.0_dp)
136 2337774 : ELSE IF ((m == 0) .AND. (MODULO(lx, 2) == 0)) THEN
137 570436 : expo = k - lx/2
138 570436 : s = (-1.0_dp)**expo
139 : ELSE
140 : s = 0.0_dp
141 : END IF
142 5782228 : s2 = s2 + binomial(j, k)*binomial(ma, lx - 2*k)*s
143 : END DO
144 : s1 = s1 + binomial(l, i)*binomial(i, j)* &
145 3006704 : (-1.0_dp)**i*fac(2*l - 2*i)/fac(l - ma - 2*i)*s2
146 : END DO
147 : orbtramat(l)%c2s(is, ic) = &
148 : SQRT((fac(2*lx)*fac(2*ly)*fac(2*lz)*fac(l)*fac(l - ma))/ &
149 : (fac(lx)*fac(ly)*fac(lz)*fac(2*l)*fac(l + ma)))*s1/ &
150 1017772 : (2.0_dp**l*fac(l))
151 : ELSE
152 1461996 : orbtramat(l)%c2s(is, ic) = 0.0_dp
153 : END IF
154 : END DO
155 : END DO
156 : END DO
157 :
158 : ! Build the corresponding transformation matrix for the transformation from
159 : ! spherical to Cartesian orbitals (s2c = s*TRANSPOSE(c2s), formulas 18 and 19)
160 :
161 142325 : DO lx1 = 0, l
162 429099 : DO ly1 = 0, l - lx1
163 286774 : lz1 = l - lx1 - ly1
164 286774 : ic1 = co(lx1, ly1, lz1)
165 : s1 = SQRT((fac(lx1)*fac(ly1)*fac(lz1))/ &
166 286774 : (fac(2*lx1)*fac(2*ly1)*fac(2*lz1)))
167 1779362 : DO lx2 = 0, l
168 6128025 : DO ly2 = 0, l - lx2
169 4457980 : lz2 = l - lx2 - ly2
170 4457980 : ic2 = co(lx2, ly2, lz2)
171 4457980 : lx = lx1 + lx2
172 4457980 : ly = ly1 + ly2
173 4457980 : lz = lz1 + lz2
174 : IF ((MODULO(lx, 2) == 0) .AND. &
175 4457980 : (MODULO(ly, 2) == 0) .AND. &
176 1383271 : (MODULO(lz, 2) == 0)) THEN
177 : s2 = SQRT((fac(lx2)*fac(ly2)*fac(lz2))/ &
178 1222042 : (fac(2*lx2)*fac(2*ly2)*fac(2*lz2)))
179 : s = fac(lx)*fac(ly)*fac(lz)*s1*s2/ &
180 1222042 : (fac(lx/2)*fac(ly/2)*fac(lz/2))
181 14237078 : DO is = 1, nso(l)
182 : orbtramat(l)%s2c(is, ic1) = orbtramat(l)%s2c(is, ic1) + &
183 14237078 : s*orbtramat(l)%c2s(is, ic2)
184 : END DO
185 : END IF
186 : END DO
187 : END DO
188 : END DO
189 : END DO
190 :
191 : ! Build up the real spherical harmonics
192 :
193 33008 : s = SQRT(0.25_dp*dfac(2*l + 1)/pi)
194 :
195 148527 : DO lx = 0, l
196 429099 : DO ly = 0, l - lx
197 286774 : lz = l - lx - ly
198 286774 : ic = co(lx, ly, lz)
199 2875859 : DO m = -l, l
200 2479768 : is = l + m + 1
201 2479768 : s1 = SQRT(dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1))
202 2479768 : orbtramat(l)%slm(is, ic) = s*orbtramat(l)%c2s(is, ic)/s1
203 : !MK s2 = (-1.0_dp)**m*s ! alternative S(lm) definition
204 : !MK orbtramat(l)%slm(is, ic) = s2*orbtramat(l)%c2s(is,ic)/s1
205 2766542 : orbtramat(l)%slm_inv(is, ic) = s1*orbtramat(l)%s2c(is, ic)/s
206 : END DO
207 : END DO
208 : END DO
209 :
210 : END DO
211 :
212 : ! Save initialization status
213 :
214 6202 : current_maxl = maxl
215 :
216 6202 : END SUBROUTINE create_spherical_harmonics
217 :
218 : ! **************************************************************************************************
219 : !> \brief Deallocate the orbital transformation matrices.
220 : !> \date 20.05.2004
221 : !> \author MK
222 : !> \version 1.0
223 : ! **************************************************************************************************
224 15014 : SUBROUTINE deallocate_spherical_harmonics()
225 :
226 : INTEGER :: l
227 :
228 15014 : !$ IF (omp_get_level() > 0) &
229 0 : !$ CPABORT("deallocate_spherical_harmonics is not thread-safe")
230 :
231 15014 : IF (current_maxl > -1) THEN
232 39210 : DO l = 0, SIZE(orbtramat, 1) - 1
233 33008 : DEALLOCATE (orbtramat(l)%c2s)
234 33008 : DEALLOCATE (orbtramat(l)%s2c)
235 33008 : DEALLOCATE (orbtramat(l)%slm)
236 39210 : DEALLOCATE (orbtramat(l)%slm_inv)
237 : END DO
238 6202 : DEALLOCATE (orbtramat)
239 6202 : current_maxl = -1
240 : END IF
241 :
242 15014 : END SUBROUTINE deallocate_spherical_harmonics
243 :
244 : ! **************************************************************************************************
245 : !> \brief Initialize or update the orbital transformation matrices.
246 : !> \param maxl ...
247 : !> \param output_unit ...
248 : !> \date 09.07.1999
249 : !> \par Variables
250 : !> - maxl : Maximum angular momentum quantum number
251 : !> \author MK
252 : !> \version 1.0
253 : ! **************************************************************************************************
254 23230 : SUBROUTINE init_spherical_harmonics(maxl, output_unit)
255 :
256 : INTEGER, INTENT(IN) :: maxl
257 : INTEGER :: output_unit
258 :
259 : CHARACTER(LEN=78) :: headline
260 : INTEGER :: l, nc, ns
261 :
262 23230 : !$ IF (omp_get_level() > 0) &
263 0 : !$ CPABORT("init_spherical_harmonics is not thread-safe")
264 :
265 23230 : IF (maxl < 0) THEN
266 : CALL cp_abort(__LOCATION__, &
267 : "A negative maximum angular momentum quantum "// &
268 0 : "number is invalid")
269 : END IF
270 :
271 23230 : IF (maxl > current_maxl) THEN
272 :
273 6202 : CALL deallocate_spherical_harmonics()
274 6202 : CALL create_spherical_harmonics(maxl)
275 :
276 : ! Print the spherical harmonics and the orbital transformation matrices
277 :
278 6202 : IF (output_unit > 0) THEN
279 :
280 4 : DO l = 0, maxl
281 :
282 3 : nc = nco(l)
283 3 : ns = nso(l)
284 :
285 : headline = "CARTESIAN ORBITAL TO SPHERICAL ORBITAL "// &
286 3 : "TRANSFORMATION MATRIX"
287 3 : CALL write_matrix(orbtramat(l)%c2s, l, output_unit, headline)
288 :
289 : headline = "SPHERICAL ORBITAL TO CARTESIAN ORBITAL "// &
290 3 : "TRANSFORMATION MATRIX"
291 3 : CALL write_matrix(orbtramat(l)%s2c, l, output_unit, headline)
292 :
293 3 : headline = "SPHERICAL HARMONICS"
294 3 : CALL write_matrix(orbtramat(l)%slm, l, output_unit, headline)
295 :
296 3 : headline = "INVERSE SPHERICAL HARMONICS"
297 4 : CALL write_matrix(orbtramat(l)%slm_inv, l, output_unit, headline)
298 :
299 : END DO
300 :
301 1 : WRITE (UNIT=output_unit, FMT="(A)") ""
302 :
303 : END IF
304 :
305 : END IF
306 :
307 23230 : END SUBROUTINE init_spherical_harmonics
308 :
309 : ! **************************************************************************************************
310 : !> \brief Calculate rotation matrices for spherical harmonics up to value l
311 : !> Joseph Ivanic and Klaus Ruedenberg
312 : !> Rotation Matrices for Real Spherical Harmonics. Direct Determination by Recursion
313 : !> J. Phys. Chem. 1996, 100, 6342-6347
314 : !> \param orbrotmat ...
315 : !> \param rotmat ...
316 : !> \param lval ...
317 : ! **************************************************************************************************
318 0 : SUBROUTINE calculate_rotmat(orbrotmat, rotmat, lval)
319 : TYPE(orbrotmat_type), DIMENSION(:), POINTER :: orbrotmat
320 : REAL(KIND=dp), DIMENSION(3, 3) :: rotmat
321 : INTEGER, INTENT(IN) :: lval
322 :
323 : INTEGER :: l, m1, m2, ns
324 : REAL(KIND=dp) :: s3, u, uf, v, vf, w, wf
325 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r
326 : REAL(KIND=dp), DIMENSION(-1:1, -1:1) :: r1
327 : REAL(KIND=dp), DIMENSION(-2:2, -2:2) :: r2
328 : REAL(KIND=dp), DIMENSION(3, 3) :: t
329 :
330 : MARK_USED(rotmat)
331 :
332 0 : CALL release_rotmat(orbrotmat)
333 :
334 0 : ALLOCATE (orbrotmat(0:lval))
335 0 : DO l = 0, lval
336 0 : ns = nso(l)
337 0 : ALLOCATE (orbrotmat(l)%mat(ns, ns))
338 : END DO
339 :
340 0 : IF (lval >= 0) THEN
341 0 : orbrotmat(0)%mat = 1.0_dp
342 : END IF
343 0 : IF (lval >= 1) THEN
344 0 : t(1, 1:3) = rotmat(2, 1:3)
345 0 : t(2, 1:3) = rotmat(3, 1:3)
346 0 : t(3, 1:3) = rotmat(1, 1:3)
347 0 : r1(-1:1, -1) = t(1:3, 2)
348 0 : r1(-1:1, 0) = t(1:3, 3)
349 0 : r1(-1:1, 1) = t(1:3, 1)
350 0 : orbrotmat(1)%mat(1:3, 1:3) = r1(-1:1, -1:1)
351 : END IF
352 0 : IF (lval >= 2) THEN
353 0 : s3 = SQRT(3.0_dp)
354 : ! Table 4
355 0 : r2(0, 0) = (3.0_dp*r1(0, 0)**2 - 1.0_dp)*0.5_dp
356 0 : r2(0, 1) = s3*r1(0, 1)*r1(0, 0)
357 0 : r2(0, -1) = s3*r1(0, -1)*r1(0, 0)
358 : r2(0, -2) = s3*(r1(0, 1)**2 - r1(0, -1)**2)*0.5_dp
359 0 : r2(0, -2) = s3*r1(0, 1)*r1(0, -1)
360 : !
361 0 : r2(1, 0) = s3*r1(1, 0)*r1(0, 0)
362 0 : r2(1, 1) = r1(1, 1)*r1(0, 0) + r1(1, 0)*r1(0, 1)
363 0 : r2(1, -1) = r1(1, -1)*r1(0, 0) + r1(1, 0)*r1(0, -1)
364 0 : r2(1, 2) = r1(1, 1)*r1(0, 1) - r1(1, -1)*r1(0, -1)
365 0 : r2(1, -2) = r1(1, 1)*r1(0, -1) + r1(1, -1)*r1(0, 1)
366 : !
367 0 : r2(-1, 0) = s3*r1(-1, 0)*r1(0, 0)
368 0 : r2(-1, 1) = r1(-1, 1)*r1(0, 0) + r1(-1, 0)*r1(0, 1)
369 0 : r2(-1, -1) = r1(-1, -1)*r1(0, 0) + r1(-1, 0)*r1(0, -1)
370 0 : r2(-1, 2) = r1(-1, 1)*r1(0, 1) - r1(-1, -1)*r1(0, -1)
371 0 : r2(-1, -2) = r1(-1, 1)*r1(0, -1) + r1(-1, -1)*r1(0, 1)
372 : !
373 0 : r2(2, 0) = s3*(r1(1, 0)**2 - r1(-1, 0)**2)*0.5_dp
374 : r2(2, 1) = r1(1, 1)*r1(1, 0) - r1(-1, 1)*r1(-1, 0)
375 : r2(2, -1) = r1(1, -1)*r1(1, 0) - r1(-1, -1)*r1(-1, 0)
376 : r2(2, 2) = (r1(1, 1)**2 - r1(1, -1)**2 - r1(-1, 1)**2 + r1(-1, -1)**2)*0.5_dp
377 : r2(2, -2) = r1(1, 1)*r1(1, -1) - r1(-1, 1)*r1(-1, -1)
378 : !
379 0 : r2(-2, 0) = s3*r1(1, 0)*r1(-1, 0)
380 0 : r2(2, 1) = r1(1, 1)*r1(-1, 0) + r1(1, 0)*r1(-1, 1)
381 0 : r2(2, -1) = r1(1, -1)*r1(-1, 0) + r1(1, 0)*r1(-1, -1)
382 0 : r2(2, 2) = r1(1, 1)*r1(-1, 1) - r1(1, -1)*r1(-1, -1)
383 0 : r2(2, -2) = r1(1, 1)*r1(-1, -1) + r1(1, -1)*r1(-1, 1)
384 : !
385 0 : orbrotmat(2)%mat(1:5, 1:5) = r2(-2:2, -2:2)
386 : END IF
387 0 : IF (lval >= 3) THEN
388 : ! General recursion
389 0 : ALLOCATE (r(0:lval, -lval:lval, -lval:lval))
390 0 : r = 0.0_dp
391 0 : r(0, 0, 0) = 1.0_dp
392 0 : r(1, -1:1, -1:1) = r1(-1:1, -1:1)
393 0 : r(2, -2:2, -2:2) = r2(-2:2, -2:2)
394 0 : DO l = 3, lval
395 0 : DO m1 = -l, l
396 0 : DO m2 = -l, l
397 0 : u = u_func(l, m1, m2)
398 0 : v = v_func(l, m1, m2)
399 0 : w = w_func(l, m1, m2)
400 0 : CALL r_recursion(l, m1, m2, r1, r, lval, uf, vf, wf)
401 0 : r(l, m1, m2) = u*uf + v*vf + w*wf
402 : END DO
403 : END DO
404 : END DO
405 0 : DO l = 3, lval
406 0 : ns = nso(l)
407 0 : orbrotmat(l)%mat(1:ns, 1:ns) = r(l, -l:l, -l:l)
408 : END DO
409 0 : DEALLOCATE (r)
410 : END IF
411 :
412 0 : END SUBROUTINE calculate_rotmat
413 :
414 : ! **************************************************************************************************
415 : !> \brief ...
416 : !> \param l ...
417 : !> \param ma ...
418 : !> \param mb ...
419 : !> \return ...
420 : ! **************************************************************************************************
421 0 : FUNCTION u_func(l, ma, mb) RESULT(u)
422 : INTEGER :: l, ma, mb
423 : REAL(KIND=dp) :: u
424 :
425 0 : IF (ABS(mb) == l) THEN
426 0 : u = REAL((l + ma)*(l - ma), KIND=dp)/REAL(2*l*(2*l - 1), KIND=dp)
427 0 : u = SQRT(u)
428 0 : ELSE IF (ABS(mb) < l) THEN
429 0 : u = REAL((l + ma)*(l - ma), KIND=dp)/REAL((l + mb)*(l - mb), KIND=dp)
430 0 : u = SQRT(u)
431 : ELSE
432 0 : CPABORT("Illegal m value")
433 : END IF
434 0 : END FUNCTION u_func
435 :
436 : ! **************************************************************************************************
437 : !> \brief ...
438 : !> \param l ...
439 : !> \param ma ...
440 : !> \param mb ...
441 : !> \return ...
442 : ! **************************************************************************************************
443 0 : FUNCTION v_func(l, ma, mb) RESULT(v)
444 : INTEGER :: l, ma, mb
445 : REAL(KIND=dp) :: v
446 :
447 : INTEGER :: a1, a2, dm0
448 :
449 0 : dm0 = 0
450 0 : IF (ma == 0) dm0 = 1
451 0 : IF (ABS(mb) == l) THEN
452 0 : a1 = (1 + dm0)*(l + ABS(ma) - 1)*(l + ABS(ma))
453 0 : a2 = 2*l*(2*l - 1)
454 0 : ELSE IF (ABS(mb) < l) THEN
455 0 : a1 = (1 + dm0)*(l + ABS(ma) - 1)*(l + ABS(ma))
456 0 : a2 = (l + mb)*(l - mb)
457 : ELSE
458 0 : CPABORT("Illegal m value")
459 : END IF
460 0 : v = 0.5_dp*SQRT(REAL(a1, KIND=dp)/REAL(a2, KIND=dp))*REAL(1 - 2*dm0, KIND=dp)
461 0 : END FUNCTION v_func
462 :
463 : ! **************************************************************************************************
464 : !> \brief ...
465 : !> \param l ...
466 : !> \param ma ...
467 : !> \param mb ...
468 : !> \return ...
469 : ! **************************************************************************************************
470 0 : FUNCTION w_func(l, ma, mb) RESULT(w)
471 : INTEGER :: l, ma, mb
472 : REAL(KIND=dp) :: w
473 :
474 : INTEGER :: a1, a2, dm0
475 :
476 0 : dm0 = 0
477 0 : IF (ma == 0) dm0 = 1
478 0 : IF (ABS(mb) == l) THEN
479 0 : a1 = (l - ABS(ma) - 1)*(l - ABS(ma))
480 0 : a2 = 2*l*(2*l - 1)
481 0 : ELSE IF (ABS(mb) < l) THEN
482 0 : a1 = (l - ABS(ma) - 1)*(l - ABS(ma))
483 0 : a2 = (l + mb)*(l - mb)
484 : ELSE
485 0 : CPABORT("Illegal m value")
486 : END IF
487 0 : w = -0.5_dp*SQRT(REAL(a1, KIND=dp)/REAL(a2, KIND=dp))*REAL(1 - dm0, KIND=dp)
488 0 : END FUNCTION w_func
489 :
490 : ! **************************************************************************************************
491 : !> \brief ...
492 : !> \param i ...
493 : !> \param l ...
494 : !> \param mu ...
495 : !> \param m ...
496 : !> \param r1 ...
497 : !> \param r ...
498 : !> \param lr ...
499 : !> \return ...
500 : ! **************************************************************************************************
501 0 : FUNCTION p_func(i, l, mu, m, r1, r, lr) RESULT(p)
502 : INTEGER :: i, l, mu, m
503 : REAL(KIND=dp), DIMENSION(-1:1, -1:1) :: r1
504 : INTEGER :: lr
505 : REAL(KIND=dp), DIMENSION(0:lr, -lr:lr, -lr:lr) :: r
506 : REAL(KIND=dp) :: p
507 :
508 0 : IF (m == l) THEN
509 0 : p = r1(i, 1)*r(l - 1, mu, m) - r1(i, -1)*r(l - 1, mu, -m)
510 0 : ELSE IF (m == -l) THEN
511 0 : p = r1(i, 1)*r(l - 1, mu, m) + r1(i, -1)*r(l - 1, mu, -m)
512 0 : ELSE IF (ABS(m) < l) THEN
513 0 : p = r1(i, 0)*r(l - 1, mu, m)
514 : ELSE
515 0 : CPABORT("Illegal m value")
516 : END IF
517 0 : END FUNCTION p_func
518 :
519 : ! **************************************************************************************************
520 : !> \brief ...
521 : !> \param l ...
522 : !> \param ma ...
523 : !> \param mb ...
524 : !> \param r1 ...
525 : !> \param r ...
526 : !> \param lr ...
527 : !> \param u ...
528 : !> \param v ...
529 : !> \param w ...
530 : ! **************************************************************************************************
531 0 : SUBROUTINE r_recursion(l, ma, mb, r1, r, lr, u, v, w)
532 : INTEGER :: l, ma, mb
533 : REAL(KIND=dp), DIMENSION(-1:1, -1:1) :: r1
534 : INTEGER :: lr
535 : REAL(KIND=dp), DIMENSION(0:lr, -lr:lr, -lr:lr) :: r
536 : REAL(KIND=dp) :: u, v, w
537 :
538 : REAL(KIND=dp) :: dd
539 :
540 0 : IF (ma == 0) THEN
541 0 : u = p_func(0, l, 0, mb, r1, r, lr)
542 0 : v = p_func(1, l, 1, mb, r1, r, lr) + p_func(-1, l, -1, mb, r1, r, lr)
543 0 : w = 0.0_dp
544 0 : ELSE IF (ma > 0) THEN
545 : dd = 1.0_dp
546 : IF (ma == 1) dd = 1.0_dp
547 0 : u = p_func(0, l, ma, mb, r1, r, lr)
548 0 : v = p_func(1, l, ma - 1, mb, r1, r, lr)*SQRT(1.0_dp + dd) - p_func(-1, l, -ma + 1, mb, r1, r, lr)*(1.0_dp - dd)
549 0 : w = p_func(1, l, ma + 1, mb, r1, r, lr) + p_func(-1, l, -1 - ma, mb, r1, r, lr)
550 0 : ELSE IF (ma < 0) THEN
551 : dd = 1.0_dp
552 : IF (ma == -1) dd = 1.0_dp
553 0 : u = p_func(0, l, ma, mb, r1, r, lr)
554 0 : v = p_func(1, l, ma + 1, mb, r1, r, lr)*(1.0_dp - dd) + p_func(-1, l, -ma - 1, mb, r1, r, lr)*SQRT(1.0_dp + dd)
555 0 : w = p_func(1, l, ma - 1, mb, r1, r, lr) - p_func(-1, l, -ma + 1, mb, r1, r, lr)
556 : END IF
557 0 : END SUBROUTINE
558 :
559 : ! **************************************************************************************************
560 : !> \brief Release rotation matrices
561 : !> \param orbrotmat ...
562 : ! **************************************************************************************************
563 48 : SUBROUTINE release_rotmat(orbrotmat)
564 : TYPE(orbrotmat_type), DIMENSION(:), POINTER :: orbrotmat
565 :
566 : INTEGER :: i
567 :
568 48 : IF (ASSOCIATED(orbrotmat)) THEN
569 0 : DO i = LBOUND(orbrotmat, 1), UBOUND(orbrotmat, 1)
570 0 : IF (ASSOCIATED(orbrotmat(i)%mat)) DEALLOCATE (orbrotmat(i)%mat)
571 : END DO
572 0 : DEALLOCATE (orbrotmat)
573 : END IF
574 :
575 48 : END SUBROUTINE release_rotmat
576 :
577 : ! **************************************************************************************************
578 : !> \brief Print a matrix to the logical unit "lunit".
579 : !> \param matrix ...
580 : !> \param l ...
581 : !> \param lunit ...
582 : !> \param headline ...
583 : !> \date 07.06.2000
584 : !> \par Variables
585 : !> - matrix : Matrix to be printed.
586 : !> - l : Angular momentum quantum number
587 : !> - lunit : Logical unit number.
588 : !> - headline: Headline of the matrix.
589 : !> \author MK
590 : !> \version 1.0
591 : ! **************************************************************************************************
592 12 : SUBROUTINE write_matrix(matrix, l, lunit, headline)
593 :
594 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
595 : INTEGER, INTENT(IN) :: l, lunit
596 : CHARACTER(LEN=*), INTENT(IN) :: headline
597 :
598 : CHARACTER(12) :: symbol
599 : CHARACTER(LEN=78) :: string
600 : INTEGER :: from, i, ic, is, jc, lx, ly, lz, m, nc, &
601 : to
602 :
603 : ! Write headline
604 :
605 12 : WRITE (UNIT=lunit, FMT="(/,/,T2,A)") TRIM(headline)
606 :
607 : ! Write the matrix in the defined format
608 :
609 12 : nc = nco(l)
610 :
611 28 : DO ic = 1, nc, 3
612 16 : from = ic
613 16 : to = MIN(nc, from + 2)
614 16 : i = 1
615 52 : DO lx = l, 0, -1
616 116 : DO ly = l - lx, 0, -1
617 64 : lz = l - lx - ly
618 64 : jc = co(lx, ly, lz)
619 100 : IF ((jc >= from) .AND. (jc <= to)) THEN
620 160 : symbol = cgf_symbol(1, (/lx, ly, lz/))
621 40 : WRITE (UNIT=string(i:), FMT="(A18)") TRIM(symbol(3:12))
622 40 : i = i + 18
623 : END IF
624 : END DO
625 : END DO
626 16 : WRITE (UNIT=lunit, FMT="(/,T8,A)") TRIM(string)
627 16 : symbol = ""
628 84 : DO m = -l, l
629 56 : is = l + m + 1
630 56 : symbol = sgf_symbol(1, l, m)
631 : WRITE (UNIT=lunit, FMT="(T4,A4,3(1X,F17.12))") &
632 72 : symbol(3:6), (matrix(is, jc), jc=from, to)
633 : END DO
634 : END DO
635 :
636 12 : END SUBROUTINE write_matrix
637 :
638 0 : END MODULE orbital_transformation_matrices
|
