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

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

Generated by: LCOV version 2.0-1