(git:ccc2433)
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 ! **************************************************************************************************
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  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 
60 CONTAINS
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 
371 END 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::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