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             :    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        6182 :    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        6182 : !$    IF (omp_get_level() > 0) &
      81           0 : !$       CPABORT("create_spherical_harmonics is not thread-safe")
      82             : 
      83        6182 :       IF (current_maxl > -1) THEN
      84             :          CALL cp_abort(__LOCATION__, &
      85             :                        "Spherical harmonics are already allocated. "// &
      86           0 :                        "Use the init routine for an update")
      87             :       END IF
      88             : 
      89        6182 :       IF (maxl < 0) THEN
      90             :          CALL cp_abort(__LOCATION__, &
      91             :                        "A negative maximum angular momentum quantum "// &
      92           0 :                        "number is invalid")
      93             :       END IF
      94             : 
      95       51446 :       ALLOCATE (orbtramat(0:maxl))
      96        6182 :       nc = ncoset(maxl)
      97        6182 :       ns = nsoset(maxl)
      98             : 
      99       39082 :       DO l = 0, maxl
     100       32900 :          nc = nco(l)
     101       32900 :          ns = nso(l)
     102      131600 :          ALLOCATE (orbtramat(l)%c2s(ns, nc))
     103     2791214 :          orbtramat(l)%c2s = 0.0_dp
     104      131600 :          ALLOCATE (orbtramat(l)%s2c(ns, nc))
     105     2791214 :          orbtramat(l)%s2c = 0.0_dp
     106      131600 :          ALLOCATE (orbtramat(l)%slm(ns, nc))
     107     2791214 :          orbtramat(l)%slm = 0.0_dp
     108      131600 :          ALLOCATE (orbtramat(l)%slm_inv(ns, nc))
     109     2797396 :          orbtramat(l)%slm_inv = 0.0_dp
     110             :       END DO
     111             : 
     112             :       ! Loop over all angular momentum quantum numbers l
     113             : 
     114       39082 :       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      141863 :          DO lx = 0, l
     120      427733 :             DO ly = 0, l - lx
     121      285870 :                lz = l - lx - ly
     122      285870 :                ic = co(lx, ly, lz)
     123     2867277 :                DO m = -l, l
     124     2472444 :                   is = l + m + 1
     125     2472444 :                   ma = ABS(m)
     126     2472444 :                   j = lx + ly - ma
     127     2758314 :                   IF ((j >= 0) .AND. (MODULO(j, 2) == 0)) THEN
     128     1014728 :                      j = j/2
     129     1014728 :                      s1 = 0.0_dp
     130     2998144 :                      DO i = 0, (l - ma)/2
     131     1983416 :                         s2 = 0.0_dp
     132     5767080 :                         DO k = 0, j
     133     3103091 :                            IF (((m < 0) .AND. (MODULO(ABS(ma - lx), 2) == 1)) .OR. &
     134     3783664 :                                ((m > 0) .AND. (MODULO(ABS(ma - lx), 2) == 0))) THEN
     135     1451930 :                               expo = (ma - lx + 2*k)/2
     136     1451930 :                               s = (-1.0_dp)**expo*SQRT(2.0_dp)
     137     2331734 :                            ELSE IF ((m == 0) .AND. (MODULO(lx, 2) == 0)) THEN
     138      568820 :                               expo = k - lx/2
     139      568820 :                               s = (-1.0_dp)**expo
     140             :                            ELSE
     141             :                               s = 0.0_dp
     142             :                            END IF
     143     5767080 :                            s2 = s2 + binomial(j, k)*binomial(ma, lx - 2*k)*s
     144             :                         END DO
     145             :                         s1 = s1 + binomial(l, i)*binomial(i, j)* &
     146     2998144 :                              (-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     1014728 :                         (2.0_dp**l*fac(l))
     152             :                   ELSE
     153     1457716 :                      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      141863 :          DO lx1 = 0, l
     163      427733 :             DO ly1 = 0, l - lx1
     164      285870 :                lz1 = l - lx1 - ly1
     165      285870 :                ic1 = co(lx1, ly1, lz1)
     166             :                s1 = SQRT((fac(lx1)*fac(ly1)*fac(lz1))/ &
     167      285870 :                          (fac(2*lx1)*fac(2*ly1)*fac(2*lz1)))
     168     1773990 :                DO lx2 = 0, l
     169     6110763 :                   DO ly2 = 0, l - lx2
     170     4445736 :                      lz2 = l - lx2 - ly2
     171     4445736 :                      ic2 = co(lx2, ly2, lz2)
     172     4445736 :                      lx = lx1 + lx2
     173     4445736 :                      ly = ly1 + ly2
     174     4445736 :                      lz = lz1 + lz2
     175             :                      IF ((MODULO(lx, 2) == 0) .AND. &
     176     4445736 :                          (MODULO(ly, 2) == 0) .AND. &
     177     1379157 :                          (MODULO(lz, 2) == 0)) THEN
     178             :                         s2 = SQRT((fac(lx2)*fac(ly2)*fac(lz2))/ &
     179     1218642 :                                   (fac(2*lx2)*fac(2*ly2)*fac(2*lz2)))
     180             :                         s = fac(lx)*fac(ly)*fac(lz)*s1*s2/ &
     181     1218642 :                             (fac(lx/2)*fac(ly/2)*fac(lz/2))
     182    14202402 :                         DO is = 1, nso(l)
     183             :                            orbtramat(l)%s2c(is, ic1) = orbtramat(l)%s2c(is, ic1) + &
     184    14202402 :                                                        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       32900 :          s = SQRT(0.25_dp*dfac(2*l + 1)/pi)
     195             : 
     196      148045 :          DO lx = 0, l
     197      427733 :             DO ly = 0, l - lx
     198      285870 :                lz = l - lx - ly
     199      285870 :                ic = co(lx, ly, lz)
     200     2867277 :                DO m = -l, l
     201     2472444 :                   is = l + m + 1
     202     2472444 :                   s1 = SQRT(dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1))
     203     2472444 :                   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     2758314 :                   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        6182 :       current_maxl = maxl
     216             : 
     217        6182 :    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       14974 :    SUBROUTINE deallocate_spherical_harmonics()
     226             : 
     227             :       INTEGER                                            :: l
     228             : 
     229       14974 : !$    IF (omp_get_level() > 0) &
     230           0 : !$       CPABORT("deallocate_spherical_harmonics is not thread-safe")
     231             : 
     232       14974 :       IF (current_maxl > -1) THEN
     233       39082 :          DO l = 0, SIZE(orbtramat, 1) - 1
     234       32900 :             DEALLOCATE (orbtramat(l)%c2s)
     235       32900 :             DEALLOCATE (orbtramat(l)%s2c)
     236       32900 :             DEALLOCATE (orbtramat(l)%slm)
     237       39082 :             DEALLOCATE (orbtramat(l)%slm_inv)
     238             :          END DO
     239        6182 :          DEALLOCATE (orbtramat)
     240        6182 :          current_maxl = -1
     241             :       END IF
     242             : 
     243       14974 :    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        6619 :    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        6619 : !$    IF (omp_get_level() > 0) &
     264           0 : !$       CPABORT("init_spherical_harmonics is not thread-safe")
     265             : 
     266        6619 :       IF (maxl < 0) THEN
     267             :          CALL cp_abort(__LOCATION__, &
     268             :                        "A negative maximum angular momentum quantum "// &
     269           0 :                        "number is invalid")
     270             :       END IF
     271             : 
     272        6619 :       IF (maxl > current_maxl) THEN
     273             : 
     274        6182 :          CALL deallocate_spherical_harmonics()
     275        6182 :          CALL create_spherical_harmonics(maxl)
     276             : 
     277             :          ! Print the spherical harmonics and the orbital transformation matrices
     278             : 
     279        6182 :          IF (output_unit > 0) THEN
     280             : 
     281           4 :             DO l = 0, maxl
     282             : 
     283           3 :                nc = nco(l)
     284           3 :                ns = nso(l)
     285             : 
     286             :                headline = "CARTESIAN ORBITAL TO SPHERICAL ORBITAL "// &
     287           3 :                           "TRANSFORMATION MATRIX"
     288           3 :                CALL write_matrix(orbtramat(l)%c2s, l, output_unit, headline)
     289             : 
     290             :                headline = "SPHERICAL ORBITAL TO CARTESIAN ORBITAL "// &
     291           3 :                           "TRANSFORMATION MATRIX"
     292           3 :                CALL write_matrix(orbtramat(l)%s2c, l, output_unit, headline)
     293             : 
     294           3 :                headline = "SPHERICAL HARMONICS"
     295           3 :                CALL write_matrix(orbtramat(l)%slm, l, output_unit, headline)
     296             : 
     297           3 :                headline = "INVERSE SPHERICAL HARMONICS"
     298           4 :                CALL write_matrix(orbtramat(l)%slm_inv, l, output_unit, headline)
     299             : 
     300             :             END DO
     301             : 
     302           1 :             WRITE (UNIT=output_unit, FMT="(A)") ""
     303             : 
     304             :          END IF
     305             : 
     306             :       END IF
     307             : 
     308        6619 :    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          12 :    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          12 :       WRITE (UNIT=lunit, FMT="(/,/,T2,A)") TRIM(headline)
     339             : 
     340             :       ! Write the matrix in the defined format
     341             : 
     342          12 :       nc = nco(l)
     343             : 
     344          28 :       DO ic = 1, nc, 3
     345          16 :          from = ic
     346          16 :          to = MIN(nc, from + 2)
     347          16 :          i = 1
     348          52 :          DO lx = l, 0, -1
     349         116 :             DO ly = l - lx, 0, -1
     350          64 :                lz = l - lx - ly
     351          64 :                jc = co(lx, ly, lz)
     352         100 :                IF ((jc >= from) .AND. (jc <= to)) THEN
     353         160 :                   symbol = cgf_symbol(1, (/lx, ly, lz/))
     354          40 :                   WRITE (UNIT=string(i:), FMT="(A18)") TRIM(symbol(3:12))
     355          40 :                   i = i + 18
     356             :                END IF
     357             :             END DO
     358             :          END DO
     359          16 :          WRITE (UNIT=lunit, FMT="(/,T8,A)") TRIM(string)
     360          16 :          symbol = ""
     361          84 :          DO m = -l, l
     362          56 :             is = l + m + 1
     363          56 :             symbol = sgf_symbol(1, l, m)
     364             :             WRITE (UNIT=lunit, FMT="(T4,A4,3(1X,F17.12))") &
     365          72 :                symbol(3:6), (matrix(is, jc), jc=from, to)
     366             :          END DO
     367             :       END DO
     368             : 
     369          12 :    END SUBROUTINE write_matrix
     370             : 
     371           0 : END MODULE orbital_transformation_matrices