LCOV - code coverage report
Current view: top level - src - qs_tddfpt2_eigensolver.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 96.1 % 306 294
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 7 7

            Line data    Source code
       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              : MODULE qs_tddfpt2_eigensolver
       9              :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      10              :    USE cp_control_types,                ONLY: tddfpt2_control_type
      11              :    USE cp_dbcsr_api,                    ONLY: dbcsr_get_info,&
      12              :                                               dbcsr_p_type,&
      13              :                                               dbcsr_type
      14              :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply
      15              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_contracted_trace,&
      16              :                                               cp_fm_scale,&
      17              :                                               cp_fm_scale_and_add,&
      18              :                                               cp_fm_trace
      19              :    USE cp_fm_diag,                      ONLY: choose_eigv_solver
      20              :    USE cp_fm_pool_types,                ONLY: fm_pool_create_fm,&
      21              :                                               fm_pool_give_back_fm
      22              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      23              :                                               cp_fm_struct_release,&
      24              :                                               cp_fm_struct_type
      25              :    USE cp_fm_types,                     ONLY: &
      26              :         cp_fm_copy_general, cp_fm_create, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_maxabsval, &
      27              :         cp_fm_release, cp_fm_set_all, cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_type
      28              :    USE cp_log_handling,                 ONLY: cp_logger_type
      29              :    USE cp_output_handling,              ONLY: cp_iterate
      30              :    USE input_constants,                 ONLY: tddfpt_kernel_full,&
      31              :                                               tddfpt_kernel_none,&
      32              :                                               tddfpt_kernel_stda
      33              :    USE input_section_types,             ONLY: section_vals_type
      34              :    USE kinds,                           ONLY: dp,&
      35              :                                               int_8
      36              :    USE machine,                         ONLY: m_flush,&
      37              :                                               m_walltime
      38              :    USE memory_utilities,                ONLY: reallocate
      39              :    USE message_passing,                 ONLY: mp_para_env_type
      40              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      41              :    USE physcon,                         ONLY: evolt
      42              :    USE qs_environment_types,            ONLY: get_qs_env,&
      43              :                                               qs_environment_type
      44              :    USE qs_kernel_types,                 ONLY: full_kernel_env_type,&
      45              :                                               kernel_env_type
      46              :    USE qs_scf_methods,                  ONLY: eigensolver
      47              :    USE qs_tddfpt2_fhxc,                 ONLY: fhxc_kernel,&
      48              :                                               stda_kernel
      49              :    USE qs_tddfpt2_operators,            ONLY: tddfpt_apply_energy_diff,&
      50              :                                               tddfpt_apply_hfx,&
      51              :                                               tddfpt_apply_hfxlr_kernel,&
      52              :                                               tddfpt_apply_hfxsr_kernel
      53              :    USE qs_tddfpt2_restart,              ONLY: tddfpt_write_restart
      54              :    USE qs_tddfpt2_smearing_methods,     ONLY: add_smearing_aterm,&
      55              :                                               compute_fermib,&
      56              :                                               orthogonalize_smeared_occupation
      57              :    USE qs_tddfpt2_subgroups,            ONLY: tddfpt_subgroup_env_type
      58              :    USE qs_tddfpt2_types,                ONLY: tddfpt_ground_state_mos,&
      59              :                                               tddfpt_work_matrices
      60              :    USE qs_tddfpt2_utils,                ONLY: tddfpt_total_number_of_states
      61              : #include "./base/base_uses.f90"
      62              : 
      63              :    IMPLICIT NONE
      64              : 
      65              :    PRIVATE
      66              : 
      67              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_eigensolver'
      68              : 
      69              :    LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .FALSE.
      70              :    ! number of first derivative components (3: d/dx, d/dy, d/dz)
      71              :    INTEGER, PARAMETER, PRIVATE          :: nderivs = 3
      72              :    INTEGER, PARAMETER, PRIVATE          :: maxspins = 2
      73              : 
      74              :    PUBLIC :: tddfpt_davidson_solver, tddfpt_orthogonalize_psi1_psi0, &
      75              :              tddfpt_orthonormalize_psi1_psi1
      76              : 
      77              : ! **************************************************************************************************
      78              : 
      79              : CONTAINS
      80              : 
      81              : ! **************************************************************************************************
      82              : !> \brief Make TDDFPT trial vectors orthogonal to all occupied molecular orbitals.
      83              : !> \param evects            trial vectors distributed across all processors (modified on exit)
      84              : !> \param S_C0_C0T          matrix product S * C_0 * C_0^T, where C_0 is the ground state
      85              : !>                          wave function for each spin expressed in atomic basis set,
      86              : !>                          and S is the corresponding overlap matrix
      87              : !> \param qs_env ...
      88              : !> \param gs_mos ...
      89              : !> \param evals ...
      90              : !> \param tddfpt_control ...
      91              : !> \param S_C0 ...
      92              : !> \par History
      93              : !>    * 05.2016 created [Sergey Chulkov]
      94              : !>    * 05.2019 use a temporary work matrix [JHU]
      95              : !> \note  Based on the subroutine p_preortho() which was created by Thomas Chassaing on 09.2002.
      96              : !>        Should be useless when ground state MOs are computed with extremely high accuracy,
      97              : !>        as all virtual orbitals are already orthogonal to the occupied ones by design.
      98              : !>        However, when the norm of residual vectors is relatively small (e.g. less then SCF_EPS),
      99              : !>        new Krylov's vectors seem to be random and should be orthogonalised even with respect to
     100              : !>        the occupied MOs.
     101              : ! **************************************************************************************************
     102         6344 :    SUBROUTINE tddfpt_orthogonalize_psi1_psi0(evects, S_C0_C0T, qs_env, gs_mos, evals, &
     103         6344 :                                              tddfpt_control, S_C0)
     104              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in)      :: evects
     105              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(in)         :: S_C0_C0T
     106              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     107              :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     108              :          INTENT(in)                                      :: gs_mos
     109              :       REAL(kind=dp), DIMENSION(:), INTENT(in)            :: evals
     110              :       TYPE(tddfpt2_control_type), INTENT(in), POINTER    :: tddfpt_control
     111              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(in)         :: S_C0
     112              : 
     113              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_orthogonalize_psi1_psi0'
     114              : 
     115              :       INTEGER                                            :: handle, ispin, ivect, nactive, nao, &
     116              :                                                             nspins, nvects
     117              :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct
     118              :       TYPE(cp_fm_type)                                   :: evortho
     119         6344 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos
     120              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     121              : 
     122         6344 :       CALL timeset(routineN, handle)
     123              : 
     124         6344 :       nspins = SIZE(evects, 1)
     125         6344 :       nvects = SIZE(evects, 2)
     126              : 
     127         6344 :       IF (nvects > 0) THEN
     128         6344 :          IF (.NOT. tddfpt_control%do_smearing) THEN
     129        13462 :             DO ispin = 1, nspins
     130              :                CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
     131         7132 :                                    nrow_global=nao, ncol_global=nactive)
     132         7132 :                CALL cp_fm_create(evortho, matrix_struct)
     133        25904 :                DO ivect = 1, nvects
     134              :                   ! evortho: C0 * C0^T * S * C1 == (S * C0 * C0^T)^T * C1
     135              :                   CALL parallel_gemm('T', 'N', nao, nactive, nao, 1.0_dp, S_C0_C0T(ispin), &
     136        18772 :                                      evects(ispin, ivect), 0.0_dp, evortho)
     137        25904 :                   CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, ivect), -1.0_dp, evortho)
     138              :                END DO
     139        20594 :                CALL cp_fm_release(evortho)
     140              :             END DO
     141              :          ELSE
     142           14 :             NULLIFY (para_env)
     143           14 :             CALL get_qs_env(qs_env, para_env=para_env)
     144           14 :             NULLIFY (mos)
     145           56 :             ALLOCATE (mos(nspins))
     146           28 :             DO ispin = 1, nspins
     147              :                CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
     148           14 :                                    nrow_global=nao, ncol_global=nactive)
     149           14 :                CALL cp_fm_create(mos(ispin), matrix_struct)
     150           42 :                CALL cp_fm_copy_general(gs_mos(ispin)%mos_occ, mos(ispin), para_env)
     151              :             END DO
     152           28 :             DO ivect = 1, nvects
     153           14 :                CALL compute_fermib(qs_env, gs_mos, evals(ivect))
     154           28 :                CALL orthogonalize_smeared_occupation(evects(:, ivect), qs_env, mos, S_C0)
     155              :             END DO
     156           28 :             DO ispin = 1, nspins
     157           28 :                CALL cp_fm_release(mos(ispin))
     158              :             END DO
     159           14 :             DEALLOCATE (mos)
     160              :          END IF
     161              :       END IF
     162              : 
     163         6344 :       CALL timestop(handle)
     164              : 
     165         6344 :    END SUBROUTINE tddfpt_orthogonalize_psi1_psi0
     166              : 
     167              : ! **************************************************************************************************
     168              : !> \brief Check that orthogonalised TDDFPT trial vectors remain orthogonal to
     169              : !>        occupied molecular orbitals.
     170              : !> \param evects    trial vectors
     171              : !> \param S_C0      matrix product S * C_0, where C_0 is the ground state wave function
     172              : !>                  for each spin in atomic basis set, and S is the corresponding overlap matrix
     173              : !> \param max_norm  the largest possible overlap between the ground state and
     174              : !>                  excited state wave functions
     175              : !> \return true if trial vectors are non-orthogonal to occupied molecular orbitals
     176              : !> \par History
     177              : !>    * 07.2016 created [Sergey Chulkov]
     178              : !>    * 05.2019 use temporary work matrices [JHU]
     179              : ! **************************************************************************************************
     180         4094 :    FUNCTION tddfpt_is_nonorthogonal_psi1_psi0(evects, S_C0, max_norm) &
     181              :       RESULT(is_nonortho)
     182              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in)      :: evects
     183              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(in)         :: S_C0
     184              :       REAL(kind=dp), INTENT(in)                          :: max_norm
     185              :       LOGICAL                                            :: is_nonortho
     186              : 
     187              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_is_nonorthogonal_psi1_psi0'
     188              : 
     189              :       INTEGER                                            :: handle, ispin, ivect, nactive, nao, &
     190              :                                                             nocc, nspins, nvects
     191              :       REAL(kind=dp)                                      :: maxabs_val
     192              :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct, matrix_struct_tmp
     193              :       TYPE(cp_fm_type)                                   :: aortho
     194              : 
     195         4094 :       CALL timeset(routineN, handle)
     196              : 
     197         4094 :       nspins = SIZE(evects, 1)
     198         4094 :       nvects = SIZE(evects, 2)
     199              : 
     200         4094 :       is_nonortho = .FALSE.
     201              : 
     202         8738 :       loop: DO ispin = 1, nspins
     203         4648 :          CALL cp_fm_get_info(matrix=S_C0(ispin), ncol_global=nocc)
     204              :          CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
     205         4648 :                              nrow_global=nao, ncol_global=nactive)
     206              :          CALL cp_fm_struct_create(matrix_struct_tmp, nrow_global=nocc, &
     207         4648 :                                   ncol_global=nactive, template_fmstruct=matrix_struct)
     208         4648 :          CALL cp_fm_create(aortho, matrix_struct_tmp)
     209         4648 :          CALL cp_fm_struct_release(matrix_struct_tmp)
     210        16390 :          DO ivect = 1, nvects
     211              :             ! aortho = S_0^T * S * C_1
     212              :             CALL parallel_gemm('T', 'N', nocc, nactive, nao, 1.0_dp, S_C0(ispin), &
     213        11746 :                                evects(ispin, ivect), 0.0_dp, aortho)
     214        11746 :             CALL cp_fm_maxabsval(aortho, maxabs_val)
     215        11746 :             is_nonortho = maxabs_val > max_norm
     216        16390 :             IF (is_nonortho) THEN
     217            4 :                CALL cp_fm_release(aortho)
     218            4 :                EXIT loop
     219              :             END IF
     220              :          END DO
     221        18030 :          CALL cp_fm_release(aortho)
     222              :       END DO loop
     223              : 
     224         4094 :       CALL timestop(handle)
     225              : 
     226         4094 :    END FUNCTION tddfpt_is_nonorthogonal_psi1_psi0
     227              : 
     228              : ! **************************************************************************************************
     229              : !> \brief Make new TDDFPT trial vectors orthonormal to all previous trial vectors.
     230              : !> \param evects      trial vectors (modified on exit)
     231              : !> \param nvects_new  number of new trial vectors to orthogonalise
     232              : !> \param S_evects    set of matrices to store matrix product S * evects (modified on exit)
     233              : !> \param matrix_s    overlap matrix
     234              : !> \par History
     235              : !>    * 05.2016 created [Sergey Chulkov]
     236              : !>    * 02.2017 caching the matrix product S * evects [Sergey Chulkov]
     237              : !> \note \parblock
     238              : !>       Based on the subroutines reorthogonalize() and normalize() which were originally created
     239              : !>       by Thomas Chassaing on 03.2003.
     240              : !>
     241              : !>       In order to orthogonalise a trial vector C3 = evects(:,3) with respect to previously
     242              : !>       orthogonalised vectors C1 = evects(:,1) and C2 = evects(:,2) we need to compute the
     243              : !>       quantity C3'' using the following formulae:
     244              : !>          C3'  = C3  - Tr(C3^T  * S * C1) * C1,
     245              : !>          C3'' = C3' - Tr(C3'^T * S * C2) * C2,
     246              : !>       which can be expanded as:
     247              : !>          C3'' = C3 - Tr(C3^T  * S * C1) * C1 - Tr(C3^T * S * C2) * C2 +
     248              : !>                 Tr(C3^T * S * C1) * Tr(C2^T * S * C1) * C2 .
     249              : !>       In case of unlimited float-point precision, the last term in above expression is exactly 0,
     250              : !>       due to orthogonality condition between C1 and C2. In this case the expression could be
     251              : !>       simplified as (taking into account the identity: Tr(A * S * B) = Tr(B * S * A)):
     252              : !>          C3'' = C3 - Tr(C1^T  * S * C3) * C1 - Tr(C2^T * S * C3) * C2 ,
     253              : !>       which means we do not need the variable S_evects to keep the matrix products S * Ci .
     254              : !>
     255              : !>       In reality, however, we deal with limited float-point precision arithmetic meaning that
     256              : !>       the trace Tr(C2^T * S * C1) is close to 0 but does not equal to 0 exactly. The term
     257              : !>          Tr(C3^T * S * C1) * Tr(C2^T * S * C1) * C2
     258              : !>       can not be ignored anymore. Ignorance of this term will lead to numerical instability
     259              : !>       when the trace Tr(C3^T * S * C1) is large enough.
     260              : !>       \endparblock
     261              : ! **************************************************************************************************
     262         6344 :    SUBROUTINE tddfpt_orthonormalize_psi1_psi1(evects, nvects_new, S_evects, matrix_s)
     263              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in)      :: evects
     264              :       INTEGER, INTENT(in)                                :: nvects_new
     265              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT)   :: S_evects
     266              :       TYPE(dbcsr_type), POINTER                          :: matrix_s
     267              : 
     268              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_orthonormalize_psi1_psi1'
     269              : 
     270              :       INTEGER                                            :: handle, ispin, ivect, jvect, nspins, &
     271              :                                                             nvects_old, nvects_total
     272              :       INTEGER, DIMENSION(maxspins)                       :: nactive
     273              :       REAL(kind=dp)                                      :: norm
     274              :       REAL(kind=dp), DIMENSION(maxspins)                 :: weights
     275              : 
     276         6344 :       CALL timeset(routineN, handle)
     277              : 
     278         6344 :       nspins = SIZE(evects, 1)
     279         6344 :       nvects_total = SIZE(evects, 2)
     280         6344 :       nvects_old = nvects_total - nvects_new
     281              : 
     282              :       IF (debug_this_module) THEN
     283              :          CPASSERT(SIZE(S_evects, 1) == nspins)
     284              :          CPASSERT(SIZE(S_evects, 2) == nvects_total)
     285              :          CPASSERT(nvects_old >= 0)
     286              :       END IF
     287              : 
     288        13490 :       DO ispin = 1, nspins
     289        13490 :          CALL cp_fm_get_info(matrix=evects(ispin, 1), ncol_global=nactive(ispin))
     290              :       END DO
     291              : 
     292        22724 :       DO jvect = nvects_old + 1, nvects_total
     293              :          ! <psi1_i | psi1_j>
     294       140026 :          DO ivect = 1, jvect - 1
     295       123646 :             CALL cp_fm_trace(evects(:, jvect), S_evects(:, ivect), weights(1:nspins), accurate=.FALSE.)
     296       268388 :             norm = SUM(weights(1:nspins))
     297              : 
     298       284768 :             DO ispin = 1, nspins
     299       268388 :                CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, jvect), -norm, evects(ispin, ivect))
     300              :             END DO
     301              :          END DO
     302              : 
     303              :          ! <psi1_j | psi1_j>
     304        35166 :          DO ispin = 1, nspins
     305              :             CALL cp_dbcsr_sm_fm_multiply(matrix_s, evects(ispin, jvect), S_evects(ispin, jvect), &
     306        35166 :                                          ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
     307              :          END DO
     308              : 
     309        16380 :          CALL cp_fm_trace(evects(:, jvect), S_evects(:, jvect), weights(1:nspins), accurate=.FALSE.)
     310              : 
     311        35166 :          norm = SUM(weights(1:nspins))
     312        16380 :          norm = 1.0_dp/SQRT(norm)
     313              : 
     314        41510 :          DO ispin = 1, nspins
     315        18786 :             CALL cp_fm_scale(norm, evects(ispin, jvect))
     316        35166 :             CALL cp_fm_scale(norm, S_evects(ispin, jvect))
     317              :          END DO
     318              :       END DO
     319              : 
     320         6344 :       CALL timestop(handle)
     321              : 
     322         6344 :    END SUBROUTINE tddfpt_orthonormalize_psi1_psi1
     323              : 
     324              : ! **************************************************************************************************
     325              : !> \brief Compute action matrix-vector products.
     326              : !> \param Aop_evects            action of TDDFPT operator on trial vectors (modified on exit)
     327              : !> \param evects                TDDFPT trial vectors
     328              : !> \param S_evects              cached matrix product S * evects where S is the overlap matrix
     329              : !>                              in primary basis set
     330              : !> \param gs_mos                molecular orbitals optimised for the ground state
     331              : !> \param tddfpt_control        control section for tddfpt
     332              : !> \param matrix_ks             Kohn-Sham matrix
     333              : !> \param qs_env                Quickstep environment
     334              : !> \param kernel_env            kernel environment
     335              : !> \param sub_env               parallel (sub)group environment
     336              : !> \param work_matrices         collection of work matrices (modified on exit)
     337              : !> \param matrix_s ...
     338              : !> \par History
     339              : !>    * 06.2016 created [Sergey Chulkov]
     340              : !>    * 03.2017 refactored [Sergey Chulkov]
     341              : ! **************************************************************************************************
     342         5256 :    SUBROUTINE tddfpt_compute_Aop_evects(Aop_evects, evects, S_evects, gs_mos, tddfpt_control, &
     343         5256 :                                         matrix_ks, qs_env, kernel_env, &
     344              :                                         sub_env, work_matrices, matrix_s)
     345              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT)   :: Aop_evects
     346              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: evects, S_evects
     347              :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     348              :          INTENT(in)                                      :: gs_mos
     349              :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     350              :       TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in)       :: matrix_ks
     351              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     352              :       TYPE(kernel_env_type), INTENT(in)                  :: kernel_env
     353              :       TYPE(tddfpt_subgroup_env_type), INTENT(in)         :: sub_env
     354              :       TYPE(tddfpt_work_matrices), INTENT(inout)          :: work_matrices
     355              :       TYPE(dbcsr_type), POINTER                          :: matrix_s
     356              : 
     357              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_compute_Aop_evects'
     358              : 
     359              :       INTEGER                                            :: handle, ispin, ivect, nspins, nvects
     360              :       INTEGER, DIMENSION(maxspins)                       :: nmo_occ
     361              :       LOGICAL                                            :: do_admm, do_hfx, do_lri_response, &
     362              :                                                             is_rks_triplets, re_int
     363              :       REAL(KIND=dp)                                      :: rcut, scale
     364              :       TYPE(cp_fm_type)                                   :: fm_dummy
     365              :       TYPE(full_kernel_env_type), POINTER                :: kernel_env_admm_aux
     366              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     367              : 
     368         5256 :       CALL timeset(routineN, handle)
     369              : 
     370         5256 :       nspins = SIZE(evects, 1)
     371         5256 :       nvects = SIZE(evects, 2)
     372         5256 :       do_hfx = tddfpt_control%do_hfx
     373         5256 :       do_admm = tddfpt_control%do_admm
     374         5256 :       IF (do_admm) THEN
     375          622 :          kernel_env_admm_aux => kernel_env%admm_kernel
     376              :       ELSE
     377         4634 :          NULLIFY (kernel_env_admm_aux)
     378              :       END IF
     379         5256 :       is_rks_triplets = tddfpt_control%rks_triplets
     380         5256 :       do_lri_response = tddfpt_control%do_lrigpw
     381              : 
     382              :       IF (debug_this_module) THEN
     383              :          CPASSERT(nspins > 0)
     384              :          CPASSERT(SIZE(Aop_evects, 1) == nspins)
     385              :          CPASSERT(SIZE(Aop_evects, 2) == nvects)
     386              :          CPASSERT(SIZE(S_evects, 1) == nspins)
     387              :          CPASSERT(SIZE(S_evects, 2) == nvects)
     388              :          CPASSERT(SIZE(gs_mos) == nspins)
     389              :       END IF
     390              : 
     391        11190 :       DO ispin = 1, nspins
     392         5256 :          nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
     393              :       END DO
     394              : 
     395         5256 :       IF (nvects > 0) THEN
     396         5256 :          CALL cp_fm_get_info(evects(1, 1), para_env=para_env)
     397         5256 :          IF (ALLOCATED(work_matrices%evects_sub)) THEN
     398           24 :             DO ivect = 1, nvects
     399           40 :                DO ispin = 1, nspins
     400           16 :                   ASSOCIATE (evect => evects(ispin, ivect), work_matrix => work_matrices%evects_sub(ispin, ivect))
     401           16 :                   IF (ASSOCIATED(evect%matrix_struct)) THEN
     402           16 :                   IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
     403            8 :                      CALL cp_fm_copy_general(evect, work_matrix, para_env)
     404              :                   ELSE
     405            8 :                      CALL cp_fm_copy_general(evect, fm_dummy, para_env)
     406              :                   END IF
     407            0 :                   ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
     408            0 :                   CALL cp_fm_copy_general(fm_dummy, work_matrix, para_env)
     409              :                   ELSE
     410            0 :                   CALL cp_fm_copy_general(fm_dummy, fm_dummy, para_env)
     411              :                   END IF
     412              :                   END ASSOCIATE
     413              :                END DO
     414              :             END DO
     415              :          END IF
     416              : 
     417         5256 :          IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
     418              :             ! full TDDFPT kernel
     419              :             CALL fhxc_kernel(Aop_evects, evects, is_rks_triplets, do_hfx, do_admm, qs_env, &
     420              :                              kernel_env%full_kernel, kernel_env_admm_aux, sub_env, work_matrices, &
     421              :                              tddfpt_control%admm_symm, tddfpt_control%admm_xc_correction, &
     422         2898 :                              do_lri_response, tddfpt_mgrid=tddfpt_control%mgrid_is_explicit)
     423         2358 :          ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
     424              :             ! sTDA kernel
     425              :             CALL stda_kernel(Aop_evects, evects, is_rks_triplets, qs_env, tddfpt_control%stda_control, &
     426         2264 :                              kernel_env%stda_kernel, sub_env, work_matrices)
     427           94 :          ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
     428              :             ! No kernel
     429          340 :             DO ivect = 1, nvects
     430          586 :                DO ispin = 1, nspins
     431          492 :                   CALL cp_fm_set_all(Aop_evects(ispin, ivect), 0.0_dp)
     432              :                END DO
     433              :             END DO
     434              :          ELSE
     435            0 :             CPABORT("Kernel type not implemented")
     436              :          END IF
     437              : 
     438         5256 :          IF (ALLOCATED(work_matrices%evects_sub)) THEN
     439           24 :             DO ivect = 1, nvects
     440           40 :                DO ispin = 1, nspins
     441              :                   ASSOCIATE (Aop_evect => Aop_evects(ispin, ivect), &
     442           16 :                              work_matrix => work_matrices%Aop_evects_sub(ispin, ivect))
     443           16 :                   IF (ASSOCIATED(Aop_evect%matrix_struct)) THEN
     444           16 :                   IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
     445            8 :                      CALL cp_fm_copy_general(work_matrix, Aop_evect, para_env)
     446              :                   ELSE
     447            8 :                      CALL cp_fm_copy_general(fm_dummy, Aop_evect, para_env)
     448              :                   END IF
     449            0 :                   ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
     450            0 :                   CALL cp_fm_copy_general(work_matrix, fm_dummy, para_env)
     451              :                   ELSE
     452            0 :                   CALL cp_fm_copy_general(fm_dummy, fm_dummy, para_env)
     453              :                   END IF
     454              :                   END ASSOCIATE
     455              :                END DO
     456              :             END DO
     457              :          END IF
     458              : 
     459              :          ! orbital energy difference term
     460              :          CALL tddfpt_apply_energy_diff(Aop_evects=Aop_evects, evects=evects, S_evects=S_evects, &
     461         5256 :                                        gs_mos=gs_mos, matrix_ks=matrix_ks)
     462              : 
     463              :          ! if smeared occupation, then add aCCSX here
     464         5256 :          IF (tddfpt_control%do_smearing) THEN
     465           24 :             DO ivect = 1, nvects
     466           36 :                DO ispin = 1, nspins
     467              :                   CALL add_smearing_aterm(qs_env, Aop_evects(ispin, ivect), evects(ispin, ivect), &
     468              :                                           S_evects(ispin, ivect), gs_mos(ispin)%mos_occ, &
     469           24 :                                           tddfpt_control%smeared_occup(ispin)%fermia, matrix_s)
     470              :                END DO
     471              :             END DO
     472              :          END IF
     473              : 
     474         5256 :          IF (do_hfx) THEN
     475         1146 :             IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
     476              :                ! full TDDFPT kernel
     477              :                CALL tddfpt_apply_hfx(Aop_evects=Aop_evects, evects=evects, gs_mos=gs_mos, do_admm=do_admm, &
     478              :                                      qs_env=qs_env, wfm_rho_orb=work_matrices%hfx_fm_ao_ao, &
     479              :                                      work_hmat_symm=work_matrices%hfx_hmat_symm, &
     480              :                                      work_hmat_asymm=work_matrices%hfx_hmat_asymm, &
     481              :                                      work_rho_ia_ao_symm=work_matrices%hfx_rho_ao_symm, &
     482         1146 :                                      work_rho_ia_ao_asymm=work_matrices%hfx_rho_ao_asymm)
     483            0 :             ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
     484              :                ! sTDA kernel
     485              :                ! special treatment of HFX term
     486            0 :             ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
     487              :                ! No kernel
     488              :                ! drop kernel contribution of HFX term
     489              :             ELSE
     490            0 :                CPABORT("Kernel type not implemented")
     491              :             END IF
     492              :          END IF
     493              :          ! short/long range HFX
     494         5256 :          IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
     495         2898 :             IF (tddfpt_control%do_hfxsr) THEN
     496           22 :                re_int = tddfpt_control%hfxsr_re_int
     497              :                ! symmetric dmat
     498              :                CALL tddfpt_apply_hfxsr_kernel(Aop_evects, evects, gs_mos, qs_env, &
     499              :                                               kernel_env%full_kernel%admm_env, &
     500              :                                               kernel_env%full_kernel%hfxsr_section, &
     501              :                                               kernel_env%full_kernel%x_data, 1, re_int, &
     502              :                                               work_rho_ia_ao=work_matrices%hfxsr_rho_ao_symm, &
     503              :                                               work_hmat=work_matrices%hfxsr_hmat_symm, &
     504           22 :                                               wfm_rho_orb=work_matrices%hfxsr_fm_ao_ao)
     505              :                ! antisymmetric dmat
     506              :                CALL tddfpt_apply_hfxsr_kernel(Aop_evects, evects, gs_mos, qs_env, &
     507              :                                               kernel_env%full_kernel%admm_env, &
     508              :                                               kernel_env%full_kernel%hfxsr_section, &
     509              :                                               kernel_env%full_kernel%x_data, -1, .FALSE., &
     510              :                                               work_rho_ia_ao=work_matrices%hfxsr_rho_ao_asymm, &
     511              :                                               work_hmat=work_matrices%hfxsr_hmat_asymm, &
     512           22 :                                               wfm_rho_orb=work_matrices%hfxsr_fm_ao_ao)
     513           22 :                tddfpt_control%hfxsr_re_int = .FALSE.
     514              :             END IF
     515         2898 :             IF (tddfpt_control%do_hfxlr) THEN
     516           36 :                rcut = tddfpt_control%hfxlr_rcut
     517           36 :                scale = tddfpt_control%hfxlr_scale
     518          108 :                DO ivect = 1, nvects
     519          108 :                   IF (ALLOCATED(work_matrices%evects_sub)) THEN
     520            0 :                      IF (ASSOCIATED(work_matrices%evects_sub(1, ivect)%matrix_struct)) THEN
     521              :                         CALL tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, scale, work_matrices, &
     522              :                                                        work_matrices%evects_sub(:, ivect), &
     523            0 :                                                        work_matrices%Aop_evects_sub(:, ivect))
     524              :                      ELSE
     525              :                         ! skip trial vectors which are assigned to different parallel groups
     526              :                         CYCLE
     527              :                      END IF
     528              :                   ELSE
     529              :                      CALL tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, scale, work_matrices, &
     530           72 :                                                     evects(:, ivect), Aop_evects(:, ivect))
     531              :                   END IF
     532              :                END DO
     533              :             END IF
     534              :          END IF
     535              : 
     536              :       END IF
     537              : 
     538         5256 :       CALL timestop(handle)
     539              : 
     540         5256 :    END SUBROUTINE tddfpt_compute_Aop_evects
     541              : 
     542              : ! **************************************************************************************************
     543              : !> \brief Solve eigenproblem for the reduced action matrix and find new Ritz eigenvectors and
     544              : !>        eigenvalues.
     545              : !> \param ritz_vects       Ritz eigenvectors (initialised on exit)
     546              : !> \param Aop_ritz         approximate action of TDDFPT operator on Ritz vectors (initialised on exit)
     547              : !> \param evals            Ritz eigenvalues (initialised on exit)
     548              : !> \param krylov_vects     Krylov's vectors
     549              : !> \param Aop_krylov       action of TDDFPT operator on Krylov's vectors
     550              : !> \param Atilde ...
     551              : !> \param nkvo ...
     552              : !> \param nkvn ...
     553              : !> \par History
     554              : !>    * 06.2016 created [Sergey Chulkov]
     555              : !>    * 03.2017 altered prototype, OpenMP parallelisation [Sergey Chulkov]
     556              : ! **************************************************************************************************
     557         5256 :    SUBROUTINE tddfpt_compute_ritz_vects(ritz_vects, Aop_ritz, evals, krylov_vects, Aop_krylov, &
     558              :                                         Atilde, nkvo, nkvn)
     559              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: ritz_vects, Aop_ritz
     560              :       REAL(kind=dp), DIMENSION(:), INTENT(out)           :: evals
     561              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: krylov_vects, Aop_krylov
     562              :       REAL(kind=dp), DIMENSION(:, :), POINTER            :: Atilde
     563              :       INTEGER, INTENT(IN)                                :: nkvo, nkvn
     564              : 
     565              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_compute_ritz_vects'
     566              : 
     567              :       INTEGER                                            :: handle, ikv, irv, ispin, nkvs, nrvs, &
     568              :                                                             nspins
     569              :       REAL(kind=dp)                                      :: act
     570         5256 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: at12, at21, at22, evects_Atilde
     571              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env_global
     572              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
     573              :       TYPE(cp_fm_type)                                   :: Atilde_fm, evects_Atilde_fm
     574              : 
     575         5256 :       CALL timeset(routineN, handle)
     576              : 
     577         5256 :       nspins = SIZE(krylov_vects, 1)
     578         5256 :       nkvs = SIZE(krylov_vects, 2)
     579         5256 :       nrvs = SIZE(ritz_vects, 2)
     580         5256 :       CPASSERT(nkvs == nkvo + nkvn)
     581              : 
     582         5256 :       CALL cp_fm_get_info(krylov_vects(1, 1), context=blacs_env_global)
     583              : 
     584         5256 :       CALL cp_fm_struct_create(fm_struct, nrow_global=nkvs, ncol_global=nkvs, context=blacs_env_global)
     585         5256 :       CALL cp_fm_create(Atilde_fm, fm_struct)
     586         5256 :       CALL cp_fm_create(evects_Atilde_fm, fm_struct)
     587         5256 :       CALL cp_fm_struct_release(fm_struct)
     588              : 
     589              :       ! *** compute upper-diagonal reduced action matrix ***
     590         5256 :       CALL reallocate(Atilde, 1, nkvs, 1, nkvs)
     591              :       ! TO DO: the subroutine 'cp_fm_contracted_trace' will compute all elements of
     592              :       ! the matrix 'Atilde', however only upper-triangular elements are actually needed
     593              :       !
     594         5256 :       IF (nkvo == 0) THEN
     595              :          CALL cp_fm_contracted_trace(Aop_krylov(:, 1:nkvs), krylov_vects(:, 1:nkvs), &
     596         1166 :                                      Atilde(1:nkvs, 1:nkvs), accurate=.FALSE.)
     597              :       ELSE
     598        36810 :          ALLOCATE (at12(nkvn, nkvo), at21(nkvo, nkvn), at22(nkvn, nkvn))
     599              :          CALL cp_fm_contracted_trace(Aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, 1:nkvo), &
     600         4090 :                                      at12, accurate=.FALSE.)
     601       136662 :          Atilde(nkvo + 1:nkvs, 1:nkvo) = at12(1:nkvn, 1:nkvo)
     602              :          CALL cp_fm_contracted_trace(Aop_krylov(:, 1:nkvo), krylov_vects(:, nkvo + 1:nkvs), &
     603         4090 :                                      at21, accurate=.FALSE.)
     604       113822 :          Atilde(1:nkvo, nkvo + 1:nkvs) = at21(1:nkvo, 1:nkvn)
     605              :          CALL cp_fm_contracted_trace(Aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, nkvo + 1:nkvs), &
     606         4090 :                                      at22, accurate=.FALSE.)
     607        51722 :          Atilde(nkvo + 1:nkvs, nkvo + 1:nkvs) = at22(1:nkvn, 1:nkvn)
     608         4090 :          DEALLOCATE (at12, at21, at22)
     609              :       END IF
     610      1597936 :       Atilde = 0.5_dp*(Atilde + TRANSPOSE(Atilde))
     611         5256 :       CALL cp_fm_set_submatrix(Atilde_fm, Atilde)
     612              : 
     613              :       ! *** solve an eigenproblem for the reduced matrix ***
     614         5256 :       CALL choose_eigv_solver(Atilde_fm, evects_Atilde_fm, evals(1:nkvs))
     615              : 
     616        21024 :       ALLOCATE (evects_Atilde(nkvs, nrvs))
     617         5256 :       CALL cp_fm_get_submatrix(evects_Atilde_fm, evects_Atilde, start_row=1, start_col=1, n_rows=nkvs, n_cols=nrvs)
     618         5256 :       CALL cp_fm_release(evects_Atilde_fm)
     619         5256 :       CALL cp_fm_release(Atilde_fm)
     620              : 
     621              : !$OMP PARALLEL DO DEFAULT(NONE), &
     622              : !$OMP             PRIVATE(act, ikv, irv, ispin), &
     623         5256 : !$OMP             SHARED(Aop_krylov, Aop_ritz, krylov_vects, evects_Atilde, nkvs, nrvs, nspins, ritz_vects)
     624              :       DO irv = 1, nrvs
     625              :          DO ispin = 1, nspins
     626              :             CALL cp_fm_set_all(ritz_vects(ispin, irv), 0.0_dp)
     627              :             CALL cp_fm_set_all(Aop_ritz(ispin, irv), 0.0_dp)
     628              :          END DO
     629              : 
     630              :          DO ikv = 1, nkvs
     631              :             act = evects_Atilde(ikv, irv)
     632              :             DO ispin = 1, nspins
     633              :                CALL cp_fm_scale_and_add(1.0_dp, ritz_vects(ispin, irv), &
     634              :                                         act, krylov_vects(ispin, ikv))
     635              :                CALL cp_fm_scale_and_add(1.0_dp, Aop_ritz(ispin, irv), &
     636              :                                         act, Aop_krylov(ispin, ikv))
     637              :             END DO
     638              :          END DO
     639              :       END DO
     640              : !$OMP END PARALLEL DO
     641              : 
     642         5256 :       DEALLOCATE (evects_Atilde)
     643              : 
     644         5256 :       CALL timestop(handle)
     645              : 
     646        15768 :    END SUBROUTINE tddfpt_compute_ritz_vects
     647              : 
     648              : ! **************************************************************************************************
     649              : !> \brief Expand Krylov space by computing residual vectors.
     650              : !> \param residual_vects          residual vectors (modified on exit)
     651              : !> \param evals                   Ritz eigenvalues (modified on exit)
     652              : !> \param ritz_vects              Ritz eigenvectors
     653              : !> \param Aop_ritz                approximate action of TDDFPT operator on Ritz vectors
     654              : !> \param gs_mos                  molecular orbitals optimised for the ground state
     655              : !> \param matrix_s                overlap matrix
     656              : !> \par History
     657              : !>    * 06.2016 created [Sergey Chulkov]
     658              : !>    * 03.2017 refactored to achieve significant performance gain [Sergey Chulkov]
     659              : ! **************************************************************************************************
     660         4094 :    SUBROUTINE tddfpt_compute_residual_vects(residual_vects, evals, ritz_vects, Aop_ritz, gs_mos, &
     661              :                                             matrix_s)
     662              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in)      :: residual_vects
     663              :       REAL(kind=dp), DIMENSION(:), INTENT(in)            :: evals
     664              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in)      :: ritz_vects, Aop_ritz
     665              :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     666              :          INTENT(in)                                      :: gs_mos
     667              :       TYPE(dbcsr_type), POINTER                          :: matrix_s
     668              : 
     669              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_compute_residual_vects'
     670              :       REAL(kind=dp), PARAMETER :: eref_scale = 0.99_dp, threshold = 16.0_dp*EPSILON(1.0_dp)
     671              : 
     672              :       INTEGER                                            :: handle, icol_local, irow_local, irv, &
     673              :                                                             ispin, nao, ncols_local, nrows_local, &
     674              :                                                             nrvs, nspins
     675         4094 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices_local, row_indices_local
     676              :       INTEGER, DIMENSION(maxspins)                       :: nactive, nmo_virt
     677              :       REAL(kind=dp)                                      :: e_occ_plus_lambda, eref, lambda
     678              :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
     679         4094 :          POINTER                                         :: weights_ldata
     680              :       TYPE(cp_fm_struct_type), POINTER                   :: ao_mo_struct, virt_mo_struct
     681         4094 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: awork, vomat
     682              : 
     683         4094 :       CALL timeset(routineN, handle)
     684              : 
     685         4094 :       nspins = SIZE(residual_vects, 1)
     686         4094 :       nrvs = SIZE(residual_vects, 2)
     687              : 
     688         4094 :       IF (nrvs > 0) THEN
     689         4094 :          CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
     690        29766 :          ALLOCATE (awork(nspins), vomat(nspins))
     691         8742 :          DO ispin = 1, nspins
     692         4648 :             nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
     693              :             !
     694              :             CALL cp_fm_get_info(matrix=ritz_vects(ispin, 1), matrix_struct=ao_mo_struct, &
     695         4648 :                                 ncol_global=nactive(ispin))
     696         4648 :             CALL cp_fm_create(awork(ispin), ao_mo_struct)
     697              :             !
     698              :             CALL cp_fm_struct_create(virt_mo_struct, nrow_global=nmo_virt(ispin), &
     699         4648 :                                      ncol_global=nactive(ispin), template_fmstruct=ao_mo_struct)
     700         4648 :             CALL cp_fm_create(vomat(ispin), virt_mo_struct)
     701         8742 :             CALL cp_fm_struct_release(virt_mo_struct)
     702              :          END DO
     703              : 
     704              :          ! *** actually compute residual vectors ***
     705        14202 :          DO irv = 1, nrvs
     706        10108 :             lambda = evals(irv)
     707              : 
     708        25960 :             DO ispin = 1, nspins
     709              :                CALL cp_fm_get_info(vomat(ispin), nrow_local=nrows_local, &
     710              :                                    ncol_local=ncols_local, row_indices=row_indices_local, &
     711        11758 :                                    col_indices=col_indices_local, local_data=weights_ldata)
     712              : 
     713              :                ! awork := Ab(ispin, irv) - evals(irv) b(ispin, irv), where 'b' is a Ritz vector
     714              :                CALL cp_dbcsr_sm_fm_multiply(matrix_s, ritz_vects(ispin, irv), awork(ispin), &
     715        11758 :                                             ncol=nactive(ispin), alpha=-lambda, beta=0.0_dp)
     716        11758 :                CALL cp_fm_scale_and_add(1.0_dp, awork(ispin), 1.0_dp, Aop_ritz(ispin, irv))
     717              :                !
     718              :                CALL parallel_gemm('T', 'N', nmo_virt(ispin), nactive(ispin), nao, 1.0_dp, gs_mos(ispin)%mos_virt, &
     719        11758 :                                   awork(ispin), 0.0_dp, vomat(ispin))
     720              : 
     721       105412 :                DO icol_local = 1, ncols_local
     722        93654 :                   e_occ_plus_lambda = gs_mos(ispin)%evals_occ(col_indices_local(icol_local)) + lambda
     723              : 
     724      3403204 :                   DO irow_local = 1, nrows_local
     725      3297792 :                      eref = gs_mos(ispin)%evals_virt(row_indices_local(irow_local)) - e_occ_plus_lambda
     726              : 
     727              :                      ! eref = e_virt - e_occ - lambda = e_virt - e_occ - (eref_scale*lambda + (1-eref_scale)*lambda);
     728              :                      ! eref_new = e_virt - e_occ - eref_scale*lambda = eref + (1 - eref_scale)*lambda
     729      3297792 :                      IF (ABS(eref) < threshold) &
     730           74 :                         eref = eref + (1.0_dp - eref_scale)*lambda
     731              : 
     732      3391446 :                      weights_ldata(irow_local, icol_local) = weights_ldata(irow_local, icol_local)/eref
     733              :                   END DO
     734              :                END DO
     735              : 
     736              :                CALL parallel_gemm('N', 'N', nao, nactive(ispin), nmo_virt(ispin), 1.0_dp, gs_mos(ispin)%mos_virt, &
     737        33624 :                                   vomat(ispin), 0.0_dp, residual_vects(ispin, irv))
     738              :             END DO
     739              :          END DO
     740              :          !
     741         4094 :          CALL cp_fm_release(awork)
     742         8188 :          CALL cp_fm_release(vomat)
     743              :       END IF
     744              : 
     745         4094 :       CALL timestop(handle)
     746              : 
     747         8188 :    END SUBROUTINE tddfpt_compute_residual_vects
     748              : 
     749              : ! **************************************************************************************************
     750              : !> \brief Perform Davidson iterations.
     751              : !> \param evects                TDDFPT trial vectors (modified on exit)
     752              : !> \param evals                 TDDFPT eigenvalues (modified on exit)
     753              : !> \param S_evects              cached matrix product S * evects (modified on exit)
     754              : !> \param gs_mos                molecular orbitals optimised for the ground state
     755              : !> \param tddfpt_control        TDDFPT control parameters
     756              : !> \param matrix_ks             Kohn-Sham matrix
     757              : !> \param qs_env                Quickstep environment
     758              : !> \param kernel_env            kernel environment
     759              : !> \param sub_env               parallel (sub)group environment
     760              : !> \param logger                CP2K logger
     761              : !> \param iter_unit             I/O unit to write basic iteration information
     762              : !> \param energy_unit           I/O unit to write detailed energy information
     763              : !> \param tddfpt_print_section  TDDFPT print input section (need to write TDDFPT restart files)
     764              : !> \param work_matrices         collection of work matrices (modified on exit)
     765              : !> \return energy convergence achieved (in Hartree)
     766              : !> \par History
     767              : !>    * 03.2017 code related to Davidson eigensolver has been moved here from the main subroutine
     768              : !>              tddfpt() [Sergey Chulkov]
     769              : !> \note Based on the subroutines apply_op() and iterative_solver() originally created by
     770              : !>       Thomas Chassaing in 2002.
     771              : ! **************************************************************************************************
     772         1166 :    FUNCTION tddfpt_davidson_solver(evects, evals, S_evects, gs_mos, tddfpt_control, &
     773              :                                    matrix_ks, qs_env, kernel_env, &
     774              :                                    sub_env, logger, iter_unit, energy_unit, &
     775              :                                    tddfpt_print_section, work_matrices) RESULT(conv)
     776              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(inout)   :: evects
     777              :       REAL(kind=dp), DIMENSION(:), INTENT(inout)         :: evals
     778              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(inout)   :: S_evects
     779              :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     780              :          INTENT(in)                                      :: gs_mos
     781              :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     782              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks
     783              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     784              :       TYPE(kernel_env_type), INTENT(in)                  :: kernel_env
     785              :       TYPE(tddfpt_subgroup_env_type), INTENT(in)         :: sub_env
     786              :       TYPE(cp_logger_type), POINTER                      :: logger
     787              :       INTEGER, INTENT(in)                                :: iter_unit, energy_unit
     788              :       TYPE(section_vals_type), POINTER                   :: tddfpt_print_section
     789              :       TYPE(tddfpt_work_matrices), INTENT(inout)          :: work_matrices
     790              :       REAL(kind=dp)                                      :: conv
     791              : 
     792              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_davidson_solver'
     793              : 
     794              :       INTEGER                                            :: handle, ispin, istate, iter, &
     795              :                                                             max_krylov_vects, nspins, nstates, &
     796              :                                                             nstates_conv, nvects_exist, nvects_new
     797              :       INTEGER(kind=int_8)                                :: nstates_total
     798              :       LOGICAL                                            :: is_nonortho
     799              :       REAL(kind=dp)                                      :: t1, t2
     800         1166 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: evals_last
     801         1166 :       REAL(kind=dp), DIMENSION(:, :), POINTER            :: Atilde
     802         1166 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: Aop_krylov, Aop_ritz, krylov_vects, &
     803         1166 :                                                             S_krylov
     804         1166 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     805              : 
     806         1166 :       CALL timeset(routineN, handle)
     807              : 
     808         1166 :       nspins = SIZE(gs_mos)
     809         1166 :       nstates = tddfpt_control%nstates
     810         1166 :       nstates_total = tddfpt_total_number_of_states(gs_mos)
     811              : 
     812              :       IF (debug_this_module) THEN
     813              :          CPASSERT(SIZE(evects, 1) == nspins)
     814              :          CPASSERT(SIZE(evects, 2) == nstates)
     815              :          CPASSERT(SIZE(evals) == nstates)
     816              :       END IF
     817              : 
     818         1166 :       CALL get_qs_env(qs_env, matrix_s=matrix_s)
     819              : 
     820              :       ! adjust the number of Krylov vectors
     821         1166 :       max_krylov_vects = tddfpt_control%nkvs
     822         1166 :       IF (max_krylov_vects < nstates) max_krylov_vects = nstates
     823         1166 :       IF (INT(max_krylov_vects, kind=int_8) > nstates_total) max_krylov_vects = INT(nstates_total)
     824              : 
     825        11718 :       ALLOCATE (Aop_ritz(nspins, nstates))
     826         4504 :       DO istate = 1, nstates
     827         8220 :          DO ispin = 1, nspins
     828         7054 :             CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, Aop_ritz(ispin, istate))
     829              :          END DO
     830              :       END DO
     831              : 
     832         3498 :       ALLOCATE (evals_last(max_krylov_vects))
     833              :       ALLOCATE (Aop_krylov(nspins, max_krylov_vects), krylov_vects(nspins, max_krylov_vects), &
     834       828898 :                 S_krylov(nspins, max_krylov_vects))
     835              : 
     836         4504 :       DO istate = 1, nstates
     837         8220 :          DO ispin = 1, nspins
     838         3716 :             CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, krylov_vects(ispin, istate))
     839         3716 :             CALL cp_fm_to_fm(evects(ispin, istate), krylov_vects(ispin, istate))
     840              : 
     841         3716 :             CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, S_krylov(ispin, istate))
     842         3716 :             CALL cp_fm_to_fm(S_evects(ispin, istate), S_krylov(ispin, istate))
     843              : 
     844         7054 :             CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, Aop_krylov(ispin, istate))
     845              :          END DO
     846              :       END DO
     847              : 
     848         1166 :       nvects_exist = 0
     849         1166 :       nvects_new = nstates
     850              : 
     851         1166 :       t1 = m_walltime()
     852              : 
     853         1166 :       ALLOCATE (Atilde(1, 1))
     854              : 
     855         5256 :       DO
     856              :          ! davidson iteration
     857         5256 :          CALL cp_iterate(logger%iter_info, iter_nr_out=iter)
     858              : 
     859              :          CALL tddfpt_compute_Aop_evects(Aop_evects=Aop_krylov(:, nvects_exist + 1:nvects_exist + nvects_new), &
     860              :                                         evects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
     861              :                                         S_evects=S_krylov(:, nvects_exist + 1:nvects_exist + nvects_new), &
     862              :                                         gs_mos=gs_mos, tddfpt_control=tddfpt_control, &
     863              :                                         matrix_ks=matrix_ks, &
     864              :                                         qs_env=qs_env, kernel_env=kernel_env, &
     865              :                                         sub_env=sub_env, &
     866              :                                         work_matrices=work_matrices, &
     867         5256 :                                         matrix_s=matrix_s(1)%matrix)
     868              : 
     869              :          CALL tddfpt_compute_ritz_vects(ritz_vects=evects, Aop_ritz=Aop_ritz, &
     870              :                                         evals=evals_last(1:nvects_exist + nvects_new), &
     871              :                                         krylov_vects=krylov_vects(:, 1:nvects_exist + nvects_new), &
     872              :                                         Aop_krylov=Aop_krylov(:, 1:nvects_exist + nvects_new), &
     873         5256 :                                         Atilde=Atilde, nkvo=nvects_exist, nkvn=nvects_new)
     874              : 
     875              :          CALL tddfpt_write_restart(evects=evects, evals=evals_last(1:nstates), gs_mos=gs_mos, &
     876         5256 :                                    logger=logger, tddfpt_print_section=tddfpt_print_section)
     877              : 
     878        24634 :          conv = MAXVAL(ABS(evals_last(1:nstates) - evals(1:nstates)))
     879              : 
     880         5256 :          nvects_exist = nvects_exist + nvects_new
     881         5256 :          IF (nvects_exist + nvects_new > max_krylov_vects) &
     882          394 :             nvects_new = max_krylov_vects - nvects_exist
     883         5256 :          IF (iter >= tddfpt_control%niters) nvects_new = 0
     884              : 
     885         5256 :          IF (conv > tddfpt_control%conv .AND. nvects_new > 0) THEN
     886              :             ! compute residual vectors for the next iteration
     887        14202 :             DO istate = 1, nvects_new
     888        25960 :                DO ispin = 1, nspins
     889              :                   CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, &
     890        11758 :                                          krylov_vects(ispin, nvects_exist + istate))
     891              :                   CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, &
     892        11758 :                                          S_krylov(ispin, nvects_exist + istate))
     893              :                   CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, &
     894        21866 :                                          Aop_krylov(ispin, nvects_exist + istate))
     895              :                END DO
     896              :             END DO
     897              : 
     898              :             CALL tddfpt_compute_residual_vects(residual_vects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
     899              :                                                evals=evals_last(1:nvects_new), &
     900              :                                                ritz_vects=evects(:, 1:nvects_new), Aop_ritz=Aop_ritz(:, 1:nvects_new), &
     901         4094 :                                                gs_mos=gs_mos, matrix_s=matrix_s(1)%matrix)
     902              : 
     903              :             CALL tddfpt_orthogonalize_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
     904              :                                                 work_matrices%S_C0_C0T, qs_env, &
     905         4094 :                                                 gs_mos, evals(1:nstates), tddfpt_control, work_matrices%S_C0)
     906              : 
     907              :             CALL tddfpt_orthonormalize_psi1_psi1(krylov_vects(:, 1:nvects_exist + nvects_new), nvects_new, &
     908         4094 :                                                  S_krylov(:, 1:nvects_exist + nvects_new), matrix_s(1)%matrix)
     909              : 
     910              :             is_nonortho = tddfpt_is_nonorthogonal_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
     911         4094 :                                                             work_matrices%S_C0, tddfpt_control%orthogonal_eps)
     912              :          ELSE
     913              :             ! convergence or the maximum number of Krylov vectors have been achieved
     914         1162 :             nvects_new = 0
     915         1162 :             is_nonortho = .FALSE.
     916              :          END IF
     917              : 
     918         5256 :          t2 = m_walltime()
     919         5256 :          IF (energy_unit > 0) THEN
     920          315 :             WRITE (energy_unit, '(/,4X,A,T14,A,T36,A)') "State", "Exc. energy (eV)", "Convergence (eV)"
     921          727 :             DO istate = 1, nstates
     922          412 :                WRITE (energy_unit, '(1X,I8,T12,F14.7,T38,ES11.4)') istate, &
     923         1139 :                   evals_last(istate)*evolt, (evals_last(istate) - evals(istate))*evolt
     924              :             END DO
     925          315 :             WRITE (energy_unit, *)
     926          315 :             CALL m_flush(energy_unit)
     927              :          END IF
     928              : 
     929         5256 :          IF (iter_unit > 0) THEN
     930         2628 :             nstates_conv = 0
     931         9689 :             DO istate = 1, nstates
     932         7061 :                IF (ABS(evals_last(istate) - evals(istate)) <= tddfpt_control%conv) &
     933         4983 :                   nstates_conv = nstates_conv + 1
     934              :             END DO
     935              : 
     936         2628 :             WRITE (iter_unit, '(T7,I8,T24,F7.1,T40,ES11.4,T66,I8)') iter, t2 - t1, conv, nstates_conv
     937         2628 :             CALL m_flush(iter_unit)
     938              :          END IF
     939              : 
     940        19378 :          t1 = t2
     941        19378 :          evals(1:nstates) = evals_last(1:nstates)
     942              : 
     943              :          ! nvects_new == 0 if iter >= tddfpt_control%niters
     944         5256 :          IF (nvects_new == 0 .OR. is_nonortho) THEN
     945              :             ! restart Davidson iterations
     946              :             CALL tddfpt_orthogonalize_psi1_psi0(evects, work_matrices%S_C0_C0T, qs_env, &
     947              :                                                 gs_mos, &
     948         1166 :                                                 evals(1:nstates), tddfpt_control, work_matrices%S_C0)
     949         1166 :             CALL tddfpt_orthonormalize_psi1_psi1(evects, nstates, S_evects, matrix_s(1)%matrix)
     950              : 
     951              :             EXIT
     952              :          END IF
     953              :       END DO
     954              : 
     955         1166 :       DEALLOCATE (Atilde)
     956              : 
     957        14612 :       DO istate = nvects_exist + nvects_new, 1, -1
     958        30086 :          DO ispin = nspins, 1, -1
     959        15474 :             CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, Aop_krylov(ispin, istate))
     960        15474 :             CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, S_krylov(ispin, istate))
     961        28920 :             CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, krylov_vects(ispin, istate))
     962              :          END DO
     963              :       END DO
     964         1166 :       DEALLOCATE (Aop_krylov, krylov_vects, S_krylov)
     965         1166 :       DEALLOCATE (evals_last)
     966              : 
     967         4504 :       DO istate = nstates, 1, -1
     968         8220 :          DO ispin = nspins, 1, -1
     969         7054 :             CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, Aop_ritz(ispin, istate))
     970              :          END DO
     971              :       END DO
     972         1166 :       DEALLOCATE (Aop_ritz)
     973              : 
     974         1166 :       CALL timestop(handle)
     975              : 
     976         1166 :    END FUNCTION tddfpt_davidson_solver
     977              : 
     978              : END MODULE qs_tddfpt2_eigensolver
        

Generated by: LCOV version 2.0-1