LCOV - code coverage report
Current view: top level - src - qs_tddfpt2_fhxc_forces.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 96.8 % 928 898
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 2 2

            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_fhxc_forces
       9              :    USE admm_methods,                    ONLY: admm_projection_derivative
      10              :    USE admm_types,                      ONLY: admm_type,&
      11              :                                               get_admm_env
      12              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      13              :                                               get_atomic_kind_set
      14              :    USE cell_types,                      ONLY: cell_type,&
      15              :                                               pbc
      16              :    USE cp_control_types,                ONLY: dft_control_type,&
      17              :                                               stda_control_type,&
      18              :                                               tddfpt2_control_type
      19              :    USE cp_dbcsr_api,                    ONLY: &
      20              :         dbcsr_add, dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, dbcsr_filter, &
      21              :         dbcsr_get_block_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
      22              :         dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
      23              :         dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_transposed, dbcsr_type, &
      24              :         dbcsr_type_antisymmetric, dbcsr_type_no_symmetry
      25              :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      26              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      27              :                                               copy_fm_to_dbcsr,&
      28              :                                               cp_dbcsr_plus_fm_fm_t,&
      29              :                                               cp_dbcsr_sm_fm_multiply,&
      30              :                                               dbcsr_allocate_matrix_set,&
      31              :                                               dbcsr_deallocate_matrix_set
      32              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_row_scale,&
      33              :                                               cp_fm_schur_product
      34              :    USE cp_fm_pool_types,                ONLY: fm_pool_create_fm,&
      35              :                                               fm_pool_give_back_fm
      36              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      37              :                                               cp_fm_struct_release,&
      38              :                                               cp_fm_struct_type
      39              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      40              :                                               cp_fm_get_info,&
      41              :                                               cp_fm_release,&
      42              :                                               cp_fm_to_fm,&
      43              :                                               cp_fm_type,&
      44              :                                               cp_fm_vectorssum
      45              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      46              :                                               cp_logger_get_default_unit_nr,&
      47              :                                               cp_logger_type
      48              :    USE ewald_environment_types,         ONLY: ewald_env_get,&
      49              :                                               ewald_environment_type
      50              :    USE ewald_methods_tb,                ONLY: tb_ewald_overlap,&
      51              :                                               tb_spme_evaluate
      52              :    USE ewald_pw_types,                  ONLY: ewald_pw_type
      53              :    USE exstates_types,                  ONLY: excited_energy_type
      54              :    USE hartree_local_methods,           ONLY: Vh_1c_gg_integrals,&
      55              :                                               init_coulomb_local
      56              :    USE hartree_local_types,             ONLY: hartree_local_create,&
      57              :                                               hartree_local_release,&
      58              :                                               hartree_local_type
      59              :    USE hfx_derivatives,                 ONLY: derivatives_four_center
      60              :    USE hfx_energy_potential,            ONLY: integrate_four_center
      61              :    USE hfx_ri,                          ONLY: hfx_ri_update_forces,&
      62              :                                               hfx_ri_update_ks
      63              :    USE hfx_types,                       ONLY: hfx_type
      64              :    USE input_constants,                 ONLY: do_admm_aux_exch_func_none,&
      65              :                                               tddfpt_kernel_full,&
      66              :                                               xc_kernel_method_analytic,&
      67              :                                               xc_kernel_method_best,&
      68              :                                               xc_kernel_method_numeric,&
      69              :                                               xc_none
      70              :    USE input_section_types,             ONLY: section_vals_get,&
      71              :                                               section_vals_get_subs_vals,&
      72              :                                               section_vals_type,&
      73              :                                               section_vals_val_get
      74              :    USE kinds,                           ONLY: default_string_length,&
      75              :                                               dp
      76              :    USE mathconstants,                   ONLY: oorootpi
      77              :    USE message_passing,                 ONLY: mp_para_env_type
      78              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      79              :    USE particle_methods,                ONLY: get_particle_set
      80              :    USE particle_types,                  ONLY: particle_type
      81              :    USE pw_env_types,                    ONLY: pw_env_get,&
      82              :                                               pw_env_type
      83              :    USE pw_methods,                      ONLY: pw_axpy,&
      84              :                                               pw_scale,&
      85              :                                               pw_transfer,&
      86              :                                               pw_zero
      87              :    USE pw_poisson_methods,              ONLY: pw_poisson_solve
      88              :    USE pw_poisson_types,                ONLY: pw_poisson_type
      89              :    USE pw_pool_types,                   ONLY: pw_pool_type
      90              :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      91              :                                               pw_r3d_rs_type
      92              :    USE qs_collocate_density,            ONLY: calculate_rho_elec
      93              :    USE qs_environment_types,            ONLY: get_qs_env,&
      94              :                                               qs_environment_type,&
      95              :                                               set_qs_env
      96              :    USE qs_force_types,                  ONLY: qs_force_type
      97              :    USE qs_fxc,                          ONLY: qs_fgxc_create,&
      98              :                                               qs_fgxc_gdiff,&
      99              :                                               qs_fgxc_release
     100              :    USE qs_gapw_densities,               ONLY: prepare_gapw_den
     101              :    USE qs_integrate_potential,          ONLY: integrate_v_rspace
     102              :    USE qs_kernel_types,                 ONLY: full_kernel_env_type
     103              :    USE qs_kind_types,                   ONLY: qs_kind_type
     104              :    USE qs_ks_atom,                      ONLY: update_ks_atom
     105              :    USE qs_ks_types,                     ONLY: qs_ks_env_type
     106              :    USE qs_local_rho_types,              ONLY: local_rho_set_create,&
     107              :                                               local_rho_set_release,&
     108              :                                               local_rho_type
     109              :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
     110              :    USE qs_oce_methods,                  ONLY: build_oce_matrices
     111              :    USE qs_oce_types,                    ONLY: allocate_oce_set,&
     112              :                                               create_oce_set,&
     113              :                                               oce_matrix_type
     114              :    USE qs_overlap,                      ONLY: build_overlap_matrix
     115              :    USE qs_rho0_ggrid,                   ONLY: integrate_vhg0_rspace,&
     116              :                                               rho0_s_grid_create
     117              :    USE qs_rho0_methods,                 ONLY: init_rho0
     118              :    USE qs_rho_atom_methods,             ONLY: allocate_rho_atom_internals,&
     119              :                                               calculate_rho_atom_coeff
     120              :    USE qs_rho_atom_types,               ONLY: rho_atom_type
     121              :    USE qs_rho_types,                    ONLY: qs_rho_create,&
     122              :                                               qs_rho_get,&
     123              :                                               qs_rho_set,&
     124              :                                               qs_rho_type
     125              :    USE qs_tddfpt2_stda_types,           ONLY: stda_env_type
     126              :    USE qs_tddfpt2_stda_utils,           ONLY: get_lowdin_x,&
     127              :                                               setup_gamma
     128              :    USE qs_tddfpt2_subgroups,            ONLY: tddfpt_subgroup_env_type
     129              :    USE qs_tddfpt2_types,                ONLY: tddfpt_ground_state_mos,&
     130              :                                               tddfpt_work_matrices
     131              :    USE qs_vxc_atom,                     ONLY: calculate_gfxc_atom,&
     132              :                                               gfxc_atom_diff
     133              :    USE task_list_types,                 ONLY: task_list_type
     134              :    USE util,                            ONLY: get_limit
     135              :    USE virial_types,                    ONLY: virial_type
     136              : #include "./base/base_uses.f90"
     137              : 
     138              :    IMPLICIT NONE
     139              : 
     140              :    PRIVATE
     141              : 
     142              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_fhxc_forces'
     143              : 
     144              :    PUBLIC :: fhxc_force, stda_force
     145              : 
     146              : ! **************************************************************************************************
     147              : 
     148              : CONTAINS
     149              : 
     150              : ! **************************************************************************************************
     151              : !> \brief Calculate direct tddft forces
     152              : !> \param qs_env ...
     153              : !> \param ex_env ...
     154              : !> \param gs_mos ...
     155              : !> \param full_kernel ...
     156              : !> \param debug_forces ...
     157              : !> \par History
     158              : !>    * 01.2020 screated [JGH]
     159              : ! **************************************************************************************************
     160          342 :    SUBROUTINE fhxc_force(qs_env, ex_env, gs_mos, full_kernel, debug_forces)
     161              : 
     162              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     163              :       TYPE(excited_energy_type), POINTER                 :: ex_env
     164              :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     165              :          POINTER                                         :: gs_mos
     166              :       TYPE(full_kernel_env_type), INTENT(IN)             :: full_kernel
     167              :       LOGICAL, INTENT(IN)                                :: debug_forces
     168              : 
     169              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'fhxc_force'
     170              : 
     171              :       CHARACTER(LEN=default_string_length)               :: basis_type
     172              :       INTEGER                                            :: handle, iounit, ispin, mspin, myfun, &
     173              :                                                             n_rep_hf, nao, nao_aux, natom, nkind, &
     174              :                                                             norb, nspins, order
     175              :       LOGICAL :: distribute_fock_matrix, do_admm, do_analytic, do_hfx, do_hfxlr, do_hfxsr, &
     176              :          do_numeric, gapw, gapw_xc, hfx_treat_lsd_in_core, is_rks_triplets, s_mstruct_changed, &
     177              :          use_virial
     178              :       REAL(KIND=dp)                                      :: eh1, eh1c, eps_delta, eps_fit, focc, &
     179              :                                                             fscal, fval, kval, xehartree
     180              :       REAL(KIND=dp), DIMENSION(3)                        :: fodeb
     181              :       TYPE(admm_type), POINTER                           :: admm_env
     182          342 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     183              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct, fm_struct_mat
     184              :       TYPE(cp_fm_type)                                   :: cvcmat, vcvec
     185          342 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: cpmos, evect
     186              :       TYPE(cp_fm_type), POINTER                          :: mos
     187              :       TYPE(cp_logger_type), POINTER                      :: logger
     188          342 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_fx, matrix_gx, matrix_hfx, &
     189          342 :          matrix_hfx_admm, matrix_hfx_admm_asymm, matrix_hfx_asymm, matrix_hx, matrix_p, &
     190          342 :          matrix_p_admm, matrix_px1, matrix_px1_admm, matrix_px1_admm_asymm, matrix_px1_asymm, &
     191          342 :          matrix_s, matrix_s_aux_fit, matrix_wx1, mdum, mfx, mgx
     192          342 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mhe, mpe, mpga
     193              :       TYPE(dbcsr_type), POINTER                          :: dbwork, dbwork_asymm
     194              :       TYPE(dft_control_type), POINTER                    :: dft_control
     195              :       TYPE(hartree_local_type), POINTER                  :: hartree_local
     196          342 :       TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
     197              :       TYPE(local_rho_type), POINTER :: local_rho_set, local_rho_set_admm, local_rho_set_f, &
     198              :          local_rho_set_f_admm, local_rho_set_g, local_rho_set_g_admm
     199              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     200              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     201          342 :          POINTER                                         :: sab, sab_aux_fit, sab_orb, sap_oce
     202              :       TYPE(oce_matrix_type), POINTER                     :: oce
     203          342 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     204              :       TYPE(pw_c1d_gs_type)                               :: rhox_tot_gspace, xv_hartree_gspace
     205          342 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g_aux, rhox_g, rhox_g_aux, rhoxx_g
     206              :       TYPE(pw_env_type), POINTER                         :: pw_env
     207              :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
     208              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     209              :       TYPE(pw_r3d_rs_type)                               :: xv_hartree_rspace
     210          342 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: fxc_rho, fxc_tau, gxc_rho, gxc_tau, &
     211          342 :                                                             rho_r_aux, rhox_r, rhox_r_aux, rhoxx_r
     212          342 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     213          342 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     214              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     215              :       TYPE(qs_rho_type), POINTER                         :: rho, rho_aux_fit, rhox, rhox_aux
     216          342 :       TYPE(rho_atom_type), DIMENSION(:), POINTER         :: rho_atom_set, rho_atom_set_f, &
     217          342 :                                                             rho_atom_set_g
     218              :       TYPE(section_vals_type), POINTER                   :: hfx_section, xc_section
     219              :       TYPE(task_list_type), POINTER                      :: task_list
     220              :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     221              : 
     222          342 :       CALL timeset(routineN, handle)
     223              : 
     224          342 :       logger => cp_get_default_logger()
     225          342 :       IF (logger%para_env%is_source()) THEN
     226          171 :          iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     227              :       ELSE
     228              :          iounit = -1
     229              :       END IF
     230              : 
     231          342 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     232          342 :       tddfpt_control => dft_control%tddfpt2_control
     233          342 :       nspins = dft_control%nspins
     234          342 :       is_rks_triplets = tddfpt_control%rks_triplets .AND. (nspins == 1)
     235          342 :       CPASSERT(tddfpt_control%kernel == tddfpt_kernel_full)
     236          342 :       do_hfx = tddfpt_control%do_hfx
     237          342 :       do_hfxsr = tddfpt_control%do_hfxsr
     238          342 :       do_hfxlr = tddfpt_control%do_hfxlr
     239          342 :       do_admm = tddfpt_control%do_admm
     240          342 :       gapw = dft_control%qs_control%gapw
     241          342 :       gapw_xc = dft_control%qs_control%gapw_xc
     242              : 
     243          342 :       evect => ex_env%evect
     244          342 :       matrix_px1 => ex_env%matrix_px1
     245          342 :       matrix_px1_admm => ex_env%matrix_px1_admm
     246          342 :       matrix_px1_asymm => ex_env%matrix_px1_asymm
     247          342 :       matrix_px1_admm_asymm => ex_env%matrix_px1_admm_asymm
     248              : 
     249          342 :       focc = 1.0_dp
     250          342 :       IF (nspins == 2) focc = 0.5_dp
     251          754 :       DO ispin = 1, nspins
     252          412 :          CALL dbcsr_set(matrix_px1(ispin)%matrix, 0.0_dp)
     253          412 :          CALL cp_fm_get_info(evect(ispin), ncol_global=norb)
     254              :          CALL cp_dbcsr_plus_fm_fm_t(matrix_px1(ispin)%matrix, &
     255              :                                     matrix_v=evect(ispin), &
     256              :                                     matrix_g=gs_mos(ispin)%mos_occ, &
     257          412 :                                     ncol=norb, alpha=2.0_dp*focc, symmetry_mode=1)
     258              : 
     259          412 :          CALL dbcsr_set(matrix_px1_asymm(ispin)%matrix, 0.0_dp)
     260          412 :          CALL cp_fm_get_info(evect(ispin), ncol_global=norb)
     261              :          CALL cp_dbcsr_plus_fm_fm_t(matrix_px1_asymm(ispin)%matrix, &
     262              :                                     matrix_v=gs_mos(ispin)%mos_occ, &
     263              :                                     matrix_g=evect(ispin), &
     264              :                                     ncol=norb, alpha=2.0_dp*focc, &
     265         1166 :                                     symmetry_mode=-1)
     266              :       END DO
     267              :       !
     268          342 :       CALL get_qs_env(qs_env, ks_env=ks_env, pw_env=pw_env, para_env=para_env)
     269              :       !
     270          342 :       NULLIFY (hartree_local, local_rho_set, local_rho_set_admm)
     271          342 :       IF (gapw .OR. gapw_xc) THEN
     272           58 :          IF (nspins == 2) THEN
     273            0 :             DO ispin = 1, nspins
     274            0 :                CALL dbcsr_scale(matrix_px1(ispin)%matrix, 2.0_dp)
     275              :             END DO
     276              :          END IF
     277              :          CALL get_qs_env(qs_env, &
     278              :                          atomic_kind_set=atomic_kind_set, &
     279           58 :                          qs_kind_set=qs_kind_set)
     280           58 :          CALL local_rho_set_create(local_rho_set)
     281              :          CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, &
     282           58 :                                           qs_kind_set, dft_control, para_env)
     283           58 :          IF (gapw) THEN
     284           48 :             CALL get_qs_env(qs_env, natom=natom)
     285              :             CALL init_rho0(local_rho_set, qs_env, dft_control%qs_control%gapw_control, &
     286           48 :                            zcore=0.0_dp)
     287           48 :             CALL rho0_s_grid_create(pw_env, local_rho_set%rho0_mpole)
     288           48 :             CALL hartree_local_create(hartree_local)
     289           48 :             CALL init_coulomb_local(hartree_local, natom)
     290              :          END IF
     291              : 
     292           58 :          CALL get_qs_env(qs_env=qs_env, oce=oce, sap_oce=sap_oce, sab_orb=sab)
     293           58 :          CALL create_oce_set(oce)
     294           58 :          CALL get_qs_env(qs_env=qs_env, nkind=nkind, particle_set=particle_set)
     295           58 :          CALL allocate_oce_set(oce, nkind)
     296           58 :          eps_fit = dft_control%qs_control%gapw_control%eps_fit
     297              :          CALL build_oce_matrices(oce%intac, .TRUE., 1, qs_kind_set, particle_set, &
     298           58 :                                  sap_oce, eps_fit)
     299           58 :          CALL set_qs_env(qs_env, oce=oce)
     300              : 
     301           58 :          mpga(1:nspins, 1:1) => matrix_px1(1:nspins)
     302              :          CALL calculate_rho_atom_coeff(qs_env, mpga, local_rho_set%rho_atom_set, &
     303           58 :                                        qs_kind_set, oce, sab, para_env)
     304           58 :          CALL prepare_gapw_den(qs_env, local_rho_set, do_rho0=gapw)
     305              :          !
     306           58 :          CALL local_rho_set_create(local_rho_set_f)
     307              :          CALL allocate_rho_atom_internals(local_rho_set_f%rho_atom_set, atomic_kind_set, &
     308           58 :                                           qs_kind_set, dft_control, para_env)
     309              :          CALL calculate_rho_atom_coeff(qs_env, mpga, local_rho_set_f%rho_atom_set, &
     310           58 :                                        qs_kind_set, oce, sab, para_env)
     311           58 :          CALL prepare_gapw_den(qs_env, local_rho_set_f, do_rho0=.FALSE.)
     312              :          !
     313           58 :          CALL local_rho_set_create(local_rho_set_g)
     314              :          CALL allocate_rho_atom_internals(local_rho_set_g%rho_atom_set, atomic_kind_set, &
     315           58 :                                           qs_kind_set, dft_control, para_env)
     316              :          CALL calculate_rho_atom_coeff(qs_env, mpga, local_rho_set_g%rho_atom_set, &
     317           58 :                                        qs_kind_set, oce, sab, para_env)
     318           58 :          CALL prepare_gapw_den(qs_env, local_rho_set_g, do_rho0=.FALSE.)
     319           58 :          IF (nspins == 2) THEN
     320            0 :             DO ispin = 1, nspins
     321            0 :                CALL dbcsr_scale(matrix_px1(ispin)%matrix, 0.5_dp)
     322              :             END DO
     323              :          END IF
     324              :       END IF
     325              :       !
     326          342 :       IF (do_admm) THEN
     327           64 :          CALL get_qs_env(qs_env, admm_env=admm_env)
     328           64 :          nao_aux = admm_env%nao_aux_fit
     329           64 :          nao = admm_env%nao_orb
     330              :          !
     331          132 :          DO ispin = 1, nspins
     332           68 :             CALL copy_dbcsr_to_fm(matrix_px1(ispin)%matrix, admm_env%work_orb_orb)
     333              :             CALL parallel_gemm('N', 'N', nao_aux, nao, nao, &
     334              :                                1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
     335           68 :                                admm_env%work_aux_orb)
     336              :             CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, &
     337              :                                1.0_dp, admm_env%work_aux_orb, admm_env%A, 0.0_dp, &
     338           68 :                                admm_env%work_aux_aux)
     339              :             CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, matrix_px1_admm(ispin)%matrix, &
     340           68 :                                   keep_sparsity=.TRUE.)
     341              : 
     342           68 :             CALL copy_dbcsr_to_fm(matrix_px1_asymm(ispin)%matrix, admm_env%work_orb_orb)
     343              :             CALL parallel_gemm('N', 'N', nao_aux, nao, nao, &
     344              :                                1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
     345           68 :                                admm_env%work_aux_orb)
     346              :             CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, &
     347              :                                1.0_dp, admm_env%work_aux_orb, admm_env%A, 0.0_dp, &
     348           68 :                                admm_env%work_aux_aux)
     349              :             CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, matrix_px1_admm_asymm(ispin)%matrix, &
     350          132 :                                   keep_sparsity=.TRUE.)
     351              :          END DO
     352              :          !
     353           64 :          IF (admm_env%do_gapw) THEN
     354           10 :             IF (do_admm .AND. tddfpt_control%admm_xc_correction) THEN
     355            8 :                IF (admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
     356              :                   ! nothing to do
     357              :                ELSE
     358            2 :                   CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
     359            2 :                   CALL local_rho_set_create(local_rho_set_admm)
     360              :                   CALL allocate_rho_atom_internals(local_rho_set_admm%rho_atom_set, atomic_kind_set, &
     361            2 :                                                    admm_env%admm_gapw_env%admm_kind_set, dft_control, para_env)
     362            2 :                   mpga(1:nspins, 1:1) => matrix_px1_admm(1:nspins)
     363            2 :                   CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
     364              :                   CALL calculate_rho_atom_coeff(qs_env, mpga, local_rho_set_admm%rho_atom_set, &
     365              :                                                 admm_env%admm_gapw_env%admm_kind_set, &
     366            2 :                                                 admm_env%admm_gapw_env%oce, sab_aux_fit, para_env)
     367              :                   CALL prepare_gapw_den(qs_env, local_rho_set_admm, do_rho0=.FALSE., &
     368            2 :                                         kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
     369              :                   !
     370            2 :                   CALL local_rho_set_create(local_rho_set_f_admm)
     371              :                   CALL allocate_rho_atom_internals(local_rho_set_f_admm%rho_atom_set, atomic_kind_set, &
     372            2 :                                                    admm_env%admm_gapw_env%admm_kind_set, dft_control, para_env)
     373              :                   CALL calculate_rho_atom_coeff(qs_env, mpga, local_rho_set_f_admm%rho_atom_set, &
     374              :                                                 admm_env%admm_gapw_env%admm_kind_set, &
     375            2 :                                                 admm_env%admm_gapw_env%oce, sab_aux_fit, para_env)
     376              :                   CALL prepare_gapw_den(qs_env, local_rho_set_f_admm, do_rho0=.FALSE., &
     377            2 :                                         kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
     378              :                   !
     379            2 :                   CALL local_rho_set_create(local_rho_set_g_admm)
     380              :                   CALL allocate_rho_atom_internals(local_rho_set_g_admm%rho_atom_set, atomic_kind_set, &
     381            2 :                                                    admm_env%admm_gapw_env%admm_kind_set, dft_control, para_env)
     382              :                   CALL calculate_rho_atom_coeff(qs_env, mpga, local_rho_set_g_admm%rho_atom_set, &
     383              :                                                 admm_env%admm_gapw_env%admm_kind_set, &
     384            2 :                                                 admm_env%admm_gapw_env%oce, sab_aux_fit, para_env)
     385              :                   CALL prepare_gapw_den(qs_env, local_rho_set_g_admm, do_rho0=.FALSE., &
     386            2 :                                         kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
     387              :                END IF
     388              :             END IF
     389              :          END IF
     390              :       END IF
     391              :       !
     392              :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
     393          342 :                       poisson_env=poisson_env)
     394              : 
     395         2534 :       ALLOCATE (rhox_r(nspins), rhox_g(nspins))
     396          754 :       DO ispin = 1, nspins
     397          412 :          CALL auxbas_pw_pool%create_pw(rhox_r(ispin))
     398          754 :          CALL auxbas_pw_pool%create_pw(rhox_g(ispin))
     399              :       END DO
     400          342 :       CALL auxbas_pw_pool%create_pw(rhox_tot_gspace)
     401              : 
     402          342 :       CALL pw_zero(rhox_tot_gspace)
     403          754 :       DO ispin = 1, nspins
     404          412 :          IF (nspins == 2) CALL dbcsr_scale(matrix_px1(ispin)%matrix, 2.0_dp)
     405              :          CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_px1(ispin)%matrix, &
     406              :                                  rho=rhox_r(ispin), rho_gspace=rhox_g(ispin), &
     407          412 :                                  soft_valid=gapw)
     408          412 :          CALL pw_axpy(rhox_g(ispin), rhox_tot_gspace)
     409          754 :          IF (nspins == 2) CALL dbcsr_scale(matrix_px1(ispin)%matrix, 0.5_dp)
     410              :       END DO
     411              : 
     412          342 :       IF (gapw_xc) THEN
     413           50 :          ALLOCATE (rhoxx_r(nspins), rhoxx_g(nspins))
     414           20 :          DO ispin = 1, nspins
     415           10 :             CALL auxbas_pw_pool%create_pw(rhoxx_r(ispin))
     416           20 :             CALL auxbas_pw_pool%create_pw(rhoxx_g(ispin))
     417              :          END DO
     418           20 :          DO ispin = 1, nspins
     419           10 :             IF (nspins == 2) CALL dbcsr_scale(matrix_px1(ispin)%matrix, 2.0_dp)
     420              :             CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_px1(ispin)%matrix, &
     421              :                                     rho=rhoxx_r(ispin), rho_gspace=rhoxx_g(ispin), &
     422           10 :                                     soft_valid=gapw_xc)
     423           20 :             IF (nspins == 2) CALL dbcsr_scale(matrix_px1(ispin)%matrix, 0.5_dp)
     424              :          END DO
     425              :       END IF
     426              : 
     427          342 :       CALL get_qs_env(qs_env, matrix_s=matrix_s, force=force)
     428              : 
     429          342 :       IF (.NOT. is_rks_triplets) THEN
     430          304 :          CALL auxbas_pw_pool%create_pw(xv_hartree_rspace)
     431          304 :          CALL auxbas_pw_pool%create_pw(xv_hartree_gspace)
     432              :          ! calculate associated hartree potential
     433          304 :          IF (gapw) THEN
     434           40 :             CALL pw_axpy(local_rho_set%rho0_mpole%rho0_s_gs, rhox_tot_gspace)
     435              :          END IF
     436              :          CALL pw_poisson_solve(poisson_env, rhox_tot_gspace, xehartree, &
     437          304 :                                xv_hartree_gspace)
     438          304 :          CALL pw_transfer(xv_hartree_gspace, xv_hartree_rspace)
     439          304 :          CALL pw_scale(xv_hartree_rspace, xv_hartree_rspace%pw_grid%dvol)
     440              :          !
     441          430 :          IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
     442          304 :          NULLIFY (matrix_hx)
     443          304 :          CALL dbcsr_allocate_matrix_set(matrix_hx, nspins)
     444          678 :          DO ispin = 1, nspins
     445          374 :             ALLOCATE (matrix_hx(ispin)%matrix)
     446          374 :             CALL dbcsr_create(matrix_hx(ispin)%matrix, template=matrix_s(1)%matrix)
     447          374 :             CALL dbcsr_copy(matrix_hx(ispin)%matrix, matrix_s(1)%matrix)
     448          374 :             CALL dbcsr_set(matrix_hx(ispin)%matrix, 0.0_dp)
     449              :             CALL integrate_v_rspace(qs_env=qs_env, v_rspace=xv_hartree_rspace, &
     450              :                                     pmat=matrix_px1(ispin), hmat=matrix_hx(ispin), &
     451          678 :                                     gapw=gapw, calculate_forces=.TRUE.)
     452              :          END DO
     453          304 :          IF (debug_forces) THEN
     454          168 :             fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
     455           42 :             CALL para_env%sum(fodeb)
     456           42 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Px*dKh[X]   ", fodeb
     457              :          END IF
     458          304 :          IF (gapw) THEN
     459          148 :             IF (debug_forces) fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1)
     460              :             CALL Vh_1c_gg_integrals(qs_env, eh1c, hartree_local%ecoul_1c, local_rho_set, para_env, tddft=.TRUE., &
     461           40 :                                     core_2nd=.TRUE.)
     462           40 :             IF (nspins == 1) THEN
     463           40 :                kval = 1.0_dp
     464              :             ELSE
     465            0 :                kval = 0.5_dp
     466              :             END IF
     467              :             CALL integrate_vhg0_rspace(qs_env, xv_hartree_rspace, para_env, calculate_forces=.TRUE., &
     468           40 :                                        local_rho_set=local_rho_set, kforce=kval)
     469           40 :             IF (debug_forces) THEN
     470          144 :                fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1) - fodeb(1:3)
     471           36 :                CALL para_env%sum(fodeb)
     472           36 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Px*dKh[X]PAWg0", fodeb
     473              :             END IF
     474          148 :             IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
     475              :             CALL update_ks_atom(qs_env, matrix_hx, matrix_px1, forces=.TRUE., &
     476           40 :                                 rho_atom_external=local_rho_set%rho_atom_set)
     477           40 :             IF (debug_forces) THEN
     478          144 :                fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
     479           36 :                CALL para_env%sum(fodeb)
     480           36 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Px*dKh[X]PAW ", fodeb
     481              :             END IF
     482              :          END IF
     483              :       END IF
     484              : 
     485              :       ! XC
     486          342 :       IF (full_kernel%do_exck) THEN
     487            0 :          CPABORT("NYA")
     488              :       END IF
     489          342 :       NULLIFY (fxc_rho, fxc_tau, gxc_rho, gxc_tau)
     490          342 :       xc_section => full_kernel%xc_section
     491              :       CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
     492          342 :                                 i_val=myfun)
     493          342 :       IF (myfun /= xc_none) THEN
     494          232 :          SELECT CASE (ex_env%xc_kernel_method)
     495              :          CASE (xc_kernel_method_best)
     496              :             do_analytic = .TRUE.
     497            0 :             do_numeric = .TRUE.
     498              :          CASE (xc_kernel_method_analytic)
     499            0 :             do_analytic = .TRUE.
     500            0 :             do_numeric = .FALSE.
     501              :          CASE (xc_kernel_method_numeric)
     502            0 :             do_analytic = .FALSE.
     503            0 :             do_numeric = .TRUE.
     504              :          CASE DEFAULT
     505          232 :             CPABORT("invalid xc_kernel_method")
     506              :          END SELECT
     507          232 :          order = ex_env%diff_order
     508          232 :          eps_delta = ex_env%eps_delta_rho
     509              : 
     510          232 :          IF (gapw_xc) THEN
     511           10 :             CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, rho_xc=rho)
     512              :          ELSE
     513          222 :             CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, rho=rho)
     514              :          END IF
     515          232 :          CALL qs_rho_get(rho, rho_ao=matrix_p)
     516              :          NULLIFY (rhox)
     517          232 :          ALLOCATE (rhox)
     518          232 :          CALL qs_rho_create(rhox)
     519          232 :          IF (gapw_xc) THEN
     520              :             CALL qs_rho_set(rho_struct=rhox, rho_ao=matrix_px1, rho_r=rhoxx_r, rho_g=rhoxx_g, &
     521           10 :                             rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
     522              :          ELSE
     523              :             CALL qs_rho_set(rho_struct=rhox, rho_ao=matrix_px1, rho_r=rhox_r, rho_g=rhox_g, &
     524          222 :                             rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
     525              :          END IF
     526          232 :          IF (do_analytic .AND. .NOT. do_numeric) THEN
     527            0 :             CPABORT("Analytic 3rd EXC derivatives not available")
     528          232 :          ELSEIF (do_numeric) THEN
     529          232 :             IF (do_analytic) THEN
     530              :                CALL qs_fgxc_gdiff(ks_env, rho, rhox, xc_section, order, eps_delta, is_rks_triplets, &
     531          232 :                                   fxc_rho, fxc_tau, gxc_rho, gxc_tau)
     532              :             ELSE
     533              :                CALL qs_fgxc_create(ks_env, rho, rhox, xc_section, order, is_rks_triplets, &
     534            0 :                                    fxc_rho, fxc_tau, gxc_rho, gxc_tau)
     535              :             END IF
     536              :          ELSE
     537            0 :             CPABORT("FHXC forces analytic/numeric")
     538              :          END IF
     539              : 
     540              :          ! Well, this is a hack :-(
     541              :          ! When qs_rho_set() was called on rhox it assumed ownership of the passed arrays.
     542              :          ! However, these arrays actually belong to ex_env. Hence, we can not call qs_rho_release()
     543              :          ! because this would release the arrays. Instead we're simply going to deallocate rhox.
     544          232 :          DEALLOCATE (rhox)
     545              : 
     546          232 :          IF (nspins == 2) THEN
     547          108 :             DO ispin = 1, nspins
     548           72 :                CALL pw_scale(gxc_rho(ispin), 0.5_dp)
     549          108 :                IF (ASSOCIATED(gxc_tau)) CALL pw_scale(gxc_tau(ispin), 0.5_dp)
     550              :             END DO
     551              :          END IF
     552          232 :          IF (gapw .OR. gapw_xc) THEN
     553           50 :             IF (do_analytic .AND. .NOT. do_numeric) THEN
     554            0 :                CPABORT("Analytic 3rd EXC derivatives not available")
     555           50 :             ELSEIF (do_numeric) THEN
     556           50 :                IF (do_analytic) THEN
     557              :                   CALL gfxc_atom_diff(qs_env, ex_env%local_rho_set%rho_atom_set, &
     558              :                                       local_rho_set_f%rho_atom_set, local_rho_set_g%rho_atom_set, &
     559           50 :                                       qs_kind_set, xc_section, is_rks_triplets, order, eps_delta)
     560              :                ELSE
     561              :                   CALL calculate_gfxc_atom(qs_env, ex_env%local_rho_set%rho_atom_set, &
     562              :                                            local_rho_set_f%rho_atom_set, local_rho_set_g%rho_atom_set, &
     563            0 :                                            qs_kind_set, xc_section, is_rks_triplets, order)
     564              :                END IF
     565              :             ELSE
     566            0 :                CPABORT("FHXC forces analytic/numeric")
     567              :             END IF
     568              :          END IF
     569              : 
     570          358 :          IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
     571          232 :          NULLIFY (matrix_fx)
     572          232 :          CALL dbcsr_allocate_matrix_set(matrix_fx, nspins)
     573          500 :          DO ispin = 1, nspins
     574          268 :             ALLOCATE (matrix_fx(ispin)%matrix)
     575          268 :             CALL dbcsr_create(matrix_fx(ispin)%matrix, template=matrix_s(1)%matrix)
     576          268 :             CALL dbcsr_copy(matrix_fx(ispin)%matrix, matrix_s(1)%matrix)
     577          268 :             CALL dbcsr_set(matrix_fx(ispin)%matrix, 0.0_dp)
     578          268 :             CALL pw_scale(fxc_rho(ispin), fxc_rho(ispin)%pw_grid%dvol)
     579              :             CALL integrate_v_rspace(qs_env=qs_env, v_rspace=fxc_rho(ispin), &
     580              :                                     pmat=matrix_px1(ispin), hmat=matrix_fx(ispin), &
     581          718 :                                     gapw=(gapw .OR. gapw_xc), calculate_forces=.TRUE.)
     582              :          END DO
     583          232 :          IF (debug_forces) THEN
     584          168 :             fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
     585           42 :             CALL para_env%sum(fodeb)
     586           42 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Px*dKf[X]   ", fodeb
     587              :          END IF
     588              : 
     589          358 :          IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
     590          232 :          NULLIFY (matrix_gx)
     591          232 :          CALL dbcsr_allocate_matrix_set(matrix_gx, nspins)
     592          500 :          DO ispin = 1, nspins
     593          268 :             ALLOCATE (matrix_gx(ispin)%matrix)
     594          268 :             CALL dbcsr_create(matrix_gx(ispin)%matrix, template=matrix_s(1)%matrix)
     595          268 :             CALL dbcsr_copy(matrix_gx(ispin)%matrix, matrix_s(1)%matrix)
     596          268 :             CALL dbcsr_set(matrix_gx(ispin)%matrix, 0.0_dp)
     597          268 :             CALL pw_scale(gxc_rho(ispin), gxc_rho(ispin)%pw_grid%dvol)
     598          268 :             CALL pw_scale(gxc_rho(ispin), 0.5_dp)
     599              :             CALL integrate_v_rspace(qs_env=qs_env, v_rspace=gxc_rho(ispin), &
     600              :                                     pmat=matrix_p(ispin), hmat=matrix_gx(ispin), &
     601          486 :                                     gapw=(gapw .OR. gapw_xc), calculate_forces=.TRUE.)
     602          500 :             CALL dbcsr_scale(matrix_gx(ispin)%matrix, 2.0_dp)
     603              :          END DO
     604          232 :          IF (debug_forces) THEN
     605          168 :             fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
     606           42 :             CALL para_env%sum(fodeb)
     607           42 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dKg[X]   ", fodeb
     608              :          END IF
     609          232 :          CALL qs_fgxc_release(ks_env, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
     610              : 
     611          232 :          IF (gapw .OR. gapw_xc) THEN
     612          176 :             IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
     613              :             CALL update_ks_atom(qs_env, matrix_fx, matrix_px1, forces=.TRUE., tddft=.TRUE., &
     614              :                                 rho_atom_external=local_rho_set_f%rho_atom_set, &
     615           50 :                                 kintegral=1.0_dp, kforce=1.0_dp)
     616           50 :             IF (debug_forces) THEN
     617          168 :                fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
     618           42 :                CALL para_env%sum(fodeb)
     619           42 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Px*dKf[X]PAW ", fodeb
     620              :             END IF
     621          176 :             IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
     622           50 :             IF (nspins == 1) THEN
     623              :                CALL update_ks_atom(qs_env, matrix_gx, matrix_p, forces=.TRUE., tddft=.TRUE., &
     624              :                                    rho_atom_external=local_rho_set_g%rho_atom_set, &
     625           50 :                                    kscale=0.5_dp)
     626              :             ELSE
     627              :                CALL update_ks_atom(qs_env, matrix_gx, matrix_p, forces=.TRUE., &
     628              :                                    rho_atom_external=local_rho_set_g%rho_atom_set, &
     629            0 :                                    kintegral=0.5_dp, kforce=0.25_dp)
     630              :             END IF
     631           50 :             IF (debug_forces) THEN
     632          168 :                fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
     633           42 :                CALL para_env%sum(fodeb)
     634           42 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Px*dKg[X]PAW ", fodeb
     635              :             END IF
     636              :          END IF
     637              :       END IF
     638              : 
     639              :       ! ADMM XC correction Exc[rho_admm]
     640          342 :       IF (do_admm .AND. tddfpt_control%admm_xc_correction) THEN
     641           52 :          IF (admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
     642              :             ! nothing to do
     643              :          ELSE
     644           32 :             IF (.NOT. tddfpt_control%admm_symm) THEN
     645            0 :                CALL cp_warn(__LOCATION__, "Forces need symmetric ADMM kernel corrections")
     646            0 :                CPABORT("ADMM KERNEL CORRECTION")
     647              :             END IF
     648           32 :             xc_section => admm_env%xc_section_aux
     649              :             CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit, matrix_s_aux_fit=matrix_s_aux_fit, &
     650           32 :                               task_list_aux_fit=task_list)
     651           32 :             basis_type = "AUX_FIT"
     652           32 :             IF (admm_env%do_gapw) THEN
     653            2 :                basis_type = "AUX_FIT_SOFT"
     654            2 :                task_list => admm_env%admm_gapw_env%task_list
     655              :             END IF
     656              :             !
     657           32 :             NULLIFY (mfx, mgx)
     658           32 :             CALL dbcsr_allocate_matrix_set(mfx, nspins)
     659           32 :             CALL dbcsr_allocate_matrix_set(mgx, nspins)
     660           64 :             DO ispin = 1, nspins
     661           32 :                ALLOCATE (mfx(ispin)%matrix, mgx(ispin)%matrix)
     662           32 :                CALL dbcsr_create(mfx(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
     663           32 :                CALL dbcsr_copy(mfx(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
     664           32 :                CALL dbcsr_set(mfx(ispin)%matrix, 0.0_dp)
     665           32 :                CALL dbcsr_create(mgx(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
     666           32 :                CALL dbcsr_copy(mgx(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
     667           64 :                CALL dbcsr_set(mgx(ispin)%matrix, 0.0_dp)
     668              :             END DO
     669              : 
     670              :             ! ADMM density and response density
     671           32 :             NULLIFY (rho_g_aux, rho_r_aux, rhox_g_aux, rhox_r_aux)
     672           32 :             CALL qs_rho_get(rho_aux_fit, rho_r=rho_r_aux, rho_g=rho_g_aux)
     673           32 :             CALL qs_rho_get(rho_aux_fit, rho_ao=matrix_p_admm)
     674              :             ! rhox_aux
     675          224 :             ALLOCATE (rhox_r_aux(nspins), rhox_g_aux(nspins))
     676           64 :             DO ispin = 1, nspins
     677           32 :                CALL auxbas_pw_pool%create_pw(rhox_r_aux(ispin))
     678           64 :                CALL auxbas_pw_pool%create_pw(rhox_g_aux(ispin))
     679              :             END DO
     680           64 :             DO ispin = 1, nspins
     681              :                CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_px1_admm(ispin)%matrix, &
     682              :                                        rho=rhox_r_aux(ispin), rho_gspace=rhox_g_aux(ispin), &
     683              :                                        basis_type=basis_type, &
     684           64 :                                        task_list_external=task_list)
     685              :             END DO
     686              :             !
     687              :             NULLIFY (rhox_aux)
     688           32 :             ALLOCATE (rhox_aux)
     689           32 :             CALL qs_rho_create(rhox_aux)
     690              :             CALL qs_rho_set(rho_struct=rhox_aux, rho_ao=matrix_px1_admm, &
     691              :                             rho_r=rhox_r_aux, rho_g=rhox_g_aux, &
     692           32 :                             rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
     693           32 :             IF (do_analytic .AND. .NOT. do_numeric) THEN
     694            0 :                CPABORT("Analytic 3rd derivatives of EXC not available")
     695           32 :             ELSEIF (do_numeric) THEN
     696           32 :                IF (do_analytic) THEN
     697              :                   CALL qs_fgxc_gdiff(ks_env, rho_aux_fit, rhox_aux, xc_section, order, eps_delta, &
     698           32 :                                      is_rks_triplets, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
     699              :                ELSE
     700              :                   CALL qs_fgxc_create(ks_env, rho_aux_fit, rhox_aux, xc_section, &
     701            0 :                                       order, is_rks_triplets, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
     702              :                END IF
     703              :             ELSE
     704            0 :                CPABORT("FHXC forces analytic/numeric")
     705              :             END IF
     706              : 
     707              :             ! Well, this is a hack :-(
     708              :             ! When qs_rho_set() was called on rhox_aux it assumed ownership of the passed arrays.
     709              :             ! However, these arrays actually belong to ex_env. Hence, we can not call qs_rho_release()
     710              :             ! because this would release the arrays. Instead we're simply going to deallocate rhox_aux.
     711           32 :             DEALLOCATE (rhox_aux)
     712              : 
     713           64 :             DO ispin = 1, nspins
     714           32 :                CALL auxbas_pw_pool%give_back_pw(rhox_r_aux(ispin))
     715           64 :                CALL auxbas_pw_pool%give_back_pw(rhox_g_aux(ispin))
     716              :             END DO
     717           32 :             DEALLOCATE (rhox_r_aux, rhox_g_aux)
     718           32 :             fscal = 1.0_dp
     719           32 :             IF (nspins == 2) fscal = 2.0_dp
     720              :             !
     721           38 :             IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
     722           64 :             DO ispin = 1, nspins
     723           32 :                CALL pw_scale(fxc_rho(ispin), fxc_rho(ispin)%pw_grid%dvol)
     724              :                CALL integrate_v_rspace(qs_env=qs_env, v_rspace=fxc_rho(ispin), &
     725              :                                        hmat=mfx(ispin), &
     726              :                                        pmat=matrix_px1_admm(ispin), &
     727              :                                        basis_type=basis_type, &
     728              :                                        calculate_forces=.TRUE., &
     729              :                                        force_adm=fscal, &
     730           64 :                                        task_list_external=task_list)
     731              :             END DO
     732           32 :             IF (debug_forces) THEN
     733            8 :                fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
     734            2 :                CALL para_env%sum(fodeb)
     735            2 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Px*dKf[X]ADMM", fodeb
     736              :             END IF
     737              : 
     738           38 :             IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
     739           64 :             DO ispin = 1, nspins
     740           32 :                CALL pw_scale(gxc_rho(ispin), gxc_rho(ispin)%pw_grid%dvol)
     741           32 :                CALL pw_scale(gxc_rho(ispin), 0.5_dp)
     742              :                CALL integrate_v_rspace(qs_env=qs_env, v_rspace=gxc_rho(ispin), &
     743              :                                        hmat=mgx(ispin), &
     744              :                                        pmat=matrix_p_admm(ispin), &
     745              :                                        basis_type=basis_type, &
     746              :                                        calculate_forces=.TRUE., &
     747              :                                        force_adm=fscal, &
     748           32 :                                        task_list_external=task_list)
     749           64 :                CALL dbcsr_scale(mgx(ispin)%matrix, 2.0_dp)
     750              :             END DO
     751           32 :             IF (debug_forces) THEN
     752            8 :                fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
     753            2 :                CALL para_env%sum(fodeb)
     754            2 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dKg[X]ADMM", fodeb
     755              :             END IF
     756           32 :             CALL qs_fgxc_release(ks_env, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
     757              :             !
     758           32 :             IF (admm_env%do_gapw) THEN
     759            2 :                CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
     760            2 :                rho_atom_set => admm_env%admm_gapw_env%local_rho_set%rho_atom_set
     761            2 :                rho_atom_set_f => local_rho_set_f_admm%rho_atom_set
     762            2 :                rho_atom_set_g => local_rho_set_g_admm%rho_atom_set
     763              : 
     764            2 :                IF (do_analytic .AND. .NOT. do_numeric) THEN
     765            0 :                   CPABORT("Analytic 3rd EXC derivatives not available")
     766            2 :                ELSEIF (do_numeric) THEN
     767            2 :                   IF (do_analytic) THEN
     768              :                      CALL gfxc_atom_diff(qs_env, rho_atom_set, &
     769              :                                          rho_atom_set_f, rho_atom_set_g, &
     770              :                                          admm_env%admm_gapw_env%admm_kind_set, xc_section, &
     771            2 :                                          is_rks_triplets, order, eps_delta)
     772              :                   ELSE
     773              :                      CALL calculate_gfxc_atom(qs_env, rho_atom_set, &
     774              :                                               rho_atom_set_f, rho_atom_set_g, &
     775              :                                               admm_env%admm_gapw_env%admm_kind_set, xc_section, &
     776            0 :                                               is_rks_triplets, order)
     777              :                   END IF
     778              :                ELSE
     779            0 :                   CPABORT("FHXC forces analytic/numeric")
     780              :                END IF
     781              : 
     782            8 :                IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
     783            2 :                IF (nspins == 1) THEN
     784              :                   CALL update_ks_atom(qs_env, mfx, matrix_px1_admm, forces=.TRUE., &
     785              :                                       rho_atom_external=rho_atom_set_f, &
     786              :                                       kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
     787              :                                       oce_external=admm_env%admm_gapw_env%oce, sab_external=sab_aux_fit, &
     788            2 :                                       kintegral=2.0_dp, kforce=0.5_dp)
     789              :                ELSE
     790              :                   CALL update_ks_atom(qs_env, mfx, matrix_px1_admm, forces=.TRUE., &
     791              :                                       rho_atom_external=rho_atom_set_f, &
     792              :                                       kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
     793              :                                       oce_external=admm_env%admm_gapw_env%oce, sab_external=sab_aux_fit, &
     794            0 :                                       kintegral=2.0_dp, kforce=1.0_dp)
     795              :                END IF
     796            2 :                IF (debug_forces) THEN
     797            8 :                   fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
     798            2 :                   CALL para_env%sum(fodeb)
     799            2 :                   IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Px*dKf[X]ADMM-PAW ", fodeb
     800              :                END IF
     801            8 :                IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
     802            2 :                IF (nspins == 1) THEN
     803              :                   CALL update_ks_atom(qs_env, mgx, matrix_p, forces=.TRUE., &
     804              :                                       rho_atom_external=rho_atom_set_g, &
     805              :                                       kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
     806              :                                       oce_external=admm_env%admm_gapw_env%oce, sab_external=sab_aux_fit, &
     807            2 :                                       kintegral=1.0_dp, kforce=0.5_dp)
     808              :                ELSE
     809              :                   CALL update_ks_atom(qs_env, mgx, matrix_p, forces=.TRUE., &
     810              :                                       rho_atom_external=rho_atom_set_g, &
     811              :                                       kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
     812              :                                       oce_external=admm_env%admm_gapw_env%oce, sab_external=sab_aux_fit, &
     813            0 :                                       kintegral=1.0_dp, kforce=1.0_dp)
     814              :                END IF
     815            2 :                IF (debug_forces) THEN
     816            8 :                   fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
     817            2 :                   CALL para_env%sum(fodeb)
     818            2 :                   IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dKg[X]ADMM-PAW ", fodeb
     819              :                END IF
     820              :             END IF
     821              :             !
     822              :             ! A' fx A - Forces
     823              :             !
     824           38 :             IF (debug_forces) fodeb(1:3) = force(1)%overlap_admm(1:3, 1)
     825           32 :             fval = 2.0_dp*REAL(nspins, KIND=dp)
     826           32 :             CALL admm_projection_derivative(qs_env, mfx, matrix_px1, fval)
     827           32 :             IF (debug_forces) THEN
     828            8 :                fodeb(1:3) = force(1)%overlap_admm(1:3, 1) - fodeb(1:3)
     829            2 :                CALL para_env%sum(fodeb)
     830            2 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dfXC(P)*S' ", fodeb
     831              :             END IF
     832           38 :             IF (debug_forces) fodeb(1:3) = force(1)%overlap_admm(1:3, 1)
     833              :             fval = 2.0_dp*REAL(nspins, KIND=dp)
     834           32 :             CALL admm_projection_derivative(qs_env, mgx, matrix_p, fval)
     835           32 :             IF (debug_forces) THEN
     836            8 :                fodeb(1:3) = force(1)%overlap_admm(1:3, 1) - fodeb(1:3)
     837            2 :                CALL para_env%sum(fodeb)
     838            2 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dgXC(P)*S' ", fodeb
     839              :             END IF
     840              :             !
     841              :             ! Add ADMM fx/gx to the full basis fx/gx
     842           32 :             fscal = 1.0_dp
     843           32 :             IF (nspins == 2) fscal = 2.0_dp
     844           32 :             nao = admm_env%nao_orb
     845           32 :             nao_aux = admm_env%nao_aux_fit
     846           32 :             ALLOCATE (dbwork)
     847           32 :             CALL dbcsr_create(dbwork, template=matrix_fx(1)%matrix)
     848           64 :             DO ispin = 1, nspins
     849              :                ! fx
     850              :                CALL cp_dbcsr_sm_fm_multiply(mfx(ispin)%matrix, admm_env%A, &
     851           32 :                                             admm_env%work_aux_orb, nao)
     852              :                CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
     853              :                                   1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
     854           32 :                                   admm_env%work_orb_orb)
     855           32 :                CALL dbcsr_copy(dbwork, matrix_fx(1)%matrix)
     856           32 :                CALL dbcsr_set(dbwork, 0.0_dp)
     857           32 :                CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.TRUE.)
     858           32 :                CALL dbcsr_add(matrix_fx(ispin)%matrix, dbwork, 1.0_dp, fscal)
     859              :                ! gx
     860              :                CALL cp_dbcsr_sm_fm_multiply(mgx(ispin)%matrix, admm_env%A, &
     861           32 :                                             admm_env%work_aux_orb, nao)
     862              :                CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
     863              :                                   1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
     864           32 :                                   admm_env%work_orb_orb)
     865           32 :                CALL dbcsr_set(dbwork, 0.0_dp)
     866           32 :                CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.TRUE.)
     867           64 :                CALL dbcsr_add(matrix_gx(ispin)%matrix, dbwork, 1.0_dp, fscal)
     868              :             END DO
     869           32 :             CALL dbcsr_release(dbwork)
     870           32 :             DEALLOCATE (dbwork)
     871           32 :             CALL dbcsr_deallocate_matrix_set(mfx)
     872           64 :             CALL dbcsr_deallocate_matrix_set(mgx)
     873              : 
     874              :          END IF
     875              :       END IF
     876              : 
     877          754 :       DO ispin = 1, nspins
     878          412 :          CALL auxbas_pw_pool%give_back_pw(rhox_r(ispin))
     879          754 :          CALL auxbas_pw_pool%give_back_pw(rhox_g(ispin))
     880              :       END DO
     881          342 :       DEALLOCATE (rhox_r, rhox_g)
     882          342 :       CALL auxbas_pw_pool%give_back_pw(rhox_tot_gspace)
     883          342 :       IF (gapw_xc) THEN
     884           20 :          DO ispin = 1, nspins
     885           10 :             CALL auxbas_pw_pool%give_back_pw(rhoxx_r(ispin))
     886           20 :             CALL auxbas_pw_pool%give_back_pw(rhoxx_g(ispin))
     887              :          END DO
     888           10 :          DEALLOCATE (rhoxx_r, rhoxx_g)
     889              :       END IF
     890          342 :       IF (.NOT. is_rks_triplets) THEN
     891          304 :          CALL auxbas_pw_pool%give_back_pw(xv_hartree_rspace)
     892          304 :          CALL auxbas_pw_pool%give_back_pw(xv_hartree_gspace)
     893              :       END IF
     894              : 
     895              :       ! HFX
     896          342 :       IF (do_hfx) THEN
     897          124 :          NULLIFY (matrix_hfx, matrix_hfx_asymm)
     898          124 :          CALL dbcsr_allocate_matrix_set(matrix_hfx, nspins)
     899          124 :          CALL dbcsr_allocate_matrix_set(matrix_hfx_asymm, nspins)
     900          260 :          DO ispin = 1, nspins
     901          136 :             ALLOCATE (matrix_hfx(ispin)%matrix)
     902          136 :             CALL dbcsr_create(matrix_hfx(ispin)%matrix, template=matrix_s(1)%matrix)
     903          136 :             CALL dbcsr_copy(matrix_hfx(ispin)%matrix, matrix_s(1)%matrix)
     904          136 :             CALL dbcsr_set(matrix_hfx(ispin)%matrix, 0.0_dp)
     905              : 
     906          136 :             ALLOCATE (matrix_hfx_asymm(ispin)%matrix)
     907              :             CALL dbcsr_create(matrix_hfx_asymm(ispin)%matrix, template=matrix_s(1)%matrix, &
     908          136 :                               matrix_type=dbcsr_type_antisymmetric)
     909          260 :             CALL dbcsr_complete_redistribute(matrix_hfx(ispin)%matrix, matrix_hfx_asymm(ispin)%matrix)
     910              :          END DO
     911              :          !
     912          124 :          xc_section => full_kernel%xc_section
     913          124 :          hfx_section => section_vals_get_subs_vals(xc_section, "HF")
     914          124 :          CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
     915          124 :          CPASSERT(n_rep_hf == 1)
     916              :          CALL section_vals_val_get(hfx_section, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
     917          124 :                                    i_rep_section=1)
     918          124 :          mspin = 1
     919          124 :          IF (hfx_treat_lsd_in_core) mspin = nspins
     920              :          !
     921          124 :          CALL get_qs_env(qs_env=qs_env, x_data=x_data, s_mstruct_changed=s_mstruct_changed)
     922          124 :          distribute_fock_matrix = .TRUE.
     923              :          !
     924          124 :          IF (do_admm) THEN
     925           64 :             CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux_fit)
     926           64 :             NULLIFY (matrix_hfx_admm, matrix_hfx_admm_asymm)
     927           64 :             CALL dbcsr_allocate_matrix_set(matrix_hfx_admm, nspins)
     928           64 :             CALL dbcsr_allocate_matrix_set(matrix_hfx_admm_asymm, nspins)
     929          132 :             DO ispin = 1, nspins
     930           68 :                ALLOCATE (matrix_hfx_admm(ispin)%matrix)
     931           68 :                CALL dbcsr_create(matrix_hfx_admm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
     932           68 :                CALL dbcsr_copy(matrix_hfx_admm(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
     933           68 :                CALL dbcsr_set(matrix_hfx_admm(ispin)%matrix, 0.0_dp)
     934              : 
     935           68 :                ALLOCATE (matrix_hfx_admm_asymm(ispin)%matrix)
     936              :                CALL dbcsr_create(matrix_hfx_admm_asymm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix, &
     937           68 :                                  matrix_type=dbcsr_type_antisymmetric)
     938          132 :                CALL dbcsr_complete_redistribute(matrix_hfx_admm(ispin)%matrix, matrix_hfx_admm_asymm(ispin)%matrix)
     939              :             END DO
     940              :             !
     941           64 :             NULLIFY (mpe, mhe)
     942          520 :             ALLOCATE (mpe(nspins, 1), mhe(nspins, 1))
     943          132 :             DO ispin = 1, nspins
     944           68 :                mhe(ispin, 1)%matrix => matrix_hfx_admm(ispin)%matrix
     945          132 :                mpe(ispin, 1)%matrix => matrix_px1_admm(ispin)%matrix
     946              :             END DO
     947           64 :             IF (x_data(1, 1)%do_hfx_ri) THEN
     948              :                eh1 = 0.0_dp
     949              :                CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhe, eh1, rho_ao=mpe, &
     950              :                                      geometry_did_change=s_mstruct_changed, nspins=nspins, &
     951            6 :                                      hf_fraction=x_data(1, 1)%general_parameter%fraction)
     952              :             ELSE
     953          116 :                DO ispin = 1, mspin
     954              :                   eh1 = 0.0
     955              :                   CALL integrate_four_center(qs_env, x_data, mhe, eh1, mpe, hfx_section, &
     956              :                                              para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
     957          116 :                                              ispin=ispin)
     958              :                END DO
     959              :             END IF
     960              :             !anti-symmetric density
     961          132 :             DO ispin = 1, nspins
     962           68 :                mhe(ispin, 1)%matrix => matrix_hfx_admm_asymm(ispin)%matrix
     963          132 :                mpe(ispin, 1)%matrix => matrix_px1_admm_asymm(ispin)%matrix
     964              :             END DO
     965           64 :             IF (x_data(1, 1)%do_hfx_ri) THEN
     966              :                eh1 = 0.0_dp
     967              :                CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhe, eh1, rho_ao=mpe, &
     968              :                                      geometry_did_change=s_mstruct_changed, nspins=nspins, &
     969            6 :                                      hf_fraction=x_data(1, 1)%general_parameter%fraction)
     970              :             ELSE
     971          116 :                DO ispin = 1, mspin
     972              :                   eh1 = 0.0
     973              :                   CALL integrate_four_center(qs_env, x_data, mhe, eh1, mpe, hfx_section, &
     974              :                                              para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
     975          116 :                                              ispin=ispin)
     976              :                END DO
     977              :             END IF
     978              :             !
     979           64 :             nao = admm_env%nao_orb
     980           64 :             nao_aux = admm_env%nao_aux_fit
     981           64 :             ALLOCATE (dbwork, dbwork_asymm)
     982           64 :             CALL dbcsr_create(dbwork, template=matrix_hfx(1)%matrix)
     983           64 :             CALL dbcsr_create(dbwork_asymm, template=matrix_hfx(1)%matrix, matrix_type=dbcsr_type_antisymmetric)
     984          132 :             DO ispin = 1, nspins
     985              :                CALL cp_dbcsr_sm_fm_multiply(matrix_hfx_admm(ispin)%matrix, admm_env%A, &
     986           68 :                                             admm_env%work_aux_orb, nao)
     987              :                CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
     988              :                                   1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
     989           68 :                                   admm_env%work_orb_orb)
     990           68 :                CALL dbcsr_copy(dbwork, matrix_hfx(1)%matrix)
     991           68 :                CALL dbcsr_set(dbwork, 0.0_dp)
     992           68 :                CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.TRUE.)
     993           68 :                CALL dbcsr_add(matrix_hfx(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
     994              :                !anti-symmetric case
     995              :                CALL cp_dbcsr_sm_fm_multiply(matrix_hfx_admm_asymm(ispin)%matrix, admm_env%A, &
     996           68 :                                             admm_env%work_aux_orb, nao)
     997              :                CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
     998              :                                   1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
     999           68 :                                   admm_env%work_orb_orb)
    1000           68 :                CALL dbcsr_copy(dbwork_asymm, matrix_hfx_asymm(1)%matrix)
    1001           68 :                CALL dbcsr_set(dbwork_asymm, 0.0_dp)
    1002           68 :                CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork_asymm, keep_sparsity=.TRUE.)
    1003          132 :                CALL dbcsr_add(matrix_hfx_asymm(ispin)%matrix, dbwork_asymm, 1.0_dp, 1.0_dp)
    1004              :             END DO
    1005           64 :             CALL dbcsr_release(dbwork)
    1006           64 :             CALL dbcsr_release(dbwork_asymm)
    1007           64 :             DEALLOCATE (dbwork, dbwork_asymm)
    1008              :             ! forces
    1009              :             ! ADMM Projection force
    1010           82 :             IF (debug_forces) fodeb(1:3) = force(1)%overlap_admm(1:3, 1)
    1011           64 :             fval = 4.0_dp*REAL(nspins, KIND=dp)*0.5_dp !0.5 for symm/anti-symm
    1012           64 :             CALL admm_projection_derivative(qs_env, matrix_hfx_admm, matrix_px1, fval)
    1013           64 :             CALL admm_projection_derivative(qs_env, matrix_hfx_admm_asymm, matrix_px1_asymm, fval)
    1014           64 :             IF (debug_forces) THEN
    1015           24 :                fodeb(1:3) = force(1)%overlap_admm(1:3, 1) - fodeb(1:3)
    1016            6 :                CALL para_env%sum(fodeb)
    1017            6 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*Hx(P)*S' ", fodeb
    1018              :             END IF
    1019              :             !
    1020           64 :             use_virial = .FALSE.
    1021           64 :             NULLIFY (mdum)
    1022           64 :             fval = 2.0_dp*REAL(nspins, KIND=dp)*0.5_dp !0.5 factor because of symemtry/anti-symmetry
    1023           82 :             IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
    1024          132 :             DO ispin = 1, nspins
    1025          132 :                mpe(ispin, 1)%matrix => matrix_px1_admm(ispin)%matrix
    1026              :             END DO
    1027           64 :             IF (x_data(1, 1)%do_hfx_ri) THEN
    1028              :                CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
    1029              :                                          x_data(1, 1)%general_parameter%fraction, &
    1030              :                                          rho_ao=mpe, rho_ao_resp=mdum, &
    1031            6 :                                          use_virial=use_virial, rescale_factor=fval)
    1032              :             ELSE
    1033              :                CALL derivatives_four_center(qs_env, mpe, mdum, hfx_section, para_env, 1, use_virial, &
    1034           58 :                                             adiabatic_rescale_factor=fval)
    1035              :             END IF
    1036          132 :             DO ispin = 1, nspins
    1037          132 :                mpe(ispin, 1)%matrix => matrix_px1_admm_asymm(ispin)%matrix
    1038              :             END DO
    1039           64 :             IF (x_data(1, 1)%do_hfx_ri) THEN
    1040              :                CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
    1041              :                                          x_data(1, 1)%general_parameter%fraction, &
    1042              :                                          rho_ao=mpe, rho_ao_resp=mdum, &
    1043            6 :                                          use_virial=use_virial, rescale_factor=fval)
    1044              :             ELSE
    1045              :                CALL derivatives_four_center(qs_env, mpe, mdum, hfx_section, para_env, 1, use_virial, &
    1046           58 :                                             adiabatic_rescale_factor=fval)
    1047              :             END IF
    1048           64 :             IF (debug_forces) THEN
    1049           24 :                fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
    1050            6 :                CALL para_env%sum(fodeb)
    1051            6 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Px*hfx'*Px ", fodeb
    1052              :             END IF
    1053              :             !
    1054           64 :             DEALLOCATE (mpe, mhe)
    1055              :             !
    1056           64 :             CALL dbcsr_deallocate_matrix_set(matrix_hfx_admm)
    1057           64 :             CALL dbcsr_deallocate_matrix_set(matrix_hfx_admm_asymm)
    1058              :          ELSE
    1059           60 :             NULLIFY (mpe, mhe)
    1060          496 :             ALLOCATE (mpe(nspins, 1), mhe(nspins, 1))
    1061          128 :             DO ispin = 1, nspins
    1062           68 :                mhe(ispin, 1)%matrix => matrix_hfx(ispin)%matrix
    1063          128 :                mpe(ispin, 1)%matrix => matrix_px1(ispin)%matrix
    1064              :             END DO
    1065           60 :             IF (x_data(1, 1)%do_hfx_ri) THEN
    1066              :                eh1 = 0.0_dp
    1067              :                CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhe, eh1, rho_ao=mpe, &
    1068              :                                      geometry_did_change=s_mstruct_changed, nspins=nspins, &
    1069           18 :                                      hf_fraction=x_data(1, 1)%general_parameter%fraction)
    1070              :             ELSE
    1071           84 :                DO ispin = 1, mspin
    1072              :                   eh1 = 0.0
    1073              :                   CALL integrate_four_center(qs_env, x_data, mhe, eh1, mpe, hfx_section, &
    1074              :                                              para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
    1075           84 :                                              ispin=ispin)
    1076              :                END DO
    1077              :             END IF
    1078              : 
    1079              :             !anti-symmetric density matrix
    1080          128 :             DO ispin = 1, nspins
    1081           68 :                mhe(ispin, 1)%matrix => matrix_hfx_asymm(ispin)%matrix
    1082          128 :                mpe(ispin, 1)%matrix => matrix_px1_asymm(ispin)%matrix
    1083              :             END DO
    1084           60 :             IF (x_data(1, 1)%do_hfx_ri) THEN
    1085              :                eh1 = 0.0_dp
    1086              :                CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhe, eh1, rho_ao=mpe, &
    1087              :                                      geometry_did_change=s_mstruct_changed, nspins=nspins, &
    1088           18 :                                      hf_fraction=x_data(1, 1)%general_parameter%fraction)
    1089              :             ELSE
    1090           84 :                DO ispin = 1, mspin
    1091              :                   eh1 = 0.0
    1092              :                   CALL integrate_four_center(qs_env, x_data, mhe, eh1, mpe, hfx_section, &
    1093              :                                              para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
    1094           84 :                                              ispin=ispin)
    1095              :                END DO
    1096              :             END IF
    1097              :             ! forces
    1098           60 :             use_virial = .FALSE.
    1099           60 :             NULLIFY (mdum)
    1100           60 :             fval = 2.0_dp*REAL(nspins, KIND=dp)*0.5_dp !extra 0.5 factor because of symmetry/antisymemtry
    1101           84 :             IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
    1102          128 :             DO ispin = 1, nspins
    1103          128 :                mpe(ispin, 1)%matrix => matrix_px1(ispin)%matrix
    1104              :             END DO
    1105           60 :             IF (x_data(1, 1)%do_hfx_ri) THEN
    1106              :                CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
    1107              :                                          x_data(1, 1)%general_parameter%fraction, &
    1108              :                                          rho_ao=mpe, rho_ao_resp=mdum, &
    1109           18 :                                          use_virial=use_virial, rescale_factor=fval)
    1110              :             ELSE
    1111              :                CALL derivatives_four_center(qs_env, mpe, mdum, hfx_section, para_env, 1, use_virial, &
    1112           42 :                                             adiabatic_rescale_factor=fval)
    1113              :             END IF
    1114          128 :             DO ispin = 1, nspins
    1115          128 :                mpe(ispin, 1)%matrix => matrix_px1_asymm(ispin)%matrix
    1116              :             END DO
    1117           60 :             IF (x_data(1, 1)%do_hfx_ri) THEN
    1118              :                CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
    1119              :                                          x_data(1, 1)%general_parameter%fraction, &
    1120              :                                          rho_ao=mpe, rho_ao_resp=mdum, &
    1121           18 :                                          use_virial=use_virial, rescale_factor=fval)
    1122              :             ELSE
    1123              :                CALL derivatives_four_center(qs_env, mpe, mdum, hfx_section, para_env, 1, use_virial, &
    1124           42 :                                             adiabatic_rescale_factor=fval)
    1125              :             END IF
    1126           60 :             IF (debug_forces) THEN
    1127           32 :                fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
    1128            8 :                CALL para_env%sum(fodeb)
    1129            8 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Px*hfx'*Px ", fodeb
    1130              :             END IF
    1131              :             !
    1132           60 :             DEALLOCATE (mpe, mhe)
    1133              :          END IF
    1134          124 :          fval = 2.0_dp*REAL(nspins, KIND=dp)*0.5_dp !extra 0.5 because of symm/antisymm
    1135          260 :          DO ispin = 1, nspins
    1136          136 :             CALL dbcsr_scale(matrix_hfx(ispin)%matrix, fval)
    1137          260 :             CALL dbcsr_scale(matrix_hfx_asymm(ispin)%matrix, fval)
    1138              :          END DO
    1139              :       END IF
    1140              : 
    1141          342 :       IF (gapw .OR. gapw_xc) THEN
    1142           58 :          CALL local_rho_set_release(local_rho_set)
    1143           58 :          CALL local_rho_set_release(local_rho_set_f)
    1144           58 :          CALL local_rho_set_release(local_rho_set_g)
    1145           58 :          IF (gapw) THEN
    1146           48 :             CALL hartree_local_release(hartree_local)
    1147              :          END IF
    1148              :       END IF
    1149          342 :       IF (do_admm) THEN
    1150           64 :          IF (admm_env%do_gapw) THEN
    1151           10 :             IF (tddfpt_control%admm_xc_correction) THEN
    1152            8 :                IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
    1153            2 :                   CALL local_rho_set_release(local_rho_set_admm)
    1154            2 :                   CALL local_rho_set_release(local_rho_set_f_admm)
    1155            2 :                   CALL local_rho_set_release(local_rho_set_g_admm)
    1156              :                END IF
    1157              :             END IF
    1158              :          END IF
    1159              :       END IF
    1160              : 
    1161              :       ! HFX short range
    1162          342 :       IF (do_hfxsr) THEN
    1163            0 :          CPABORT("HFXSR not implemented")
    1164              :       END IF
    1165              :       ! HFX long range
    1166          342 :       IF (do_hfxlr) THEN
    1167            0 :          CPABORT("HFXLR not implemented")
    1168              :       END IF
    1169              : 
    1170          342 :       CALL get_qs_env(qs_env, sab_orb=sab_orb)
    1171          342 :       NULLIFY (matrix_wx1)
    1172          342 :       CALL dbcsr_allocate_matrix_set(matrix_wx1, nspins)
    1173          342 :       cpmos => ex_env%cpmos
    1174          342 :       focc = 2.0_dp
    1175          342 :       IF (nspins == 2) focc = 1.0_dp
    1176          754 :       DO ispin = 1, nspins
    1177          412 :          mos => gs_mos(ispin)%mos_occ
    1178          412 :          CALL cp_fm_get_info(evect(ispin), ncol_global=norb)
    1179          412 :          CALL cp_fm_create(vcvec, mos%matrix_struct, "vcvec")
    1180          412 :          CALL cp_fm_get_info(vcvec, matrix_struct=fm_struct, nrow_global=nao)
    1181              :          CALL cp_fm_struct_create(fm_struct_mat, context=fm_struct%context, nrow_global=norb, &
    1182          412 :                                   ncol_global=norb, para_env=fm_struct%para_env)
    1183          412 :          CALL cp_fm_create(cvcmat, fm_struct_mat)
    1184          412 :          CALL cp_fm_struct_release(fm_struct_mat)
    1185              :          !
    1186          412 :          ALLOCATE (matrix_wx1(ispin)%matrix)
    1187          412 :          CALL dbcsr_create(matrix=matrix_wx1(ispin)%matrix, template=matrix_s(1)%matrix)
    1188          412 :          CALL cp_dbcsr_alloc_block_from_nbl(matrix_wx1(ispin)%matrix, sab_orb)
    1189          412 :          CALL dbcsr_set(matrix_wx1(ispin)%matrix, 0.0_dp)
    1190              :          !
    1191          412 :          IF (.NOT. is_rks_triplets) THEN
    1192              :             CALL cp_dbcsr_sm_fm_multiply(matrix_hx(ispin)%matrix, evect(ispin), &
    1193          374 :                                          cpmos(ispin), norb, alpha=focc, beta=1.0_dp)
    1194          374 :             CALL cp_dbcsr_sm_fm_multiply(matrix_hx(ispin)%matrix, mos, vcvec, norb, alpha=1.0_dp, beta=0.0_dp)
    1195          374 :             CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, mos, vcvec, 0.0_dp, cvcmat)
    1196          374 :             CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, evect(ispin), cvcmat, 0.0_dp, vcvec)
    1197              :             CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), &
    1198          374 :                                          norb, alpha=-focc, beta=1.0_dp)
    1199              :             !
    1200              :             CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=mos, matrix_g=vcvec, &
    1201          374 :                                        ncol=norb, alpha=2.0_dp, symmetry_mode=1)
    1202              :          END IF
    1203              :          !
    1204          412 :          IF (myfun /= xc_none) THEN
    1205              :             CALL cp_dbcsr_sm_fm_multiply(matrix_fx(ispin)%matrix, evect(ispin), &
    1206          268 :                                          cpmos(ispin), norb, alpha=focc, beta=1.0_dp)
    1207          268 :             CALL cp_dbcsr_sm_fm_multiply(matrix_fx(ispin)%matrix, mos, vcvec, norb, alpha=1.0_dp, beta=0.0_dp)
    1208          268 :             CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, mos, vcvec, 0.0_dp, cvcmat)
    1209          268 :             CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, evect(ispin), cvcmat, 0.0_dp, vcvec)
    1210              :             CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), &
    1211          268 :                                          norb, alpha=-focc, beta=1.0_dp)
    1212              :             !
    1213              :             CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=mos, matrix_g=vcvec, &
    1214          268 :                                        ncol=norb, alpha=2.0_dp, symmetry_mode=1)
    1215              :             !
    1216              :             CALL cp_dbcsr_sm_fm_multiply(matrix_gx(ispin)%matrix, mos, &
    1217          268 :                                          cpmos(ispin), norb, alpha=focc, beta=1.0_dp)
    1218              :             !
    1219          268 :             CALL cp_dbcsr_sm_fm_multiply(matrix_gx(ispin)%matrix, mos, vcvec, norb, alpha=1.0_dp, beta=0.0_dp)
    1220          268 :             CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, mos, vcvec, 0.0_dp, cvcmat)
    1221          268 :             CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, mos, cvcmat, 0.0_dp, vcvec)
    1222              :             CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=mos, matrix_g=vcvec, &
    1223          268 :                                        ncol=norb, alpha=1.0_dp, symmetry_mode=1)
    1224              :          END IF
    1225              :          !
    1226          412 :          IF (do_hfx) THEN
    1227              :             CALL cp_dbcsr_sm_fm_multiply(matrix_hfx(ispin)%matrix, evect(ispin), &
    1228          136 :                                          cpmos(ispin), norb, alpha=focc, beta=1.0_dp)
    1229          136 :             CALL cp_dbcsr_sm_fm_multiply(matrix_hfx(ispin)%matrix, mos, vcvec, norb, alpha=1.0_dp, beta=0.0_dp)
    1230              :             CALL cp_dbcsr_sm_fm_multiply(matrix_hfx_asymm(ispin)%matrix, evect(ispin), &
    1231          136 :                                          cpmos(ispin), norb, alpha=focc, beta=1.0_dp)
    1232          136 :             CALL cp_dbcsr_sm_fm_multiply(matrix_hfx_asymm(ispin)%matrix, mos, vcvec, norb, alpha=1.0_dp, beta=1.0_dp)
    1233              :             !
    1234          136 :             CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, mos, vcvec, 0.0_dp, cvcmat)
    1235          136 :             CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, evect(ispin), cvcmat, 0.0_dp, vcvec)
    1236              :             CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), &
    1237          136 :                                          norb, alpha=-focc, beta=1.0_dp)
    1238              :             !
    1239              :             CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=mos, matrix_g=vcvec, &
    1240          136 :                                        ncol=norb, alpha=2.0_dp, symmetry_mode=1)
    1241              :          END IF
    1242              :          !
    1243          412 :          CALL cp_fm_release(vcvec)
    1244         1578 :          CALL cp_fm_release(cvcmat)
    1245              :       END DO
    1246              : 
    1247          342 :       IF (.NOT. is_rks_triplets) THEN
    1248          304 :          CALL dbcsr_deallocate_matrix_set(matrix_hx)
    1249              :       END IF
    1250          342 :       IF (ASSOCIATED(ex_env%matrix_wx1)) CALL dbcsr_deallocate_matrix_set(ex_env%matrix_wx1)
    1251          342 :       ex_env%matrix_wx1 => matrix_wx1
    1252          342 :       IF (myfun /= xc_none) THEN
    1253          232 :          CALL dbcsr_deallocate_matrix_set(matrix_fx)
    1254          232 :          CALL dbcsr_deallocate_matrix_set(matrix_gx)
    1255              :       END IF
    1256          342 :       IF (do_hfx) THEN
    1257          124 :          CALL dbcsr_deallocate_matrix_set(matrix_hfx)
    1258          124 :          CALL dbcsr_deallocate_matrix_set(matrix_hfx_asymm)
    1259              :       END IF
    1260              : 
    1261          342 :       CALL timestop(handle)
    1262              : 
    1263          684 :    END SUBROUTINE fhxc_force
    1264              : 
    1265              : ! **************************************************************************************************
    1266              : !> \brief Simplified Tamm Dancoff approach (sTDA). Kernel contribution to forces
    1267              : !> \param qs_env ...
    1268              : !> \param ex_env ...
    1269              : !> \param gs_mos ...
    1270              : !> \param stda_env ...
    1271              : !> \param sub_env ...
    1272              : !> \param work ...
    1273              : !> \param debug_forces ...
    1274              : ! **************************************************************************************************
    1275          158 :    SUBROUTINE stda_force(qs_env, ex_env, gs_mos, stda_env, sub_env, work, debug_forces)
    1276              : 
    1277              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1278              :       TYPE(excited_energy_type), POINTER                 :: ex_env
    1279              :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
    1280              :          POINTER                                         :: gs_mos
    1281              :       TYPE(stda_env_type), POINTER                       :: stda_env
    1282              :       TYPE(tddfpt_subgroup_env_type)                     :: sub_env
    1283              :       TYPE(tddfpt_work_matrices)                         :: work
    1284              :       LOGICAL, INTENT(IN)                                :: debug_forces
    1285              : 
    1286              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'stda_force'
    1287              : 
    1288              :       INTEGER                                            :: atom_i, atom_j, ewald_type, handle, i, &
    1289              :                                                             ia, iatom, idimk, ikind, iounit, is, &
    1290              :                                                             ispin, jatom, jkind, jspin, nao, &
    1291              :                                                             natom, norb, nsgf, nspins
    1292          158 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, first_sgf, kind_of, &
    1293          158 :                                                             last_sgf
    1294              :       INTEGER, DIMENSION(2)                              :: nactive, nlim
    1295              :       LOGICAL                                            :: calculate_forces, do_coulomb, do_ewald, &
    1296              :                                                             found, is_rks_triplets, use_virial
    1297              :       REAL(KIND=dp)                                      :: alpha, bp, dgabr, dr, eta, fdim, gabr, &
    1298              :                                                             hfx, rbeta, spinfac, xfac
    1299          158 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: tcharge, tv
    1300          158 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: gtcharge
    1301              :       REAL(KIND=dp), DIMENSION(3)                        :: fij, fodeb, rij
    1302          158 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: gab, pblock
    1303          158 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1304              :       TYPE(cell_type), POINTER                           :: cell
    1305              :       TYPE(cp_fm_struct_type), POINTER                   :: fmstruct, fmstruct_mat, fmstructjspin
    1306              :       TYPE(cp_fm_type)                                   :: cvcmat, cvec, cvecjspin, t0matrix, &
    1307              :                                                             t1matrix, vcvec, xvec
    1308          158 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: xtransformed
    1309          158 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: cpmos, X
    1310              :       TYPE(cp_fm_type), POINTER                          :: ct, ctjspin, ucmatrix, uxmatrix
    1311              :       TYPE(cp_logger_type), POINTER                      :: logger
    1312              :       TYPE(dbcsr_iterator_type)                          :: iter
    1313          158 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: gamma_matrix, matrix_plo, matrix_s, &
    1314          158 :                                                             matrix_wx1, scrm
    1315              :       TYPE(dbcsr_type)                                   :: pdens, ptrans
    1316              :       TYPE(dbcsr_type), POINTER                          :: tempmat
    1317              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1318              :       TYPE(ewald_environment_type), POINTER              :: ewald_env
    1319              :       TYPE(ewald_pw_type), POINTER                       :: ewald_pw
    1320              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1321              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1322          158 :          POINTER                                         :: n_list, sab_orb
    1323          158 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1324          158 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
    1325          158 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1326              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
    1327              :       TYPE(stda_control_type), POINTER                   :: stda_control
    1328              :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
    1329              :       TYPE(virial_type), POINTER                         :: virial
    1330              : 
    1331          158 :       CALL timeset(routineN, handle)
    1332              : 
    1333          158 :       CPASSERT(ASSOCIATED(ex_env))
    1334          158 :       CPASSERT(ASSOCIATED(gs_mos))
    1335              : 
    1336          158 :       logger => cp_get_default_logger()
    1337          158 :       IF (logger%para_env%is_source()) THEN
    1338           79 :          iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
    1339              :       ELSE
    1340              :          iounit = -1
    1341              :       END IF
    1342              : 
    1343          158 :       CALL get_qs_env(qs_env, dft_control=dft_control)
    1344          158 :       tddfpt_control => dft_control%tddfpt2_control
    1345          158 :       stda_control => tddfpt_control%stda_control
    1346          158 :       nspins = dft_control%nspins
    1347          158 :       is_rks_triplets = tddfpt_control%rks_triplets .AND. (nspins == 1)
    1348              : 
    1349          158 :       X => ex_env%evect
    1350              : 
    1351          474 :       nactive(:) = stda_env%nactive(:)
    1352          158 :       xfac = 2.0_dp
    1353          158 :       spinfac = 2.0_dp
    1354          158 :       IF (nspins == 2) spinfac = 1.0_dp
    1355          158 :       NULLIFY (matrix_wx1)
    1356          158 :       CALL dbcsr_allocate_matrix_set(matrix_wx1, nspins)
    1357          158 :       NULLIFY (matrix_plo)
    1358          158 :       CALL dbcsr_allocate_matrix_set(matrix_plo, nspins)
    1359              : 
    1360          158 :       IF (nspins == 1 .AND. is_rks_triplets) THEN
    1361              :          do_coulomb = .FALSE.
    1362              :       ELSE
    1363          142 :          do_coulomb = .TRUE.
    1364              :       END IF
    1365          158 :       do_ewald = stda_control%do_ewald
    1366              : 
    1367          158 :       CALL get_qs_env(qs_env, para_env=para_env, force=force)
    1368              : 
    1369              :       CALL get_qs_env(qs_env, natom=natom, cell=cell, &
    1370          158 :                       particle_set=particle_set, qs_kind_set=qs_kind_set)
    1371          474 :       ALLOCATE (first_sgf(natom))
    1372          316 :       ALLOCATE (last_sgf(natom))
    1373          158 :       CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf, last_sgf=last_sgf)
    1374              : 
    1375          158 :       CALL get_qs_env(qs_env, ks_env=ks_env, matrix_s=matrix_s, sab_orb=sab_orb, atomic_kind_set=atomic_kind_set)
    1376          158 :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of, atom_of_kind=atom_of_kind)
    1377              : 
    1378              :       ! calculate Loewdin transformed Davidson trial vector tilde(X)=S^1/2*X
    1379              :       ! and tilde(tilde(X))=S^1/2_A*tilde(X)_A
    1380          658 :       ALLOCATE (xtransformed(nspins))
    1381          342 :       DO ispin = 1, nspins
    1382          184 :          NULLIFY (fmstruct)
    1383          184 :          ct => work%ctransformed(ispin)
    1384          184 :          CALL cp_fm_get_info(ct, matrix_struct=fmstruct)
    1385          342 :          CALL cp_fm_create(matrix=xtransformed(ispin), matrix_struct=fmstruct, name="XTRANSFORMED")
    1386              :       END DO
    1387          158 :       CALL get_lowdin_x(work%shalf, X, xtransformed)
    1388              : 
    1389          790 :       ALLOCATE (tcharge(natom), gtcharge(natom, 4))
    1390              : 
    1391          158 :       cpmos => ex_env%cpmos
    1392              : 
    1393          342 :       DO ispin = 1, nspins
    1394          184 :          ct => work%ctransformed(ispin)
    1395          184 :          CALL cp_fm_get_info(ct, matrix_struct=fmstruct, nrow_global=nsgf)
    1396          552 :          ALLOCATE (tv(nsgf))
    1397          184 :          CALL cp_fm_create(cvec, fmstruct)
    1398          184 :          CALL cp_fm_create(xvec, fmstruct)
    1399              :          !
    1400          184 :          ALLOCATE (matrix_wx1(ispin)%matrix)
    1401          184 :          CALL dbcsr_create(matrix=matrix_wx1(ispin)%matrix, template=matrix_s(1)%matrix)
    1402          184 :          CALL cp_dbcsr_alloc_block_from_nbl(matrix_wx1(ispin)%matrix, sab_orb)
    1403          184 :          CALL dbcsr_set(matrix_wx1(ispin)%matrix, 0.0_dp)
    1404          184 :          ALLOCATE (matrix_plo(ispin)%matrix)
    1405          184 :          CALL dbcsr_create(matrix=matrix_plo(ispin)%matrix, template=matrix_s(1)%matrix)
    1406          184 :          CALL cp_dbcsr_alloc_block_from_nbl(matrix_plo(ispin)%matrix, sab_orb)
    1407          184 :          CALL dbcsr_set(matrix_plo(ispin)%matrix, 0.0_dp)
    1408              :          !
    1409              :          ! *** Coulomb contribution
    1410              :          !
    1411          184 :          IF (do_coulomb) THEN
    1412              :             !
    1413          174 :             IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
    1414              :             !
    1415          870 :             tcharge(:) = 0.0_dp
    1416          388 :             DO jspin = 1, nspins
    1417          220 :                ctjspin => work%ctransformed(jspin)
    1418          220 :                CALL cp_fm_get_info(ctjspin, matrix_struct=fmstructjspin)
    1419          220 :                CALL cp_fm_get_info(ctjspin, matrix_struct=fmstructjspin, nrow_global=nsgf)
    1420          220 :                CALL cp_fm_create(cvecjspin, fmstructjspin)
    1421              :                ! CV(mu,j) = CT(mu,j)*XT(mu,j)
    1422          220 :                CALL cp_fm_schur_product(ctjspin, xtransformed(jspin), cvecjspin)
    1423              :                ! TV(mu) = SUM_j CV(mu,j)
    1424          220 :                CALL cp_fm_vectorssum(cvecjspin, tv, "R")
    1425              :                ! contract charges
    1426              :                ! TC(a) = SUM_(mu of a) TV(mu)
    1427         1078 :                DO ia = 1, natom
    1428         5692 :                   DO is = first_sgf(ia), last_sgf(ia)
    1429         5472 :                      tcharge(ia) = tcharge(ia) + tv(is)
    1430              :                   END DO
    1431              :                END DO
    1432          608 :                CALL cp_fm_release(cvecjspin)
    1433              :             END DO !jspin
    1434              :             ! Apply tcharge*gab -> gtcharge
    1435              :             ! gT(b) = SUM_a g(a,b)*TC(a)
    1436              :             ! gab = work%gamma_exchange(1)%matrix
    1437         3648 :             gtcharge = 0.0_dp
    1438              :             ! short range contribution
    1439          168 :             NULLIFY (gamma_matrix)
    1440          168 :             CALL setup_gamma(qs_env, stda_env, sub_env, gamma_matrix, ndim=4)
    1441          168 :             tempmat => gamma_matrix(1)%matrix
    1442          168 :             CALL dbcsr_iterator_start(iter, tempmat)
    1443         5323 :             DO WHILE (dbcsr_iterator_blocks_left(iter))
    1444         5155 :                CALL dbcsr_iterator_next_block(iter, iatom, jatom, gab)
    1445         5155 :                gtcharge(iatom, 1) = gtcharge(iatom, 1) + gab(1, 1)*tcharge(jatom)
    1446         5155 :                IF (iatom /= jatom) THEN
    1447         4804 :                   gtcharge(jatom, 1) = gtcharge(jatom, 1) + gab(1, 1)*tcharge(iatom)
    1448              :                END IF
    1449        20788 :                DO idimk = 2, 4
    1450        15465 :                   fdim = -1.0_dp
    1451              :                   CALL dbcsr_get_block_p(matrix=gamma_matrix(idimk)%matrix, &
    1452        15465 :                                          row=iatom, col=jatom, block=gab, found=found)
    1453        20620 :                   IF (found) THEN
    1454        15465 :                      gtcharge(iatom, idimk) = gtcharge(iatom, idimk) + gab(1, 1)*tcharge(jatom)
    1455        15465 :                      IF (iatom /= jatom) THEN
    1456        14412 :                         gtcharge(jatom, idimk) = gtcharge(jatom, idimk) + fdim*gab(1, 1)*tcharge(iatom)
    1457              :                      END IF
    1458              :                   END IF
    1459              :                END DO
    1460              :             END DO
    1461          168 :             CALL dbcsr_iterator_stop(iter)
    1462          168 :             CALL dbcsr_deallocate_matrix_set(gamma_matrix)
    1463              :             ! Ewald long range contribution
    1464          168 :             IF (do_ewald) THEN
    1465           40 :                ewald_env => work%ewald_env
    1466           40 :                ewald_pw => work%ewald_pw
    1467           40 :                CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
    1468           40 :                CALL get_qs_env(qs_env=qs_env, sab_orb=n_list, virial=virial)
    1469           40 :                use_virial = .FALSE.
    1470           40 :                calculate_forces = .FALSE.
    1471           40 :                CALL tb_ewald_overlap(gtcharge, tcharge, alpha, n_list, virial, use_virial)
    1472              :                CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, cell, &
    1473           40 :                                      gtcharge, tcharge, calculate_forces, virial, use_virial)
    1474              :                ! add self charge interaction contribution
    1475           40 :                IF (para_env%is_source()) THEN
    1476          173 :                   gtcharge(:, 1) = gtcharge(:, 1) - 2._dp*alpha*oorootpi*tcharge(:)
    1477              :                END IF
    1478              :             ELSE
    1479          128 :                nlim = get_limit(natom, para_env%num_pe, para_env%mepos)
    1480          326 :                DO iatom = nlim(1), nlim(2)
    1481          536 :                   DO jatom = 1, iatom - 1
    1482          840 :                      rij = particle_set(iatom)%r - particle_set(jatom)%r
    1483          840 :                      rij = pbc(rij, cell)
    1484          840 :                      dr = SQRT(SUM(rij(:)**2))
    1485          408 :                      IF (dr > 1.e-6_dp) THEN
    1486          210 :                         gtcharge(iatom, 1) = gtcharge(iatom, 1) + tcharge(jatom)/dr
    1487          210 :                         gtcharge(jatom, 1) = gtcharge(jatom, 1) + tcharge(iatom)/dr
    1488          840 :                         DO idimk = 2, 4
    1489          630 :                            gtcharge(iatom, idimk) = gtcharge(iatom, idimk) + rij(idimk - 1)*tcharge(jatom)/dr**3
    1490          840 :                            gtcharge(jatom, idimk) = gtcharge(jatom, idimk) - rij(idimk - 1)*tcharge(iatom)/dr**3
    1491              :                         END DO
    1492              :                      END IF
    1493              :                   END DO
    1494              :                END DO
    1495              :             END IF
    1496          168 :             CALL para_env%sum(gtcharge(:, 1))
    1497              :             ! expand charges
    1498              :             ! TV(mu) = TC(a of mu)
    1499         3898 :             tv(1:nsgf) = 0.0_dp
    1500          870 :             DO ia = 1, natom
    1501         4600 :                DO is = first_sgf(ia), last_sgf(ia)
    1502         4432 :                   tv(is) = gtcharge(ia, 1)
    1503              :                END DO
    1504              :             END DO
    1505              :             !
    1506          870 :             DO iatom = 1, natom
    1507          702 :                ikind = kind_of(iatom)
    1508          702 :                atom_i = atom_of_kind(iatom)
    1509         2808 :                DO i = 1, 3
    1510         2808 :                   fij(i) = spinfac*spinfac*gtcharge(iatom, i + 1)*tcharge(iatom)
    1511              :                END DO
    1512          702 :                force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
    1513          702 :                force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
    1514          870 :                force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
    1515              :             END DO
    1516              :             !
    1517          168 :             IF (debug_forces) THEN
    1518            8 :                fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
    1519            2 :                CALL para_env%sum(fodeb)
    1520            2 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Coul[X]   ", fodeb
    1521              :             END IF
    1522          168 :             norb = nactive(ispin)
    1523              :             ! forces from Lowdin charge derivative
    1524          168 :             CALL cp_fm_get_info(work%S_C0_C0T(ispin), matrix_struct=fmstruct)
    1525          168 :             CALL cp_fm_create(t0matrix, matrix_struct=fmstruct, name="T0 SCRATCH")
    1526          168 :             CALL cp_fm_create(t1matrix, matrix_struct=fmstruct, name="T1 SCRATCH")
    1527          168 :             ALLOCATE (ucmatrix)
    1528          168 :             CALL fm_pool_create_fm(work%fm_pool_ao_mo_occ(ispin)%pool, ucmatrix)
    1529          168 :             ALLOCATE (uxmatrix)
    1530          168 :             CALL fm_pool_create_fm(work%fm_pool_ao_mo_occ(ispin)%pool, uxmatrix)
    1531          168 :             ct => work%ctransformed(ispin)
    1532          168 :             CALL cp_fm_to_fm(ct, cvec)
    1533          168 :             CALL cp_fm_row_scale(cvec, tv)
    1534              :             CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
    1535          168 :                                cvec, 0.0_dp, ucmatrix)
    1536              :             CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
    1537          168 :                                X(ispin), 0.0_dp, uxmatrix)
    1538          168 :             CALL parallel_gemm('N', 'T', nsgf, nsgf, norb, 1.0_dp, uxmatrix, ucmatrix, 0.0_dp, t0matrix)
    1539          168 :             CALL cp_fm_to_fm(xtransformed(ispin), cvec)
    1540          168 :             CALL cp_fm_row_scale(cvec, tv)
    1541              :             CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
    1542          168 :                                cvec, 0.0_dp, uxmatrix)
    1543              :             CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
    1544          168 :                                gs_mos(ispin)%mos_occ, 0.0_dp, ucmatrix)
    1545          168 :             CALL parallel_gemm('N', 'T', nsgf, nsgf, norb, 1.0_dp, ucmatrix, uxmatrix, 1.0_dp, t0matrix)
    1546          168 :             CALL cp_fm_schur_product(work%slambda, t0matrix, t1matrix)
    1547              :             !
    1548              :             CALL parallel_gemm('N', 'N', nsgf, nsgf, nsgf, spinfac, work%S_eigenvectors, t1matrix, &
    1549          168 :                                0.0_dp, t0matrix)
    1550              :             CALL cp_dbcsr_plus_fm_fm_t(matrix_plo(ispin)%matrix, matrix_v=t0matrix, &
    1551          168 :                                        matrix_g=work%S_eigenvectors, ncol=nsgf, alpha=2.0_dp, symmetry_mode=1)
    1552          168 :             CALL fm_pool_give_back_fm(work%fm_pool_ao_mo_occ(ispin)%pool, ucmatrix)
    1553          168 :             DEALLOCATE (ucmatrix)
    1554          168 :             CALL fm_pool_give_back_fm(work%fm_pool_ao_mo_occ(ispin)%pool, uxmatrix)
    1555          168 :             DEALLOCATE (uxmatrix)
    1556          168 :             CALL cp_fm_release(t0matrix)
    1557          168 :             CALL cp_fm_release(t1matrix)
    1558              :             !
    1559              :             ! CV(mu,i) = TV(mu)*XT(mu,i)
    1560          168 :             CALL cp_fm_to_fm(xtransformed(ispin), cvec)
    1561          168 :             CALL cp_fm_row_scale(cvec, tv)
    1562          168 :             CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, cpmos(ispin), norb, 2.0_dp*spinfac, 1.0_dp)
    1563              :             ! CV(mu,i) = TV(mu)*CT(mu,i)
    1564          168 :             ct => work%ctransformed(ispin)
    1565          168 :             CALL cp_fm_to_fm(ct, cvec)
    1566          168 :             CALL cp_fm_row_scale(cvec, tv)
    1567              :             ! Shalf(nu,mu)*CV(mu,i)
    1568          168 :             CALL cp_fm_get_info(cvec, matrix_struct=fmstruct, nrow_global=nao)
    1569          168 :             CALL cp_fm_create(vcvec, fmstruct)
    1570          168 :             CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, vcvec, norb, 1.0_dp, 0.0_dp)
    1571              :             CALL cp_fm_struct_create(fmstruct_mat, context=fmstruct%context, nrow_global=norb, &
    1572          168 :                                      ncol_global=norb, para_env=fmstruct%para_env)
    1573          168 :             CALL cp_fm_create(cvcmat, fmstruct_mat)
    1574          168 :             CALL cp_fm_struct_release(fmstruct_mat)
    1575          168 :             CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, gs_mos(ispin)%mos_occ, vcvec, 0.0_dp, cvcmat)
    1576          168 :             CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, X(ispin), cvcmat, 0.0_dp, vcvec)
    1577              :             CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), &
    1578          168 :                                          nactive(ispin), alpha=-2.0_dp*spinfac, beta=1.0_dp)
    1579              :             ! wx1
    1580          168 :             alpha = 2.0_dp
    1581              :             CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=gs_mos(ispin)%mos_occ, &
    1582          168 :                                        matrix_g=vcvec, ncol=norb, alpha=2.0_dp*alpha, symmetry_mode=1)
    1583          168 :             CALL cp_fm_release(vcvec)
    1584          168 :             CALL cp_fm_release(cvcmat)
    1585              :          END IF
    1586              :          !
    1587              :          ! *** Exchange contribution
    1588              :          !
    1589          184 :          IF (stda_env%do_exchange) THEN
    1590              :             !
    1591          166 :             IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
    1592              :             !
    1593          160 :             norb = nactive(ispin)
    1594              :             !
    1595          160 :             tempmat => work%shalf
    1596          160 :             CALL dbcsr_create(pdens, template=tempmat, matrix_type=dbcsr_type_no_symmetry)
    1597              :             ! P(nu,mu) = SUM_j XT(nu,j)*CT(mu,j)
    1598          160 :             ct => work%ctransformed(ispin)
    1599          160 :             CALL dbcsr_set(pdens, 0.0_dp)
    1600              :             CALL cp_dbcsr_plus_fm_fm_t(pdens, xtransformed(ispin), ct, nactive(ispin), &
    1601          160 :                                        1.0_dp, keep_sparsity=.FALSE.)
    1602          160 :             CALL dbcsr_filter(pdens, stda_env%eps_td_filter)
    1603              :             ! Apply PP*gab -> PP; gab = gamma_coulomb
    1604              :             ! P(nu,mu) = P(nu,mu)*g(a of nu,b of mu)
    1605          160 :             bp = stda_env%beta_param
    1606          160 :             hfx = stda_env%hfx_fraction
    1607          160 :             CALL dbcsr_iterator_start(iter, pdens)
    1608        10062 :             DO WHILE (dbcsr_iterator_blocks_left(iter))
    1609         9902 :                CALL dbcsr_iterator_next_block(iter, iatom, jatom, pblock)
    1610        39608 :                rij = particle_set(iatom)%r - particle_set(jatom)%r
    1611        39608 :                rij = pbc(rij, cell)
    1612        39608 :                dr = SQRT(SUM(rij(:)**2))
    1613         9902 :                ikind = kind_of(iatom)
    1614         9902 :                jkind = kind_of(jatom)
    1615              :                eta = (stda_env%kind_param_set(ikind)%kind_param%hardness_param + &
    1616         9902 :                       stda_env%kind_param_set(jkind)%kind_param%hardness_param)/2.0_dp
    1617         9902 :                rbeta = dr**bp
    1618         9902 :                IF (hfx > 0.0_dp) THEN
    1619         9843 :                   gabr = (1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp)
    1620              :                ELSE
    1621           59 :                   IF (dr < 1.0e-6_dp) THEN
    1622              :                      gabr = 0.0_dp
    1623              :                   ELSE
    1624           42 :                      gabr = 1._dp/dr
    1625              :                   END IF
    1626              :                END IF
    1627              :                !      gabr = (1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp)
    1628              :                ! forces
    1629         9885 :                IF (dr > 1.0e-6_dp) THEN
    1630         9583 :                   IF (hfx > 0.0_dp) THEN
    1631         9541 :                      dgabr = -(1._dp/bp)*(1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp + 1._dp)
    1632         9541 :                      dgabr = bp*rbeta/dr**2*dgabr
    1633       112529 :                      dgabr = SUM(pblock**2)*dgabr
    1634              :                   ELSE
    1635           42 :                      dgabr = -1.0_dp/dr**3
    1636         3142 :                      dgabr = SUM(pblock**2)*dgabr
    1637              :                   END IF
    1638         9583 :                   atom_i = atom_of_kind(iatom)
    1639         9583 :                   atom_j = atom_of_kind(jatom)
    1640        38332 :                   DO i = 1, 3
    1641        38332 :                      fij(i) = dgabr*rij(i)
    1642              :                   END DO
    1643         9583 :                   force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
    1644         9583 :                   force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
    1645         9583 :                   force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
    1646         9583 :                   force(jkind)%rho_elec(1, atom_j) = force(jkind)%rho_elec(1, atom_j) + fij(1)
    1647         9583 :                   force(jkind)%rho_elec(2, atom_j) = force(jkind)%rho_elec(2, atom_j) + fij(2)
    1648         9583 :                   force(jkind)%rho_elec(3, atom_j) = force(jkind)%rho_elec(3, atom_j) + fij(3)
    1649              :                END IF
    1650              :                !
    1651       132828 :                pblock = gabr*pblock
    1652              :             END DO
    1653          160 :             CALL dbcsr_iterator_stop(iter)
    1654              :             !
    1655              :             ! Transpose pdens matrix
    1656          160 :             CALL dbcsr_create(ptrans, template=pdens)
    1657          160 :             CALL dbcsr_transposed(ptrans, pdens)
    1658              :             !
    1659              :             ! forces from Lowdin charge derivative
    1660          160 :             CALL cp_fm_get_info(work%S_C0_C0T(ispin), matrix_struct=fmstruct)
    1661          160 :             CALL cp_fm_create(t0matrix, matrix_struct=fmstruct, name="T0 SCRATCH")
    1662          160 :             CALL cp_fm_create(t1matrix, matrix_struct=fmstruct, name="T1 SCRATCH")
    1663          160 :             ALLOCATE (ucmatrix)
    1664          160 :             CALL fm_pool_create_fm(work%fm_pool_ao_mo_occ(ispin)%pool, ucmatrix)
    1665          160 :             ALLOCATE (uxmatrix)
    1666          160 :             CALL fm_pool_create_fm(work%fm_pool_ao_mo_occ(ispin)%pool, uxmatrix)
    1667          160 :             ct => work%ctransformed(ispin)
    1668          160 :             CALL cp_dbcsr_sm_fm_multiply(pdens, ct, cvec, norb, 1.0_dp, 0.0_dp)
    1669              :             CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
    1670          160 :                                cvec, 0.0_dp, ucmatrix)
    1671              :             CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
    1672          160 :                                X(ispin), 0.0_dp, uxmatrix)
    1673          160 :             CALL parallel_gemm('N', 'T', nsgf, nsgf, norb, 1.0_dp, uxmatrix, ucmatrix, 0.0_dp, t0matrix)
    1674          160 :             CALL cp_dbcsr_sm_fm_multiply(ptrans, xtransformed(ispin), cvec, norb, 1.0_dp, 0.0_dp)
    1675              :             CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
    1676          160 :                                cvec, 0.0_dp, uxmatrix)
    1677              :             CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
    1678          160 :                                gs_mos(ispin)%mos_occ, 0.0_dp, ucmatrix)
    1679          160 :             CALL parallel_gemm('N', 'T', nsgf, nsgf, norb, 1.0_dp, ucmatrix, uxmatrix, 1.0_dp, t0matrix)
    1680          160 :             CALL cp_fm_schur_product(work%slambda, t0matrix, t1matrix)
    1681              :             CALL parallel_gemm('N', 'N', nsgf, nsgf, nsgf, -1.0_dp, work%S_eigenvectors, t1matrix, &
    1682          160 :                                0.0_dp, t0matrix)
    1683              :             CALL cp_dbcsr_plus_fm_fm_t(matrix_plo(ispin)%matrix, matrix_v=t0matrix, &
    1684          160 :                                        matrix_g=work%S_eigenvectors, ncol=nsgf, alpha=2.0_dp, symmetry_mode=1)
    1685          160 :             CALL fm_pool_give_back_fm(work%fm_pool_ao_mo_occ(ispin)%pool, ucmatrix)
    1686          160 :             DEALLOCATE (ucmatrix)
    1687          160 :             CALL fm_pool_give_back_fm(work%fm_pool_ao_mo_occ(ispin)%pool, uxmatrix)
    1688          160 :             DEALLOCATE (uxmatrix)
    1689          160 :             CALL cp_fm_release(t0matrix)
    1690          160 :             CALL cp_fm_release(t1matrix)
    1691              : 
    1692              :             ! RHS contribution to response matrix
    1693              :             ! CV(nu,i) = P(nu,mu)*XT(mu,i)
    1694          160 :             CALL cp_dbcsr_sm_fm_multiply(ptrans, xtransformed(ispin), cvec, norb, 1.0_dp, 0.0_dp)
    1695              :             CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, cpmos(ispin), norb, &
    1696          160 :                                          alpha=-xfac, beta=1.0_dp)
    1697              :             !
    1698          160 :             CALL cp_fm_get_info(cvec, matrix_struct=fmstruct, nrow_global=nao)
    1699          160 :             CALL cp_fm_create(vcvec, fmstruct)
    1700              :             ! CV(nu,i) = P(nu,mu)*CT(mu,i)
    1701          160 :             CALL cp_dbcsr_sm_fm_multiply(ptrans, ct, cvec, norb, 1.0_dp, 0.0_dp)
    1702          160 :             CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, vcvec, norb, 1.0_dp, 0.0_dp)
    1703              :             CALL cp_fm_struct_create(fmstruct_mat, context=fmstruct%context, nrow_global=norb, &
    1704          160 :                                      ncol_global=norb, para_env=fmstruct%para_env)
    1705          160 :             CALL cp_fm_create(cvcmat, fmstruct_mat)
    1706          160 :             CALL cp_fm_struct_release(fmstruct_mat)
    1707          160 :             CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, gs_mos(ispin)%mos_occ, vcvec, 0.0_dp, cvcmat)
    1708          160 :             CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, X(ispin), cvcmat, 0.0_dp, vcvec)
    1709              :             CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), &
    1710          160 :                                          norb, alpha=xfac, beta=1.0_dp)
    1711              :             ! wx1
    1712          160 :             IF (nspins == 2) THEN
    1713           44 :                alpha = -2.0_dp
    1714              :             ELSE
    1715          116 :                alpha = -1.0_dp
    1716              :             END IF
    1717              :             CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=gs_mos(ispin)%mos_occ, &
    1718              :                                        matrix_g=vcvec, &
    1719          160 :                                        ncol=norb, alpha=2.0_dp*alpha, symmetry_mode=1)
    1720          160 :             CALL cp_fm_release(vcvec)
    1721          160 :             CALL cp_fm_release(cvcmat)
    1722              :             !
    1723          160 :             CALL dbcsr_release(pdens)
    1724          160 :             CALL dbcsr_release(ptrans)
    1725              :             !
    1726          160 :             IF (debug_forces) THEN
    1727            8 :                fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
    1728            2 :                CALL para_env%sum(fodeb)
    1729            2 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Exch[X]   ", fodeb
    1730              :             END IF
    1731              :          END IF
    1732              :          !
    1733          184 :          CALL cp_fm_release(cvec)
    1734          184 :          CALL cp_fm_release(xvec)
    1735          710 :          DEALLOCATE (tv)
    1736              :       END DO
    1737              : 
    1738          158 :       CALL cp_fm_release(xtransformed)
    1739          158 :       DEALLOCATE (tcharge, gtcharge)
    1740          158 :       DEALLOCATE (first_sgf, last_sgf)
    1741              : 
    1742              :       ! Lowdin forces
    1743          158 :       IF (nspins == 2) THEN
    1744              :          CALL dbcsr_add(matrix_plo(1)%matrix, matrix_plo(2)%matrix, &
    1745           26 :                         alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
    1746              :       END IF
    1747          158 :       CALL dbcsr_scale(matrix_plo(1)%matrix, -1.0_dp)
    1748          158 :       NULLIFY (scrm)
    1749          164 :       IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
    1750              :       CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
    1751              :                                 matrix_name="OVERLAP MATRIX", &
    1752              :                                 basis_type_a="ORB", basis_type_b="ORB", &
    1753              :                                 sab_nl=sab_orb, calculate_forces=.TRUE., &
    1754          158 :                                 matrix_p=matrix_plo(1)%matrix)
    1755          158 :       CALL dbcsr_deallocate_matrix_set(scrm)
    1756          158 :       CALL dbcsr_deallocate_matrix_set(matrix_plo)
    1757          158 :       IF (debug_forces) THEN
    1758            8 :          fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
    1759            2 :          CALL para_env%sum(fodeb)
    1760            2 :          IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Lowdin ", fodeb
    1761              :       END IF
    1762              : 
    1763          158 :       IF (ASSOCIATED(ex_env%matrix_wx1)) CALL dbcsr_deallocate_matrix_set(ex_env%matrix_wx1)
    1764          158 :       ex_env%matrix_wx1 => matrix_wx1
    1765              : 
    1766          158 :       CALL timestop(handle)
    1767              : 
    1768          316 :    END SUBROUTINE stda_force
    1769              : 
    1770              : ! **************************************************************************************************
    1771              : 
    1772              : END MODULE qs_tddfpt2_fhxc_forces
        

Generated by: LCOV version 2.0-1