LCOV - code coverage report
Current view: top level - src - qs_tddfpt2_forces.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 576 600 96.0 %
Date: 2024-03-28 07:31:50 Functions: 9 9 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : MODULE qs_tddfpt2_forces
       9             :    USE admm_types,                      ONLY: admm_type,&
      10             :                                               get_admm_env
      11             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      12             :                                               get_atomic_kind,&
      13             :                                               get_atomic_kind_set
      14             :    USE cp_control_types,                ONLY: dft_control_type,&
      15             :                                               tddfpt2_control_type
      16             :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      17             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      18             :                                               copy_fm_to_dbcsr,&
      19             :                                               cp_dbcsr_plus_fm_fm_t,&
      20             :                                               cp_dbcsr_sm_fm_multiply,&
      21             :                                               dbcsr_allocate_matrix_set,&
      22             :                                               dbcsr_deallocate_matrix_set
      23             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      24             :                                               cp_fm_struct_release,&
      25             :                                               cp_fm_struct_type
      26             :    USE cp_fm_types,                     ONLY: cp_fm_copy_general,&
      27             :                                               cp_fm_create,&
      28             :                                               cp_fm_get_info,&
      29             :                                               cp_fm_release,&
      30             :                                               cp_fm_set_all,&
      31             :                                               cp_fm_type
      32             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      33             :                                               cp_logger_get_default_unit_nr,&
      34             :                                               cp_logger_type
      35             :    USE dbcsr_api,                       ONLY: &
      36             :         dbcsr_add, dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, dbcsr_dot, dbcsr_p_type, &
      37             :         dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric
      38             :    USE exstates_types,                  ONLY: excited_energy_type,&
      39             :                                               exstate_potential_release
      40             :    USE hartree_local_methods,           ONLY: Vh_1c_gg_integrals,&
      41             :                                               init_coulomb_local
      42             :    USE hartree_local_types,             ONLY: hartree_local_create,&
      43             :                                               hartree_local_release,&
      44             :                                               hartree_local_type
      45             :    USE hfx_energy_potential,            ONLY: integrate_four_center
      46             :    USE hfx_ri,                          ONLY: hfx_ri_update_ks
      47             :    USE hfx_types,                       ONLY: hfx_type
      48             :    USE input_constants,                 ONLY: do_admm_aux_exch_func_none,&
      49             :                                               oe_shift,&
      50             :                                               tddfpt_kernel_full,&
      51             :                                               tddfpt_kernel_none,&
      52             :                                               tddfpt_kernel_stda
      53             :    USE input_section_types,             ONLY: section_get_lval,&
      54             :                                               section_vals_get,&
      55             :                                               section_vals_get_subs_vals,&
      56             :                                               section_vals_type,&
      57             :                                               section_vals_val_get
      58             :    USE kinds,                           ONLY: default_string_length,&
      59             :                                               dp
      60             :    USE message_passing,                 ONLY: mp_para_env_type
      61             :    USE mulliken,                        ONLY: ao_charges
      62             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      63             :    USE particle_types,                  ONLY: particle_type
      64             :    USE pw_env_types,                    ONLY: pw_env_get,&
      65             :                                               pw_env_type
      66             :    USE pw_methods,                      ONLY: pw_axpy,&
      67             :                                               pw_scale,&
      68             :                                               pw_transfer,&
      69             :                                               pw_zero
      70             :    USE pw_poisson_methods,              ONLY: pw_poisson_solve
      71             :    USE pw_poisson_types,                ONLY: pw_poisson_type
      72             :    USE pw_pool_types,                   ONLY: pw_pool_type
      73             :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      74             :                                               pw_r3d_rs_type
      75             :    USE qs_collocate_density,            ONLY: calculate_rho_elec
      76             :    USE qs_density_matrices,             ONLY: calculate_wx_matrix,&
      77             :                                               calculate_xwx_matrix
      78             :    USE qs_environment_types,            ONLY: get_qs_env,&
      79             :                                               qs_environment_type,&
      80             :                                               set_qs_env
      81             :    USE qs_force_types,                  ONLY: allocate_qs_force,&
      82             :                                               deallocate_qs_force,&
      83             :                                               qs_force_type,&
      84             :                                               sum_qs_force,&
      85             :                                               total_qs_force,&
      86             :                                               zero_qs_force
      87             :    USE qs_fxc,                          ONLY: qs_fxc_analytic,&
      88             :                                               qs_fxc_fdiff
      89             :    USE qs_gapw_densities,               ONLY: prepare_gapw_den
      90             :    USE qs_integrate_potential,          ONLY: integrate_v_rspace
      91             :    USE qs_kernel_types,                 ONLY: kernel_env_type
      92             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      93             :                                               get_qs_kind_set,&
      94             :                                               qs_kind_type
      95             :    USE qs_ks_atom,                      ONLY: update_ks_atom
      96             :    USE qs_ks_reference,                 ONLY: ks_ref_potential,&
      97             :                                               ks_ref_potential_atom
      98             :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      99             :    USE qs_local_rho_types,              ONLY: local_rho_set_create,&
     100             :                                               local_rho_set_release,&
     101             :                                               local_rho_type
     102             :    USE qs_mo_types,                     ONLY: get_mo_set,&
     103             :                                               mo_set_type
     104             :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
     105             :    USE qs_oce_types,                    ONLY: oce_matrix_type
     106             :    USE qs_overlap,                      ONLY: build_overlap_matrix
     107             :    USE qs_rho0_ggrid,                   ONLY: integrate_vhg0_rspace,&
     108             :                                               rho0_s_grid_create
     109             :    USE qs_rho0_methods,                 ONLY: init_rho0
     110             :    USE qs_rho0_types,                   ONLY: get_rho0_mpole
     111             :    USE qs_rho_atom_methods,             ONLY: allocate_rho_atom_internals,&
     112             :                                               calculate_rho_atom_coeff
     113             :    USE qs_rho_atom_types,               ONLY: rho_atom_type
     114             :    USE qs_rho_types,                    ONLY: qs_rho_create,&
     115             :                                               qs_rho_get,&
     116             :                                               qs_rho_set,&
     117             :                                               qs_rho_type
     118             :    USE qs_tddfpt2_fhxc_forces,          ONLY: fhxc_force,&
     119             :                                               stda_force
     120             :    USE qs_tddfpt2_subgroups,            ONLY: tddfpt_subgroup_env_type
     121             :    USE qs_tddfpt2_types,                ONLY: tddfpt_ground_state_mos,&
     122             :                                               tddfpt_work_matrices
     123             :    USE qs_vxc_atom,                     ONLY: calculate_xc_2nd_deriv_atom
     124             :    USE task_list_types,                 ONLY: task_list_type
     125             :    USE xtb_ehess,                       ONLY: xtb_coulomb_hessian
     126             :    USE xtb_types,                       ONLY: get_xtb_atom_param,&
     127             :                                               xtb_atom_type
     128             : #include "./base/base_uses.f90"
     129             : 
     130             :    IMPLICIT NONE
     131             : 
     132             :    PRIVATE
     133             : 
     134             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_forces'
     135             : 
     136             :    PUBLIC :: tddfpt_forces_main
     137             : 
     138             : ! **************************************************************************************************
     139             : 
     140             : CONTAINS
     141             : 
     142             : ! **************************************************************************************************
     143             : !> \brief Perform TDDFPT gradient calculation.
     144             : !> \param qs_env  Quickstep environment
     145             : !> \param gs_mos ...
     146             : !> \param ex_env ...
     147             : !> \param kernel_env ...
     148             : !> \param sub_env ...
     149             : !> \param work_matrices ...
     150             : !> \par History
     151             : !>    * 10.2022 created JHU
     152             : ! **************************************************************************************************
     153         566 :    SUBROUTINE tddfpt_forces_main(qs_env, gs_mos, ex_env, kernel_env, sub_env, work_matrices)
     154             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     155             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     156             :          POINTER                                         :: gs_mos
     157             :       TYPE(excited_energy_type), POINTER                 :: ex_env
     158             :       TYPE(kernel_env_type)                              :: kernel_env
     159             :       TYPE(tddfpt_subgroup_env_type)                     :: sub_env
     160             :       TYPE(tddfpt_work_matrices)                         :: work_matrices
     161             : 
     162             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_forces_main'
     163             : 
     164             :       INTEGER                                            :: handle, ispin, nspins
     165             :       TYPE(admm_type), POINTER                           :: admm_env
     166             :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct
     167         566 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_pe_asymm, matrix_pe_symm, &
     168         566 :                                                             matrix_s, matrix_s_aux_fit
     169             :       TYPE(dft_control_type), POINTER                    :: dft_control
     170             :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     171             : 
     172         566 :       CALL timeset(routineN, handle)
     173             : 
     174         566 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     175         566 :       nspins = dft_control%nspins
     176         566 :       tddfpt_control => dft_control%tddfpt2_control
     177             :       ! rhs of linres equation
     178         566 :       IF (ASSOCIATED(ex_env%cpmos)) THEN
     179         504 :          DO ispin = 1, SIZE(ex_env%cpmos)
     180         504 :             CALL cp_fm_release(ex_env%cpmos(ispin))
     181             :          END DO
     182         230 :          DEALLOCATE (ex_env%cpmos)
     183             :       END IF
     184        2360 :       ALLOCATE (ex_env%cpmos(nspins))
     185        1228 :       DO ispin = 1, nspins
     186         662 :          CALL cp_fm_get_info(matrix=ex_env%evect(ispin), matrix_struct=matrix_struct)
     187         662 :          CALL cp_fm_create(ex_env%cpmos(ispin), matrix_struct)
     188        1228 :          CALL cp_fm_set_all(ex_env%cpmos(ispin), 0.0_dp)
     189             :       END DO
     190         566 :       CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
     191         566 :       NULLIFY (matrix_pe_asymm, matrix_pe_symm)
     192         566 :       CALL dbcsr_allocate_matrix_set(ex_env%matrix_pe, nspins)
     193         566 :       CALL dbcsr_allocate_matrix_set(matrix_pe_symm, nspins)
     194         566 :       CALL dbcsr_allocate_matrix_set(matrix_pe_asymm, nspins)
     195        1228 :       DO ispin = 1, nspins
     196         662 :          ALLOCATE (ex_env%matrix_pe(ispin)%matrix)
     197         662 :          CALL dbcsr_create(ex_env%matrix_pe(ispin)%matrix, template=matrix_s(1)%matrix)
     198         662 :          CALL dbcsr_copy(ex_env%matrix_pe(ispin)%matrix, matrix_s(1)%matrix)
     199         662 :          CALL dbcsr_set(ex_env%matrix_pe(ispin)%matrix, 0.0_dp)
     200             : 
     201         662 :          ALLOCATE (matrix_pe_symm(ispin)%matrix)
     202         662 :          CALL dbcsr_create(matrix_pe_symm(ispin)%matrix, template=matrix_s(1)%matrix)
     203         662 :          CALL dbcsr_copy(matrix_pe_symm(ispin)%matrix, ex_env%matrix_pe(ispin)%matrix)
     204             : 
     205         662 :          ALLOCATE (matrix_pe_asymm(ispin)%matrix)
     206             :          CALL dbcsr_create(matrix_pe_asymm(ispin)%matrix, template=matrix_s(1)%matrix, &
     207         662 :                            matrix_type=dbcsr_type_antisymmetric)
     208         662 :          CALL dbcsr_complete_redistribute(ex_env%matrix_pe(ispin)%matrix, matrix_pe_asymm(ispin)%matrix)
     209             : 
     210             :          CALL tddfpt_resvec1(ex_env%evect(ispin), gs_mos(ispin)%mos_occ, &
     211        1228 :                              matrix_s(1)%matrix, ex_env%matrix_pe(ispin)%matrix)
     212             :       END DO
     213             :       !
     214             :       ! ground state ADMM!
     215         566 :       IF (dft_control%do_admm) THEN
     216         128 :          CALL get_qs_env(qs_env, admm_env=admm_env)
     217         128 :          CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit)
     218         128 :          CALL dbcsr_allocate_matrix_set(ex_env%matrix_pe_admm, nspins)
     219         276 :          DO ispin = 1, nspins
     220         148 :             ALLOCATE (ex_env%matrix_pe_admm(ispin)%matrix)
     221         148 :             CALL dbcsr_create(ex_env%matrix_pe_admm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
     222         148 :             CALL dbcsr_copy(ex_env%matrix_pe_admm(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
     223         148 :             CALL dbcsr_set(ex_env%matrix_pe_admm(ispin)%matrix, 0.0_dp)
     224             :             CALL tddfpt_resvec1_admm(ex_env%matrix_pe(ispin)%matrix, &
     225         276 :                                      admm_env, ex_env%matrix_pe_admm(ispin)%matrix)
     226             :          END DO
     227             :       END IF
     228             :       !
     229         566 :       CALL dbcsr_allocate_matrix_set(ex_env%matrix_hz, nspins)
     230        1228 :       DO ispin = 1, nspins
     231         662 :          ALLOCATE (ex_env%matrix_hz(ispin)%matrix)
     232         662 :          CALL dbcsr_create(ex_env%matrix_hz(ispin)%matrix, template=matrix_s(1)%matrix)
     233         662 :          CALL dbcsr_copy(ex_env%matrix_hz(ispin)%matrix, matrix_s(1)%matrix)
     234        1228 :          CALL dbcsr_set(ex_env%matrix_hz(ispin)%matrix, 0.0_dp)
     235             :       END DO
     236         566 :       IF (dft_control%qs_control%xtb) THEN
     237          16 :          CALL tddfpt_resvec2_xtb(qs_env, ex_env%matrix_pe, gs_mos, ex_env%matrix_hz, ex_env%cpmos)
     238             :       ELSE
     239             :          CALL tddfpt_resvec2(qs_env, ex_env%matrix_pe, ex_env%matrix_pe_admm, &
     240         550 :                              gs_mos, ex_env%matrix_hz, ex_env%cpmos)
     241             :       END IF
     242             :       !
     243         566 :       CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1, nspins)
     244         566 :       CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_asymm, nspins)
     245        1228 :       DO ispin = 1, nspins
     246         662 :          ALLOCATE (ex_env%matrix_px1(ispin)%matrix)
     247         662 :          CALL dbcsr_create(ex_env%matrix_px1(ispin)%matrix, template=matrix_s(1)%matrix)
     248         662 :          CALL dbcsr_copy(ex_env%matrix_px1(ispin)%matrix, matrix_s(1)%matrix)
     249         662 :          CALL dbcsr_set(ex_env%matrix_px1(ispin)%matrix, 0.0_dp)
     250             : 
     251         662 :          ALLOCATE (ex_env%matrix_px1_asymm(ispin)%matrix)
     252             :          CALL dbcsr_create(ex_env%matrix_px1_asymm(ispin)%matrix, template=matrix_s(1)%matrix, &
     253         662 :                            matrix_type=dbcsr_type_antisymmetric)
     254        1228 :          CALL dbcsr_complete_redistribute(ex_env%matrix_px1(ispin)%matrix, ex_env%matrix_px1_asymm(ispin)%matrix)
     255             :       END DO
     256             :       ! Kernel ADMM
     257         566 :       IF (tddfpt_control%do_admm) THEN
     258          64 :          CALL get_qs_env(qs_env, admm_env=admm_env)
     259          64 :          CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit)
     260          64 :          CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_admm, nspins)
     261          64 :          CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_admm_asymm, nspins)
     262         132 :          DO ispin = 1, nspins
     263          68 :             ALLOCATE (ex_env%matrix_px1_admm(ispin)%matrix)
     264          68 :             CALL dbcsr_create(ex_env%matrix_px1_admm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
     265          68 :             CALL dbcsr_copy(ex_env%matrix_px1_admm(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
     266          68 :             CALL dbcsr_set(ex_env%matrix_px1_admm(ispin)%matrix, 0.0_dp)
     267             : 
     268          68 :             ALLOCATE (ex_env%matrix_px1_admm_asymm(ispin)%matrix)
     269             :             CALL dbcsr_create(ex_env%matrix_px1_admm_asymm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix, &
     270          68 :                               matrix_type=dbcsr_type_antisymmetric)
     271             :             CALL dbcsr_complete_redistribute(ex_env%matrix_px1_admm(ispin)%matrix, &
     272         132 :                                              ex_env%matrix_px1_admm_asymm(ispin)%matrix)
     273             :          END DO
     274             :       END IF
     275             :       ! TDA forces
     276         566 :       CALL tddfpt_forces(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices)
     277             :       ! Rotate res vector cpmos into original frame of occupied orbitals
     278         566 :       CALL tddfpt_resvec3(qs_env, ex_env%cpmos, work_matrices)
     279             : 
     280         566 :       CALL dbcsr_deallocate_matrix_set(matrix_pe_symm)
     281         566 :       CALL dbcsr_deallocate_matrix_set(matrix_pe_asymm)
     282             : 
     283         566 :       CALL timestop(handle)
     284             : 
     285         566 :    END SUBROUTINE tddfpt_forces_main
     286             : 
     287             : ! **************************************************************************************************
     288             : !> \brief Calculate direct tddft forces
     289             : !> \param qs_env ...
     290             : !> \param ex_env ...
     291             : !> \param gs_mos ...
     292             : !> \param kernel_env ...
     293             : !> \param sub_env ...
     294             : !> \param work_matrices ...
     295             : !> \par History
     296             : !>    * 01.2020 screated [JGH]
     297             : ! **************************************************************************************************
     298         566 :    SUBROUTINE tddfpt_forces(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices)
     299             : 
     300             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     301             :       TYPE(excited_energy_type), POINTER                 :: ex_env
     302             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     303             :          POINTER                                         :: gs_mos
     304             :       TYPE(kernel_env_type), INTENT(IN)                  :: kernel_env
     305             :       TYPE(tddfpt_subgroup_env_type)                     :: sub_env
     306             :       TYPE(tddfpt_work_matrices)                         :: work_matrices
     307             : 
     308             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'tddfpt_forces'
     309             : 
     310             :       INTEGER                                            :: handle
     311         566 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: natom_of_kind
     312             :       LOGICAL                                            :: debug_forces
     313             :       REAL(KIND=dp)                                      :: ehartree, exc
     314         566 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     315             :       TYPE(dft_control_type), POINTER                    :: dft_control
     316         566 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: ks_force, td_force
     317             : 
     318         566 :       CALL timeset(routineN, handle)
     319             : 
     320             :       ! for extended debug output
     321         566 :       debug_forces = ex_env%debug_forces
     322             :       ! prepare force array
     323             :       CALL get_qs_env(qs_env, dft_control=dft_control, force=ks_force, &
     324         566 :                       atomic_kind_set=atomic_kind_set)
     325         566 :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind)
     326         566 :       NULLIFY (td_force)
     327         566 :       CALL allocate_qs_force(td_force, natom_of_kind)
     328         566 :       DEALLOCATE (natom_of_kind)
     329         566 :       CALL zero_qs_force(td_force)
     330         566 :       CALL set_qs_env(qs_env, force=td_force)
     331             :       !
     332         566 :       IF (dft_control%qs_control%xtb) THEN
     333             :          CALL tddfpt_force_direct(qs_env, ex_env, gs_mos, kernel_env, sub_env, &
     334          16 :                                   work_matrices, debug_forces)
     335             :       ELSE
     336             :          !
     337         550 :          CALL exstate_potential_release(ex_env)
     338             :          CALL ks_ref_potential(qs_env, ex_env%vh_rspace, ex_env%vxc_rspace, &
     339         550 :                                ex_env%vtau_rspace, ex_env%vadmm_rspace, ehartree, exc)
     340             :          CALL ks_ref_potential_atom(qs_env, ex_env%local_rho_set, ex_env%local_rho_set_admm, &
     341         550 :                                     ex_env%vh_rspace)
     342             :          CALL tddfpt_force_direct(qs_env, ex_env, gs_mos, kernel_env, sub_env, &
     343         550 :                                   work_matrices, debug_forces)
     344             :       END IF
     345             :       !
     346             :       ! add TD and KS forces
     347         566 :       CALL get_qs_env(qs_env, force=td_force)
     348         566 :       CALL sum_qs_force(ks_force, td_force)
     349         566 :       CALL set_qs_env(qs_env, force=ks_force)
     350         566 :       CALL deallocate_qs_force(td_force)
     351             :       !
     352         566 :       CALL timestop(handle)
     353             : 
     354         566 :    END SUBROUTINE tddfpt_forces
     355             : 
     356             : ! **************************************************************************************************
     357             : !> \brief Calculate direct tddft forces
     358             : !> \param qs_env ...
     359             : !> \param ex_env ...
     360             : !> \param gs_mos ...
     361             : !> \param kernel_env ...
     362             : !> \param sub_env ...
     363             : !> \param work_matrices ...
     364             : !> \param debug_forces ...
     365             : !> \par History
     366             : !>    * 01.2020 screated [JGH]
     367             : ! **************************************************************************************************
     368         566 :    SUBROUTINE tddfpt_force_direct(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices, &
     369             :                                   debug_forces)
     370             : 
     371             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     372             :       TYPE(excited_energy_type), POINTER                 :: ex_env
     373             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     374             :          POINTER                                         :: gs_mos
     375             :       TYPE(kernel_env_type), INTENT(IN)                  :: kernel_env
     376             :       TYPE(tddfpt_subgroup_env_type)                     :: sub_env
     377             :       TYPE(tddfpt_work_matrices)                         :: work_matrices
     378             :       LOGICAL                                            :: debug_forces
     379             : 
     380             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_force_direct'
     381             : 
     382             :       INTEGER                                            :: handle, iounit, ispin, natom, norb, &
     383             :                                                             nspins
     384             :       REAL(KIND=dp)                                      :: evalue
     385         566 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: ftot1, ftot2
     386             :       REAL(KIND=dp), DIMENSION(3)                        :: fodeb
     387         566 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     388         566 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: evect
     389             :       TYPE(cp_logger_type), POINTER                      :: logger
     390         566 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_s, matrix_wx1, &
     391         566 :                                                             matrix_wz, scrm
     392             :       TYPE(dft_control_type), POINTER                    :: dft_control
     393             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     394             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     395         566 :          POINTER                                         :: sab_orb
     396         566 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     397             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     398             :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     399             : 
     400         566 :       CALL timeset(routineN, handle)
     401             : 
     402         566 :       logger => cp_get_default_logger()
     403         566 :       IF (logger%para_env%is_source()) THEN
     404         283 :          iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     405             :       ELSE
     406             :          iounit = -1
     407             :       END IF
     408             : 
     409         566 :       evect => ex_env%evect
     410             : 
     411             :       CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, para_env=para_env, &
     412         566 :                       sab_orb=sab_orb, dft_control=dft_control, force=force)
     413         566 :       NULLIFY (tddfpt_control)
     414         566 :       tddfpt_control => dft_control%tddfpt2_control
     415         566 :       nspins = dft_control%nspins
     416             : 
     417         566 :       IF (debug_forces) THEN
     418          56 :          CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
     419         168 :          ALLOCATE (ftot1(3, natom))
     420          56 :          CALL total_qs_force(ftot1, force, atomic_kind_set)
     421             :       END IF
     422             : 
     423         566 :       CALL tddfpt_kernel_force(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices, debug_forces)
     424             : 
     425             :       ! Overlap matrix
     426         566 :       matrix_wx1 => ex_env%matrix_wx1
     427         566 :       CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, matrix_ks=matrix_ks)
     428         566 :       NULLIFY (matrix_wz)
     429         566 :       CALL dbcsr_allocate_matrix_set(matrix_wz, nspins)
     430        1228 :       DO ispin = 1, nspins
     431         662 :          ALLOCATE (matrix_wz(ispin)%matrix)
     432         662 :          CALL dbcsr_create(matrix=matrix_wz(ispin)%matrix, template=matrix_s(1)%matrix)
     433         662 :          CALL cp_dbcsr_alloc_block_from_nbl(matrix_wz(ispin)%matrix, sab_orb)
     434         662 :          CALL dbcsr_set(matrix_wz(ispin)%matrix, 0.0_dp)
     435         662 :          CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=norb)
     436         662 :          CALL cp_dbcsr_plus_fm_fm_t(matrix_wz(ispin)%matrix, matrix_v=evect(ispin), ncol=norb)
     437         662 :          evalue = ex_env%evalue
     438         662 :          IF (tddfpt_control%oe_corr == oe_shift) THEN
     439           4 :             evalue = ex_env%evalue - tddfpt_control%ev_shift
     440             :          END IF
     441         662 :          CALL dbcsr_scale(matrix_wz(ispin)%matrix, evalue)
     442             :          CALL calculate_wx_matrix(gs_mos(ispin)%mos_occ, evect(ispin), matrix_ks(ispin)%matrix, &
     443        1890 :                                   matrix_wz(ispin)%matrix)
     444             :       END DO
     445         566 :       IF (nspins == 2) THEN
     446             :          CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
     447          96 :                         alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
     448             :       END IF
     449         566 :       NULLIFY (scrm)
     450         734 :       IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
     451             :       CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
     452             :                                 matrix_name="OVERLAP MATRIX", &
     453             :                                 basis_type_a="ORB", basis_type_b="ORB", &
     454             :                                 sab_nl=sab_orb, calculate_forces=.TRUE., &
     455         566 :                                 matrix_p=matrix_wz(1)%matrix)
     456         566 :       CALL dbcsr_deallocate_matrix_set(scrm)
     457         566 :       CALL dbcsr_deallocate_matrix_set(matrix_wz)
     458         566 :       IF (debug_forces) THEN
     459         224 :          fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
     460          56 :          CALL para_env%sum(fodeb)
     461          56 :          IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wx*dS ", fodeb
     462             :       END IF
     463             : 
     464             :       ! Overlap matrix
     465         566 :       CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, matrix_ks=matrix_ks)
     466         566 :       NULLIFY (matrix_wz)
     467         566 :       CALL dbcsr_allocate_matrix_set(matrix_wz, nspins)
     468        1228 :       DO ispin = 1, nspins
     469         662 :          ALLOCATE (matrix_wz(ispin)%matrix)
     470         662 :          CALL dbcsr_create(matrix=matrix_wz(ispin)%matrix, template=matrix_s(1)%matrix)
     471         662 :          CALL cp_dbcsr_alloc_block_from_nbl(matrix_wz(ispin)%matrix, sab_orb)
     472         662 :          CALL dbcsr_set(matrix_wz(ispin)%matrix, 0.0_dp)
     473         662 :          CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=norb)
     474         662 :          evalue = ex_env%evalue
     475         662 :          IF (tddfpt_control%oe_corr == oe_shift) THEN
     476           4 :             evalue = ex_env%evalue - tddfpt_control%ev_shift
     477             :          END IF
     478             :          CALL calculate_xwx_matrix(gs_mos(ispin)%mos_occ, evect(ispin), matrix_s(1)%matrix, &
     479        1890 :                                    matrix_ks(ispin)%matrix, matrix_wz(ispin)%matrix, evalue)
     480             :       END DO
     481         566 :       IF (nspins == 2) THEN
     482             :          CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
     483          96 :                         alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
     484             :       END IF
     485         566 :       NULLIFY (scrm)
     486         734 :       IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
     487             :       CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
     488             :                                 matrix_name="OVERLAP MATRIX", &
     489             :                                 basis_type_a="ORB", basis_type_b="ORB", &
     490             :                                 sab_nl=sab_orb, calculate_forces=.TRUE., &
     491         566 :                                 matrix_p=matrix_wz(1)%matrix)
     492         566 :       CALL dbcsr_deallocate_matrix_set(scrm)
     493         566 :       CALL dbcsr_deallocate_matrix_set(matrix_wz)
     494         566 :       IF (debug_forces) THEN
     495         224 :          fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
     496          56 :          CALL para_env%sum(fodeb)
     497          56 :          IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: xWx*dS ", fodeb
     498             :       END IF
     499             : 
     500             :       ! Overlap matrix
     501         566 :       IF (ASSOCIATED(matrix_wx1)) THEN
     502         498 :          IF (nspins == 2) THEN
     503             :             CALL dbcsr_add(matrix_wx1(1)%matrix, matrix_wx1(2)%matrix, &
     504          96 :                            alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
     505             :          END IF
     506         498 :          NULLIFY (scrm)
     507         648 :          IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
     508             :          CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
     509             :                                    matrix_name="OVERLAP MATRIX", &
     510             :                                    basis_type_a="ORB", basis_type_b="ORB", &
     511             :                                    sab_nl=sab_orb, calculate_forces=.TRUE., &
     512         498 :                                    matrix_p=matrix_wx1(1)%matrix)
     513         498 :          CALL dbcsr_deallocate_matrix_set(scrm)
     514         498 :          IF (debug_forces) THEN
     515         200 :             fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
     516          50 :             CALL para_env%sum(fodeb)
     517          50 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: WK*dS ", fodeb
     518             :          END IF
     519             :       END IF
     520             : 
     521         566 :       IF (debug_forces) THEN
     522         168 :          ALLOCATE (ftot2(3, natom))
     523          56 :          CALL total_qs_force(ftot2, force, atomic_kind_set)
     524         224 :          fodeb(1:3) = ftot2(1:3, 1) - ftot1(1:3, 1)
     525          56 :          CALL para_env%sum(fodeb)
     526          56 :          IF (iounit > 0) WRITE (iounit, "(T3,A,T30,3F16.8)") "DEBUG:: Excitation Force", fodeb
     527          56 :          DEALLOCATE (ftot1, ftot2)
     528             :       END IF
     529             : 
     530         566 :       CALL timestop(handle)
     531             : 
     532        1132 :    END SUBROUTINE tddfpt_force_direct
     533             : 
     534             : ! **************************************************************************************************
     535             : !> \brief ...
     536             : !> \param evect ...
     537             : !> \param mos_occ ...
     538             : !> \param matrix_s ...
     539             : !> \param matrix_pe ...
     540             : ! **************************************************************************************************
     541        1324 :    SUBROUTINE tddfpt_resvec1(evect, mos_occ, matrix_s, matrix_pe)
     542             : 
     543             :       TYPE(cp_fm_type), INTENT(IN)                       :: evect, mos_occ
     544             :       TYPE(dbcsr_type), POINTER                          :: matrix_s, matrix_pe
     545             : 
     546             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'tddfpt_resvec1'
     547             : 
     548             :       INTEGER                                            :: handle, iounit, nao, norb
     549             :       REAL(KIND=dp)                                      :: tmp
     550             :       TYPE(cp_fm_struct_type), POINTER                   :: fmstruct, fmstruct2
     551             :       TYPE(cp_fm_type)                                   :: cxmat, xxmat
     552             :       TYPE(cp_logger_type), POINTER                      :: logger
     553             : 
     554         662 :       CALL timeset(routineN, handle)
     555             :       ! X*X^T
     556         662 :       CALL cp_fm_get_info(mos_occ, nrow_global=nao, ncol_global=norb)
     557         662 :       CALL cp_dbcsr_plus_fm_fm_t(matrix_pe, matrix_v=evect, ncol=norb)
     558             :       ! X^T*S*X
     559         662 :       CALL cp_fm_get_info(evect, matrix_struct=fmstruct)
     560         662 :       NULLIFY (fmstruct2)
     561             :       CALL cp_fm_struct_create(fmstruct=fmstruct2, template_fmstruct=fmstruct, &
     562         662 :                                nrow_global=norb, ncol_global=norb)
     563         662 :       CALL cp_fm_create(xxmat, matrix_struct=fmstruct2)
     564         662 :       CALL cp_fm_struct_release(fmstruct2)
     565         662 :       CALL cp_fm_create(cxmat, matrix_struct=fmstruct)
     566         662 :       CALL cp_dbcsr_sm_fm_multiply(matrix_s, evect, cxmat, norb, alpha=1.0_dp, beta=0.0_dp)
     567         662 :       CALL parallel_gemm('T', 'N', norb, norb, nao, 1.0_dp, cxmat, evect, 0.0_dp, xxmat)
     568         662 :       CALL parallel_gemm('N', 'N', nao, norb, norb, 1.0_dp, mos_occ, xxmat, 0.0_dp, cxmat)
     569         662 :       CALL cp_fm_release(xxmat)
     570             :       ! C*C^T*XX
     571             :       CALL cp_dbcsr_plus_fm_fm_t(matrix_pe, matrix_v=mos_occ, matrix_g=cxmat, &
     572         662 :                                  ncol=norb, alpha=-1.0_dp, symmetry_mode=1)
     573         662 :       CALL cp_fm_release(cxmat)
     574             :       !
     575             :       ! Test for Tr(Pe*S)=0
     576         662 :       CALL dbcsr_dot(matrix_pe, matrix_s, tmp)
     577         662 :       IF (ABS(tmp) > 1.e-08_dp) THEN
     578           0 :          logger => cp_get_default_logger()
     579           0 :          IF (logger%para_env%is_source()) THEN
     580           0 :             iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     581             :          ELSE
     582             :             iounit = -1
     583             :          END IF
     584           0 :          CPWARN("Electron count of excitation density matrix is non-zero.")
     585           0 :          IF (iounit > 0) THEN
     586           0 :             WRITE (iounit, "(T2,A,T61,G20.10)") "Measured electron count is ", tmp
     587           0 :             WRITE (iounit, "(T2,A,/)") REPEAT("*", 79)
     588             :          END IF
     589             :       END IF
     590             :       !
     591             : 
     592         662 :       CALL timestop(handle)
     593             : 
     594         662 :    END SUBROUTINE tddfpt_resvec1
     595             : 
     596             : ! **************************************************************************************************
     597             : !> \brief PA = A * P * A(T)
     598             : !> \param matrix_pe ...
     599             : !> \param admm_env ...
     600             : !> \param matrix_pe_admm ...
     601             : ! **************************************************************************************************
     602         148 :    SUBROUTINE tddfpt_resvec1_admm(matrix_pe, admm_env, matrix_pe_admm)
     603             : 
     604             :       TYPE(dbcsr_type), POINTER                          :: matrix_pe
     605             :       TYPE(admm_type), POINTER                           :: admm_env
     606             :       TYPE(dbcsr_type), POINTER                          :: matrix_pe_admm
     607             : 
     608             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_resvec1_admm'
     609             : 
     610             :       INTEGER                                            :: handle, nao, nao_aux
     611             : 
     612         148 :       CALL timeset(routineN, handle)
     613             :       !
     614         148 :       nao_aux = admm_env%nao_aux_fit
     615         148 :       nao = admm_env%nao_orb
     616             :       !
     617         148 :       CALL copy_dbcsr_to_fm(matrix_pe, admm_env%work_orb_orb)
     618             :       CALL parallel_gemm('N', 'N', nao_aux, nao, nao, &
     619             :                          1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
     620         148 :                          admm_env%work_aux_orb)
     621             :       CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, &
     622             :                          1.0_dp, admm_env%work_aux_orb, admm_env%A, 0.0_dp, &
     623         148 :                          admm_env%work_aux_aux)
     624         148 :       CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, matrix_pe_admm, keep_sparsity=.TRUE.)
     625             :       !
     626         148 :       CALL timestop(handle)
     627             : 
     628         148 :    END SUBROUTINE tddfpt_resvec1_admm
     629             : 
     630             : ! **************************************************************************************************
     631             : !> \brief ...
     632             : !> \param qs_env ...
     633             : !> \param matrix_pe ...
     634             : !> \param matrix_pe_admm ...
     635             : !> \param gs_mos ...
     636             : !> \param matrix_hz ...
     637             : !> \param cpmos ...
     638             : ! **************************************************************************************************
     639         550 :    SUBROUTINE tddfpt_resvec2(qs_env, matrix_pe, matrix_pe_admm, gs_mos, matrix_hz, cpmos)
     640             : 
     641             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     642             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_pe, matrix_pe_admm
     643             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     644             :          POINTER                                         :: gs_mos
     645             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_hz
     646             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: cpmos
     647             : 
     648             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'tddfpt_resvec2'
     649             : 
     650             :       CHARACTER(LEN=default_string_length)               :: basis_type
     651             :       INTEGER                                            :: handle, iounit, ispin, mspin, n_rep_hf, &
     652             :                                                             nao, nao_aux, natom, norb, nspins
     653             :       LOGICAL                                            :: deriv2_analytic, distribute_fock_matrix, &
     654             :                                                             do_hfx, gapw, gapw_xc, &
     655             :                                                             hfx_treat_lsd_in_core, &
     656             :                                                             s_mstruct_changed
     657             :       REAL(KIND=dp)                                      :: eh1, focc, rhotot, thartree
     658             :       REAL(KIND=dp), DIMENSION(2)                        :: total_rho
     659         550 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: Qlm_tot
     660             :       TYPE(admm_type), POINTER                           :: admm_env
     661         550 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     662             :       TYPE(cp_fm_type), POINTER                          :: mos
     663             :       TYPE(cp_logger_type), POINTER                      :: logger
     664         550 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: msaux
     665         550 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mhz, mpe
     666             :       TYPE(dbcsr_type), POINTER                          :: dbwork
     667             :       TYPE(dft_control_type), POINTER                    :: dft_control
     668             :       TYPE(hartree_local_type), POINTER                  :: hartree_local
     669         550 :       TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
     670             :       TYPE(local_rho_type), POINTER                      :: local_rho_set, local_rho_set_admm
     671             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     672             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     673         550 :          POINTER                                         :: sab, sab_aux_fit
     674             :       TYPE(oce_matrix_type), POINTER                     :: oce
     675             :       TYPE(pw_c1d_gs_type)                               :: rho_tot_gspace, v_hartree_gspace
     676         550 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g, rho_g_aux, rhoz_g_aux, trho_g, &
     677         550 :                                                             trho_xc_g
     678             :       TYPE(pw_env_type), POINTER                         :: pw_env
     679             :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
     680             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     681             :       TYPE(pw_r3d_rs_type)                               :: v_hartree_rspace
     682         550 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r, rho_r_aux, rhoz_r_aux, tau_r, &
     683         550 :                                                             trho_r, trho_xc_r, v_xc, v_xc_tau
     684         550 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     685             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     686             :       TYPE(qs_rho_type), POINTER                         :: rho, rho_aux_fit, rho_xc, rhoz_aux, trho
     687         550 :       TYPE(rho_atom_type), DIMENSION(:), POINTER         :: rho1_atom_set, rho_atom_set
     688             :       TYPE(section_vals_type), POINTER                   :: hfx_section, input, xc_section
     689             :       TYPE(task_list_type), POINTER                      :: task_list
     690             : 
     691         550 :       CALL timeset(routineN, handle)
     692             : 
     693         550 :       NULLIFY (pw_env)
     694             :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, ks_env=ks_env, &
     695         550 :                       dft_control=dft_control, para_env=para_env)
     696         550 :       CPASSERT(ASSOCIATED(pw_env))
     697         550 :       nspins = dft_control%nspins
     698         550 :       gapw = dft_control%qs_control%gapw
     699         550 :       gapw_xc = dft_control%qs_control%gapw_xc
     700             : 
     701         550 :       CPASSERT(.NOT. dft_control%tddfpt2_control%do_exck)
     702         550 :       CPASSERT(.NOT. dft_control%tddfpt2_control%do_hfxsr)
     703         550 :       CPASSERT(.NOT. dft_control%tddfpt2_control%do_hfxlr)
     704             : 
     705         550 :       NULLIFY (auxbas_pw_pool, poisson_env)
     706             :       ! gets the tmp grids
     707             :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
     708         550 :                       poisson_env=poisson_env)
     709             : 
     710         550 :       CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
     711         550 :       CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
     712         550 :       CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
     713             : 
     714        4042 :       ALLOCATE (trho_r(nspins), trho_g(nspins))
     715        1196 :       DO ispin = 1, nspins
     716         646 :          CALL auxbas_pw_pool%create_pw(trho_r(ispin))
     717        1196 :          CALL auxbas_pw_pool%create_pw(trho_g(ispin))
     718             :       END DO
     719         550 :       IF (gapw_xc) THEN
     720          98 :          ALLOCATE (trho_xc_r(nspins), trho_xc_g(nspins))
     721          28 :          DO ispin = 1, nspins
     722          14 :             CALL auxbas_pw_pool%create_pw(trho_xc_r(ispin))
     723          28 :             CALL auxbas_pw_pool%create_pw(trho_xc_g(ispin))
     724             :          END DO
     725             :       END IF
     726             : 
     727             :       ! GAPW/GAPW_XC initializations
     728         550 :       NULLIFY (hartree_local, local_rho_set)
     729         550 :       IF (gapw) THEN
     730             :          CALL get_qs_env(qs_env, &
     731             :                          atomic_kind_set=atomic_kind_set, &
     732             :                          natom=natom, &
     733          62 :                          qs_kind_set=qs_kind_set)
     734          62 :          CALL local_rho_set_create(local_rho_set)
     735             :          CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, &
     736          62 :                                           qs_kind_set, dft_control, para_env)
     737             :          CALL init_rho0(local_rho_set, qs_env, dft_control%qs_control%gapw_control, &
     738          62 :                         zcore=0.0_dp)
     739          62 :          CALL rho0_s_grid_create(pw_env, local_rho_set%rho0_mpole)
     740          62 :          CALL hartree_local_create(hartree_local)
     741          62 :          CALL init_coulomb_local(hartree_local, natom)
     742         488 :       ELSEIF (gapw_xc) THEN
     743             :          CALL get_qs_env(qs_env, &
     744             :                          atomic_kind_set=atomic_kind_set, &
     745          14 :                          qs_kind_set=qs_kind_set)
     746          14 :          CALL local_rho_set_create(local_rho_set)
     747             :          CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, &
     748          14 :                                           qs_kind_set, dft_control, para_env)
     749             :       END IF
     750             : 
     751         550 :       total_rho = 0.0_dp
     752         550 :       CALL pw_zero(rho_tot_gspace)
     753        1196 :       DO ispin = 1, nspins
     754             :          CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_pe(ispin)%matrix, &
     755             :                                  rho=trho_r(ispin), &
     756             :                                  rho_gspace=trho_g(ispin), &
     757             :                                  soft_valid=gapw, &
     758         646 :                                  total_rho=total_rho(ispin))
     759         646 :          CALL pw_axpy(trho_g(ispin), rho_tot_gspace)
     760        1196 :          IF (gapw_xc) THEN
     761             :             CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_pe(ispin)%matrix, &
     762             :                                     rho=trho_xc_r(ispin), &
     763             :                                     rho_gspace=trho_xc_g(ispin), &
     764             :                                     soft_valid=gapw_xc, &
     765          14 :                                     total_rho=rhotot)
     766             :          END IF
     767             :       END DO
     768             : 
     769             :       ! GAPW o GAPW_XC require the calculation of hard and soft local densities
     770         550 :       IF (gapw .OR. gapw_xc) THEN
     771          76 :          CALL get_qs_env(qs_env=qs_env, oce=oce, sab_orb=sab)
     772             :          CALL calculate_rho_atom_coeff(qs_env, matrix_pe, local_rho_set%rho_atom_set, &
     773          76 :                                        qs_kind_set, oce, sab, para_env)
     774          76 :          CALL prepare_gapw_den(qs_env, local_rho_set, do_rho0=gapw)
     775             :       END IF
     776        1650 :       rhotot = SUM(total_rho)
     777         550 :       IF (gapw) THEN
     778          62 :          CALL get_rho0_mpole(local_rho_set%rho0_mpole, Qlm_tot=Qlm_tot)
     779          62 :          rhotot = rhotot + local_rho_set%rho0_mpole%total_rho0_h
     780          62 :          CALL pw_axpy(local_rho_set%rho0_mpole%rho0_s_gs, rho_tot_gspace)
     781             :       END IF
     782             : 
     783         550 :       IF (ABS(rhotot) > 1.e-05_dp) THEN
     784          20 :          logger => cp_get_default_logger()
     785          20 :          IF (logger%para_env%is_source()) THEN
     786          10 :             iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     787             :          ELSE
     788             :             iounit = -1
     789             :          END IF
     790          20 :          CPWARN("Real space electron count of excitation density is non-zero.")
     791          20 :          IF (iounit > 0) THEN
     792          10 :             WRITE (iounit, "(T2,A,T61,G20.10)") "Measured electron count is ", rhotot
     793          10 :             WRITE (iounit, "(T2,A,/)") REPEAT("*", 79)
     794             :          END IF
     795             :       END IF
     796             : 
     797             :       ! calculate associated hartree potential
     798             :       CALL pw_poisson_solve(poisson_env, rho_tot_gspace, thartree, &
     799         550 :                             v_hartree_gspace)
     800         550 :       CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
     801         550 :       CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
     802         550 :       IF (gapw) THEN
     803             :          CALL Vh_1c_gg_integrals(qs_env, thartree, hartree_local%ecoul_1c, &
     804          62 :                                  local_rho_set, para_env, tddft=.TRUE.)
     805             :          CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, &
     806             :                                     calculate_forces=.FALSE., &
     807          62 :                                     local_rho_set=local_rho_set)
     808             :       END IF
     809             : 
     810             :       ! Fxc*drho term
     811         550 :       CALL get_qs_env(qs_env, rho=rho)
     812         550 :       CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g)
     813             :       !
     814         550 :       CALL get_qs_env(qs_env, input=input)
     815         550 :       IF (dft_control%do_admm) THEN
     816         128 :          CALL get_qs_env(qs_env, admm_env=admm_env)
     817         128 :          xc_section => admm_env%xc_section_primary
     818             :       ELSE
     819         422 :          xc_section => section_vals_get_subs_vals(input, "DFT%XC")
     820             :       END IF
     821             :       !
     822         550 :       deriv2_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
     823         550 :       IF (deriv2_analytic) THEN
     824         550 :          NULLIFY (v_xc, v_xc_tau, tau_r)
     825         550 :          IF (gapw_xc) THEN
     826          14 :             CALL get_qs_env(qs_env=qs_env, rho_xc=rho_xc)
     827          14 :             CALL qs_fxc_analytic(rho_xc, trho_xc_r, tau_r, xc_section, auxbas_pw_pool, .FALSE., v_xc, v_xc_tau)
     828             :          ELSE
     829         536 :             CALL qs_fxc_analytic(rho, trho_r, tau_r, xc_section, auxbas_pw_pool, .FALSE., v_xc, v_xc_tau)
     830             :          END IF
     831         550 :          IF (gapw .OR. gapw_xc) THEN
     832          76 :             CALL get_qs_env(qs_env, rho_atom_set=rho_atom_set)
     833          76 :             rho1_atom_set => local_rho_set%rho_atom_set
     834             :             CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
     835          76 :                                              do_tddft=.TRUE., do_triplet=.FALSE.)
     836             :          END IF
     837             :       ELSE
     838           0 :          CPABORT("NYA 00006")
     839           0 :          NULLIFY (v_xc, trho)
     840           0 :          ALLOCATE (trho)
     841           0 :          CALL qs_rho_create(trho)
     842           0 :          CALL qs_rho_set(trho, rho_r=trho_r, rho_g=trho_g)
     843           0 :          CALL qs_fxc_fdiff(ks_env, rho, trho, xc_section, 6, .FALSE., v_xc, v_xc_tau)
     844           0 :          DEALLOCATE (trho)
     845             :       END IF
     846             : 
     847        1196 :       DO ispin = 1, nspins
     848         646 :          CALL dbcsr_set(matrix_hz(ispin)%matrix, 0.0_dp)
     849        1196 :          CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
     850             :       END DO
     851         550 :       IF (gapw_xc) THEN
     852          28 :          DO ispin = 1, nspins
     853             :             CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_hartree_rspace, &
     854             :                                     hmat=matrix_hz(ispin), &
     855          14 :                                     calculate_forces=.FALSE.)
     856             :             CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), &
     857             :                                     hmat=matrix_hz(ispin), &
     858          28 :                                     gapw=gapw_xc, calculate_forces=.FALSE.)
     859             :          END DO
     860             :       ELSE
     861             :          ! vtot = v_xc(ispin) + v_hartree
     862        1168 :          DO ispin = 1, nspins
     863             :             CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), &
     864             :                                     hmat=matrix_hz(ispin), &
     865         632 :                                     gapw=gapw, calculate_forces=.FALSE.)
     866             :             CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_hartree_rspace, &
     867             :                                     hmat=matrix_hz(ispin), &
     868        1168 :                                     gapw=gapw, calculate_forces=.FALSE.)
     869             :          END DO
     870             :       END IF
     871         550 :       IF (gapw .OR. gapw_xc) THEN
     872          76 :          mhz(1:nspins, 1:1) => matrix_hz(1:nspins)
     873          76 :          mpe(1:nspins, 1:1) => matrix_pe(1:nspins)
     874             :          CALL update_ks_atom(qs_env, mhz, mpe, forces=.FALSE., &
     875          76 :                              rho_atom_external=local_rho_set%rho_atom_set)
     876             :       END IF
     877             : 
     878         550 :       CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
     879         550 :       CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
     880         550 :       CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
     881        1196 :       DO ispin = 1, nspins
     882         646 :          CALL auxbas_pw_pool%give_back_pw(trho_r(ispin))
     883         646 :          CALL auxbas_pw_pool%give_back_pw(trho_g(ispin))
     884        1196 :          CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
     885             :       END DO
     886         550 :       DEALLOCATE (trho_r, trho_g, v_xc)
     887         550 :       IF (gapw_xc) THEN
     888          28 :          DO ispin = 1, nspins
     889          14 :             CALL auxbas_pw_pool%give_back_pw(trho_xc_r(ispin))
     890          28 :             CALL auxbas_pw_pool%give_back_pw(trho_xc_g(ispin))
     891             :          END DO
     892          14 :          DEALLOCATE (trho_xc_r, trho_xc_g)
     893             :       END IF
     894         550 :       IF (ASSOCIATED(v_xc_tau)) THEN
     895           0 :          DO ispin = 1, nspins
     896           0 :             CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
     897             :          END DO
     898           0 :          DEALLOCATE (v_xc_tau)
     899             :       END IF
     900         550 :       IF (dft_control%do_admm) THEN
     901         128 :          IF (qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
     902             :             ! nothing to do
     903             :          ELSE
     904             :             ! add ADMM xc_section_aux terms: f_x[rhoz_ADMM]
     905          78 :             CALL get_qs_env(qs_env, admm_env=admm_env)
     906             :             CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit, matrix_s_aux_fit=msaux, &
     907          78 :                               task_list_aux_fit=task_list)
     908          78 :             basis_type = "AUX_FIT"
     909             :             !
     910          78 :             NULLIFY (mpe, mhz)
     911         398 :             ALLOCATE (mpe(nspins, 1))
     912          78 :             CALL dbcsr_allocate_matrix_set(mhz, nspins, 1)
     913         164 :             DO ispin = 1, nspins
     914          86 :                ALLOCATE (mhz(ispin, 1)%matrix)
     915          86 :                CALL dbcsr_create(mhz(ispin, 1)%matrix, template=msaux(1)%matrix)
     916          86 :                CALL dbcsr_copy(mhz(ispin, 1)%matrix, msaux(1)%matrix)
     917          86 :                CALL dbcsr_set(mhz(ispin, 1)%matrix, 0.0_dp)
     918         164 :                mpe(ispin, 1)%matrix => matrix_pe_admm(ispin)%matrix
     919             :             END DO
     920             :             !
     921             :             ! GAPW/GAPW_XC initializations
     922          78 :             NULLIFY (local_rho_set_admm)
     923          78 :             IF (admm_env%do_gapw) THEN
     924           4 :                basis_type = "AUX_FIT_SOFT"
     925           4 :                task_list => admm_env%admm_gapw_env%task_list
     926           4 :                CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
     927           4 :                CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
     928           4 :                CALL local_rho_set_create(local_rho_set_admm)
     929             :                CALL allocate_rho_atom_internals(local_rho_set_admm%rho_atom_set, atomic_kind_set, &
     930           4 :                                                 admm_env%admm_gapw_env%admm_kind_set, dft_control, para_env)
     931             :                CALL calculate_rho_atom_coeff(qs_env, matrix_pe_admm, &
     932             :                                              rho_atom_set=local_rho_set_admm%rho_atom_set, &
     933             :                                              qs_kind_set=admm_env%admm_gapw_env%admm_kind_set, &
     934           4 :                                              oce=admm_env%admm_gapw_env%oce, sab=sab_aux_fit, para_env=para_env)
     935             :                CALL prepare_gapw_den(qs_env, local_rho_set=local_rho_set_admm, &
     936           4 :                                      do_rho0=.FALSE., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
     937             :             END IF
     938             :             !
     939          78 :             xc_section => admm_env%xc_section_aux
     940             :             !
     941          78 :             NULLIFY (rho_g_aux, rho_r_aux, rhoz_g_aux, rhoz_r_aux)
     942          78 :             CALL qs_rho_get(rho_aux_fit, rho_r=rho_r_aux, rho_g=rho_g_aux)
     943             :             ! rhoz_aux
     944         562 :             ALLOCATE (rhoz_r_aux(nspins), rhoz_g_aux(nspins))
     945         164 :             DO ispin = 1, nspins
     946          86 :                CALL auxbas_pw_pool%create_pw(rhoz_r_aux(ispin))
     947         164 :                CALL auxbas_pw_pool%create_pw(rhoz_g_aux(ispin))
     948             :             END DO
     949         164 :             DO ispin = 1, nspins
     950             :                CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpe(ispin, 1)%matrix, &
     951             :                                        rho=rhoz_r_aux(ispin), rho_gspace=rhoz_g_aux(ispin), &
     952             :                                        basis_type=basis_type, &
     953         164 :                                        task_list_external=task_list)
     954             :             END DO
     955             :             !
     956          78 :             NULLIFY (v_xc)
     957          78 :             deriv2_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
     958          78 :             IF (deriv2_analytic) THEN
     959          78 :                NULLIFY (tau_r)
     960          78 :                CALL qs_fxc_analytic(rho_aux_fit, rhoz_r_aux, tau_r, xc_section, auxbas_pw_pool, .FALSE., v_xc, v_xc_tau)
     961             :             ELSE
     962           0 :                CPABORT("NYA 00007")
     963             :                NULLIFY (rhoz_aux)
     964           0 :                ALLOCATE (rhoz_aux)
     965           0 :                CALL qs_rho_create(rhoz_aux)
     966           0 :                CALL qs_rho_set(rhoz_aux, rho_r=rhoz_r_aux, rho_g=rhoz_g_aux)
     967           0 :                CALL qs_fxc_fdiff(ks_env, rho_aux_fit, rhoz_aux, xc_section, 6, .FALSE., v_xc, v_xc_tau)
     968           0 :                DEALLOCATE (rhoz_aux)
     969             :             END IF
     970             :             !
     971         164 :             DO ispin = 1, nspins
     972          86 :                CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
     973             :                CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), &
     974             :                                        hmat=mhz(ispin, 1), basis_type=basis_type, &
     975             :                                        calculate_forces=.FALSE., &
     976         164 :                                        task_list_external=task_list)
     977             :             END DO
     978         164 :             DO ispin = 1, nspins
     979          86 :                CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
     980          86 :                CALL auxbas_pw_pool%give_back_pw(rhoz_r_aux(ispin))
     981         164 :                CALL auxbas_pw_pool%give_back_pw(rhoz_g_aux(ispin))
     982             :             END DO
     983          78 :             DEALLOCATE (v_xc, rhoz_r_aux, rhoz_g_aux)
     984             :             !
     985          78 :             IF (admm_env%do_gapw) THEN
     986           4 :                rho_atom_set => admm_env%admm_gapw_env%local_rho_set%rho_atom_set
     987           4 :                rho1_atom_set => local_rho_set_admm%rho_atom_set
     988             :                CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, &
     989           4 :                                                 para_env, kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
     990             :                CALL update_ks_atom(qs_env, mhz(:, 1), matrix_pe_admm, forces=.FALSE., tddft=.FALSE., &
     991             :                                    rho_atom_external=rho1_atom_set, &
     992             :                                    kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
     993             :                                    oce_external=admm_env%admm_gapw_env%oce, &
     994           4 :                                    sab_external=sab_aux_fit)
     995             :             END IF
     996             :             !
     997          78 :             nao = admm_env%nao_orb
     998          78 :             nao_aux = admm_env%nao_aux_fit
     999          78 :             ALLOCATE (dbwork)
    1000          78 :             CALL dbcsr_create(dbwork, template=matrix_hz(1)%matrix)
    1001         164 :             DO ispin = 1, nspins
    1002             :                CALL cp_dbcsr_sm_fm_multiply(mhz(ispin, 1)%matrix, admm_env%A, &
    1003          86 :                                             admm_env%work_aux_orb, nao)
    1004             :                CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
    1005             :                                   1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
    1006          86 :                                   admm_env%work_orb_orb)
    1007          86 :                CALL dbcsr_copy(dbwork, matrix_hz(1)%matrix)
    1008          86 :                CALL dbcsr_set(dbwork, 0.0_dp)
    1009          86 :                CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.TRUE.)
    1010         164 :                CALL dbcsr_add(matrix_hz(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
    1011             :             END DO
    1012          78 :             CALL dbcsr_release(dbwork)
    1013          78 :             DEALLOCATE (dbwork)
    1014          78 :             CALL dbcsr_deallocate_matrix_set(mhz)
    1015          78 :             DEALLOCATE (mpe)
    1016          78 :             IF (admm_env%do_gapw) THEN
    1017           4 :                IF (ASSOCIATED(local_rho_set_admm)) CALL local_rho_set_release(local_rho_set_admm)
    1018             :             END IF
    1019             :          END IF
    1020             :       END IF
    1021         550 :       IF (gapw .OR. gapw_xc) THEN
    1022          76 :          IF (ASSOCIATED(local_rho_set)) CALL local_rho_set_release(local_rho_set)
    1023          76 :          IF (ASSOCIATED(hartree_local)) CALL hartree_local_release(hartree_local)
    1024             :       END IF
    1025             : 
    1026             :       ! HFX
    1027         550 :       hfx_section => section_vals_get_subs_vals(xc_section, "HF")
    1028         550 :       CALL section_vals_get(hfx_section, explicit=do_hfx)
    1029         550 :       IF (do_hfx) THEN
    1030         248 :          CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
    1031         248 :          CPASSERT(n_rep_hf == 1)
    1032             :          CALL section_vals_val_get(hfx_section, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
    1033         248 :                                    i_rep_section=1)
    1034         248 :          mspin = 1
    1035         248 :          IF (hfx_treat_lsd_in_core) mspin = nspins
    1036             :          !
    1037             :          CALL get_qs_env(qs_env=qs_env, rho=rho, x_data=x_data, para_env=para_env, &
    1038         248 :                          s_mstruct_changed=s_mstruct_changed)
    1039         248 :          distribute_fock_matrix = .TRUE.
    1040         248 :          IF (dft_control%do_admm) THEN
    1041         128 :             CALL get_qs_env(qs_env, admm_env=admm_env)
    1042         128 :             CALL get_admm_env(admm_env, matrix_s_aux_fit=msaux)
    1043         128 :             NULLIFY (mpe, mhz)
    1044         660 :             ALLOCATE (mpe(nspins, 1))
    1045         128 :             CALL dbcsr_allocate_matrix_set(mhz, nspins, 1)
    1046         276 :             DO ispin = 1, nspins
    1047         148 :                ALLOCATE (mhz(ispin, 1)%matrix)
    1048         148 :                CALL dbcsr_create(mhz(ispin, 1)%matrix, template=msaux(1)%matrix)
    1049         148 :                CALL dbcsr_copy(mhz(ispin, 1)%matrix, msaux(1)%matrix)
    1050         148 :                CALL dbcsr_set(mhz(ispin, 1)%matrix, 0.0_dp)
    1051         276 :                mpe(ispin, 1)%matrix => matrix_pe_admm(ispin)%matrix
    1052             :             END DO
    1053         128 :             IF (x_data(1, 1)%do_hfx_ri) THEN
    1054             :                eh1 = 0.0_dp
    1055             :                CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhz, eh1, rho_ao=mpe, &
    1056             :                                      geometry_did_change=s_mstruct_changed, nspins=nspins, &
    1057           6 :                                      hf_fraction=x_data(1, 1)%general_parameter%fraction)
    1058             :             ELSE
    1059         244 :                DO ispin = 1, mspin
    1060             :                   eh1 = 0.0
    1061             :                   CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpe, hfx_section, &
    1062             :                                              para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
    1063         244 :                                              ispin=ispin)
    1064             :                END DO
    1065             :             END IF
    1066             :             !
    1067         128 :             CPASSERT(ASSOCIATED(admm_env%work_aux_orb))
    1068         128 :             CPASSERT(ASSOCIATED(admm_env%work_orb_orb))
    1069         128 :             nao = admm_env%nao_orb
    1070         128 :             nao_aux = admm_env%nao_aux_fit
    1071         128 :             ALLOCATE (dbwork)
    1072         128 :             CALL dbcsr_create(dbwork, template=matrix_hz(1)%matrix)
    1073         276 :             DO ispin = 1, nspins
    1074             :                CALL cp_dbcsr_sm_fm_multiply(mhz(ispin, 1)%matrix, admm_env%A, &
    1075         148 :                                             admm_env%work_aux_orb, nao)
    1076             :                CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
    1077             :                                   1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
    1078         148 :                                   admm_env%work_orb_orb)
    1079         148 :                CALL dbcsr_copy(dbwork, matrix_hz(ispin)%matrix)
    1080         148 :                CALL dbcsr_set(dbwork, 0.0_dp)
    1081         148 :                CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.TRUE.)
    1082         276 :                CALL dbcsr_add(matrix_hz(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
    1083             :             END DO
    1084         128 :             CALL dbcsr_release(dbwork)
    1085         128 :             DEALLOCATE (dbwork)
    1086         128 :             CALL dbcsr_deallocate_matrix_set(mhz)
    1087         128 :             DEALLOCATE (mpe)
    1088             :          ELSE
    1089         120 :             NULLIFY (mpe, mhz)
    1090        1120 :             ALLOCATE (mpe(nspins, 1), mhz(nspins, 1))
    1091         260 :             DO ispin = 1, nspins
    1092         140 :                mhz(ispin, 1)%matrix => matrix_hz(ispin)%matrix
    1093         260 :                mpe(ispin, 1)%matrix => matrix_pe(ispin)%matrix
    1094             :             END DO
    1095         120 :             IF (x_data(1, 1)%do_hfx_ri) THEN
    1096             :                eh1 = 0.0_dp
    1097             :                CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhz, eh1, rho_ao=mpe, &
    1098             :                                      geometry_did_change=s_mstruct_changed, nspins=nspins, &
    1099          18 :                                      hf_fraction=x_data(1, 1)%general_parameter%fraction)
    1100             :             ELSE
    1101         204 :                DO ispin = 1, mspin
    1102             :                   eh1 = 0.0
    1103             :                   CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpe, hfx_section, &
    1104             :                                              para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
    1105         204 :                                              ispin=ispin)
    1106             :                END DO
    1107             :             END IF
    1108         120 :             DEALLOCATE (mpe, mhz)
    1109             :          END IF
    1110             :       END IF
    1111             : 
    1112         550 :       focc = 4.0_dp
    1113         550 :       IF (nspins == 2) focc = 2.0_dp
    1114        1196 :       DO ispin = 1, nspins
    1115         646 :          mos => gs_mos(ispin)%mos_occ
    1116         646 :          CALL cp_fm_get_info(mos, ncol_global=norb)
    1117             :          CALL cp_dbcsr_sm_fm_multiply(matrix_hz(ispin)%matrix, mos, cpmos(ispin), &
    1118        1196 :                                       norb, alpha=focc, beta=0.0_dp)
    1119             :       END DO
    1120             : 
    1121         550 :       CALL timestop(handle)
    1122             : 
    1123        1650 :    END SUBROUTINE tddfpt_resvec2
    1124             : 
    1125             : ! **************************************************************************************************
    1126             : !> \brief ...
    1127             : !> \param qs_env ...
    1128             : !> \param matrix_pe ...
    1129             : !> \param gs_mos ...
    1130             : !> \param matrix_hz ...
    1131             : !> \param cpmos ...
    1132             : ! **************************************************************************************************
    1133          16 :    SUBROUTINE tddfpt_resvec2_xtb(qs_env, matrix_pe, gs_mos, matrix_hz, cpmos)
    1134             : 
    1135             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1136             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_pe
    1137             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
    1138             :          POINTER                                         :: gs_mos
    1139             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_hz
    1140             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: cpmos
    1141             : 
    1142             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_resvec2_xtb'
    1143             : 
    1144             :       INTEGER                                            :: atom_a, handle, iatom, ikind, is, ispin, &
    1145             :                                                             na, natom, natorb, nkind, norb, ns, &
    1146             :                                                             nsgf, nspins
    1147             :       INTEGER, DIMENSION(25)                             :: lao
    1148             :       INTEGER, DIMENSION(5)                              :: occ
    1149          16 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: mcharge, mcharge1
    1150          16 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: aocg, aocg1, charges, charges1
    1151             :       REAL(KIND=dp)                                      :: focc
    1152          16 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1153             :       TYPE(cp_fm_type), POINTER                          :: mos
    1154          16 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: p_matrix
    1155          16 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_p, matrix_s
    1156             :       TYPE(dbcsr_type), POINTER                          :: s_matrix
    1157             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1158             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1159          16 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1160          16 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1161             :       TYPE(qs_rho_type), POINTER                         :: rho
    1162             :       TYPE(xtb_atom_type), POINTER                       :: xtb_kind
    1163             : 
    1164          16 :       CALL timeset(routineN, handle)
    1165             : 
    1166          16 :       CPASSERT(ASSOCIATED(matrix_pe))
    1167             : 
    1168          16 :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
    1169          16 :       nspins = dft_control%nspins
    1170             : 
    1171          32 :       DO ispin = 1, nspins
    1172          32 :          CALL dbcsr_set(matrix_hz(ispin)%matrix, 0.0_dp)
    1173             :       END DO
    1174             : 
    1175          16 :       IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
    1176             :          ! Mulliken charges
    1177             :          CALL get_qs_env(qs_env, rho=rho, particle_set=particle_set, &
    1178          14 :                          matrix_s_kp=matrix_s, para_env=para_env)
    1179          14 :          natom = SIZE(particle_set)
    1180          14 :          CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
    1181          70 :          ALLOCATE (mcharge(natom), charges(natom, 5))
    1182          70 :          ALLOCATE (mcharge1(natom), charges1(natom, 5))
    1183        1254 :          charges = 0.0_dp
    1184        1254 :          charges1 = 0.0_dp
    1185          14 :          CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
    1186          14 :          nkind = SIZE(atomic_kind_set)
    1187          14 :          CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
    1188          56 :          ALLOCATE (aocg(nsgf, natom))
    1189        1184 :          aocg = 0.0_dp
    1190          56 :          ALLOCATE (aocg1(nsgf, natom))
    1191        1184 :          aocg1 = 0.0_dp
    1192          14 :          p_matrix => matrix_p(:, 1)
    1193          14 :          s_matrix => matrix_s(1, 1)%matrix
    1194          14 :          CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
    1195          14 :          CALL ao_charges(matrix_pe, s_matrix, aocg1, para_env)
    1196          48 :          DO ikind = 1, nkind
    1197          34 :             CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
    1198          34 :             CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
    1199          34 :             CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
    1200         316 :             DO iatom = 1, na
    1201         234 :                atom_a = atomic_kind_set(ikind)%atom_list(iatom)
    1202        1404 :                charges(atom_a, :) = REAL(occ(:), KIND=dp)
    1203         900 :                DO is = 1, natorb
    1204         632 :                   ns = lao(is) + 1
    1205         632 :                   charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
    1206         866 :                   charges1(atom_a, ns) = charges1(atom_a, ns) - aocg1(is, atom_a)
    1207             :                END DO
    1208             :             END DO
    1209             :          END DO
    1210          14 :          DEALLOCATE (aocg, aocg1)
    1211         248 :          DO iatom = 1, natom
    1212        1404 :             mcharge(iatom) = SUM(charges(iatom, :))
    1213        1418 :             mcharge1(iatom) = SUM(charges1(iatom, :))
    1214             :          END DO
    1215             :          ! Coulomb Kernel
    1216          14 :          CALL xtb_coulomb_hessian(qs_env, matrix_hz, charges1, mcharge1, mcharge)
    1217             :          !
    1218          28 :          DEALLOCATE (charges, mcharge, charges1, mcharge1)
    1219             :       END IF
    1220             : 
    1221          16 :       focc = 2.0_dp
    1222          16 :       IF (nspins == 2) focc = 1.0_dp
    1223          32 :       DO ispin = 1, nspins
    1224          16 :          mos => gs_mos(ispin)%mos_occ
    1225          16 :          CALL cp_fm_get_info(mos, ncol_global=norb)
    1226             :          CALL cp_dbcsr_sm_fm_multiply(matrix_hz(ispin)%matrix, mos, cpmos(ispin), &
    1227          32 :                                       norb, alpha=focc, beta=0.0_dp)
    1228             :       END DO
    1229             : 
    1230          16 :       CALL timestop(handle)
    1231             : 
    1232          32 :    END SUBROUTINE tddfpt_resvec2_xtb
    1233             : 
    1234             : ! **************************************************************************************************
    1235             : !> \brief ...
    1236             : !> \param qs_env ...
    1237             : !> \param cpmos ...
    1238             : !> \param work ...
    1239             : ! **************************************************************************************************
    1240         566 :    SUBROUTINE tddfpt_resvec3(qs_env, cpmos, work)
    1241             : 
    1242             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1243             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: cpmos
    1244             :       TYPE(tddfpt_work_matrices)                         :: work
    1245             : 
    1246             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'tddfpt_resvec3'
    1247             : 
    1248             :       INTEGER                                            :: handle, ispin, nao, norb, nspins
    1249             :       TYPE(cp_fm_struct_type), POINTER                   :: fmstruct
    1250             :       TYPE(cp_fm_type)                                   :: cvec, umat
    1251             :       TYPE(cp_fm_type), POINTER                          :: omos
    1252             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1253         566 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    1254             : 
    1255         566 :       CALL timeset(routineN, handle)
    1256             : 
    1257         566 :       CALL get_qs_env(qs_env, mos=mos, dft_control=dft_control)
    1258         566 :       nspins = dft_control%nspins
    1259             : 
    1260        1228 :       DO ispin = 1, nspins
    1261         662 :          CALL get_mo_set(mos(ispin), mo_coeff=omos)
    1262             :          ASSOCIATE (rvecs => cpmos(ispin))
    1263         662 :             CALL cp_fm_get_info(rvecs, nrow_global=nao, ncol_global=norb)
    1264         662 :             CALL cp_fm_create(cvec, rvecs%matrix_struct, "cvec")
    1265             :             CALL cp_fm_struct_create(fmstruct, context=rvecs%matrix_struct%context, nrow_global=norb, &
    1266         662 :                                      ncol_global=norb, para_env=rvecs%matrix_struct%para_env)
    1267         662 :             CALL cp_fm_create(umat, fmstruct, "umat")
    1268         662 :             CALL cp_fm_struct_release(fmstruct)
    1269             :             !
    1270         662 :             CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, omos, work%S_C0(ispin), 0.0_dp, umat)
    1271         662 :             CALL cp_fm_copy_general(rvecs, cvec, rvecs%matrix_struct%para_env)
    1272         662 :             CALL parallel_gemm("N", "T", nao, norb, norb, 1.0_dp, cvec, umat, 0.0_dp, rvecs)
    1273             :          END ASSOCIATE
    1274         662 :          CALL cp_fm_release(cvec)
    1275        1890 :          CALL cp_fm_release(umat)
    1276             :       END DO
    1277             : 
    1278         566 :       CALL timestop(handle)
    1279             : 
    1280         566 :    END SUBROUTINE tddfpt_resvec3
    1281             : 
    1282             : ! **************************************************************************************************
    1283             : !> \brief Calculate direct tddft forces
    1284             : !> \param qs_env ...
    1285             : !> \param ex_env ...
    1286             : !> \param gs_mos ...
    1287             : !> \param kernel_env ...
    1288             : !> \param sub_env ...
    1289             : !> \param work_matrices ...
    1290             : !> \param debug_forces ...
    1291             : !> \par History
    1292             : !>    * 01.2020 screated [JGH]
    1293             : ! **************************************************************************************************
    1294         566 :    SUBROUTINE tddfpt_kernel_force(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices, debug_forces)
    1295             : 
    1296             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1297             :       TYPE(excited_energy_type), POINTER                 :: ex_env
    1298             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
    1299             :          POINTER                                         :: gs_mos
    1300             :       TYPE(kernel_env_type), INTENT(IN)                  :: kernel_env
    1301             :       TYPE(tddfpt_subgroup_env_type)                     :: sub_env
    1302             :       TYPE(tddfpt_work_matrices)                         :: work_matrices
    1303             :       LOGICAL, INTENT(IN)                                :: debug_forces
    1304             : 
    1305             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_kernel_force'
    1306             : 
    1307             :       INTEGER                                            :: handle
    1308             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1309             :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
    1310             : 
    1311             :       MARK_USED(work_matrices)
    1312             : 
    1313         566 :       CALL timeset(routineN, handle)
    1314             : 
    1315         566 :       CALL get_qs_env(qs_env, dft_control=dft_control)
    1316         566 :       tddfpt_control => dft_control%tddfpt2_control
    1317             : 
    1318         566 :       IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
    1319             :          ! full Kernel
    1320         340 :          CALL fhxc_force(qs_env, ex_env, gs_mos, kernel_env%full_kernel, debug_forces)
    1321         226 :       ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
    1322             :          ! sTDA Kernel
    1323         158 :          CALL stda_force(qs_env, ex_env, gs_mos, kernel_env%stda_kernel, sub_env, work_matrices, debug_forces)
    1324          68 :       ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
    1325             :          ! nothing to be done here
    1326          68 :          ex_env%matrix_wx1 => NULL()
    1327             :       ELSE
    1328           0 :          CPABORT('Unknown kernel type')
    1329             :       END IF
    1330             : 
    1331         566 :       CALL timestop(handle)
    1332             : 
    1333         566 :    END SUBROUTINE tddfpt_kernel_force
    1334             : 
    1335             : END MODULE qs_tddfpt2_forces

Generated by: LCOV version 1.15