LCOV - code coverage report
Current view: top level - src - qs_tddfpt2_operators.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 95.3 % 258 246
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 9 9

            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_operators
       9              :    USE admm_types,                      ONLY: admm_type
      10              :    USE cell_types,                      ONLY: cell_type,&
      11              :                                               pbc
      12              :    USE cp_dbcsr_api,                    ONLY: &
      13              :         dbcsr_create, dbcsr_filter, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
      14              :         dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
      15              :         dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
      16              :    USE cp_dbcsr_operations,             ONLY: copy_fm_to_dbcsr,&
      17              :                                               cp_dbcsr_plus_fm_fm_t,&
      18              :                                               cp_dbcsr_sm_fm_multiply
      19              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_column_scale,&
      20              :                                               cp_fm_scale_and_add
      21              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_type
      22              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      23              :                                               cp_fm_get_info,&
      24              :                                               cp_fm_release,&
      25              :                                               cp_fm_to_fm,&
      26              :                                               cp_fm_type
      27              :    USE hartree_local_methods,           ONLY: Vh_1c_gg_integrals
      28              :    USE hartree_local_types,             ONLY: hartree_local_type
      29              :    USE hfx_admm_utils,                  ONLY: tddft_hfx_matrix
      30              :    USE hfx_types,                       ONLY: hfx_type
      31              :    USE input_section_types,             ONLY: section_vals_get,&
      32              :                                               section_vals_get_subs_vals,&
      33              :                                               section_vals_type
      34              :    USE kinds,                           ONLY: dp
      35              :    USE message_passing,                 ONLY: mp_para_env_type
      36              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      37              :    USE particle_types,                  ONLY: particle_type
      38              :    USE pw_env_types,                    ONLY: pw_env_get,&
      39              :                                               pw_env_type
      40              :    USE pw_methods,                      ONLY: pw_axpy,&
      41              :                                               pw_multiply,&
      42              :                                               pw_scale,&
      43              :                                               pw_transfer,&
      44              :                                               pw_zero
      45              :    USE pw_poisson_methods,              ONLY: pw_poisson_solve
      46              :    USE pw_poisson_types,                ONLY: pw_poisson_type
      47              :    USE pw_pool_types,                   ONLY: pw_pool_p_type,&
      48              :                                               pw_pool_type
      49              :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      50              :                                               pw_r3d_rs_type
      51              :    USE qs_environment_types,            ONLY: get_qs_env,&
      52              :                                               qs_environment_type
      53              :    USE qs_kernel_types,                 ONLY: full_kernel_env_type
      54              :    USE qs_local_rho_types,              ONLY: local_rho_type
      55              :    USE qs_rho0_ggrid,                   ONLY: integrate_vhg0_rspace
      56              :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      57              :                                               qs_rho_type
      58              :    USE qs_tddfpt2_stda_utils,           ONLY: get_lowdin_x
      59              :    USE qs_tddfpt2_subgroups,            ONLY: tddfpt_subgroup_env_type
      60              :    USE qs_tddfpt2_types,                ONLY: tddfpt_ground_state_mos,&
      61              :                                               tddfpt_work_matrices
      62              :    USE realspace_grid_types,            ONLY: realspace_grid_desc_p_type,&
      63              :                                               realspace_grid_type
      64              :    USE xc,                              ONLY: xc_calc_2nd_deriv_analytical,&
      65              :                                               xc_calc_2nd_deriv_numerical
      66              :    USE xc_rho_set_types,                ONLY: xc_rho_set_type,&
      67              :                                               xc_rho_set_update
      68              : #include "./base/base_uses.f90"
      69              : 
      70              :    IMPLICIT NONE
      71              : 
      72              :    PRIVATE
      73              : 
      74              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_operators'
      75              : 
      76              :    LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .FALSE.
      77              :    ! number of first derivative components (3: d/dx, d/dy, d/dz)
      78              :    INTEGER, PARAMETER, PRIVATE          :: nderivs = 3
      79              :    INTEGER, PARAMETER, PRIVATE          :: maxspins = 2
      80              : 
      81              :    PUBLIC :: tddfpt_apply_energy_diff, tddfpt_apply_coulomb, tddfpt_apply_xc, tddfpt_apply_hfx, &
      82              :              tddfpt_apply_xc_potential, tddfpt_apply_hfxlr_kernel, tddfpt_apply_hfxsr_kernel
      83              : 
      84              : ! **************************************************************************************************
      85              : 
      86              : CONTAINS
      87              : 
      88              : ! **************************************************************************************************
      89              : !> \brief Apply orbital energy difference term:
      90              : !>        Aop_evects(spin,state) += KS(spin) * evects(spin,state) -
      91              : !>                                  S * evects(spin,state) * diag(evals_occ(spin))
      92              : !> \param Aop_evects  action of TDDFPT operator on trial vectors (modified on exit)
      93              : !> \param evects      trial vectors C_{1,i}
      94              : !> \param S_evects    S * C_{1,i}
      95              : !> \param gs_mos      molecular orbitals optimised for the ground state (only occupied orbital
      96              : !>                    energies [component %evals_occ] are needed)
      97              : !> \param matrix_ks   Kohn-Sham matrix
      98              : !> \par History
      99              : !>    * 05.2016 initialise all matrix elements in one go [Sergey Chulkov]
     100              : !>    * 03.2017 renamed from tddfpt_init_energy_diff(), altered prototype [Sergey Chulkov]
     101              : !> \note Based on the subroutine p_op_l1() which was originally created by
     102              : !>       Thomas Chassaing on 08.2002.
     103              : ! **************************************************************************************************
     104         5256 :    SUBROUTINE tddfpt_apply_energy_diff(Aop_evects, evects, S_evects, gs_mos, matrix_ks)
     105              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT)   :: Aop_evects
     106              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: evects, S_evects
     107              :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     108              :          INTENT(in)                                      :: gs_mos
     109              :       TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in)       :: matrix_ks
     110              : 
     111              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_energy_diff'
     112              : 
     113              :       INTEGER                                            :: handle, ispin, ivect, nactive, nao, &
     114              :                                                             nspins, nvects
     115              :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct
     116              :       TYPE(cp_fm_type)                                   :: hevec
     117              : 
     118         5256 :       CALL timeset(routineN, handle)
     119              : 
     120         5256 :       nspins = SIZE(evects, 1)
     121         5256 :       nvects = SIZE(evects, 2)
     122              : 
     123        11190 :       DO ispin = 1, nspins
     124              :          CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
     125         5934 :                              nrow_global=nao, ncol_global=nactive)
     126         5934 :          CALL cp_fm_create(hevec, matrix_struct)
     127              : 
     128        21388 :          DO ivect = 1, nvects
     129              :             CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, evects(ispin, ivect), &
     130              :                                          Aop_evects(ispin, ivect), ncol=nactive, &
     131        15454 :                                          alpha=1.0_dp, beta=1.0_dp)
     132              : 
     133        15454 :             IF (ASSOCIATED(gs_mos(ispin)%evals_occ_matrix)) THEN
     134              :                ! orbital energy correction: evals_occ_matrix is not a diagonal matrix
     135              :                CALL parallel_gemm('N', 'N', nao, nactive, nactive, 1.0_dp, &
     136              :                                   S_evects(ispin, ivect), gs_mos(ispin)%evals_occ_matrix, &
     137          756 :                                   0.0_dp, hevec)
     138              :             ELSE
     139        14698 :                CALL cp_fm_to_fm(S_evects(ispin, ivect), hevec)
     140        14698 :                CALL cp_fm_column_scale(hevec, gs_mos(ispin)%evals_occ)
     141              :             END IF
     142              : 
     143              :             ! KS * C1 - S * C1 * occupied_orbital_energies
     144        21388 :             CALL cp_fm_scale_and_add(1.0_dp, Aop_evects(ispin, ivect), -1.0_dp, hevec)
     145              :          END DO
     146        17124 :          CALL cp_fm_release(hevec)
     147              :       END DO
     148              : 
     149         5256 :       CALL timestop(handle)
     150              : 
     151         5256 :    END SUBROUTINE tddfpt_apply_energy_diff
     152              : 
     153              : ! **************************************************************************************************
     154              : !> \brief Update v_rspace by adding coulomb term.
     155              : !> \param A_ia_rspace    action of TDDFPT operator on the trial vector expressed in a plane wave
     156              : !>                       representation (modified on exit)
     157              : !> \param rho_ia_g       response density in reciprocal space for the given trial vector
     158              : !> \param local_rho_set ...
     159              : !> \param hartree_local ...
     160              : !> \param qs_env ...
     161              : !> \param sub_env        the full sub_environment needed for calculation
     162              : !> \param gapw           Flag indicating GAPW cacluation
     163              : !> \param work_v_gspace  work reciprocal-space grid to store Coulomb potential (modified on exit)
     164              : !> \param work_v_rspace  work real-space grid to store Coulomb potential (modified on exit)
     165              : !> \param tddfpt_mgrid ...
     166              : !> \par History
     167              : !>    * 05.2016 compute all coulomb terms in one go [Sergey Chulkov]
     168              : !>    * 03.2017 proceed excited states sequentially; minimise the number of conversions between
     169              : !>              DBCSR and FM matrices [Sergey Chulkov]
     170              : !>    * 06.2018 return the action expressed in the plane wave representation instead of the one
     171              : !>              in the atomic basis set representation
     172              : !> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by
     173              : !>       Mohamed Fawzi on 10.2002.
     174              : ! **************************************************************************************************
     175         6186 :    SUBROUTINE tddfpt_apply_coulomb(A_ia_rspace, rho_ia_g, local_rho_set, hartree_local, &
     176              :                                    qs_env, sub_env, gapw, work_v_gspace, work_v_rspace, tddfpt_mgrid)
     177              :       TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(INOUT)  :: A_ia_rspace
     178              :       TYPE(pw_c1d_gs_type), INTENT(INOUT)                :: rho_ia_g
     179              :       TYPE(local_rho_type), POINTER                      :: local_rho_set
     180              :       TYPE(hartree_local_type), POINTER                  :: hartree_local
     181              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     182              :       TYPE(tddfpt_subgroup_env_type), INTENT(in)         :: sub_env
     183              :       LOGICAL, INTENT(IN)                                :: gapw
     184              :       TYPE(pw_c1d_gs_type), INTENT(INOUT)                :: work_v_gspace
     185              :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                :: work_v_rspace
     186              :       LOGICAL, INTENT(IN)                                :: tddfpt_mgrid
     187              : 
     188              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_coulomb'
     189              : 
     190              :       INTEGER                                            :: handle, ispin, nspins
     191              :       REAL(kind=dp)                                      :: alpha, pair_energy
     192              :       TYPE(pw_env_type), POINTER                         :: pw_env
     193              :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
     194         6186 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: my_pools
     195              :       TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
     196         6186 :          POINTER                                         :: my_rs_descs
     197         6186 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: my_rs_grids
     198              : 
     199         6186 :       CALL timeset(routineN, handle)
     200              : 
     201         6186 :       nspins = SIZE(A_ia_rspace)
     202         6186 :       pw_env => sub_env%pw_env
     203         6186 :       IF (tddfpt_mgrid) THEN
     204              :          CALL pw_env_get(pw_env, poisson_env=poisson_env, rs_grids=my_rs_grids, &
     205           86 :                          rs_descs=my_rs_descs, pw_pools=my_pools)
     206              :       ELSE
     207         6100 :          CALL pw_env_get(pw_env, poisson_env=poisson_env)
     208              :       END IF
     209              : 
     210         6186 :       IF (nspins > 1) THEN
     211         1358 :          alpha = 1.0_dp
     212              :       ELSE
     213              :          ! spin-restricted case: alpha == 2 due to singlet state.
     214              :          ! In case of triplet states alpha == 0, so we should not call this subroutine at all.
     215         4828 :          alpha = 2.0_dp
     216              :       END IF
     217              : 
     218         6186 :       IF (gapw) THEN
     219          924 :          CPASSERT(ASSOCIATED(local_rho_set))
     220          924 :          CALL pw_axpy(local_rho_set%rho0_mpole%rho0_s_gs, rho_ia_g)
     221              :       END IF
     222              : 
     223         6186 :       CALL pw_poisson_solve(poisson_env, rho_ia_g, pair_energy, work_v_gspace)
     224         6186 :       CALL pw_transfer(work_v_gspace, work_v_rspace)
     225              : 
     226              :       ! (i a || j b) = ( i_alpha a_alpha + i_beta a_beta || j_alpha b_alpha + j_beta b_beta) =
     227              :       !                tr (Cj_alpha^T * [J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu] * Cb_alpha) +
     228              :       !                tr (Cj_beta^T * [J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu] * Cb_beta)
     229        13730 :       DO ispin = 1, nspins
     230        13730 :          CALL pw_axpy(work_v_rspace, A_ia_rspace(ispin), alpha)
     231              :       END DO
     232              : 
     233         6186 :       IF (gapw) THEN
     234              :          CALL Vh_1c_gg_integrals(qs_env, pair_energy, &
     235              :                                  hartree_local%ecoul_1c, &
     236              :                                  local_rho_set, &
     237          924 :                                  sub_env%para_env, tddft=.TRUE., core_2nd=.TRUE.)
     238          924 :          CALL pw_scale(work_v_rspace, work_v_rspace%pw_grid%dvol)
     239          924 :          IF (tddfpt_mgrid) THEN
     240              :             CALL integrate_vhg0_rspace(qs_env, work_v_rspace, sub_env%para_env, &
     241              :                                        calculate_forces=.FALSE., &
     242              :                                        local_rho_set=local_rho_set, my_pools=my_pools, &
     243           50 :                                        my_rs_descs=my_rs_descs)
     244              :          ELSE
     245              :             CALL integrate_vhg0_rspace(qs_env, work_v_rspace, sub_env%para_env, &
     246              :                                        calculate_forces=.FALSE., &
     247          874 :                                        local_rho_set=local_rho_set)
     248              :          END IF
     249              :       END IF
     250              : 
     251         6186 :       CALL timestop(handle)
     252              : 
     253         6186 :    END SUBROUTINE tddfpt_apply_coulomb
     254              : 
     255              : ! **************************************************************************************************
     256              : !> \brief Driver routine for applying fxc (analyic vs. finite difference for testing
     257              : !> \param A_ia_rspace      action of TDDFPT operator on trial vectors expressed in a plane wave
     258              : !>                         representation (modified on exit)
     259              : !> \param kernel_env       kernel environment
     260              : !> \param rho_ia_struct    response density for the given trial vector
     261              : !> \param is_rks_triplets  indicates that the triplet excited states calculation using
     262              : !>                         spin-unpolarised molecular orbitals has been requested
     263              : !> \param pw_env           plain wave environment
     264              : !> \param work_v_xc        work real-space grid to store the gradient of the exchange-correlation
     265              : !>                         potential with respect to the response density (modified on exit)
     266              : !> \param work_v_xc_tau ...
     267              : ! **************************************************************************************************
     268         7360 :    SUBROUTINE tddfpt_apply_xc(A_ia_rspace, kernel_env, rho_ia_struct, is_rks_triplets, &
     269              :                               pw_env, work_v_xc, work_v_xc_tau)
     270              : 
     271              :       TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(INOUT)  :: A_ia_rspace
     272              :       TYPE(full_kernel_env_type), INTENT(IN)             :: kernel_env
     273              :       TYPE(qs_rho_type), POINTER                         :: rho_ia_struct
     274              :       LOGICAL, INTENT(in)                                :: is_rks_triplets
     275              :       TYPE(pw_env_type), POINTER                         :: pw_env
     276              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: work_v_xc, work_v_xc_tau
     277              : 
     278              :       INTEGER                                            :: ispin, nspins
     279              : 
     280         7360 :       nspins = SIZE(A_ia_rspace)
     281              : 
     282         7360 :       IF (kernel_env%deriv2_analytic) THEN
     283              :          CALL tddfpt_apply_xc_analytic(kernel_env, rho_ia_struct, is_rks_triplets, nspins, &
     284         7320 :                                        pw_env, work_v_xc, work_v_xc_tau)
     285              :       ELSE
     286              :          CALL tddfpt_apply_xc_fd(kernel_env, rho_ia_struct, is_rks_triplets, nspins, &
     287           40 :                                  pw_env, work_v_xc, work_v_xc_tau)
     288              :       END IF
     289              : 
     290        16078 :       DO ispin = 1, nspins
     291              :          ! pw2 = pw2 + alpha * pw1
     292        16078 :          CALL pw_axpy(work_v_xc(ispin), A_ia_rspace(ispin), kernel_env%alpha)
     293              :       END DO
     294              : 
     295         7360 :    END SUBROUTINE tddfpt_apply_xc
     296              : 
     297              : ! **************************************************************************************************
     298              : !> \brief Routine for applying fxc potential
     299              : !> \param A_ia_rspace      action of TDDFPT operator on trial vectors expressed in a plane wave
     300              : !>                         representation (modified on exit)
     301              : !> \param fxc_rspace ...
     302              : !> \param rho_ia_struct    response density for the given trial vector
     303              : !> \param is_rks_triplets ...
     304              : ! **************************************************************************************************
     305          156 :    SUBROUTINE tddfpt_apply_xc_potential(A_ia_rspace, fxc_rspace, rho_ia_struct, is_rks_triplets)
     306              : 
     307              :       TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(INOUT)  :: A_ia_rspace
     308              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: fxc_rspace
     309              :       TYPE(qs_rho_type), POINTER                         :: rho_ia_struct
     310              :       LOGICAL, INTENT(in)                                :: is_rks_triplets
     311              : 
     312              :       INTEGER                                            :: nspins
     313              :       REAL(KIND=dp)                                      :: alpha
     314          156 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho1_r
     315              : 
     316          156 :       nspins = SIZE(A_ia_rspace)
     317              : 
     318          156 :       alpha = 1.0_dp
     319              : 
     320          156 :       CALL qs_rho_get(rho_ia_struct, rho_r=rho1_r)
     321              : 
     322          156 :       IF (nspins == 2) THEN
     323            0 :          CALL pw_multiply(A_ia_rspace(1), fxc_rspace(1), rho1_r(1), alpha)
     324            0 :          CALL pw_multiply(A_ia_rspace(1), fxc_rspace(2), rho1_r(2), alpha)
     325            0 :          CALL pw_multiply(A_ia_rspace(2), fxc_rspace(3), rho1_r(2), alpha)
     326            0 :          CALL pw_multiply(A_ia_rspace(2), fxc_rspace(2), rho1_r(1), alpha)
     327          156 :       ELSE IF (is_rks_triplets) THEN
     328            0 :          CALL pw_multiply(A_ia_rspace(1), fxc_rspace(1), rho1_r(1), alpha)
     329            0 :          CALL pw_multiply(A_ia_rspace(1), fxc_rspace(2), rho1_r(1), -alpha)
     330              :       ELSE
     331          156 :          CALL pw_multiply(A_ia_rspace(1), fxc_rspace(1), rho1_r(1), alpha)
     332          156 :          CALL pw_multiply(A_ia_rspace(1), fxc_rspace(2), rho1_r(1), alpha)
     333              :       END IF
     334              : 
     335          156 :    END SUBROUTINE tddfpt_apply_xc_potential
     336              : 
     337              : ! **************************************************************************************************
     338              : !> \brief Update A_ia_munu by adding exchange-correlation term.
     339              : !> \param kernel_env       kernel environment
     340              : !> \param rho_ia_struct    response density for the given trial vector
     341              : !> \param is_rks_triplets  indicates that the triplet excited states calculation using
     342              : !>                         spin-unpolarised molecular orbitals has been requested
     343              : !> \param nspins ...
     344              : !> \param pw_env           plain wave environment
     345              : !> \param work_v_xc        work real-space grid to store the gradient of the exchange-correlation
     346              : !>                         potential with respect to the response density (modified on exit)
     347              : !> \param work_v_xc_tau ...
     348              : !> \par History
     349              : !>    * 05.2016 compute all kernel terms in one go [Sergey Chulkov]
     350              : !>    * 03.2017 proceed excited states sequentially; minimise the number of conversions between
     351              : !>              DBCSR and FM matrices [Sergey Chulkov]
     352              : !>    * 06.2018 return the action expressed in the plane wave representation instead of the one
     353              : !>              in the atomic basis set representation
     354              : !> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by
     355              : !>       Mohamed Fawzi on 10.2002.
     356              : ! **************************************************************************************************
     357         7320 :    SUBROUTINE tddfpt_apply_xc_analytic(kernel_env, rho_ia_struct, is_rks_triplets, nspins, &
     358              :                                        pw_env, work_v_xc, work_v_xc_tau)
     359              :       TYPE(full_kernel_env_type), INTENT(in)             :: kernel_env
     360              :       TYPE(qs_rho_type), POINTER                         :: rho_ia_struct
     361              :       LOGICAL, INTENT(in)                                :: is_rks_triplets
     362              :       INTEGER, INTENT(in)                                :: nspins
     363              :       TYPE(pw_env_type), POINTER                         :: pw_env
     364              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: work_v_xc, work_v_xc_tau
     365              : 
     366              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_xc_analytic'
     367              : 
     368              :       INTEGER                                            :: handle, ispin
     369         7320 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_ia_g, rho_ia_g2
     370              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     371         7320 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_ia_r, rho_ia_r2, tau_ia_r, tau_ia_r2
     372              : 
     373         7320 :       CALL timeset(routineN, handle)
     374              : 
     375         7320 :       CALL qs_rho_get(rho_ia_struct, rho_g=rho_ia_g, rho_r=rho_ia_r, tau_r=tau_ia_r)
     376         7320 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     377              : 
     378              :       IF (debug_this_module) THEN
     379              :          CPASSERT(SIZE(rho_ia_g) == nspins)
     380              :          CPASSERT(SIZE(rho_ia_r) == nspins)
     381              :          CPASSERT((.NOT. ASSOCIATED(tau_ia_r)) .OR. SIZE(tau_ia_r) == nspins)
     382              :          CPASSERT((.NOT. is_rks_triplets) .OR. nspins == 1)
     383              :       END IF
     384              : 
     385         7320 :       NULLIFY (tau_ia_r2)
     386         7320 :       IF (is_rks_triplets) THEN
     387         1512 :          ALLOCATE (rho_ia_r2(2))
     388         1512 :          ALLOCATE (rho_ia_g2(2))
     389          504 :          rho_ia_r2(1) = rho_ia_r(1)
     390          504 :          rho_ia_r2(2) = rho_ia_r(1)
     391          504 :          rho_ia_g2(1) = rho_ia_g(1)
     392          504 :          rho_ia_g2(2) = rho_ia_g(1)
     393              : 
     394          504 :          IF (ASSOCIATED(tau_ia_r)) THEN
     395            0 :             ALLOCATE (tau_ia_r2(2))
     396            0 :             tau_ia_r2(1) = tau_ia_r(1)
     397            0 :             tau_ia_r2(2) = tau_ia_r(1)
     398              :          END IF
     399              :       ELSE
     400         6816 :          rho_ia_r2 => rho_ia_r
     401         6816 :          rho_ia_g2 => rho_ia_g
     402              : 
     403         6816 :          tau_ia_r2 => tau_ia_r
     404              :       END IF
     405              : 
     406        15978 :       DO ispin = 1, nspins
     407         8658 :          CALL pw_zero(work_v_xc(ispin))
     408        15978 :          IF (ASSOCIATED(work_v_xc_tau)) CALL pw_zero(work_v_xc_tau(ispin))
     409              :       END DO
     410              : 
     411              :       CALL xc_rho_set_update(rho_set=kernel_env%xc_rho1_set, rho_r=rho_ia_r2, rho_g=rho_ia_g2, tau=tau_ia_r2, &
     412              :                              needs=kernel_env%xc_rho1_cflags, xc_deriv_method_id=kernel_env%deriv_method_id, &
     413         7320 :                              xc_rho_smooth_id=kernel_env%rho_smooth_id, pw_pool=auxbas_pw_pool)
     414              : 
     415              :       CALL xc_calc_2nd_deriv_analytical(v_xc=work_v_xc, v_xc_tau=work_v_xc_tau, deriv_set=kernel_env%xc_deriv_set, &
     416              :                                         rho_set=kernel_env%xc_rho_set, &
     417              :                                         rho1_set=kernel_env%xc_rho1_set, pw_pool=auxbas_pw_pool, &
     418         7320 :                                         xc_section=kernel_env%xc_section, gapw=.FALSE., tddfpt_fac=kernel_env%beta)
     419              : 
     420         7320 :       IF (is_rks_triplets) THEN
     421          504 :          DEALLOCATE (rho_ia_r2)
     422          504 :          DEALLOCATE (rho_ia_g2)
     423          504 :          IF (ASSOCIATED(tau_ia_r2)) DEALLOCATE (tau_ia_r2)
     424              :       END IF
     425              : 
     426         7320 :       CALL timestop(handle)
     427              : 
     428         7320 :    END SUBROUTINE tddfpt_apply_xc_analytic
     429              : 
     430              : ! **************************************************************************************************
     431              : !> \brief Update A_ia_munu by adding exchange-correlation term using finite difference methods.
     432              : !> \param kernel_env       kernel environment
     433              : !> \param rho_ia_struct    response density for the given trial vector
     434              : !> \param is_rks_triplets  indicates that the triplet excited states calculation using
     435              : !>                         spin-unpolarised molecular orbitals has been requested
     436              : !> \param nspins ...
     437              : !> \param pw_env           plain wave environment
     438              : !> \param work_v_xc        work real-space grid to store the gradient of the exchange-correlation
     439              : !>                         potential with respect to the response density (modified on exit)
     440              : !> \param work_v_xc_tau ...
     441              : ! **************************************************************************************************
     442           40 :    SUBROUTINE tddfpt_apply_xc_fd(kernel_env, rho_ia_struct, is_rks_triplets, nspins, &
     443              :                                  pw_env, work_v_xc, work_v_xc_tau)
     444              :       TYPE(full_kernel_env_type), INTENT(in)             :: kernel_env
     445              :       TYPE(qs_rho_type), POINTER                         :: rho_ia_struct
     446              :       LOGICAL, INTENT(in)                                :: is_rks_triplets
     447              :       INTEGER, INTENT(in)                                :: nspins
     448              :       TYPE(pw_env_type), POINTER                         :: pw_env
     449              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: work_v_xc, work_v_xc_tau
     450              : 
     451              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_xc_fd'
     452              : 
     453              :       INTEGER                                            :: handle, ispin
     454              :       LOGICAL                                            :: lsd, singlet, triplet
     455           40 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho1_g
     456              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     457           40 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho1_r, tau1_r
     458              :       TYPE(xc_rho_set_type), POINTER                     :: rho_set
     459              : 
     460           40 :       CALL timeset(routineN, handle)
     461              : 
     462           40 :       CALL qs_rho_get(rho_ia_struct, rho_r=rho1_r, rho_g=rho1_g, tau_r=tau1_r)
     463           40 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     464          100 :       DO ispin = 1, nspins
     465          100 :          CALL pw_zero(work_v_xc(ispin))
     466              :       END DO
     467           40 :       rho_set => kernel_env%xc_rho_set
     468              : 
     469           40 :       singlet = .FALSE.
     470           40 :       triplet = .FALSE.
     471           40 :       lsd = .FALSE.
     472           40 :       IF (nspins == 1 .AND. .NOT. is_rks_triplets) THEN
     473           40 :          singlet = .TRUE.
     474           40 :       ELSE IF (nspins == 1 .AND. is_rks_triplets) THEN
     475           40 :          triplet = .TRUE.
     476           20 :       ELSE IF (nspins == 2) THEN
     477           40 :          lsd = .TRUE.
     478              :       ELSE
     479            0 :          CPABORT("illegal options")
     480              :       END IF
     481              : 
     482           40 :       IF (ASSOCIATED(tau1_r)) THEN
     483            0 :          DO ispin = 1, nspins
     484            0 :             CALL pw_zero(work_v_xc_tau(ispin))
     485              :          END DO
     486              :       END IF
     487              : 
     488              :       CALL xc_calc_2nd_deriv_numerical(work_v_xc, work_v_xc_tau, rho_set, rho1_r, rho1_g, tau1_r, &
     489              :                                        auxbas_pw_pool, kernel_env%xc_section, &
     490           40 :                                        is_rks_triplets)
     491              : 
     492           40 :       CALL timestop(handle)
     493              : 
     494           40 :    END SUBROUTINE tddfpt_apply_xc_fd
     495              : 
     496              : ! **************************************************************************************************
     497              : !> \brief Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
     498              : !> \param Aop_evects      action of TDDFPT operator on trial vectors (modified on exit)
     499              : !> \param evects          trial vectors
     500              : !> \param gs_mos          molecular orbitals optimised for the ground state (only occupied
     501              : !>                        molecular orbitals [component %mos_occ] are needed)
     502              : !> \param do_admm         perform auxiliary density matrix method calculations
     503              : !> \param qs_env          Quickstep environment
     504              : !> \param work_rho_ia_ao_symm ...
     505              : !> \param work_hmat_symm ...
     506              : !> \param work_rho_ia_ao_asymm ...
     507              : !> \param work_hmat_asymm ...
     508              : !> \param wfm_rho_orb ...
     509              : !> \par History
     510              : !>    * 05.2016 compute all exact-exchange terms in one go [Sergey Chulkov]
     511              : !>    * 03.2017 code related to ADMM correction is now moved to tddfpt_apply_admm_correction()
     512              : !>              in order to compute this correction within parallel groups [Sergey Chulkov]
     513              : !> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by
     514              : !>       Mohamed Fawzi on 10.2002.
     515              : ! **************************************************************************************************
     516         1146 :    SUBROUTINE tddfpt_apply_hfx(Aop_evects, evects, gs_mos, do_admm, qs_env, &
     517         1146 :                                work_rho_ia_ao_symm, work_hmat_symm, work_rho_ia_ao_asymm, &
     518         1146 :                                work_hmat_asymm, wfm_rho_orb)
     519              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT)   :: Aop_evects
     520              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: evects
     521              :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     522              :          INTENT(in)                                      :: gs_mos
     523              :       LOGICAL, INTENT(in)                                :: do_admm
     524              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     525              :       TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT)    :: work_rho_ia_ao_symm
     526              :       TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
     527              :          TARGET                                          :: work_hmat_symm
     528              :       TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT)    :: work_rho_ia_ao_asymm
     529              :       TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
     530              :          TARGET                                          :: work_hmat_asymm
     531              :       TYPE(cp_fm_type), INTENT(IN)                       :: wfm_rho_orb
     532              : 
     533              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'tddfpt_apply_hfx'
     534              : 
     535              :       INTEGER                                            :: handle, ispin, ivect, nao, nao_aux, &
     536              :                                                             nspins, nvects
     537              :       INTEGER, DIMENSION(maxspins)                       :: nactive
     538              :       LOGICAL                                            :: do_hfx
     539              :       REAL(kind=dp)                                      :: alpha
     540              :       TYPE(admm_type), POINTER                           :: admm_env
     541              :       TYPE(section_vals_type), POINTER                   :: hfx_section, input
     542              : 
     543         1146 :       CALL timeset(routineN, handle)
     544              : 
     545              :       ! Check for hfx section
     546         1146 :       CALL get_qs_env(qs_env, input=input)
     547         1146 :       hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
     548         1146 :       CALL section_vals_get(hfx_section, explicit=do_hfx)
     549              : 
     550         1146 :       IF (do_hfx) THEN
     551         1146 :          nspins = SIZE(evects, 1)
     552         1146 :          nvects = SIZE(evects, 2)
     553              : 
     554         1146 :          IF (nspins > 1) THEN
     555           78 :             alpha = 1.0_dp
     556              :          ELSE
     557         1068 :             alpha = 2.0_dp
     558              :          END IF
     559              : 
     560         1146 :          CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao)
     561         2370 :          DO ispin = 1, nspins
     562         2370 :             CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive(ispin))
     563              :          END DO
     564              : 
     565         1146 :          IF (do_admm) THEN
     566          622 :             CALL get_qs_env(qs_env, admm_env=admm_env)
     567          622 :             CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux)
     568              :          END IF
     569              : 
     570              :          !Note: the symmetrized transition density matrix is P = 0.5*(C*evect^T + evect*C^T)
     571              :          !      in the end, we only want evect*C^T for consistency with the MO formulation of TDDFT
     572              :          !      therefore, we go in 2 steps: with the symmetric 0.5*(C*evect^T + evect*C^T) and
     573              :          !      the antisymemtric 0.5*(C*evect^T - evect*C^T)
     574              : 
     575              :          ! some stuff from qs_ks_build_kohn_sham_matrix
     576              :          ! TO DO: add SIC support
     577         3528 :          DO ivect = 1, nvects
     578         4892 :             DO ispin = 1, nspins
     579              : 
     580              :                !The symmetric density matrix
     581              :                CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, evects(ispin, ivect), &
     582         2510 :                                   gs_mos(ispin)%mos_occ, 0.0_dp, wfm_rho_orb)
     583              :                CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, gs_mos(ispin)%mos_occ, &
     584         2510 :                                   evects(ispin, ivect), 1.0_dp, wfm_rho_orb)
     585              : 
     586         2510 :                CALL dbcsr_set(work_hmat_symm(ispin)%matrix, 0.0_dp)
     587         4892 :                IF (do_admm) THEN
     588              :                   CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, &
     589         1320 :                                      wfm_rho_orb, 0.0_dp, admm_env%work_aux_orb)
     590              :                   CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
     591         1320 :                                      0.0_dp, admm_env%work_aux_aux)
     592         1320 :                   CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_rho_ia_ao_symm(ispin)%matrix, keep_sparsity=.TRUE.)
     593              :                ELSE
     594         1190 :                   CALL copy_fm_to_dbcsr(wfm_rho_orb, work_rho_ia_ao_symm(ispin)%matrix, keep_sparsity=.TRUE.)
     595              :                END IF
     596              :             END DO
     597              : 
     598         2382 :             CALL tddft_hfx_matrix(work_hmat_symm, work_rho_ia_ao_symm, qs_env)
     599              : 
     600         2382 :             IF (do_admm) THEN
     601         2612 :                DO ispin = 1, nspins
     602              :                   CALL cp_dbcsr_sm_fm_multiply(work_hmat_symm(ispin)%matrix, admm_env%A, admm_env%work_aux_orb, &
     603         1320 :                                                ncol=nao, alpha=1.0_dp, beta=0.0_dp)
     604              : 
     605              :                   CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, &
     606         1320 :                                      admm_env%work_aux_orb, 0.0_dp, wfm_rho_orb)
     607              : 
     608              :                   CALL parallel_gemm('N', 'N', nao, nactive(ispin), nao, alpha, wfm_rho_orb, &
     609         2612 :                                      gs_mos(ispin)%mos_occ, 1.0_dp, Aop_evects(ispin, ivect))
     610              :                END DO
     611              :             ELSE
     612         2280 :                DO ispin = 1, nspins
     613              :                   CALL cp_dbcsr_sm_fm_multiply(work_hmat_symm(ispin)%matrix, gs_mos(ispin)%mos_occ, &
     614              :                                                Aop_evects(ispin, ivect), ncol=nactive(ispin), &
     615         2280 :                                                alpha=alpha, beta=1.0_dp)
     616              :                END DO
     617              :             END IF
     618              : 
     619              :             !The anti-symmetric density matrix
     620         4892 :             DO ispin = 1, nspins
     621              : 
     622              :                !The symmetric density matrix
     623              :                CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, evects(ispin, ivect), &
     624         2510 :                                   gs_mos(ispin)%mos_occ, 0.0_dp, wfm_rho_orb)
     625              :                CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), -0.5_dp, gs_mos(ispin)%mos_occ, &
     626         2510 :                                   evects(ispin, ivect), 1.0_dp, wfm_rho_orb)
     627              : 
     628         2510 :                CALL dbcsr_set(work_hmat_asymm(ispin)%matrix, 0.0_dp)
     629         4892 :                IF (do_admm) THEN
     630              :                   CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, &
     631         1320 :                                      wfm_rho_orb, 0.0_dp, admm_env%work_aux_orb)
     632              :                   CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
     633         1320 :                                      0.0_dp, admm_env%work_aux_aux)
     634         1320 :                   CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_rho_ia_ao_asymm(ispin)%matrix, keep_sparsity=.TRUE.)
     635              :                ELSE
     636         1190 :                   CALL copy_fm_to_dbcsr(wfm_rho_orb, work_rho_ia_ao_asymm(ispin)%matrix, keep_sparsity=.TRUE.)
     637              :                END IF
     638              :             END DO
     639              : 
     640         2382 :             CALL tddft_hfx_matrix(work_hmat_asymm, work_rho_ia_ao_asymm, qs_env)
     641              : 
     642         3528 :             IF (do_admm) THEN
     643         2612 :                DO ispin = 1, nspins
     644              :                   CALL cp_dbcsr_sm_fm_multiply(work_hmat_asymm(ispin)%matrix, admm_env%A, admm_env%work_aux_orb, &
     645         1320 :                                                ncol=nao, alpha=1.0_dp, beta=0.0_dp)
     646              : 
     647              :                   CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, &
     648         1320 :                                      admm_env%work_aux_orb, 0.0_dp, wfm_rho_orb)
     649              : 
     650              :                   CALL parallel_gemm('N', 'N', nao, nactive(ispin), nao, alpha, wfm_rho_orb, &
     651         2612 :                                      gs_mos(ispin)%mos_occ, 1.0_dp, Aop_evects(ispin, ivect))
     652              :                END DO
     653              :             ELSE
     654         2280 :                DO ispin = 1, nspins
     655              :                   CALL cp_dbcsr_sm_fm_multiply(work_hmat_asymm(ispin)%matrix, gs_mos(ispin)%mos_occ, &
     656              :                                                Aop_evects(ispin, ivect), ncol=nactive(ispin), &
     657         2280 :                                                alpha=alpha, beta=1.0_dp)
     658              :                END DO
     659              :             END IF
     660              :          END DO
     661              :       END IF
     662              : 
     663         1146 :       CALL timestop(handle)
     664              : 
     665         1146 :    END SUBROUTINE tddfpt_apply_hfx
     666              : 
     667              : ! **************************************************************************************************
     668              : !> \brief Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
     669              : !> \param Aop_evects      action of TDDFPT operator on trial vectors (modified on exit)
     670              : !> \param evects          trial vectors
     671              : !> \param gs_mos          molecular orbitals optimised for the ground state (only occupied
     672              : !>                        molecular orbitals [component %mos_occ] are needed)
     673              : !> \param qs_env          Quickstep environment
     674              : !> \param admm_env ...
     675              : !> \param hfx_section ...
     676              : !> \param x_data ...
     677              : !> \param symmetry ...
     678              : !> \param recalc_integrals ...
     679              : !> \param work_rho_ia_ao ...
     680              : !> \param work_hmat ...
     681              : !> \param wfm_rho_orb ...
     682              : ! **************************************************************************************************
     683           44 :    SUBROUTINE tddfpt_apply_hfxsr_kernel(Aop_evects, evects, gs_mos, qs_env, admm_env, &
     684              :                                         hfx_section, x_data, symmetry, recalc_integrals, &
     685           44 :                                         work_rho_ia_ao, work_hmat, wfm_rho_orb)
     686              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in)      :: Aop_evects, evects
     687              :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     688              :          INTENT(in)                                      :: gs_mos
     689              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     690              :       TYPE(admm_type), POINTER                           :: admm_env
     691              :       TYPE(section_vals_type), POINTER                   :: hfx_section
     692              :       TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
     693              :       INTEGER, INTENT(IN)                                :: symmetry
     694              :       LOGICAL, INTENT(IN)                                :: recalc_integrals
     695              :       TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT)    :: work_rho_ia_ao
     696              :       TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
     697              :          TARGET                                          :: work_hmat
     698              :       TYPE(cp_fm_type), INTENT(IN)                       :: wfm_rho_orb
     699              : 
     700              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_hfxsr_kernel'
     701              : 
     702              :       INTEGER                                            :: handle, ispin, ivect, nao, nao_aux, &
     703              :                                                             nspins, nvects
     704              :       INTEGER, DIMENSION(maxspins)                       :: nactive
     705              :       LOGICAL                                            :: reint
     706              :       REAL(kind=dp)                                      :: alpha
     707              : 
     708           44 :       CALL timeset(routineN, handle)
     709              : 
     710           44 :       nspins = SIZE(evects, 1)
     711           44 :       nvects = SIZE(evects, 2)
     712              : 
     713           44 :       alpha = 2.0_dp
     714           44 :       IF (nspins > 1) alpha = 1.0_dp
     715              : 
     716           44 :       CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao)
     717           44 :       CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux)
     718           88 :       DO ispin = 1, nspins
     719           88 :          CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive(ispin))
     720              :       END DO
     721              : 
     722           44 :       reint = recalc_integrals
     723              : 
     724          132 :       DO ivect = 1, nvects
     725          176 :          DO ispin = 1, nspins
     726              :             CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, evects(ispin, ivect), &
     727           88 :                                gs_mos(ispin)%mos_occ, 0.0_dp, wfm_rho_orb)
     728              :             CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp*symmetry, gs_mos(ispin)%mos_occ, &
     729           88 :                                evects(ispin, ivect), 1.0_dp, wfm_rho_orb)
     730           88 :             CALL dbcsr_set(work_hmat(ispin)%matrix, 0.0_dp)
     731              :             CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, &
     732           88 :                                wfm_rho_orb, 0.0_dp, admm_env%work_aux_orb)
     733              :             CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
     734           88 :                                0.0_dp, admm_env%work_aux_aux)
     735          176 :             CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_rho_ia_ao(ispin)%matrix, keep_sparsity=.TRUE.)
     736              :          END DO
     737              : 
     738           88 :          CALL tddft_hfx_matrix(work_hmat, work_rho_ia_ao, qs_env, .FALSE., reint, hfx_section, x_data)
     739           88 :          reint = .FALSE.
     740              : 
     741          220 :          DO ispin = 1, nspins
     742              :             CALL cp_dbcsr_sm_fm_multiply(work_hmat(ispin)%matrix, admm_env%A, admm_env%work_aux_orb, &
     743           88 :                                          ncol=nao, alpha=1.0_dp, beta=0.0_dp)
     744              :             CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, &
     745           88 :                                admm_env%work_aux_orb, 0.0_dp, wfm_rho_orb)
     746              :             CALL parallel_gemm('N', 'N', nao, nactive(ispin), nao, alpha, wfm_rho_orb, &
     747          176 :                                gs_mos(ispin)%mos_occ, 1.0_dp, Aop_evects(ispin, ivect))
     748              :          END DO
     749              :       END DO
     750              : 
     751           44 :       CALL timestop(handle)
     752              : 
     753           44 :    END SUBROUTINE tddfpt_apply_hfxsr_kernel
     754              : 
     755              : ! **************************************************************************************************
     756              : !> \brief ...Calculate the HFXLR kernel contribution by contracting the Lowdin MO coefficients --
     757              : !>           transition charges with the exchange-type integrals using the sTDA approximation
     758              : !> \param qs_env ...
     759              : !> \param sub_env ...
     760              : !> \param rcut ...
     761              : !> \param hfx_scale ...
     762              : !> \param work ...
     763              : !> \param X ...
     764              : !> \param res ... vector AX with A being the sTDA matrix and X the Davidson trial vector of the
     765              : !>                eigenvalue problem A*X = omega*X
     766              : ! **************************************************************************************************
     767           72 :    SUBROUTINE tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, hfx_scale, work, X, res)
     768              : 
     769              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     770              :       TYPE(tddfpt_subgroup_env_type)                     :: sub_env
     771              :       REAL(KIND=dp), INTENT(IN)                          :: rcut, hfx_scale
     772              :       TYPE(tddfpt_work_matrices)                         :: work
     773              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: X
     774              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT)      :: res
     775              : 
     776              :       CHARACTER(len=*), PARAMETER :: routineN = 'tddfpt_apply_hfxlr_kernel'
     777              : 
     778              :       INTEGER                                            :: handle, iatom, ispin, jatom, natom, &
     779              :                                                             nsgf, nspins
     780              :       INTEGER, DIMENSION(2)                              :: nactive
     781              :       REAL(KIND=dp)                                      :: dr, eps_filter, fcut, gabr
     782              :       REAL(KIND=dp), DIMENSION(3)                        :: rij
     783           72 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pblock
     784              :       TYPE(cell_type), POINTER                           :: cell
     785              :       TYPE(cp_fm_struct_type), POINTER                   :: fmstruct
     786              :       TYPE(cp_fm_type)                                   :: cvec
     787           72 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: xtransformed
     788              :       TYPE(cp_fm_type), POINTER                          :: ct
     789              :       TYPE(dbcsr_iterator_type)                          :: iter
     790              :       TYPE(dbcsr_type)                                   :: pdens
     791              :       TYPE(dbcsr_type), POINTER                          :: tempmat
     792              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     793           72 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     794              : 
     795           72 :       CALL timeset(routineN, handle)
     796              : 
     797              :       ! parameters
     798           72 :       eps_filter = 1.E-08_dp
     799              : 
     800           72 :       nspins = SIZE(X)
     801          144 :       DO ispin = 1, nspins
     802          144 :          CALL cp_fm_get_info(X(ispin), ncol_global=nactive(ispin))
     803              :       END DO
     804              : 
     805           72 :       para_env => sub_env%para_env
     806              : 
     807           72 :       CALL get_qs_env(qs_env, natom=natom, cell=cell, particle_set=particle_set)
     808              : 
     809              :       ! calculate Loewdin transformed Davidson trial vector tilde(X)=S^1/2*X
     810              :       ! and tilde(tilde(X))=S^1/2_A*tilde(X)_A
     811          288 :       ALLOCATE (xtransformed(nspins))
     812          144 :       DO ispin = 1, nspins
     813           72 :          NULLIFY (fmstruct)
     814           72 :          ct => work%ctransformed(ispin)
     815           72 :          CALL cp_fm_get_info(ct, matrix_struct=fmstruct)
     816          144 :          CALL cp_fm_create(matrix=xtransformed(ispin), matrix_struct=fmstruct, name="XTRANSFORMED")
     817              :       END DO
     818           72 :       CALL get_lowdin_x(work%shalf, X, xtransformed)
     819              : 
     820          144 :       DO ispin = 1, nspins
     821           72 :          ct => work%ctransformed(ispin)
     822           72 :          CALL cp_fm_get_info(ct, matrix_struct=fmstruct, nrow_global=nsgf)
     823           72 :          CALL cp_fm_create(cvec, fmstruct)
     824              :          !
     825           72 :          tempmat => work%shalf
     826           72 :          CALL dbcsr_create(pdens, template=tempmat, matrix_type=dbcsr_type_no_symmetry)
     827              :          ! P(nu,mu) = SUM_j XT(nu,j)*CT(mu,j)
     828           72 :          ct => work%ctransformed(ispin)
     829           72 :          CALL dbcsr_set(pdens, 0.0_dp)
     830              :          CALL cp_dbcsr_plus_fm_fm_t(pdens, xtransformed(ispin), ct, nactive(ispin), &
     831           72 :                                     1.0_dp, keep_sparsity=.FALSE.)
     832           72 :          CALL dbcsr_filter(pdens, eps_filter)
     833              :          ! Apply PP*gab -> PP; gab = gamma_coulomb
     834              :          ! P(nu,mu) = P(nu,mu)*g(a of nu,b of mu)
     835           72 :          CALL dbcsr_iterator_start(iter, pdens)
     836          396 :          DO WHILE (dbcsr_iterator_blocks_left(iter))
     837          324 :             CALL dbcsr_iterator_next_block(iter, iatom, jatom, pblock)
     838         1296 :             rij = particle_set(iatom)%r - particle_set(jatom)%r
     839         1296 :             rij = pbc(rij, cell)
     840         1296 :             dr = SQRT(SUM(rij(:)**2))
     841          324 :             gabr = 1._dp/rcut
     842          324 :             IF (dr < 1.e-6) THEN
     843          108 :                gabr = 2._dp*gabr/SQRT(3.1415926_dp)
     844              :             ELSE
     845          216 :                gabr = ERF(gabr*dr)/dr
     846              :                fcut = EXP(dr - 4._dp*rcut)
     847          216 :                fcut = fcut/(fcut + 1._dp)
     848              :             END IF
     849        21924 :             pblock = hfx_scale*gabr*pblock
     850              :          END DO
     851           72 :          CALL dbcsr_iterator_stop(iter)
     852              :          ! CV(mu,i) = P(nu,mu)*CT(mu,i)
     853           72 :          CALL cp_dbcsr_sm_fm_multiply(pdens, ct, cvec, nactive(ispin), 1.0_dp, 0.0_dp)
     854              :          ! rho(nu,i) = rho(nu,i) + ShalfP(nu,mu)*CV(mu,i)
     855              :          CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, res(ispin), nactive(ispin), &
     856           72 :                                       -1.0_dp, 1.0_dp)
     857              :          !
     858           72 :          CALL dbcsr_release(pdens)
     859              :          !
     860          288 :          CALL cp_fm_release(cvec)
     861              :       END DO
     862              : 
     863           72 :       CALL cp_fm_release(xtransformed)
     864              : 
     865           72 :       CALL timestop(handle)
     866              : 
     867          144 :    END SUBROUTINE tddfpt_apply_hfxlr_kernel
     868              : 
     869              : ! **************************************************************************************************
     870              : 
     871              : END MODULE qs_tddfpt2_operators
        

Generated by: LCOV version 2.0-1