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 6324 : 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 6324 : !$ IF (omp_get_level() > 0) &
80 0 : !$ CPABORT("create_spherical_harmonics is not thread-safe")
81 :
82 6324 : 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 6324 : 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 52550 : ALLOCATE (orbtramat(0:maxl))
95 6324 : nc = ncoset(maxl)
96 6324 : ns = nsoset(maxl)
97 :
98 39902 : DO l = 0, maxl
99 33578 : nc = nco(l)
100 33578 : ns = nso(l)
101 134312 : ALLOCATE (orbtramat(l)%c2s(ns, nc))
102 2835312 : orbtramat(l)%c2s = 0.0_dp
103 100734 : ALLOCATE (orbtramat(l)%s2c(ns, nc))
104 2835312 : orbtramat(l)%s2c = 0.0_dp
105 100734 : ALLOCATE (orbtramat(l)%slm(ns, nc))
106 2835312 : orbtramat(l)%slm = 0.0_dp
107 100734 : ALLOCATE (orbtramat(l)%slm_inv(ns, nc))
108 2841636 : orbtramat(l)%slm_inv = 0.0_dp
109 : END DO
110 :
111 : ! Loop over all angular momentum quantum numbers l
112 :
113 39902 : 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 144611 : DO lx = 0, l
119 435511 : DO ly = 0, l - lx
120 290900 : lz = l - lx - ly
121 290900 : ic = co(lx, ly, lz)
122 2912767 : DO m = -l, l
123 2510834 : is = l + m + 1
124 2510834 : ma = ABS(m)
125 2510834 : j = lx + ly - ma
126 2801734 : IF ((j >= 0) .AND. (MODULO(j, 2) == 0)) THEN
127 1030878 : j = j/2
128 1030878 : s1 = 0.0_dp
129 3042636 : DO i = 0, (l - ma)/2
130 2011758 : s2 = 0.0_dp
131 5843828 : DO k = 0, j
132 3143113 : IF (((m < 0) .AND. (MODULO(ABS(ma - lx), 2) == 1)) .OR. &
133 3832070 : ((m > 0) .AND. (MODULO(ABS(ma - lx), 2) == 0))) THEN
134 1469866 : expo = (ma - lx + 2*k)/2
135 1469866 : s = (-1.0_dp)**expo*SQRT(2.0_dp)
136 2362204 : ELSE IF ((m == 0) .AND. (MODULO(lx, 2) == 0)) THEN
137 577194 : expo = k - lx/2
138 577194 : s = (-1.0_dp)**expo
139 : ELSE
140 : s = 0.0_dp
141 : END IF
142 5843828 : s2 = s2 + binomial(j, k)*binomial(ma, lx - 2*k)*s
143 : END DO
144 : s1 = s1 + binomial(l, i)*binomial(i, j)* &
145 3042636 : (-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 1030878 : (2.0_dp**l*fac(l))
151 : ELSE
152 1479956 : 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 144611 : DO lx1 = 0, l
162 435511 : DO ly1 = 0, l - lx1
163 290900 : lz1 = l - lx1 - ly1
164 290900 : ic1 = co(lx1, ly1, lz1)
165 : s1 = SQRT((fac(lx1)*fac(ly1)*fac(lz1))/ &
166 290900 : (fac(2*lx1)*fac(2*ly1)*fac(2*lz1)))
167 1802800 : DO lx2 = 0, l
168 6199989 : DO ly2 = 0, l - lx2
169 4508222 : lz2 = l - lx2 - ly2
170 4508222 : ic2 = co(lx2, ly2, lz2)
171 4508222 : lx = lx1 + lx2
172 4508222 : ly = ly1 + ly2
173 4508222 : lz = lz1 + lz2
174 : IF ((MODULO(lx, 2) == 0) .AND. &
175 4508222 : (MODULO(ly, 2) == 0) .AND. &
176 1400867 : (MODULO(lz, 2) == 0)) THEN
177 : s2 = SQRT((fac(lx2)*fac(ly2)*fac(lz2))/ &
178 1236152 : (fac(2*lx2)*fac(2*ly2)*fac(2*lz2)))
179 : s = fac(lx)*fac(ly)*fac(lz)*s1*s2/ &
180 1236152 : (fac(lx/2)*fac(ly/2)*fac(lz/2))
181 14374750 : DO is = 1, nso(l)
182 : orbtramat(l)%s2c(is, ic1) = orbtramat(l)%s2c(is, ic1) + &
183 14374750 : 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 33578 : s = SQRT(0.25_dp*dfac(2*l + 1)/pi)
194 :
195 150935 : DO lx = 0, l
196 435511 : DO ly = 0, l - lx
197 290900 : lz = l - lx - ly
198 290900 : ic = co(lx, ly, lz)
199 2912767 : DO m = -l, l
200 2510834 : is = l + m + 1
201 2510834 : s1 = SQRT(dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1))
202 2510834 : 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 2801734 : 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 6324 : current_maxl = maxl
215 :
216 6324 : 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 15254 : SUBROUTINE deallocate_spherical_harmonics()
225 :
226 : INTEGER :: l
227 :
228 15254 : !$ IF (omp_get_level() > 0) &
229 0 : !$ CPABORT("deallocate_spherical_harmonics is not thread-safe")
230 :
231 15254 : IF (current_maxl > -1) THEN
232 39902 : DO l = 0, SIZE(orbtramat, 1) - 1
233 33578 : DEALLOCATE (orbtramat(l)%c2s)
234 33578 : DEALLOCATE (orbtramat(l)%s2c)
235 33578 : DEALLOCATE (orbtramat(l)%slm)
236 39902 : DEALLOCATE (orbtramat(l)%slm_inv)
237 : END DO
238 6324 : DEALLOCATE (orbtramat)
239 6324 : current_maxl = -1
240 : END IF
241 :
242 15254 : 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 23718 : 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 23718 : !$ IF (omp_get_level() > 0) &
263 0 : !$ CPABORT("init_spherical_harmonics is not thread-safe")
264 :
265 23718 : 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 23718 : IF (maxl > current_maxl) THEN
272 :
273 6324 : CALL deallocate_spherical_harmonics()
274 6324 : CALL create_spherical_harmonics(maxl)
275 :
276 : ! Print the spherical harmonics and the orbital transformation matrices
277 :
278 6324 : 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 23718 : 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 0 : CALL release_rotmat(orbrotmat)
331 :
332 0 : ALLOCATE (orbrotmat(0:lval))
333 0 : DO l = 0, lval
334 0 : ns = nso(l)
335 0 : ALLOCATE (orbrotmat(l)%mat(ns, ns))
336 : END DO
337 :
338 0 : IF (lval >= 0) THEN
339 0 : orbrotmat(0)%mat = 1.0_dp
340 : END IF
341 0 : IF (lval >= 1) THEN
342 0 : t(1, 1:3) = rotmat(2, 1:3)
343 0 : t(2, 1:3) = rotmat(3, 1:3)
344 0 : t(3, 1:3) = rotmat(1, 1:3)
345 0 : r1(-1:1, -1) = t(1:3, 2)
346 0 : r1(-1:1, 0) = t(1:3, 3)
347 0 : r1(-1:1, 1) = t(1:3, 1)
348 0 : orbrotmat(1)%mat(1:3, 1:3) = r1(-1:1, -1:1)
349 : END IF
350 0 : IF (lval >= 2) THEN
351 0 : s3 = SQRT(3.0_dp)
352 : ! Table 4
353 0 : r2(0, 0) = (3.0_dp*r1(0, 0)**2 - 1.0_dp)*0.5_dp
354 0 : r2(0, 1) = s3*r1(0, 1)*r1(0, 0)
355 0 : r2(0, -1) = s3*r1(0, -1)*r1(0, 0)
356 : r2(0, -2) = s3*(r1(0, 1)**2 - r1(0, -1)**2)*0.5_dp
357 0 : r2(0, -2) = s3*r1(0, 1)*r1(0, -1)
358 : !
359 0 : r2(1, 0) = s3*r1(1, 0)*r1(0, 0)
360 0 : r2(1, 1) = r1(1, 1)*r1(0, 0) + r1(1, 0)*r1(0, 1)
361 0 : r2(1, -1) = r1(1, -1)*r1(0, 0) + r1(1, 0)*r1(0, -1)
362 0 : r2(1, 2) = r1(1, 1)*r1(0, 1) - r1(1, -1)*r1(0, -1)
363 0 : r2(1, -2) = r1(1, 1)*r1(0, -1) + r1(1, -1)*r1(0, 1)
364 : !
365 0 : r2(-1, 0) = s3*r1(-1, 0)*r1(0, 0)
366 0 : r2(-1, 1) = r1(-1, 1)*r1(0, 0) + r1(-1, 0)*r1(0, 1)
367 0 : r2(-1, -1) = r1(-1, -1)*r1(0, 0) + r1(-1, 0)*r1(0, -1)
368 0 : r2(-1, 2) = r1(-1, 1)*r1(0, 1) - r1(-1, -1)*r1(0, -1)
369 0 : r2(-1, -2) = r1(-1, 1)*r1(0, -1) + r1(-1, -1)*r1(0, 1)
370 : !
371 0 : r2(2, 0) = s3*(r1(1, 0)**2 - r1(-1, 0)**2)*0.5_dp
372 : r2(2, 1) = r1(1, 1)*r1(1, 0) - r1(-1, 1)*r1(-1, 0)
373 : r2(2, -1) = r1(1, -1)*r1(1, 0) - r1(-1, -1)*r1(-1, 0)
374 : r2(2, 2) = (r1(1, 1)**2 - r1(1, -1)**2 - r1(-1, 1)**2 + r1(-1, -1)**2)*0.5_dp
375 : r2(2, -2) = r1(1, 1)*r1(1, -1) - r1(-1, 1)*r1(-1, -1)
376 : !
377 0 : r2(-2, 0) = s3*r1(1, 0)*r1(-1, 0)
378 0 : r2(2, 1) = r1(1, 1)*r1(-1, 0) + r1(1, 0)*r1(-1, 1)
379 0 : r2(2, -1) = r1(1, -1)*r1(-1, 0) + r1(1, 0)*r1(-1, -1)
380 0 : r2(2, 2) = r1(1, 1)*r1(-1, 1) - r1(1, -1)*r1(-1, -1)
381 0 : r2(2, -2) = r1(1, 1)*r1(-1, -1) + r1(1, -1)*r1(-1, 1)
382 : !
383 0 : orbrotmat(2)%mat(1:5, 1:5) = r2(-2:2, -2:2)
384 : END IF
385 0 : IF (lval >= 3) THEN
386 : ! General recursion
387 0 : ALLOCATE (r(0:lval, -lval:lval, -lval:lval))
388 0 : r = 0.0_dp
389 0 : r(0, 0, 0) = 1.0_dp
390 0 : r(1, -1:1, -1:1) = r1(-1:1, -1:1)
391 0 : r(2, -2:2, -2:2) = r2(-2:2, -2:2)
392 0 : DO l = 3, lval
393 0 : DO m1 = -l, l
394 0 : DO m2 = -l, l
395 0 : u = u_func(l, m1, m2)
396 0 : v = v_func(l, m1, m2)
397 0 : w = w_func(l, m1, m2)
398 0 : CALL r_recursion(l, m1, m2, r1, r, lval, uf, vf, wf)
399 0 : r(l, m1, m2) = u*uf + v*vf + w*wf
400 : END DO
401 : END DO
402 : END DO
403 0 : DO l = 3, lval
404 0 : ns = nso(l)
405 0 : orbrotmat(l)%mat(1:ns, 1:ns) = r(l, -l:l, -l:l)
406 : END DO
407 0 : DEALLOCATE (r)
408 : END IF
409 :
410 0 : END SUBROUTINE calculate_rotmat
411 :
412 : ! **************************************************************************************************
413 : !> \brief ...
414 : !> \param l ...
415 : !> \param ma ...
416 : !> \param mb ...
417 : !> \return ...
418 : ! **************************************************************************************************
419 0 : FUNCTION u_func(l, ma, mb) RESULT(u)
420 : INTEGER :: l, ma, mb
421 : REAL(KIND=dp) :: u
422 :
423 0 : IF (ABS(mb) == l) THEN
424 0 : u = REAL((l + ma)*(l - ma), KIND=dp)/REAL(2*l*(2*l - 1), KIND=dp)
425 0 : u = SQRT(u)
426 0 : ELSE IF (ABS(mb) < l) THEN
427 0 : u = REAL((l + ma)*(l - ma), KIND=dp)/REAL((l + mb)*(l - mb), KIND=dp)
428 0 : u = SQRT(u)
429 : ELSE
430 0 : CPABORT("Illegal m value")
431 : END IF
432 0 : END FUNCTION u_func
433 :
434 : ! **************************************************************************************************
435 : !> \brief ...
436 : !> \param l ...
437 : !> \param ma ...
438 : !> \param mb ...
439 : !> \return ...
440 : ! **************************************************************************************************
441 0 : 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 0 : dm0 = 0
448 0 : IF (ma == 0) dm0 = 1
449 0 : IF (ABS(mb) == l) THEN
450 0 : a1 = (1 + dm0)*(l + ABS(ma) - 1)*(l + ABS(ma))
451 0 : a2 = 2*l*(2*l - 1)
452 0 : ELSE IF (ABS(mb) < l) THEN
453 0 : a1 = (1 + dm0)*(l + ABS(ma) - 1)*(l + ABS(ma))
454 0 : a2 = (l + mb)*(l - mb)
455 : ELSE
456 0 : CPABORT("Illegal m value")
457 : END IF
458 0 : v = 0.5_dp*SQRT(REAL(a1, KIND=dp)/REAL(a2, KIND=dp))*REAL(1 - 2*dm0, KIND=dp)
459 0 : END FUNCTION v_func
460 :
461 : ! **************************************************************************************************
462 : !> \brief ...
463 : !> \param l ...
464 : !> \param ma ...
465 : !> \param mb ...
466 : !> \return ...
467 : ! **************************************************************************************************
468 0 : 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 0 : dm0 = 0
475 0 : IF (ma == 0) dm0 = 1
476 0 : IF (ABS(mb) == l) THEN
477 0 : a1 = (l - ABS(ma) - 1)*(l - ABS(ma))
478 0 : a2 = 2*l*(2*l - 1)
479 0 : ELSE IF (ABS(mb) < l) THEN
480 0 : a1 = (l - ABS(ma) - 1)*(l - ABS(ma))
481 0 : a2 = (l + mb)*(l - mb)
482 : ELSE
483 0 : CPABORT("Illegal m value")
484 : END IF
485 0 : w = -0.5_dp*SQRT(REAL(a1, KIND=dp)/REAL(a2, KIND=dp))*REAL(1 - dm0, KIND=dp)
486 0 : 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 0 : 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 0 : IF (m == l) THEN
507 0 : p = r1(i, 1)*r(l - 1, mu, m) - r1(i, -1)*r(l - 1, mu, -m)
508 0 : ELSE 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 (ABS(m) < l) THEN
511 0 : p = r1(i, 0)*r(l - 1, mu, m)
512 : ELSE
513 0 : CPABORT("Illegal m value")
514 : END IF
515 0 : END FUNCTION p_func
516 :
517 : ! **************************************************************************************************
518 : !> \brief ...
519 : !> \param l ...
520 : !> \param ma ...
521 : !> \param mb ...
522 : !> \param r1 ...
523 : !> \param r ...
524 : !> \param lr ...
525 : !> \param u ...
526 : !> \param v ...
527 : !> \param w ...
528 : ! **************************************************************************************************
529 0 : SUBROUTINE r_recursion(l, ma, mb, r1, r, lr, u, v, w)
530 : INTEGER :: l, ma, mb
531 : REAL(KIND=dp), DIMENSION(-1:1, -1:1) :: r1
532 : INTEGER :: lr
533 : REAL(KIND=dp), DIMENSION(0:lr, -lr:lr, -lr:lr) :: r
534 : REAL(KIND=dp) :: u, v, w
535 :
536 : REAL(KIND=dp) :: dd
537 :
538 0 : IF (ma == 0) THEN
539 0 : u = p_func(0, l, 0, mb, r1, r, lr)
540 0 : v = p_func(1, l, 1, mb, r1, r, lr) + p_func(-1, l, -1, mb, r1, r, lr)
541 0 : w = 0.0_dp
542 0 : ELSE IF (ma > 0) THEN
543 : dd = 1.0_dp
544 : IF (ma == 1) dd = 1.0_dp
545 0 : u = p_func(0, l, ma, mb, r1, r, lr)
546 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)
547 0 : w = p_func(1, l, ma + 1, mb, r1, r, lr) + p_func(-1, l, -1 - ma, mb, r1, r, lr)
548 0 : ELSE IF (ma < 0) THEN
549 : dd = 1.0_dp
550 : IF (ma == -1) dd = 1.0_dp
551 0 : u = p_func(0, l, ma, mb, r1, r, lr)
552 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)
553 0 : w = p_func(1, l, ma - 1, mb, r1, r, lr) - p_func(-1, l, -ma + 1, mb, r1, r, lr)
554 : END IF
555 0 : END SUBROUTINE
556 :
557 : ! **************************************************************************************************
558 : !> \brief Release rotation matrices
559 : !> \param orbrotmat ...
560 : ! **************************************************************************************************
561 48 : SUBROUTINE release_rotmat(orbrotmat)
562 : TYPE(orbrotmat_type), DIMENSION(:), POINTER :: orbrotmat
563 :
564 : INTEGER :: i
565 :
566 48 : IF (ASSOCIATED(orbrotmat)) THEN
567 0 : DO i = LBOUND(orbrotmat, 1), UBOUND(orbrotmat, 1)
568 0 : IF (ASSOCIATED(orbrotmat(i)%mat)) DEALLOCATE (orbrotmat(i)%mat)
569 : END DO
570 0 : DEALLOCATE (orbrotmat)
571 : END IF
572 :
573 48 : END SUBROUTINE release_rotmat
574 :
575 : ! **************************************************************************************************
576 : !> \brief Print a matrix to the logical unit "lunit".
577 : !> \param matrix ...
578 : !> \param l ...
579 : !> \param lunit ...
580 : !> \param headline ...
581 : !> \date 07.06.2000
582 : !> \par Variables
583 : !> - matrix : Matrix to be printed.
584 : !> - l : Angular momentum quantum number
585 : !> - lunit : Logical unit number.
586 : !> - headline: Headline of the matrix.
587 : !> \author MK
588 : !> \version 1.0
589 : ! **************************************************************************************************
590 12 : SUBROUTINE write_matrix(matrix, l, lunit, headline)
591 :
592 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
593 : INTEGER, INTENT(IN) :: l, lunit
594 : CHARACTER(LEN=*), INTENT(IN) :: headline
595 :
596 : CHARACTER(12) :: symbol
597 : CHARACTER(LEN=78) :: string
598 : INTEGER :: from, i, ic, is, jc, lx, ly, lz, m, nc, &
599 : to
600 :
601 : ! Write headline
602 :
603 12 : WRITE (UNIT=lunit, FMT="(/,/,T2,A)") TRIM(headline)
604 :
605 : ! Write the matrix in the defined format
606 :
607 12 : nc = nco(l)
608 :
609 28 : DO ic = 1, nc, 3
610 16 : from = ic
611 16 : to = MIN(nc, from + 2)
612 16 : i = 1
613 52 : DO lx = l, 0, -1
614 116 : DO ly = l - lx, 0, -1
615 64 : lz = l - lx - ly
616 64 : jc = co(lx, ly, lz)
617 100 : IF ((jc >= from) .AND. (jc <= to)) THEN
618 160 : symbol = cgf_symbol(1, (/lx, ly, lz/))
619 40 : WRITE (UNIT=string(i:), FMT="(A18)") TRIM(symbol(3:12))
620 40 : i = i + 18
621 : END IF
622 : END DO
623 : END DO
624 16 : WRITE (UNIT=lunit, FMT="(/,T8,A)") TRIM(string)
625 16 : symbol = ""
626 84 : DO m = -l, l
627 56 : is = l + m + 1
628 56 : symbol = sgf_symbol(1, l, m)
629 : WRITE (UNIT=lunit, FMT="(T4,A4,3(1X,F17.12))") &
630 72 : symbol(3:6), (matrix(is, jc), jc=from, to)
631 : END DO
632 : END DO
633 :
634 12 : END SUBROUTINE write_matrix
635 :
636 0 : END MODULE orbital_transformation_matrices
|
