LCOV - code coverage report
Current view: top level - src - qs_tddfpt_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 148 153 96.7 %
Date: 2024-03-28 07:31:50 Functions: 7 9 77.8 %

          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             : !> \par History
      10             : !>      09-JUL-2002, TCH, development started
      11             : ! **************************************************************************************************
      12             : MODULE qs_tddfpt_utils
      13             :    USE cp_control_types,                ONLY: dft_control_type,&
      14             :                                               tddfpt_control_type
      15             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      16             :                                               cp_dbcsr_sm_fm_multiply
      17             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_scale,&
      18             :                                               cp_fm_scale_and_add,&
      19             :                                               cp_fm_trace
      20             :    USE cp_fm_cholesky,                  ONLY: cp_fm_cholesky_decompose,&
      21             :                                               cp_fm_cholesky_invert
      22             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      23             :                                               cp_fm_get_submatrix,&
      24             :                                               cp_fm_init_random,&
      25             :                                               cp_fm_release,&
      26             :                                               cp_fm_set_all,&
      27             :                                               cp_fm_set_submatrix,&
      28             :                                               cp_fm_to_fm,&
      29             :                                               cp_fm_type
      30             :    USE cp_log_handling,                 ONLY: cp_logger_get_default_io_unit
      31             :    USE dbcsr_api,                       ONLY: dbcsr_p_type
      32             :    USE kinds,                           ONLY: dp
      33             :    USE physcon,                         ONLY: evolt
      34             :    USE qs_environment_types,            ONLY: get_qs_env,&
      35             :                                               qs_environment_type
      36             :    USE qs_mo_types,                     ONLY: get_mo_set
      37             :    USE qs_p_env_methods,                ONLY: p_env_create,&
      38             :                                               p_env_psi0_changed
      39             :    USE qs_p_env_types,                  ONLY: p_env_release,&
      40             :                                               qs_p_env_type
      41             :    USE qs_tddfpt_types,                 ONLY: tddfpt_env_allocate,&
      42             :                                               tddfpt_env_deallocate,&
      43             :                                               tddfpt_env_type
      44             : #include "./base/base_uses.f90"
      45             : 
      46             :    IMPLICIT NONE
      47             : 
      48             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt_utils'
      49             :    LOGICAL, PARAMETER          :: DEBUG_THIS_MODULE = .TRUE.
      50             : 
      51             : ! **************************************************************************************************
      52             :    TYPE simple_solution_sorter
      53             :       INTEGER                               :: orbit
      54             :       INTEGER                               :: lumo
      55             :       REAL(KIND=DP)                        :: value
      56             :       TYPE(simple_solution_sorter), POINTER :: next
      57             :    END TYPE simple_solution_sorter
      58             : 
      59             :    PRIVATE
      60             : 
      61             :    ! METHODS
      62             :    PUBLIC :: tddfpt_cleanup, &
      63             :              tddfpt_init, &
      64             :              co_initial_guess, &
      65             :              find_contributions, &
      66             :              normalize, &
      67             :              reorthogonalize
      68             : 
      69             : CONTAINS
      70             : 
      71             : ! **************************************************************************************************
      72             : !> \brief Initialize some necessary structures for a tddfpt calculation.
      73             : !> \param p_env perturbation environment to be initialized
      74             : !> \param t_env tddfpt environment to be initialized
      75             : !> \param qs_env Quickstep environment with the results of a
      76             : !>                   ground state calcualtion
      77             : ! **************************************************************************************************
      78          12 :    SUBROUTINE tddfpt_init(p_env, t_env, qs_env)
      79             : 
      80             :       TYPE(qs_p_env_type), INTENT(INOUT)                 :: p_env
      81             :       TYPE(tddfpt_env_type), INTENT(out)                 :: t_env
      82             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      83             : 
      84             : !------------------!
      85             : ! create the p_env !
      86             : !------------------!
      87             : 
      88          12 :       CALL p_env_create(p_env, qs_env, orthogonal_orbitals=.TRUE.)
      89          12 :       CALL p_env_psi0_changed(p_env, qs_env) ! update the m_epsilon matrix
      90             : 
      91             :       !------------------!
      92             :       ! create the t_env !
      93             :       !------------------!
      94          12 :       CALL tddfpt_env_allocate(t_env, p_env, qs_env)
      95          12 :       CALL tddfpt_env_init(t_env, qs_env)
      96             : 
      97          12 :    END SUBROUTINE tddfpt_init
      98             : 
      99             : ! **************************************************************************************************
     100             : !> \brief Initialize t_env with meaningfull values.
     101             : !> \param t_env ...
     102             : !> \param qs_env ...
     103             : ! **************************************************************************************************
     104          12 :    SUBROUTINE tddfpt_env_init(t_env, qs_env)
     105             : 
     106             :       TYPE(tddfpt_env_type), INTENT(inout)               :: t_env
     107             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     108             : 
     109             :       INTEGER                                            :: n_spins, spin
     110          12 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     111             :       TYPE(dft_control_type), POINTER                    :: dft_control
     112             : 
     113          12 :       NULLIFY (matrix_s, dft_control)
     114             : 
     115          12 :       CALL get_qs_env(qs_env, matrix_s=matrix_s, dft_control=dft_control)
     116          12 :       n_spins = dft_control%nspins
     117          12 :       IF (dft_control%tddfpt_control%invert_S) THEN
     118          28 :          DO spin = 1, n_spins
     119          16 :             CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, t_env%invS(spin))
     120          16 :             CALL cp_fm_cholesky_decompose(t_env%invS(spin))
     121          28 :             CALL cp_fm_cholesky_invert(t_env%invS(spin))
     122             :          END DO
     123             :       END IF
     124             : 
     125          12 :    END SUBROUTINE tddfpt_env_init
     126             : 
     127             : ! **************************************************************************************************
     128             : !> \brief ...
     129             : !> \param t_env ...
     130             : !> \param p_env ...
     131             : ! **************************************************************************************************
     132          12 :    SUBROUTINE tddfpt_cleanup(t_env, p_env)
     133             : 
     134             :       TYPE(tddfpt_env_type), INTENT(inout)               :: t_env
     135             :       TYPE(qs_p_env_type), INTENT(INOUT)                 :: p_env
     136             : 
     137          12 :       CALL tddfpt_env_deallocate(t_env)
     138          12 :       CALL p_env_release(p_env)
     139             : 
     140          12 :    END SUBROUTINE tddfpt_cleanup
     141             : 
     142             : ! **************************************************************************************************
     143             : !> \brief ...
     144             : !> \param matrices ...
     145             : !> \param energies ...
     146             : !> \param n_v ...
     147             : !> \param qs_env ...
     148             : ! **************************************************************************************************
     149          12 :    SUBROUTINE co_initial_guess(matrices, energies, n_v, qs_env)
     150             : 
     151             :       TYPE(cp_fm_type), DIMENSION(:, :), POINTER         :: matrices
     152             :       REAL(kind=DP), DIMENSION(:), INTENT(OUT)           :: energies
     153             :       INTEGER, INTENT(IN)                                :: n_v
     154             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     155             : 
     156             :       INTEGER                                            :: i, n_cols, n_lumos, n_orbits, n_rows, &
     157             :                                                             n_spins, oo, spin, vo
     158             :       REAL(KIND=DP)                                      :: evd
     159          12 :       REAL(KIND=DP), ALLOCATABLE, DIMENSION(:, :)        :: guess, lumos
     160          12 :       REAL(KIND=DP), DIMENSION(:), POINTER               :: orbital_eigenvalues
     161             :       TYPE(dft_control_type), POINTER                    :: dft_control
     162             :       TYPE(simple_solution_sorter), POINTER              :: sorter_iterator, sorter_pointer, &
     163             :                                                             sorter_start
     164             :       TYPE(tddfpt_control_type), POINTER                 :: tddfpt_control
     165             : 
     166             : ! number of vectors to initialize
     167             : 
     168          12 :       NULLIFY (dft_control)
     169             : 
     170          12 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     171          12 :       tddfpt_control => dft_control%tddfpt_control
     172          12 :       n_spins = dft_control%nspins
     173          44 :       energies = 0.0_dp
     174             : 
     175          12 :       IF (.NOT. ASSOCIATED(tddfpt_control%lumos)) THEN
     176           0 :          CPABORT("LUMOS missing")
     177             :       END IF
     178             : 
     179          28 :       DO spin = 1, n_spins
     180             : 
     181          16 :          n_cols = matrices(1, spin)%matrix_struct%ncol_global
     182          16 :          n_rows = matrices(1, spin)%matrix_struct%nrow_global
     183             : 
     184          60 :          DO i = 1, n_v
     185          60 :             CALL cp_fm_set_all(matrices(i, spin), 0.0_dp)
     186             :          END DO
     187             : 
     188          16 :          CALL get_mo_set(qs_env%mos(spin), eigenvalues=orbital_eigenvalues)
     189             : 
     190          16 :          n_lumos = tddfpt_control%lumos(spin)%matrix_struct%ncol_global
     191             : 
     192          16 :          n_orbits = SIZE(orbital_eigenvalues)
     193             : 
     194             :          !-----------------------------------------!
     195             :          ! create a SORTED list of initial guesses !
     196             :          !-----------------------------------------!
     197             :          ! first element
     198          16 :          evd = tddfpt_control%lumos_eigenvalues(1, spin) - orbital_eigenvalues(n_orbits)
     199          16 :          ALLOCATE (sorter_start)
     200          16 :          sorter_start%orbit = n_orbits
     201          16 :          sorter_start%lumo = 1
     202          16 :          sorter_start%value = evd
     203          16 :          NULLIFY (sorter_start%next)
     204             :          ! rest of the elements
     205          80 :          DO oo = n_orbits, 1, -1
     206         376 :             DO vo = 1, n_lumos
     207             : 
     208         296 :                IF (oo == n_orbits .AND. vo == 1) CYCLE ! already in list
     209             : 
     210         280 :                evd = tddfpt_control%lumos_eigenvalues(vo, spin) - orbital_eigenvalues(oo)
     211             : 
     212         280 :                sorter_iterator => sorter_start
     213         280 :                NULLIFY (sorter_pointer)
     214        2314 :                DO WHILE (ASSOCIATED(sorter_iterator%next))
     215        2134 :                   IF (sorter_iterator%next%value > evd) THEN
     216             :                      sorter_pointer => sorter_iterator%next
     217             :                      EXIT
     218             :                   END IF
     219         180 :                   sorter_iterator => sorter_iterator%next
     220             :                END DO
     221             : 
     222         280 :                ALLOCATE (sorter_iterator%next)
     223         280 :                sorter_iterator%next%orbit = oo
     224         280 :                sorter_iterator%next%lumo = vo
     225         280 :                sorter_iterator%next%value = evd
     226         360 :                sorter_iterator%next%next => sorter_pointer
     227             : 
     228             :             END DO
     229             :          END DO
     230             : 
     231         112 :          ALLOCATE (lumos(n_rows, n_lumos), guess(n_rows, n_orbits))
     232             :          CALL cp_fm_get_submatrix(tddfpt_control%lumos(spin), lumos, &
     233          16 :                                   start_col=1, n_cols=n_lumos)
     234             : 
     235             :          !-------------------!
     236             :          ! fill the matrices !
     237             :          !-------------------!
     238          16 :          sorter_iterator => sorter_start
     239          60 :          DO i = 1, MIN(n_v, n_orbits*n_lumos)
     240        3996 :             guess(:, :) = 0.0_dp
     241             :             CALL dcopy(n_rows, lumos(:, sorter_iterator%lumo), 1, &
     242          44 :                        guess(:, sorter_iterator%orbit), 1)
     243          44 :             CALL cp_fm_set_submatrix(matrices(i, spin), guess)
     244          44 :             energies(i) = energies(i) + sorter_iterator%value/REAL(n_spins, dp)
     245          60 :             sorter_iterator => sorter_iterator%next
     246             :          END DO
     247          16 :          IF (n_v > n_orbits*n_lumos) THEN
     248           0 :             DO i = n_orbits*n_lumos + 1, n_v
     249           0 :                CALL cp_fm_init_random(matrices(i, spin), n_orbits)
     250           0 :                energies(i) = 1.0E38_dp
     251             :             END DO
     252             :          END IF
     253             : 
     254             :          !--------------!
     255             :          ! some cleanup !
     256             :          !--------------!
     257          16 :          DEALLOCATE (lumos, guess)
     258          16 :          sorter_iterator => sorter_start
     259         324 :          DO WHILE (ASSOCIATED(sorter_iterator))
     260         296 :             sorter_pointer => sorter_iterator
     261         296 :             sorter_iterator => sorter_iterator%next
     262         296 :             DEALLOCATE (sorter_pointer)
     263             :          END DO
     264             : 
     265             :       END DO
     266             : 
     267          24 :    END SUBROUTINE co_initial_guess
     268             : 
     269             : ! **************************************************************************************************
     270             : !> \brief ...
     271             : !> \param qs_env ...
     272             : !> \param t_env ...
     273             : ! **************************************************************************************************
     274          12 :    SUBROUTINE find_contributions(qs_env, t_env)
     275             : 
     276             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     277             :       TYPE(tddfpt_env_type), INTENT(IN)                  :: t_env
     278             : 
     279             :       INTEGER                                            :: i, j, n_ev, n_spins, occ, output_unit, &
     280             :                                                             spin, virt
     281             :       INTEGER, DIMENSION(2)                              :: nhomos, nlumos, nrows
     282             :       REAL(KIND=dp)                                      :: contribution, summed_contributions
     283          12 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: homo_coeff_col, lumo_coeff_col
     284          12 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: S_lumos
     285          12 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     286             :       TYPE(dft_control_type), POINTER                    :: dft_control
     287             :       TYPE(tddfpt_control_type)                          :: t_control
     288             : 
     289          12 :       NULLIFY (matrix_s, dft_control)
     290          24 :       output_unit = cp_logger_get_default_io_unit()
     291          12 :       CALL get_qs_env(qs_env, matrix_s=matrix_s, dft_control=dft_control)
     292             : 
     293          12 :       IF (output_unit > 0) WRITE (output_unit, *)
     294          12 :       IF (output_unit > 0) WRITE (output_unit, *)
     295             : 
     296          12 :       t_control = dft_control%tddfpt_control
     297          12 :       n_ev = t_control%n_ev
     298          12 :       n_spins = dft_control%nspins
     299             : 
     300          52 :       ALLOCATE (S_lumos(n_spins))
     301             : 
     302          28 :       DO spin = 1, n_spins
     303          16 :          nrows(spin) = t_control%lumos(spin)%matrix_struct%nrow_global
     304          16 :          nhomos(spin) = t_env%evecs(1, spin)%matrix_struct%ncol_global
     305          16 :          nlumos(spin) = t_control%lumos(spin)%matrix_struct%ncol_global
     306             :          CALL cp_fm_create(S_lumos(spin), t_control%lumos(spin)%matrix_struct, &
     307          16 :                            "S times lumos")
     308             :          CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, t_control%lumos(spin), &
     309          28 :                                       S_lumos(spin), nlumos(spin), 1.0_dp, 0.0_dp)
     310             :       END DO
     311             : 
     312             :       ALLOCATE (homo_coeff_col(MAXVAL(nrows(1:n_spins)), 1), &
     313          92 :                 lumo_coeff_col(MAXVAL(nrows(1:n_spins)), 1))
     314          44 :       DO i = 1, n_ev
     315          32 :          IF (output_unit > 0) THEN
     316          16 :             WRITE (output_unit, '(A,I3,5X,F15.6)') "  excited state : ", i, t_env%evals(i)*evolt
     317          16 :             WRITE (output_unit, *)
     318             :          END IF
     319          32 :          summed_contributions = 0.0_dp
     320          76 :          DO spin = 1, n_spins
     321          44 :             IF (n_spins == 2) THEN
     322          24 :                IF (spin == 1) THEN
     323          12 :                   IF (output_unit > 0) WRITE (output_unit, *) 'alpha:'
     324             :                ELSE
     325          12 :                   IF (output_unit > 0) WRITE (output_unit, *) 'beta:'
     326             :                END IF
     327             :             END IF
     328         252 :             searchloop: DO occ = nhomos(spin), 1, -1
     329             :                CALL cp_fm_get_submatrix(t_env%evecs(i, spin), homo_coeff_col, &
     330         176 :                                         1, occ, nrows(spin), 1)
     331        1052 :                DO virt = 1, nlumos(spin)
     332             :                   CALL cp_fm_get_submatrix(S_lumos(spin), lumo_coeff_col, &
     333         832 :                                            1, virt, nrows(spin), 1)
     334         832 :                   contribution = 0.0_dp
     335       19424 :                   DO j = 1, nrows(spin)
     336       19424 :                      contribution = contribution + homo_coeff_col(j, 1)*lumo_coeff_col(j, 1)
     337             :                   END DO
     338         832 :                   summed_contributions = summed_contributions + (contribution)**2
     339         832 :                   IF (ABS(contribution) > 5.0e-2_dp) THEN
     340         234 :                      IF (output_unit > 0) WRITE (output_unit, '(14X,I5,A,I5,10X,F8.3,5X,F8.3)') &
     341         117 :                         occ, " ->", nhomos(spin) + virt, ABS(contribution), summed_contributions
     342             :                   END IF
     343        1008 :                   IF (ABS(summed_contributions - 1.0_dp) < 1.0e-3_dp) CYCLE searchloop
     344             :                END DO
     345             :             END DO searchloop
     346             :          END DO
     347          44 :          IF (output_unit > 0) WRITE (output_unit, *)
     348             :       END DO
     349             : 
     350             :       !
     351             :       ! punch a checksum for the regs
     352          12 :       IF (output_unit > 0) THEN
     353          22 :          WRITE (output_unit, '(T2,A,E14.6)') ' TDDFPT : CheckSum  =', SQRT(SUM(t_env%evals**2))
     354             :       END IF
     355             : 
     356          12 :       CALL cp_fm_release(S_lumos)
     357             : 
     358          12 :       DEALLOCATE (homo_coeff_col, lumo_coeff_col)
     359             : 
     360          24 :    END SUBROUTINE find_contributions
     361             : 
     362             : ! **************************************************************************************************
     363             : !> \brief ...
     364             : !> \param X ...
     365             : !> \param tmp_vec ...
     366             : !> \param metric ...
     367             : ! **************************************************************************************************
     368         546 :    SUBROUTINE normalize(X, tmp_vec, metric)
     369             : 
     370             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: x, tmp_vec
     371             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: metric
     372             : 
     373             :       INTEGER                                            :: n_spins, spin
     374             :       REAL(KIND=dp)                                      :: norm, tmp
     375             : 
     376         546 :       n_spins = SIZE(x)
     377         546 :       norm = 0.0_dp
     378             : 
     379        1188 :       DO spin = 1, n_spins
     380             :          tmp = 0.0_dp
     381             :          CALL cp_dbcsr_sm_fm_multiply(metric(1)%matrix, X(spin), &
     382             :                                       tmp_vec(spin), &
     383             :                                       X(spin)%matrix_struct%ncol_global, &
     384         642 :                                       1.0_dp, 0.0_dp)
     385         642 :          CALL cp_fm_trace(X(spin), tmp_vec(spin), tmp)
     386        1188 :          norm = norm + tmp
     387             :       END DO
     388             : 
     389         546 :       norm = SQRT(norm)
     390        1188 :       DO spin = 1, n_spins
     391        1188 :          CALL cp_fm_scale((1.0_dp/norm), X(spin))
     392             :       END DO
     393             : 
     394         546 :    END SUBROUTINE normalize
     395             : 
     396             :    !---------------------------------------!
     397             :    ! x must not be changed in this routine !
     398             :    ! tmp_vec may be changed                !
     399             :    !---------------------------------------!
     400             : ! **************************************************************************************************
     401             : !> \brief ...
     402             : !> \param X ...
     403             : !> \param V_set ...
     404             : !> \param SV_set ...
     405             : !> \param work ...
     406             : !> \param n ...
     407             : ! **************************************************************************************************
     408        1092 :    SUBROUTINE reorthogonalize(X, V_set, SV_set, work, n)
     409             : 
     410             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: X
     411             :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: V_set, SV_set
     412             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: work
     413             :       INTEGER, INTENT(IN)                                :: n
     414             : 
     415             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'reorthogonalize'
     416             : 
     417             :       INTEGER                                            :: handle, i, n_spins, spin
     418             :       REAL(DP)                                           :: dot_product, tmp
     419             : 
     420        1092 :       CALL timeset(routineN, handle)
     421             : 
     422        1092 :       IF (n > 0) THEN
     423             : 
     424        1020 :          n_spins = SIZE(X)
     425        2208 :          DO spin = 1, n_spins
     426        2208 :             CALL cp_fm_to_fm(X(spin), work(spin))
     427             :          END DO
     428             : 
     429       18844 :          DO i = 1, n
     430             :             dot_product = 0.0_dp
     431       36344 :             DO spin = 1, n_spins
     432       18520 :                CALL cp_fm_trace(SV_set(i, spin), work(spin), tmp)
     433       36344 :                dot_product = dot_product + tmp
     434             :             END DO
     435       37364 :             DO spin = 1, n_spins
     436             :                CALL cp_fm_scale_and_add(1.0_dp, X(spin), &
     437       36344 :                                         -1.0_dp*dot_product, V_set(i, spin))
     438             :             END DO
     439             :          END DO
     440             : 
     441             :       END IF
     442             : 
     443        1092 :       CALL timestop(handle)
     444             : 
     445        1092 :    END SUBROUTINE reorthogonalize
     446             : 
     447             : ! **************************************************************************************************
     448             : 
     449           0 : END MODULE qs_tddfpt_utils

Generated by: LCOV version 1.15