Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 8114 : 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 8114 : !$ IF (omp_get_level() > 0) &
80 0 : !$ CPABORT("create_spherical_harmonics is not thread-safe")
81 :
82 8114 : 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 8114 : 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 67322 : ALLOCATE (orbtramat(0:maxl))
95 8114 : nc = ncoset(maxl)
96 8114 : ns = nsoset(maxl)
97 :
98 51094 : DO l = 0, maxl
99 42980 : nc = nco(l)
100 42980 : ns = nso(l)
101 171920 : ALLOCATE (orbtramat(l)%c2s(ns, nc))
102 3592642 : orbtramat(l)%c2s = 0.0_dp
103 128940 : ALLOCATE (orbtramat(l)%s2c(ns, nc))
104 3592642 : orbtramat(l)%s2c = 0.0_dp
105 128940 : ALLOCATE (orbtramat(l)%slm(ns, nc))
106 3592642 : orbtramat(l)%slm = 0.0_dp
107 128940 : ALLOCATE (orbtramat(l)%slm_inv(ns, nc))
108 3600756 : orbtramat(l)%slm_inv = 0.0_dp
109 : END DO
110 :
111 : ! Loop over all angular momentum quantum numbers l
112 :
113 51094 : 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 184665 : DO lx = 0, l
119 554695 : DO ly = 0, l - lx
120 370030 : lz = l - lx - ly
121 370030 : ic = co(lx, ly, lz)
122 3691347 : DO m = -l, l
123 3179632 : is = l + m + 1
124 3179632 : ma = ABS(m)
125 3179632 : j = lx + ly - ma
126 3549662 : IF ((j >= 0) .AND. (MODULO(j, 2) == 0)) THEN
127 1306564 : j = j/2
128 1306564 : s1 = 0.0_dp
129 3850288 : DO i = 0, (l - ma)/2
130 2543724 : s2 = 0.0_dp
131 7386496 : DO k = 0, j
132 3972139 : IF (((m < 0) .AND. (MODULO(ABS(ma - lx), 2) == 1)) .OR. &
133 4842772 : ((m > 0) .AND. (MODULO(ABS(ma - lx), 2) == 0))) THEN
134 1857282 : expo = (ma - lx + 2*k)/2
135 1857282 : s = (-1.0_dp)**expo*SQRT(2.0_dp)
136 2985490 : ELSE IF ((m == 0) .AND. (MODULO(lx, 2) == 0)) THEN
137 730328 : expo = k - lx/2
138 730328 : s = (-1.0_dp)**expo
139 : ELSE
140 : s = 0.0_dp
141 : END IF
142 7386496 : s2 = s2 + binomial(j, k)*binomial(ma, lx - 2*k)*s
143 : END DO
144 : s1 = s1 + binomial(l, i)*binomial(i, j)* &
145 3850288 : (-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 1306564 : (2.0_dp**l*fac(l))
151 : ELSE
152 1873068 : 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 184665 : DO lx1 = 0, l
162 554695 : DO ly1 = 0, l - lx1
163 370030 : lz1 = l - lx1 - ly1
164 370030 : ic1 = co(lx1, ly1, lz1)
165 : s1 = SQRT((fac(lx1)*fac(ly1)*fac(lz1))/ &
166 370030 : (fac(2*lx1)*fac(2*ly1)*fac(2*lz1)))
167 2286546 : DO lx2 = 0, l
168 7842753 : DO ly2 = 0, l - lx2
169 5697892 : lz2 = l - lx2 - ly2
170 5697892 : ic2 = co(lx2, ly2, lz2)
171 5697892 : lx = lx1 + lx2
172 5697892 : ly = ly1 + ly2
173 5697892 : lz = lz1 + lz2
174 : IF ((MODULO(lx, 2) == 0) .AND. &
175 5697892 : (MODULO(ly, 2) == 0) .AND. &
176 1774831 : (MODULO(lz, 2) == 0)) THEN
177 : s2 = SQRT((fac(lx2)*fac(ly2)*fac(lz2))/ &
178 1563250 : (fac(2*lx2)*fac(2*ly2)*fac(2*lz2)))
179 : s = fac(lx)*fac(ly)*fac(lz)*s1*s2/ &
180 1563250 : (fac(lx/2)*fac(ly/2)*fac(lz/2))
181 18162158 : DO is = 1, nso(l)
182 : orbtramat(l)%s2c(is, ic1) = orbtramat(l)%s2c(is, ic1) + &
183 18162158 : 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 42980 : s = SQRT(0.25_dp*dfac(2*l + 1)/pi)
194 :
195 192779 : DO lx = 0, l
196 554695 : DO ly = 0, l - lx
197 370030 : lz = l - lx - ly
198 370030 : ic = co(lx, ly, lz)
199 3691347 : DO m = -l, l
200 3179632 : is = l + m + 1
201 3179632 : s1 = SQRT(dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1))
202 3179632 : 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 3549662 : 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 8114 : current_maxl = maxl
215 :
216 8114 : 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 18842 : SUBROUTINE deallocate_spherical_harmonics()
225 :
226 : INTEGER :: l
227 :
228 18842 : !$ IF (omp_get_level() > 0) &
229 0 : !$ CPABORT("deallocate_spherical_harmonics is not thread-safe")
230 :
231 18842 : IF (current_maxl > -1) THEN
232 51094 : DO l = 0, SIZE(orbtramat, 1) - 1
233 42980 : DEALLOCATE (orbtramat(l)%c2s)
234 42980 : DEALLOCATE (orbtramat(l)%s2c)
235 42980 : DEALLOCATE (orbtramat(l)%slm)
236 51094 : DEALLOCATE (orbtramat(l)%slm_inv)
237 : END DO
238 8114 : DEALLOCATE (orbtramat)
239 8114 : current_maxl = -1
240 : END IF
241 :
242 18842 : 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 30136 : 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 30136 : !$ IF (omp_get_level() > 0) &
263 0 : !$ CPABORT("init_spherical_harmonics is not thread-safe")
264 :
265 30136 : 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 30136 : IF (maxl > current_maxl) THEN
272 :
273 8114 : CALL deallocate_spherical_harmonics()
274 8114 : CALL create_spherical_harmonics(maxl)
275 :
276 : ! Print the spherical harmonics and the orbital transformation matrices
277 :
278 8114 : 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 30136 : 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 21188 : 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 21188 : 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 21188 : CALL release_rotmat(orbrotmat)
331 :
332 114694 : ALLOCATE (orbrotmat(0:lval))
333 72318 : DO l = 0, lval
334 51130 : ns = nso(l)
335 225708 : ALLOCATE (orbrotmat(l)%mat(ns, ns))
336 : END DO
337 :
338 21188 : IF (lval >= 0) THEN
339 63564 : orbrotmat(0)%mat = 1.0_dp
340 : END IF
341 21188 : IF (lval >= 1) THEN
342 77920 : t(1, 1:3) = rotmat(2, 1:3)
343 77920 : t(2, 1:3) = rotmat(3, 1:3)
344 77920 : t(3, 1:3) = rotmat(1, 1:3)
345 77920 : r1(-1:1, -1) = t(1:3, 2)
346 77920 : r1(-1:1, 0) = t(1:3, 3)
347 77920 : r1(-1:1, 1) = t(1:3, 1)
348 253240 : orbrotmat(1)%mat(1:3, 1:3) = r1(-1:1, -1:1)
349 : END IF
350 21188 : IF (lval >= 2) THEN
351 10444 : s3 = SQRT(3.0_dp)
352 : ! Table 4
353 10444 : r2(0, 0) = (3.0_dp*r1(0, 0)**2 - 1.0_dp)*0.5_dp
354 10444 : r2(0, 1) = s3*r1(0, 1)*r1(0, 0)
355 10444 : r2(0, -1) = s3*r1(0, -1)*r1(0, 0)
356 10444 : r2(0, 2) = s3*(r1(0, 1)**2 - r1(0, -1)**2)*0.5_dp
357 10444 : r2(0, -2) = s3*r1(0, 1)*r1(0, -1)
358 : !
359 10444 : r2(1, 0) = s3*r1(1, 0)*r1(0, 0)
360 10444 : r2(1, 1) = r1(1, 1)*r1(0, 0) + r1(1, 0)*r1(0, 1)
361 10444 : r2(1, -1) = r1(1, -1)*r1(0, 0) + r1(1, 0)*r1(0, -1)
362 10444 : r2(1, 2) = r1(1, 1)*r1(0, 1) - r1(1, -1)*r1(0, -1)
363 10444 : r2(1, -2) = r1(1, 1)*r1(0, -1) + r1(1, -1)*r1(0, 1)
364 : !
365 10444 : r2(-1, 0) = s3*r1(-1, 0)*r1(0, 0)
366 10444 : r2(-1, 1) = r1(-1, 1)*r1(0, 0) + r1(-1, 0)*r1(0, 1)
367 10444 : r2(-1, -1) = r1(-1, -1)*r1(0, 0) + r1(-1, 0)*r1(0, -1)
368 10444 : r2(-1, 2) = r1(-1, 1)*r1(0, 1) - r1(-1, -1)*r1(0, -1)
369 10444 : r2(-1, -2) = r1(-1, 1)*r1(0, -1) + r1(-1, -1)*r1(0, 1)
370 : !
371 10444 : r2(2, 0) = s3*(r1(1, 0)**2 - r1(-1, 0)**2)*0.5_dp
372 10444 : r2(2, 1) = r1(1, 1)*r1(1, 0) - r1(-1, 1)*r1(-1, 0)
373 10444 : r2(2, -1) = r1(1, -1)*r1(1, 0) - r1(-1, -1)*r1(-1, 0)
374 10444 : r2(2, 2) = (r1(1, 1)**2 - r1(1, -1)**2 - r1(-1, 1)**2 + r1(-1, -1)**2)*0.5_dp
375 10444 : r2(2, -2) = r1(1, 1)*r1(1, -1) - r1(-1, 1)*r1(-1, -1)
376 : !
377 10444 : r2(-2, 0) = s3*r1(1, 0)*r1(-1, 0)
378 10444 : r2(-2, 1) = r1(1, 1)*r1(-1, 0) + r1(1, 0)*r1(-1, 1)
379 10444 : r2(-2, -1) = r1(1, -1)*r1(-1, 0) + r1(1, 0)*r1(-1, -1)
380 10444 : r2(-2, 2) = r1(1, 1)*r1(-1, 1) - r1(1, -1)*r1(-1, -1)
381 10444 : r2(-2, -2) = r1(1, 1)*r1(-1, -1) + r1(1, -1)*r1(-1, 1)
382 : !
383 323764 : orbrotmat(2)%mat(1:5, 1:5) = r2(-2:2, -2:2)
384 : END IF
385 21188 : IF (lval >= 3) THEN
386 : ! General recursion
387 30 : ALLOCATE (r(0:lval, -lval:lval, -lval:lval))
388 6 : r = 0.0_dp
389 6 : r(0, 0, 0) = 1.0_dp
390 78 : r(1, -1:1, -1:1) = r1(-1:1, -1:1)
391 186 : r(2, -2:2, -2:2) = r2(-2:2, -2:2)
392 24 : DO l = 3, lval
393 186 : DO m1 = -l, l
394 1686 : DO m2 = -l, l
395 1506 : u = u_func(l, m1, m2)
396 1506 : v = v_func(l, m1, m2)
397 1506 : w = w_func(l, m1, m2)
398 1506 : CALL r_recursion(l, m1, m2, r1, r, lval, uf, vf, wf)
399 1668 : r(l, m1, m2) = u*uf + v*vf + w*wf
400 : END DO
401 : END DO
402 : END DO
403 24 : DO l = 3, lval
404 18 : ns = nso(l)
405 1692 : orbrotmat(l)%mat(1:ns, 1:ns) = r(l, -l:l, -l:l)
406 : END DO
407 6 : DEALLOCATE (r)
408 : END IF
409 :
410 21188 : END SUBROUTINE calculate_rotmat
411 :
412 : ! **************************************************************************************************
413 : !> \brief ...
414 : !> \param l ...
415 : !> \param ma ...
416 : !> \param mb ...
417 : !> \return ...
418 : ! **************************************************************************************************
419 1506 : FUNCTION u_func(l, ma, mb) RESULT(u)
420 : INTEGER :: l, ma, mb
421 : REAL(KIND=dp) :: u
422 :
423 1506 : IF (ABS(mb) == l) THEN
424 324 : u = REAL((l + ma)*(l - ma), KIND=dp)/REAL(2*l*(2*l - 1), KIND=dp)
425 324 : u = SQRT(u)
426 1182 : ELSE IF (ABS(mb) < l) THEN
427 1182 : u = REAL((l + ma)*(l - ma), KIND=dp)/REAL((l + mb)*(l - mb), KIND=dp)
428 1182 : u = SQRT(u)
429 : ELSE
430 0 : CPABORT("Illegal m value")
431 : END IF
432 1506 : END FUNCTION u_func
433 :
434 : ! **************************************************************************************************
435 : !> \brief ...
436 : !> \param l ...
437 : !> \param ma ...
438 : !> \param mb ...
439 : !> \return ...
440 : ! **************************************************************************************************
441 1506 : FUNCTION v_func(l, ma, mb) RESULT(v)
442 : INTEGER :: l, ma, mb
443 : REAL(KIND=dp) :: v
444 :
445 : INTEGER :: a1, a2, dm0
446 :
447 1506 : dm0 = 0
448 1506 : IF (ma == 0) dm0 = 1
449 1506 : IF (ABS(mb) == l) THEN
450 324 : a1 = (1 + dm0)*(l + ABS(ma) - 1)*(l + ABS(ma))
451 324 : a2 = 2*l*(2*l - 1)
452 1182 : ELSE IF (ABS(mb) < l) THEN
453 1182 : a1 = (1 + dm0)*(l + ABS(ma) - 1)*(l + ABS(ma))
454 1182 : a2 = (l + mb)*(l - mb)
455 : ELSE
456 0 : CPABORT("Illegal m value")
457 : END IF
458 1506 : v = 0.5_dp*SQRT(REAL(a1, KIND=dp)/REAL(a2, KIND=dp))*REAL(1 - 2*dm0, KIND=dp)
459 1506 : END FUNCTION v_func
460 :
461 : ! **************************************************************************************************
462 : !> \brief ...
463 : !> \param l ...
464 : !> \param ma ...
465 : !> \param mb ...
466 : !> \return ...
467 : ! **************************************************************************************************
468 1506 : FUNCTION w_func(l, ma, mb) RESULT(w)
469 : INTEGER :: l, ma, mb
470 : REAL(KIND=dp) :: w
471 :
472 : INTEGER :: a1, a2, dm0
473 :
474 1506 : dm0 = 0
475 1506 : IF (ma == 0) dm0 = 1
476 1506 : IF (ABS(mb) == l) THEN
477 324 : a1 = (l - ABS(ma) - 1)*(l - ABS(ma))
478 324 : a2 = 2*l*(2*l - 1)
479 1182 : ELSE IF (ABS(mb) < l) THEN
480 1182 : a1 = (l - ABS(ma) - 1)*(l - ABS(ma))
481 1182 : a2 = (l + mb)*(l - mb)
482 : ELSE
483 0 : CPABORT("Illegal m value")
484 : END IF
485 1506 : w = -0.5_dp*SQRT(REAL(a1, KIND=dp)/REAL(a2, KIND=dp))*REAL(1 - dm0, KIND=dp)
486 1506 : END FUNCTION w_func
487 :
488 : ! **************************************************************************************************
489 : !> \brief ...
490 : !> \param i ...
491 : !> \param l ...
492 : !> \param mu ...
493 : !> \param m ...
494 : !> \param r1 ...
495 : !> \param r ...
496 : !> \param lr ...
497 : !> \return ...
498 : ! **************************************************************************************************
499 7206 : FUNCTION p_func(i, l, mu, m, r1, r, lr) RESULT(p)
500 : INTEGER :: i, l, mu, m
501 : REAL(KIND=dp), DIMENSION(-1:1, -1:1) :: r1
502 : INTEGER :: lr
503 : REAL(KIND=dp), DIMENSION(0:lr, -lr:lr, -lr:lr) :: r
504 : REAL(KIND=dp) :: p
505 :
506 7206 : IF (ABS(mu) > l - 1) THEN
507 : p = 0.0_dp
508 5586 : ELSE IF (m == l) THEN
509 594 : p = r1(i, 1)*r(l - 1, mu, l - 1) - r1(i, -1)*r(l - 1, mu, -l + 1)
510 4992 : ELSE IF (m == -l) THEN
511 594 : p = r1(i, 1)*r(l - 1, mu, -l + 1) + r1(i, -1)*r(l - 1, mu, l - 1)
512 4398 : ELSE IF (ABS(m) < l) THEN
513 4398 : p = r1(i, 0)*r(l - 1, mu, m)
514 : ELSE
515 0 : CPABORT("Illegal m value")
516 : END IF
517 7206 : 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 1506 : 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 1506 : IF (ma == 0) THEN
541 162 : u = p_func(0, l, 0, mb, r1, r, lr)
542 162 : v = p_func(1, l, 1, mb, r1, r, lr) + p_func(-1, l, -1, mb, r1, r, lr)
543 162 : w = 0.0_dp
544 1344 : ELSE IF (ma > 0) THEN
545 672 : dd = 0.0_dp
546 672 : IF (ma == 1) dd = 1.0_dp
547 672 : u = p_func(0, l, ma, mb, r1, r, lr)
548 672 : 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 672 : w = p_func(1, l, ma + 1, mb, r1, r, lr) + p_func(-1, l, -1 - ma, mb, r1, r, lr)
550 672 : ELSE IF (ma < 0) THEN
551 672 : dd = 0.0_dp
552 672 : IF (ma == -1) dd = 1.0_dp
553 672 : u = p_func(0, l, ma, mb, r1, r, lr)
554 672 : 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 672 : 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 1506 : END SUBROUTINE r_recursion
558 :
559 : ! **************************************************************************************************
560 : !> \brief Release rotation matrices
561 : !> \param orbrotmat ...
562 : ! **************************************************************************************************
563 23618 : SUBROUTINE release_rotmat(orbrotmat)
564 : TYPE(orbrotmat_type), DIMENSION(:), POINTER :: orbrotmat
565 :
566 : INTEGER :: i
567 :
568 23618 : IF (ASSOCIATED(orbrotmat)) THEN
569 114694 : DO i = LBOUND(orbrotmat, 1), UBOUND(orbrotmat, 1)
570 72318 : IF (ASSOCIATED(orbrotmat(i)%mat)) DEALLOCATE (orbrotmat(i)%mat)
571 : END DO
572 21188 : DEALLOCATE (orbrotmat)
573 : END IF
574 :
575 23618 : 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
|
