(git:374b731)
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-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! **************************************************************************************************
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 INTEGER, SAVE :: current_maxl = -1
50
51 ! Public variables
52
53 PUBLIC :: orbtramat
54
55 ! Public subroutines
56
57 PUBLIC :: deallocate_spherical_harmonics, &
58 init_spherical_harmonics
59
60CONTAINS
61
62! **************************************************************************************************
63!> \brief Allocate and initialize the orbital transformation matrices for
64!> the transformation and for the backtransformation of orbitals
65!> from the Cartesian representation to the spherical representation.
66!> \param maxl ...
67!> \date 20.05.2004
68!> \author MK
69!> \version 1.0
70! **************************************************************************************************
71 SUBROUTINE create_spherical_harmonics(maxl)
72
73 INTEGER, INTENT(IN) :: maxl
74
75 INTEGER :: expo, i, ic, ic1, ic2, is, j, k, l, lx, &
76 lx1, lx2, ly, ly1, ly2, lz, lz1, lz2, &
77 m, ma, nc, ns
78 REAL(kind=dp) :: s, s1, s2
79
80!$ IF (omp_get_level() > 0) &
81!$ CPABORT("create_spherical_harmonics is not thread-safe")
82
83 IF (current_maxl > -1) THEN
84 CALL cp_abort(__location__, &
85 "Spherical harmonics are already allocated. "// &
86 "Use the init routine for an update")
87 END IF
88
89 IF (maxl < 0) THEN
90 CALL cp_abort(__location__, &
91 "A negative maximum angular momentum quantum "// &
92 "number is invalid")
93 END IF
94
95 ALLOCATE (orbtramat(0:maxl))
96 nc = ncoset(maxl)
97 ns = nsoset(maxl)
98
99 DO l = 0, maxl
100 nc = nco(l)
101 ns = nso(l)
102 ALLOCATE (orbtramat(l)%c2s(ns, nc))
103 orbtramat(l)%c2s = 0.0_dp
104 ALLOCATE (orbtramat(l)%s2c(ns, nc))
105 orbtramat(l)%s2c = 0.0_dp
106 ALLOCATE (orbtramat(l)%slm(ns, nc))
107 orbtramat(l)%slm = 0.0_dp
108 ALLOCATE (orbtramat(l)%slm_inv(ns, nc))
109 orbtramat(l)%slm_inv = 0.0_dp
110 END DO
111
112 ! Loop over all angular momentum quantum numbers l
113
114 DO l = 0, maxl
115
116 ! Build the orbital transformation matrix for the transformation
117 ! from Cartesian to spherical orbitals (c2s, formula 15)
118
119 DO lx = 0, l
120 DO ly = 0, l - lx
121 lz = l - lx - ly
122 ic = co(lx, ly, lz)
123 DO m = -l, l
124 is = l + m + 1
125 ma = abs(m)
126 j = lx + ly - ma
127 IF ((j >= 0) .AND. (modulo(j, 2) == 0)) THEN
128 j = j/2
129 s1 = 0.0_dp
130 DO i = 0, (l - ma)/2
131 s2 = 0.0_dp
132 DO k = 0, j
133 IF (((m < 0) .AND. (modulo(abs(ma - lx), 2) == 1)) .OR. &
134 ((m > 0) .AND. (modulo(abs(ma - lx), 2) == 0))) THEN
135 expo = (ma - lx + 2*k)/2
136 s = (-1.0_dp)**expo*sqrt(2.0_dp)
137 ELSE IF ((m == 0) .AND. (modulo(lx, 2) == 0)) THEN
138 expo = k - lx/2
139 s = (-1.0_dp)**expo
140 ELSE
141 s = 0.0_dp
142 END IF
143 s2 = s2 + binomial(j, k)*binomial(ma, lx - 2*k)*s
144 END DO
145 s1 = s1 + binomial(l, i)*binomial(i, j)* &
146 (-1.0_dp)**i*fac(2*l - 2*i)/fac(l - ma - 2*i)*s2
147 END DO
148 orbtramat(l)%c2s(is, ic) = &
149 sqrt((fac(2*lx)*fac(2*ly)*fac(2*lz)*fac(l)*fac(l - ma))/ &
150 (fac(lx)*fac(ly)*fac(lz)*fac(2*l)*fac(l + ma)))*s1/ &
151 (2.0_dp**l*fac(l))
152 ELSE
153 orbtramat(l)%c2s(is, ic) = 0.0_dp
154 END IF
155 END DO
156 END DO
157 END DO
158
159 ! Build the corresponding transformation matrix for the transformation from
160 ! spherical to Cartesian orbitals (s2c = s*TRANSPOSE(c2s), formulas 18 and 19)
161
162 DO lx1 = 0, l
163 DO ly1 = 0, l - lx1
164 lz1 = l - lx1 - ly1
165 ic1 = co(lx1, ly1, lz1)
166 s1 = sqrt((fac(lx1)*fac(ly1)*fac(lz1))/ &
167 (fac(2*lx1)*fac(2*ly1)*fac(2*lz1)))
168 DO lx2 = 0, l
169 DO ly2 = 0, l - lx2
170 lz2 = l - lx2 - ly2
171 ic2 = co(lx2, ly2, lz2)
172 lx = lx1 + lx2
173 ly = ly1 + ly2
174 lz = lz1 + lz2
175 IF ((modulo(lx, 2) == 0) .AND. &
176 (modulo(ly, 2) == 0) .AND. &
177 (modulo(lz, 2) == 0)) THEN
178 s2 = sqrt((fac(lx2)*fac(ly2)*fac(lz2))/ &
179 (fac(2*lx2)*fac(2*ly2)*fac(2*lz2)))
180 s = fac(lx)*fac(ly)*fac(lz)*s1*s2/ &
181 (fac(lx/2)*fac(ly/2)*fac(lz/2))
182 DO is = 1, nso(l)
183 orbtramat(l)%s2c(is, ic1) = orbtramat(l)%s2c(is, ic1) + &
184 s*orbtramat(l)%c2s(is, ic2)
185 END DO
186 END IF
187 END DO
188 END DO
189 END DO
190 END DO
191
192 ! Build up the real spherical harmonics
193
194 s = sqrt(0.25_dp*dfac(2*l + 1)/pi)
195
196 DO lx = 0, l
197 DO ly = 0, l - lx
198 lz = l - lx - ly
199 ic = co(lx, ly, lz)
200 DO m = -l, l
201 is = l + m + 1
202 s1 = sqrt(dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1))
203 orbtramat(l)%slm(is, ic) = s*orbtramat(l)%c2s(is, ic)/s1
204 !MK s2 = (-1.0_dp)**m*s ! alternative S(lm) definition
205 !MK orbtramat(l)%slm(is, ic) = s2*orbtramat(l)%c2s(is,ic)/s1
206 orbtramat(l)%slm_inv(is, ic) = s1*orbtramat(l)%s2c(is, ic)/s
207 END DO
208 END DO
209 END DO
210
211 END DO
212
213 ! Save initialization status
214
215 current_maxl = maxl
216
217 END SUBROUTINE create_spherical_harmonics
218
219! **************************************************************************************************
220!> \brief Deallocate the orbital transformation matrices.
221!> \date 20.05.2004
222!> \author MK
223!> \version 1.0
224! **************************************************************************************************
225 SUBROUTINE deallocate_spherical_harmonics()
226
227 INTEGER :: l
228
229!$ IF (omp_get_level() > 0) &
230!$ CPABORT("deallocate_spherical_harmonics is not thread-safe")
231
232 IF (current_maxl > -1) THEN
233 DO l = 0, SIZE(orbtramat, 1) - 1
234 DEALLOCATE (orbtramat(l)%c2s)
235 DEALLOCATE (orbtramat(l)%s2c)
236 DEALLOCATE (orbtramat(l)%slm)
237 DEALLOCATE (orbtramat(l)%slm_inv)
238 END DO
239 DEALLOCATE (orbtramat)
240 current_maxl = -1
241 END IF
242
243 END SUBROUTINE deallocate_spherical_harmonics
244
245! **************************************************************************************************
246!> \brief Initialize or update the orbital transformation matrices.
247!> \param maxl ...
248!> \param output_unit ...
249!> \date 09.07.1999
250!> \par Variables
251!> - maxl : Maximum angular momentum quantum number
252!> \author MK
253!> \version 1.0
254! **************************************************************************************************
255 SUBROUTINE init_spherical_harmonics(maxl, output_unit)
256
257 INTEGER, INTENT(IN) :: maxl
258 INTEGER :: output_unit
259
260 CHARACTER(LEN=78) :: headline
261 INTEGER :: l, nc, ns
262
263!$ IF (omp_get_level() > 0) &
264!$ CPABORT("init_spherical_harmonics is not thread-safe")
265
266 IF (maxl < 0) THEN
267 CALL cp_abort(__location__, &
268 "A negative maximum angular momentum quantum "// &
269 "number is invalid")
270 END IF
271
272 IF (maxl > current_maxl) THEN
273
274 CALL deallocate_spherical_harmonics()
275 CALL create_spherical_harmonics(maxl)
276
277 ! Print the spherical harmonics and the orbital transformation matrices
278
279 IF (output_unit > 0) THEN
280
281 DO l = 0, maxl
282
283 nc = nco(l)
284 ns = nso(l)
285
286 headline = "CARTESIAN ORBITAL TO SPHERICAL ORBITAL "// &
287 "TRANSFORMATION MATRIX"
288 CALL write_matrix(orbtramat(l)%c2s, l, output_unit, headline)
289
290 headline = "SPHERICAL ORBITAL TO CARTESIAN ORBITAL "// &
291 "TRANSFORMATION MATRIX"
292 CALL write_matrix(orbtramat(l)%s2c, l, output_unit, headline)
293
294 headline = "SPHERICAL HARMONICS"
295 CALL write_matrix(orbtramat(l)%slm, l, output_unit, headline)
296
297 headline = "INVERSE SPHERICAL HARMONICS"
298 CALL write_matrix(orbtramat(l)%slm_inv, l, output_unit, headline)
299
300 END DO
301
302 WRITE (unit=output_unit, fmt="(A)") ""
303
304 END IF
305
306 END IF
307
308 END SUBROUTINE init_spherical_harmonics
309
310! **************************************************************************************************
311!> \brief Print a matrix to the logical unit "lunit".
312!> \param matrix ...
313!> \param l ...
314!> \param lunit ...
315!> \param headline ...
316!> \date 07.06.2000
317!> \par Variables
318!> - matrix : Matrix to be printed.
319!> - l : Angular momentum quantum number
320!> - lunit : Logical unit number.
321!> - headline: Headline of the matrix.
322!> \author MK
323!> \version 1.0
324! **************************************************************************************************
325 SUBROUTINE write_matrix(matrix, l, lunit, headline)
326
327 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: matrix
328 INTEGER, INTENT(IN) :: l, lunit
329 CHARACTER(LEN=*), INTENT(IN) :: headline
330
331 CHARACTER(12) :: symbol
332 CHARACTER(LEN=78) :: string
333 INTEGER :: from, i, ic, is, jc, lx, ly, lz, m, nc, &
334 to
335
336 ! Write headline
337
338 WRITE (unit=lunit, fmt="(/,/,T2,A)") trim(headline)
339
340 ! Write the matrix in the defined format
341
342 nc = nco(l)
343
344 DO ic = 1, nc, 3
345 from = ic
346 to = min(nc, from + 2)
347 i = 1
348 DO lx = l, 0, -1
349 DO ly = l - lx, 0, -1
350 lz = l - lx - ly
351 jc = co(lx, ly, lz)
352 IF ((jc >= from) .AND. (jc <= to)) THEN
353 symbol = cgf_symbol(1, (/lx, ly, lz/))
354 WRITE (unit=string(i:), fmt="(A18)") trim(symbol(3:12))
355 i = i + 18
356 END IF
357 END DO
358 END DO
359 WRITE (unit=lunit, fmt="(/,T8,A)") trim(string)
360 symbol = ""
361 DO m = -l, l
362 is = l + m + 1
363 symbol = sgf_symbol(1, l, m)
364 WRITE (unit=lunit, fmt="(T4,A4,3(1X,F17.12))") &
365 symbol(3:6), (matrix(is, jc), jc=from, to)
366 END DO
367 END DO
368
369 END SUBROUTINE write_matrix
370
371END 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:117
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:205
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:256
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:226