(git:d18deda)
Loading...
Searching...
No Matches
orbital_transformation_matrices.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 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! **************************************************************************************************
17MODULE 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
59CONTAINS
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 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!$ IF (omp_get_level() > 0) &
80!$ CPABORT("create_spherical_harmonics is not thread-safe")
81
82 IF (current_maxl > -1) THEN
83 CALL cp_abort(__location__, &
84 "Spherical harmonics are already allocated. "// &
85 "Use the init routine for an update")
86 END IF
87
88 IF (maxl < 0) THEN
89 CALL cp_abort(__location__, &
90 "A negative maximum angular momentum quantum "// &
91 "number is invalid")
92 END IF
93
94 ALLOCATE (orbtramat(0:maxl))
95 nc = ncoset(maxl)
96 ns = nsoset(maxl)
97
98 DO l = 0, maxl
99 nc = nco(l)
100 ns = nso(l)
101 ALLOCATE (orbtramat(l)%c2s(ns, nc))
102 orbtramat(l)%c2s = 0.0_dp
103 ALLOCATE (orbtramat(l)%s2c(ns, nc))
104 orbtramat(l)%s2c = 0.0_dp
105 ALLOCATE (orbtramat(l)%slm(ns, nc))
106 orbtramat(l)%slm = 0.0_dp
107 ALLOCATE (orbtramat(l)%slm_inv(ns, nc))
108 orbtramat(l)%slm_inv = 0.0_dp
109 END DO
110
111 ! Loop over all angular momentum quantum numbers l
112
113 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 DO lx = 0, l
119 DO ly = 0, l - lx
120 lz = l - lx - ly
121 ic = co(lx, ly, lz)
122 DO m = -l, l
123 is = l + m + 1
124 ma = abs(m)
125 j = lx + ly - ma
126 IF ((j >= 0) .AND. (modulo(j, 2) == 0)) THEN
127 j = j/2
128 s1 = 0.0_dp
129 DO i = 0, (l - ma)/2
130 s2 = 0.0_dp
131 DO k = 0, j
132 IF (((m < 0) .AND. (modulo(abs(ma - lx), 2) == 1)) .OR. &
133 ((m > 0) .AND. (modulo(abs(ma - lx), 2) == 0))) THEN
134 expo = (ma - lx + 2*k)/2
135 s = (-1.0_dp)**expo*sqrt(2.0_dp)
136 ELSE IF ((m == 0) .AND. (modulo(lx, 2) == 0)) THEN
137 expo = k - lx/2
138 s = (-1.0_dp)**expo
139 ELSE
140 s = 0.0_dp
141 END IF
142 s2 = s2 + binomial(j, k)*binomial(ma, lx - 2*k)*s
143 END DO
144 s1 = s1 + binomial(l, i)*binomial(i, j)* &
145 (-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 (2.0_dp**l*fac(l))
151 ELSE
152 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 DO lx1 = 0, l
162 DO ly1 = 0, l - lx1
163 lz1 = l - lx1 - ly1
164 ic1 = co(lx1, ly1, lz1)
165 s1 = sqrt((fac(lx1)*fac(ly1)*fac(lz1))/ &
166 (fac(2*lx1)*fac(2*ly1)*fac(2*lz1)))
167 DO lx2 = 0, l
168 DO ly2 = 0, l - lx2
169 lz2 = l - lx2 - ly2
170 ic2 = co(lx2, ly2, lz2)
171 lx = lx1 + lx2
172 ly = ly1 + ly2
173 lz = lz1 + lz2
174 IF ((modulo(lx, 2) == 0) .AND. &
175 (modulo(ly, 2) == 0) .AND. &
176 (modulo(lz, 2) == 0)) THEN
177 s2 = sqrt((fac(lx2)*fac(ly2)*fac(lz2))/ &
178 (fac(2*lx2)*fac(2*ly2)*fac(2*lz2)))
179 s = fac(lx)*fac(ly)*fac(lz)*s1*s2/ &
180 (fac(lx/2)*fac(ly/2)*fac(lz/2))
181 DO is = 1, nso(l)
182 orbtramat(l)%s2c(is, ic1) = orbtramat(l)%s2c(is, ic1) + &
183 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 s = sqrt(0.25_dp*dfac(2*l + 1)/pi)
194
195 DO lx = 0, l
196 DO ly = 0, l - lx
197 lz = l - lx - ly
198 ic = co(lx, ly, lz)
199 DO m = -l, l
200 is = l + m + 1
201 s1 = sqrt(dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1))
202 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 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 current_maxl = maxl
215
216 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 SUBROUTINE deallocate_spherical_harmonics()
225
226 INTEGER :: l
227
228!$ IF (omp_get_level() > 0) &
229!$ CPABORT("deallocate_spherical_harmonics is not thread-safe")
230
231 IF (current_maxl > -1) THEN
232 DO l = 0, SIZE(orbtramat, 1) - 1
233 DEALLOCATE (orbtramat(l)%c2s)
234 DEALLOCATE (orbtramat(l)%s2c)
235 DEALLOCATE (orbtramat(l)%slm)
236 DEALLOCATE (orbtramat(l)%slm_inv)
237 END DO
238 DEALLOCATE (orbtramat)
239 current_maxl = -1
240 END IF
241
242 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 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!$ IF (omp_get_level() > 0) &
263!$ CPABORT("init_spherical_harmonics is not thread-safe")
264
265 IF (maxl < 0) THEN
266 CALL cp_abort(__location__, &
267 "A negative maximum angular momentum quantum "// &
268 "number is invalid")
269 END IF
270
271 IF (maxl > current_maxl) THEN
272
273 CALL deallocate_spherical_harmonics()
274 CALL create_spherical_harmonics(maxl)
275
276 ! Print the spherical harmonics and the orbital transformation matrices
277
278 IF (output_unit > 0) THEN
279
280 DO l = 0, maxl
281
282 nc = nco(l)
283 ns = nso(l)
284
285 headline = "CARTESIAN ORBITAL TO SPHERICAL ORBITAL "// &
286 "TRANSFORMATION MATRIX"
287 CALL write_matrix(orbtramat(l)%c2s, l, output_unit, headline)
288
289 headline = "SPHERICAL ORBITAL TO CARTESIAN ORBITAL "// &
290 "TRANSFORMATION MATRIX"
291 CALL write_matrix(orbtramat(l)%s2c, l, output_unit, headline)
292
293 headline = "SPHERICAL HARMONICS"
294 CALL write_matrix(orbtramat(l)%slm, l, output_unit, headline)
295
296 headline = "INVERSE SPHERICAL HARMONICS"
297 CALL write_matrix(orbtramat(l)%slm_inv, l, output_unit, headline)
298
299 END DO
300
301 WRITE (unit=output_unit, fmt="(A)") ""
302
303 END IF
304
305 END IF
306
307 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 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 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 CALL release_rotmat(orbrotmat)
331
332 ALLOCATE (orbrotmat(0:lval))
333 DO l = 0, lval
334 ns = nso(l)
335 ALLOCATE (orbrotmat(l)%mat(ns, ns))
336 END DO
337
338 IF (lval >= 0) THEN
339 orbrotmat(0)%mat = 1.0_dp
340 END IF
341 IF (lval >= 1) THEN
342 t(1, 1:3) = rotmat(2, 1:3)
343 t(2, 1:3) = rotmat(3, 1:3)
344 t(3, 1:3) = rotmat(1, 1:3)
345 r1(-1:1, -1) = t(1:3, 2)
346 r1(-1:1, 0) = t(1:3, 3)
347 r1(-1:1, 1) = t(1:3, 1)
348 orbrotmat(1)%mat(1:3, 1:3) = r1(-1:1, -1:1)
349 END IF
350 IF (lval >= 2) THEN
351 s3 = sqrt(3.0_dp)
352 ! Table 4
353 r2(0, 0) = (3.0_dp*r1(0, 0)**2 - 1.0_dp)*0.5_dp
354 r2(0, 1) = s3*r1(0, 1)*r1(0, 0)
355 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 r2(0, -2) = s3*r1(0, 1)*r1(0, -1)
358 !
359 r2(1, 0) = s3*r1(1, 0)*r1(0, 0)
360 r2(1, 1) = r1(1, 1)*r1(0, 0) + r1(1, 0)*r1(0, 1)
361 r2(1, -1) = r1(1, -1)*r1(0, 0) + r1(1, 0)*r1(0, -1)
362 r2(1, 2) = r1(1, 1)*r1(0, 1) - r1(1, -1)*r1(0, -1)
363 r2(1, -2) = r1(1, 1)*r1(0, -1) + r1(1, -1)*r1(0, 1)
364 !
365 r2(-1, 0) = s3*r1(-1, 0)*r1(0, 0)
366 r2(-1, 1) = r1(-1, 1)*r1(0, 0) + r1(-1, 0)*r1(0, 1)
367 r2(-1, -1) = r1(-1, -1)*r1(0, 0) + r1(-1, 0)*r1(0, -1)
368 r2(-1, 2) = r1(-1, 1)*r1(0, 1) - r1(-1, -1)*r1(0, -1)
369 r2(-1, -2) = r1(-1, 1)*r1(0, -1) + r1(-1, -1)*r1(0, 1)
370 !
371 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 r2(-2, 0) = s3*r1(1, 0)*r1(-1, 0)
378 r2(2, 1) = r1(1, 1)*r1(-1, 0) + r1(1, 0)*r1(-1, 1)
379 r2(2, -1) = r1(1, -1)*r1(-1, 0) + r1(1, 0)*r1(-1, -1)
380 r2(2, 2) = r1(1, 1)*r1(-1, 1) - r1(1, -1)*r1(-1, -1)
381 r2(2, -2) = r1(1, 1)*r1(-1, -1) + r1(1, -1)*r1(-1, 1)
382 !
383 orbrotmat(2)%mat(1:5, 1:5) = r2(-2:2, -2:2)
384 END IF
385 IF (lval >= 3) THEN
386 ! General recursion
387 ALLOCATE (r(0:lval, -lval:lval, -lval:lval))
388 r = 0.0_dp
389 r(0, 0, 0) = 1.0_dp
390 r(1, -1:1, -1:1) = r1(-1:1, -1:1)
391 r(2, -2:2, -2:2) = r2(-2:2, -2:2)
392 DO l = 3, lval
393 DO m1 = -l, l
394 DO m2 = -l, l
395 u = u_func(l, m1, m2)
396 v = v_func(l, m1, m2)
397 w = w_func(l, m1, m2)
398 CALL r_recursion(l, m1, m2, r1, r, lval, uf, vf, wf)
399 r(l, m1, m2) = u*uf + v*vf + w*wf
400 END DO
401 END DO
402 END DO
403 DO l = 3, lval
404 ns = nso(l)
405 orbrotmat(l)%mat(1:ns, 1:ns) = r(l, -l:l, -l:l)
406 END DO
407 DEALLOCATE (r)
408 END IF
409
410 END SUBROUTINE calculate_rotmat
411
412! **************************************************************************************************
413!> \brief ...
414!> \param l ...
415!> \param ma ...
416!> \param mb ...
417!> \return ...
418! **************************************************************************************************
419 FUNCTION u_func(l, ma, mb) RESULT(u)
420 INTEGER :: l, ma, mb
421 REAL(kind=dp) :: u
422
423 IF (abs(mb) == l) THEN
424 u = real((l + ma)*(l - ma), kind=dp)/real(2*l*(2*l - 1), kind=dp)
425 u = sqrt(u)
426 ELSE IF (abs(mb) < l) THEN
427 u = real((l + ma)*(l - ma), kind=dp)/real((l + mb)*(l - mb), kind=dp)
428 u = sqrt(u)
429 ELSE
430 cpabort("Illegal m value")
431 END IF
432 END FUNCTION u_func
433
434! **************************************************************************************************
435!> \brief ...
436!> \param l ...
437!> \param ma ...
438!> \param mb ...
439!> \return ...
440! **************************************************************************************************
441 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 dm0 = 0
448 IF (ma == 0) dm0 = 1
449 IF (abs(mb) == l) THEN
450 a1 = (1 + dm0)*(l + abs(ma) - 1)*(l + abs(ma))
451 a2 = 2*l*(2*l - 1)
452 ELSE IF (abs(mb) < l) THEN
453 a1 = (1 + dm0)*(l + abs(ma) - 1)*(l + abs(ma))
454 a2 = (l + mb)*(l - mb)
455 ELSE
456 cpabort("Illegal m value")
457 END IF
458 v = 0.5_dp*sqrt(real(a1, kind=dp)/real(a2, kind=dp))*real(1 - 2*dm0, kind=dp)
459 END FUNCTION v_func
460
461! **************************************************************************************************
462!> \brief ...
463!> \param l ...
464!> \param ma ...
465!> \param mb ...
466!> \return ...
467! **************************************************************************************************
468 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 dm0 = 0
475 IF (ma == 0) dm0 = 1
476 IF (abs(mb) == l) THEN
477 a1 = (l - abs(ma) - 1)*(l - abs(ma))
478 a2 = 2*l*(2*l - 1)
479 ELSE IF (abs(mb) < l) THEN
480 a1 = (l - abs(ma) - 1)*(l - abs(ma))
481 a2 = (l + mb)*(l - mb)
482 ELSE
483 cpabort("Illegal m value")
484 END IF
485 w = -0.5_dp*sqrt(real(a1, kind=dp)/real(a2, kind=dp))*real(1 - dm0, kind=dp)
486 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 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 IF (m == l) THEN
507 p = r1(i, 1)*r(l - 1, mu, m) - r1(i, -1)*r(l - 1, mu, -m)
508 ELSE IF (m == -l) THEN
509 p = r1(i, 1)*r(l - 1, mu, m) + r1(i, -1)*r(l - 1, mu, -m)
510 ELSE IF (abs(m) < l) THEN
511 p = r1(i, 0)*r(l - 1, mu, m)
512 ELSE
513 cpabort("Illegal m value")
514 END IF
515 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 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 IF (ma == 0) THEN
539 u = p_func(0, l, 0, mb, r1, r, lr)
540 v = p_func(1, l, 1, mb, r1, r, lr) + p_func(-1, l, -1, mb, r1, r, lr)
541 w = 0.0_dp
542 ELSE IF (ma > 0) THEN
543 dd = 1.0_dp
544 IF (ma == 1) dd = 1.0_dp
545 u = p_func(0, l, ma, mb, r1, r, lr)
546 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 w = p_func(1, l, ma + 1, mb, r1, r, lr) + p_func(-1, l, -1 - ma, mb, r1, r, lr)
548 ELSE IF (ma < 0) THEN
549 dd = 1.0_dp
550 IF (ma == -1) dd = 1.0_dp
551 u = p_func(0, l, ma, mb, r1, r, lr)
552 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 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 END SUBROUTINE
556
557! **************************************************************************************************
558!> \brief Release rotation matrices
559!> \param orbrotmat ...
560! **************************************************************************************************
561 SUBROUTINE release_rotmat(orbrotmat)
562 TYPE(orbrotmat_type), DIMENSION(:), POINTER :: orbrotmat
563
564 INTEGER :: i
565
566 IF (ASSOCIATED(orbrotmat)) THEN
567 DO i = lbound(orbrotmat, 1), ubound(orbrotmat, 1)
568 IF (ASSOCIATED(orbrotmat(i)%mat)) DEALLOCATE (orbrotmat(i)%mat)
569 END DO
570 DEALLOCATE (orbrotmat)
571 END IF
572
573 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 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 WRITE (unit=lunit, fmt="(/,/,T2,A)") trim(headline)
604
605 ! Write the matrix in the defined format
606
607 nc = nco(l)
608
609 DO ic = 1, nc, 3
610 from = ic
611 to = min(nc, from + 2)
612 i = 1
613 DO lx = l, 0, -1
614 DO ly = l - lx, 0, -1
615 lz = l - lx - ly
616 jc = co(lx, ly, lz)
617 IF ((jc >= from) .AND. (jc <= to)) THEN
618 symbol = cgf_symbol(1, (/lx, ly, lz/))
619 WRITE (unit=string(i:), fmt="(A18)") trim(symbol(3:12))
620 i = i + 18
621 END IF
622 END DO
623 END DO
624 WRITE (unit=lunit, fmt="(/,T8,A)") trim(string)
625 symbol = ""
626 DO m = -l, l
627 is = l + m + 1
628 symbol = sgf_symbol(1, l, m)
629 WRITE (unit=lunit, fmt="(T4,A4,3(1X,F17.12))") &
630 symbol(3:6), (matrix(is, jc), jc=from, to)
631 END DO
632 END DO
633
634 END SUBROUTINE write_matrix
635
636END MODULE orbital_transformation_matrices
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:120
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::dfac
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
Definition mathconstants.F:61
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::binomial
elemental real(kind=dp) function, public binomial(n, k)
The binomial coefficient n over k for 0 <= k <= n is calculated, otherwise zero is returned.
Definition mathlib.F:206
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::current_maxl
integer, save, public current_maxl
Definition orbital_pointers.F:40
orbital_pointers::co
integer, dimension(:, :, :), allocatable, public co
Definition orbital_pointers.F:45
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition orbital_pointers.F:42
orbital_pointers::nsoset
integer, dimension(:), allocatable, public nsoset
Definition orbital_pointers.F:42
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition orbital_pointers.F:42
orbital_symbols
orbital_symbols
Definition orbital_symbols.F:14
orbital_symbols::cgf_symbol
character(len=12) function, public cgf_symbol(n, lxyz)
Build a Cartesian orbital symbol (orbital labels for printing).
Definition orbital_symbols.F:44
orbital_symbols::sgf_symbol
character(len=6) function, public sgf_symbol(n, l, m)
Build a spherical orbital symbol (orbital labels for printing).
Definition orbital_symbols.F:107
orbital_transformation_matrices
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
Definition orbital_transformation_matrices.F:17
orbital_transformation_matrices::init_spherical_harmonics
subroutine, public init_spherical_harmonics(maxl, output_unit)
Initialize or update the orbital transformation matrices.
Definition orbital_transformation_matrices.F:255
orbital_transformation_matrices::calculate_rotmat
subroutine, public calculate_rotmat(orbrotmat, rotmat, lval)
Calculate rotation matrices for spherical harmonics up to value l Joseph Ivanic and Klaus Ruedenberg ...
Definition orbital_transformation_matrices.F:319
orbital_transformation_matrices::release_rotmat
subroutine, public release_rotmat(orbrotmat)
Release rotation matrices.
Definition orbital_transformation_matrices.F:562
orbital_transformation_matrices::orbtramat
type(orbtramat_type), dimension(:), pointer, public orbtramat
Definition orbital_transformation_matrices.F:47
orbital_transformation_matrices::deallocate_spherical_harmonics
subroutine, public deallocate_spherical_harmonics()
Deallocate the orbital transformation matrices.
Definition orbital_transformation_matrices.F:225
orbital_transformation_matrices::orbrotmat_type
Definition orbital_transformation_matrices.F:49