LCOV - code coverage report
Current view: top level - src - qs_linres_kernel.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:9e7f125) Lines: 362 372 97.3 %
Date: 2025-05-16 07:28:05 Functions: 6 6 100.0 %

          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             : ! **************************************************************************************************
       9             : !> \brief linres kernel functions
      10             : !> \par History
      11             : !>      created from qs_linres_methods
      12             : !> \author JGH
      13             : ! **************************************************************************************************
      14             : MODULE qs_linres_kernel
      15             :    USE admm_types,                      ONLY: admm_type,&
      16             :                                               get_admm_env
      17             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      18             :                                               get_atomic_kind
      19             :    USE cp_control_types,                ONLY: dft_control_type
      20             :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      21             :                                               dbcsr_copy,&
      22             :                                               dbcsr_create,&
      23             :                                               dbcsr_deallocate_matrix,&
      24             :                                               dbcsr_p_type,&
      25             :                                               dbcsr_set
      26             :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply,&
      27             :                                               dbcsr_allocate_matrix_set,&
      28             :                                               dbcsr_deallocate_matrix_set
      29             :    USE cp_fm_types,                     ONLY: cp_fm_get_info,&
      30             :                                               cp_fm_type
      31             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      32             :                                               cp_logger_type,&
      33             :                                               cp_to_string
      34             :    USE hartree_local_methods,           ONLY: Vh_1c_gg_integrals
      35             :    USE hfx_energy_potential,            ONLY: integrate_four_center
      36             :    USE hfx_ri,                          ONLY: hfx_ri_update_ks
      37             :    USE hfx_types,                       ONLY: hfx_type
      38             :    USE input_constants,                 ONLY: do_admm_aux_exch_func_none,&
      39             :                                               do_admm_basis_projection,&
      40             :                                               do_admm_exch_scaling_none,&
      41             :                                               do_admm_purify_none,&
      42             :                                               kg_tnadd_embed
      43             :    USE input_section_types,             ONLY: section_get_ival,&
      44             :                                               section_get_lval,&
      45             :                                               section_get_rval,&
      46             :                                               section_vals_get,&
      47             :                                               section_vals_get_subs_vals,&
      48             :                                               section_vals_type,&
      49             :                                               section_vals_val_get
      50             :    USE kg_correction,                   ONLY: kg_ekin_subset
      51             :    USE kg_environment_types,            ONLY: kg_environment_type
      52             :    USE kinds,                           ONLY: default_string_length,&
      53             :                                               dp
      54             :    USE lri_environment_types,           ONLY: lri_density_type,&
      55             :                                               lri_environment_type,&
      56             :                                               lri_kind_type
      57             :    USE lri_ks_methods,                  ONLY: calculate_lri_ks_matrix
      58             :    USE message_passing,                 ONLY: mp_para_env_type
      59             :    USE mulliken,                        ONLY: ao_charges
      60             :    USE particle_types,                  ONLY: particle_type
      61             :    USE pw_env_types,                    ONLY: pw_env_get,&
      62             :                                               pw_env_type
      63             :    USE pw_methods,                      ONLY: pw_axpy,&
      64             :                                               pw_copy,&
      65             :                                               pw_scale,&
      66             :                                               pw_transfer
      67             :    USE pw_poisson_methods,              ONLY: pw_poisson_solve
      68             :    USE pw_poisson_types,                ONLY: pw_poisson_type
      69             :    USE pw_pool_types,                   ONLY: pw_pool_type
      70             :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      71             :                                               pw_r3d_rs_type
      72             :    USE qs_environment_types,            ONLY: get_qs_env,&
      73             :                                               qs_environment_type
      74             :    USE qs_fxc,                          ONLY: qs_fxc_analytic,&
      75             :                                               qs_fxc_fdiff
      76             :    USE qs_gapw_densities,               ONLY: prepare_gapw_den
      77             :    USE qs_integrate_potential,          ONLY: integrate_v_rspace,&
      78             :                                               integrate_v_rspace_diagonal,&
      79             :                                               integrate_v_rspace_one_center
      80             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      81             :                                               get_qs_kind_set,&
      82             :                                               qs_kind_type
      83             :    USE qs_kpp1_env_types,               ONLY: qs_kpp1_env_type
      84             :    USE qs_ks_atom,                      ONLY: update_ks_atom
      85             :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      86             :    USE qs_linres_types,                 ONLY: linres_control_type
      87             :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
      88             :    USE qs_p_env_methods,                ONLY: p_env_finish_kpp1
      89             :    USE qs_p_env_types,                  ONLY: qs_p_env_type
      90             :    USE qs_rho0_ggrid,                   ONLY: integrate_vhg0_rspace
      91             :    USE qs_rho_atom_types,               ONLY: rho_atom_type
      92             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      93             :                                               qs_rho_type
      94             :    USE qs_vxc_atom,                     ONLY: calculate_xc_2nd_deriv_atom
      95             :    USE task_list_types,                 ONLY: task_list_type
      96             :    USE xc,                              ONLY: xc_calc_2nd_deriv,&
      97             :                                               xc_prep_2nd_deriv
      98             :    USE xc_derivatives,                  ONLY: xc_functionals_get_needs
      99             :    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
     100             :    USE xc_rho_set_types,                ONLY: xc_rho_set_create,&
     101             :                                               xc_rho_set_release,&
     102             :                                               xc_rho_set_type,&
     103             :                                               xc_rho_set_update
     104             :    USE xtb_ehess,                       ONLY: xtb_coulomb_hessian
     105             :    USE xtb_types,                       ONLY: get_xtb_atom_param,&
     106             :                                               xtb_atom_type
     107             : #include "./base/base_uses.f90"
     108             : 
     109             :    IMPLICIT NONE
     110             : 
     111             :    PRIVATE
     112             : 
     113             :    ! *** Public subroutines ***
     114             :    PUBLIC :: apply_xc_admm
     115             :    PUBLIC :: apply_hfx
     116             :    PUBLIC :: apply_op_2
     117             :    PUBLIC :: hfx_matrix
     118             : 
     119             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_kernel'
     120             : 
     121             : ! **************************************************************************************************
     122             : 
     123             : CONTAINS
     124             : 
     125             : ! **************************************************************************************************
     126             : !> \brief ...
     127             : !> \param qs_env ...
     128             : !> \param p_env ...
     129             : !> \param c0 ...
     130             : !> \param Av ...
     131             : ! **************************************************************************************************
     132        7730 :    SUBROUTINE apply_op_2(qs_env, p_env, c0, Av)
     133             :       !
     134             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     135             :       TYPE(qs_p_env_type)                                :: p_env
     136             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: c0
     137             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT)      :: Av
     138             : 
     139             :       INTEGER                                            :: ispin, ncol
     140             :       TYPE(dft_control_type), POINTER                    :: dft_control
     141             : 
     142        7730 :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
     143        7730 :       IF (dft_control%qs_control%semi_empirical) THEN
     144           0 :          CPABORT("Linear response not available with SE methods")
     145        7730 :       ELSEIF (dft_control%qs_control%dftb) THEN
     146           0 :          CPABORT("Linear response not available with DFTB")
     147        7730 :       ELSEIF (dft_control%qs_control%xtb) THEN
     148         110 :          CALL apply_op_2_xtb(qs_env, p_env)
     149             :       ELSE
     150        7620 :          CALL apply_op_2_dft(qs_env, p_env)
     151        7620 :          CALL apply_hfx(qs_env, p_env)
     152        7620 :          CALL apply_xc_admm(qs_env, p_env)
     153        7620 :          IF (dft_control%do_admm) CALL p_env_finish_kpp1(qs_env, p_env)
     154             :       END IF
     155             : 
     156       16316 :       DO ispin = 1, SIZE(c0)
     157        8586 :          CALL cp_fm_get_info(c0(ispin), ncol_global=ncol)
     158             :          CALL cp_dbcsr_sm_fm_multiply(p_env%kpp1(ispin)%matrix, &
     159             :                                       c0(ispin), &
     160             :                                       Av(ispin), &
     161       16316 :                                       ncol=ncol, alpha=1.0_dp, beta=1.0_dp)
     162             :       END DO
     163             : 
     164        7730 :    END SUBROUTINE apply_op_2
     165             : 
     166             : ! **************************************************************************************************
     167             : !> \brief ...
     168             : !> \param qs_env ...
     169             : !> \param p_env ...
     170             : ! **************************************************************************************************
     171        7620 :    SUBROUTINE apply_op_2_dft(qs_env, p_env)
     172             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     173             :       TYPE(qs_p_env_type)                                :: p_env
     174             : 
     175             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'apply_op_2_dft'
     176             : 
     177             :       INTEGER                                            :: handle, ikind, ispin, nkind, ns, nspins
     178             :       LOGICAL                                            :: deriv2_analytic, gapw, gapw_xc, &
     179             :                                                             lr_triplet, lrigpw
     180             :       REAL(KIND=dp)                                      :: alpha, ekin_mol, energy_hartree, &
     181             :                                                             energy_hartree_1c
     182             :       TYPE(admm_type), POINTER                           :: admm_env
     183        7620 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     184             :       TYPE(cp_logger_type), POINTER                      :: logger
     185        7620 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: k1mat, matrix_s, rho1_ao, rho_ao
     186        7620 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: ksmat, psmat
     187             :       TYPE(dft_control_type), POINTER                    :: dft_control
     188             :       TYPE(kg_environment_type), POINTER                 :: kg_env
     189             :       TYPE(linres_control_type), POINTER                 :: linres_control
     190             :       TYPE(lri_density_type), POINTER                    :: lri_density
     191             :       TYPE(lri_environment_type), POINTER                :: lri_env
     192        7620 :       TYPE(lri_kind_type), DIMENSION(:), POINTER         :: lri_v_int
     193             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     194             :       TYPE(pw_c1d_gs_type)                               :: rho1_tot_gspace, v_hartree_gspace
     195        7620 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho1_g
     196             :       TYPE(pw_env_type), POINTER                         :: pw_env
     197             :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
     198             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     199             :       TYPE(pw_r3d_rs_type)                               :: v_hartree_rspace
     200        7620 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho1_r, rho_r, tau1_r, v_rspace_new, &
     201        7620 :                                                             v_xc, v_xc_tau
     202             :       TYPE(qs_kpp1_env_type), POINTER                    :: kpp1_env
     203             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     204             :       TYPE(qs_rho_type), POINTER                         :: rho, rho0, rho1, rho1_xc, rho1a, &
     205             :                                                             rho_aux, rho_xc
     206        7620 :       TYPE(rho_atom_type), DIMENSION(:), POINTER         :: rho1_atom_set, rho_atom_set
     207             :       TYPE(section_vals_type), POINTER                   :: input, xc_section, xc_section_aux
     208             : 
     209        7620 :       CALL timeset(routineN, handle)
     210             : 
     211        7620 :       NULLIFY (auxbas_pw_pool, pw_env, v_rspace_new, para_env, rho1_r, &
     212        7620 :                v_xc, rho1_ao, rho_ao, poisson_env, input, rho, dft_control, &
     213        7620 :                logger, rho1_g, v_xc_tau)
     214        7620 :       logger => cp_get_default_logger()
     215             : 
     216        7620 :       energy_hartree = 0.0_dp
     217        7620 :       energy_hartree_1c = 0.0_dp
     218             : 
     219        7620 :       CPASSERT(ASSOCIATED(p_env%kpp1))
     220        7620 :       CPASSERT(ASSOCIATED(p_env%kpp1_env))
     221        7620 :       kpp1_env => p_env%kpp1_env
     222             : 
     223             :       CALL get_qs_env(qs_env=qs_env, &
     224             :                       ks_env=ks_env, &
     225             :                       pw_env=pw_env, &
     226             :                       input=input, &
     227             :                       admm_env=admm_env, &
     228             :                       para_env=para_env, &
     229             :                       rho=rho, &
     230             :                       rho_xc=rho_xc, &
     231             :                       linres_control=linres_control, &
     232        7620 :                       dft_control=dft_control)
     233             : 
     234        7620 :       gapw = dft_control%qs_control%gapw
     235        7620 :       gapw_xc = dft_control%qs_control%gapw_xc
     236        7620 :       lr_triplet = linres_control%lr_triplet
     237             : 
     238        7620 :       rho1 => p_env%rho1
     239        7620 :       rho1_xc => p_env%rho1_xc
     240        7620 :       CPASSERT(ASSOCIATED(rho1))
     241        7620 :       IF (gapw_xc) THEN
     242         352 :          CPASSERT(ASSOCIATED(rho1_xc))
     243             :       END IF
     244             : 
     245        7620 :       CALL qs_rho_get(rho, rho_ao=rho_ao, rho_r=rho_r)
     246        7620 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     247             : 
     248        7620 :       nspins = SIZE(p_env%kpp1)
     249        7620 :       lrigpw = dft_control%qs_control%lrigpw
     250        7620 :       IF (lrigpw) THEN
     251             :          CALL get_qs_env(qs_env, &
     252             :                          lri_env=lri_env, &
     253             :                          lri_density=lri_density, &
     254          72 :                          atomic_kind_set=atomic_kind_set)
     255             :       END IF
     256             : 
     257        7620 :       IF (.NOT. ASSOCIATED(kpp1_env%v_ao)) THEN
     258        1004 :          CALL get_qs_env(qs_env, matrix_s=matrix_s)
     259        1004 :          CALL dbcsr_allocate_matrix_set(kpp1_env%v_ao, nspins)
     260        2126 :          DO ispin = 1, nspins
     261        1122 :             ALLOCATE (kpp1_env%v_ao(ispin)%matrix)
     262             :             CALL dbcsr_copy(kpp1_env%v_ao(ispin)%matrix, matrix_s(1)%matrix, &
     263        2126 :                             name="kpp1%v_ao-"//ADJUSTL(cp_to_string(ispin)))
     264             :          END DO
     265             :       END IF
     266             : 
     267        7620 :       IF (dft_control%do_admm) THEN
     268        1728 :          xc_section => admm_env%xc_section_primary
     269             :       ELSE
     270        5892 :          xc_section => section_vals_get_subs_vals(input, "DFT%XC")
     271             :       END IF
     272             : 
     273             :       ! gets the tmp grids
     274        7620 :       CPASSERT(ASSOCIATED(pw_env))
     275             :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
     276        7620 :                       poisson_env=poisson_env)
     277        7620 :       CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
     278        7620 :       CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
     279             : 
     280        7620 :       IF (gapw .OR. gapw_xc) &
     281        1144 :          CALL prepare_gapw_den(qs_env, p_env%local_rho_set, do_rho0=(.NOT. gapw_xc))
     282             : 
     283             :       ! *** calculate the hartree potential on the total density ***
     284        7620 :       CALL auxbas_pw_pool%create_pw(rho1_tot_gspace)
     285             : 
     286        7620 :       CALL qs_rho_get(rho1, rho_g=rho1_g)
     287        7620 :       CALL pw_copy(rho1_g(1), rho1_tot_gspace)
     288        8476 :       DO ispin = 2, nspins
     289        8476 :          CALL pw_axpy(rho1_g(ispin), rho1_tot_gspace)
     290             :       END DO
     291        7620 :       IF (gapw) &
     292         792 :          CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rho0_s_gs, rho1_tot_gspace)
     293             : 
     294        7620 :       IF (.NOT. (nspins == 1 .AND. lr_triplet)) THEN
     295             :          CALL pw_poisson_solve(poisson_env, rho1_tot_gspace, &
     296             :                                energy_hartree, &
     297        7620 :                                v_hartree_gspace)
     298        7620 :          CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
     299             :       END IF
     300             : 
     301        7620 :       CALL auxbas_pw_pool%give_back_pw(rho1_tot_gspace)
     302             : 
     303             :       ! *** calculate the xc potential ***
     304        7620 :       NULLIFY (rho1a)
     305        7620 :       IF (gapw_xc) THEN
     306         352 :          rho0 => rho_xc
     307         352 :          rho1a => rho1_xc
     308             :       ELSE
     309        7268 :          rho0 => rho
     310        7268 :          rho1a => rho1
     311             :       END IF
     312             : 
     313        7620 :       deriv2_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
     314        7620 :       NULLIFY (v_xc_tau)
     315        7620 :       IF (deriv2_analytic) THEN
     316        7560 :          CALL qs_rho_get(rho1a, rho_r=rho1_r, tau_r=tau1_r)
     317        7560 :          CALL qs_fxc_analytic(rho0, rho1_r, tau1_r, xc_section, auxbas_pw_pool, lr_triplet, v_xc, v_xc_tau)
     318        7560 :          IF (gapw .OR. gapw_xc) THEN
     319        1144 :             CALL get_qs_env(qs_env, rho_atom_set=rho_atom_set)
     320        1144 :             rho1_atom_set => p_env%local_rho_set%rho_atom_set
     321             :             CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
     322        1144 :                                              do_triplet=lr_triplet)
     323             :          END IF
     324             :       ELSE
     325          60 :          CALL qs_fxc_fdiff(ks_env, rho0, rho1a, xc_section, 6, lr_triplet, v_xc, v_xc_tau)
     326          60 :          CPASSERT((.NOT. gapw) .AND. (.NOT. gapw_xc))
     327             :       END IF
     328             : 
     329        7620 :       v_rspace_new => v_xc
     330        7620 :       NULLIFY (v_xc)
     331             : 
     332        7620 :       CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
     333       16096 :       DO ispin = 1, nspins
     334        8476 :          CALL pw_scale(v_rspace_new(ispin), v_rspace_new(ispin)%pw_grid%dvol)
     335       16096 :          IF (ASSOCIATED(v_xc_tau)) CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
     336             :       END DO
     337             : 
     338             :       ! ADMM Correction
     339        7620 :       IF (dft_control%do_admm) THEN
     340        1728 :          IF (admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
     341         998 :             IF (.NOT. ASSOCIATED(kpp1_env%deriv_set_admm)) THEN
     342         154 :                CPASSERT(.NOT. lr_triplet)
     343         154 :                xc_section_aux => admm_env%xc_section_aux
     344         154 :                CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux)
     345         154 :                CALL qs_rho_get(rho_aux, rho_r=rho_r)
     346        3542 :                ALLOCATE (kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm)
     347             :                CALL xc_prep_2nd_deriv(kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm, &
     348             :                                       rho_r, auxbas_pw_pool, &
     349         154 :                                       xc_section=xc_section_aux)
     350             :             END IF
     351             :          END IF
     352             :       END IF
     353             : 
     354             :       !-------------------------------!
     355             :       ! Add both hartree and xc terms !
     356             :       !-------------------------------!
     357       16096 :       DO ispin = 1, nspins
     358        8476 :          CALL dbcsr_set(kpp1_env%v_ao(ispin)%matrix, 0.0_dp)
     359             : 
     360        8476 :          IF (gapw_xc) THEN
     361             :             ! XC and Hartree are integrated separatedly
     362             :             ! XC uses the soft basis set only
     363             : 
     364         368 :             IF (nspins == 1) THEN
     365             : 
     366         336 :                IF (.NOT. (lr_triplet)) THEN
     367         336 :                   CALL pw_scale(v_rspace_new(1), 2.0_dp)
     368         336 :                   IF (ASSOCIATED(v_xc_tau)) CALL pw_scale(v_xc_tau(1), 2.0_dp)
     369             :                END IF
     370         336 :                CALL qs_rho_get(rho1, rho_ao=rho1_ao)
     371             :                ! remove kpp1_env%v_ao and work directly on k_p_p1 ?
     372             :                CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
     373             :                                        pmat=rho1_ao(ispin), &
     374             :                                        hmat=kpp1_env%v_ao(ispin), &
     375             :                                        qs_env=qs_env, &
     376         336 :                                        calculate_forces=.FALSE., gapw=gapw_xc)
     377             : 
     378         336 :                IF (ASSOCIATED(v_xc_tau)) THEN
     379             :                   CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
     380             :                                           pmat=rho1_ao(ispin), &
     381             :                                           hmat=kpp1_env%v_ao(ispin), &
     382             :                                           qs_env=qs_env, &
     383             :                                           compute_tau=.TRUE., &
     384           0 :                                           calculate_forces=.FALSE., gapw=gapw_xc)
     385             :                END IF
     386             : 
     387             :                ! add hartree only for SINGLETS
     388         336 :                IF (.NOT. lr_triplet) THEN
     389         336 :                   CALL pw_axpy(v_hartree_rspace, v_rspace_new(1), 2.0_dp, 0.0_dp)
     390             : 
     391             :                   CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
     392             :                                           pmat=rho_ao(ispin), &
     393             :                                           hmat=kpp1_env%v_ao(ispin), &
     394             :                                           qs_env=qs_env, &
     395         336 :                                           calculate_forces=.FALSE., gapw=gapw)
     396             :                END IF
     397             :             ELSE
     398             :                ! remove kpp1_env%v_ao and work directly on k_p_p1 ?
     399             :                CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
     400             :                                        pmat=rho_ao(ispin), &
     401             :                                        hmat=kpp1_env%v_ao(ispin), &
     402             :                                        qs_env=qs_env, &
     403          32 :                                        calculate_forces=.FALSE., gapw=gapw_xc)
     404             : 
     405          32 :                IF (ASSOCIATED(v_xc_tau)) THEN
     406             :                   CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
     407             :                                           pmat=rho_ao(ispin), &
     408             :                                           hmat=kpp1_env%v_ao(ispin), &
     409             :                                           qs_env=qs_env, &
     410             :                                           compute_tau=.TRUE., &
     411           0 :                                           calculate_forces=.FALSE., gapw=gapw_xc)
     412             :                END IF
     413             : 
     414          32 :                CALL pw_copy(v_hartree_rspace, v_rspace_new(ispin))
     415             :                CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
     416             :                                        pmat=rho_ao(ispin), &
     417             :                                        hmat=kpp1_env%v_ao(ispin), &
     418             :                                        qs_env=qs_env, &
     419          32 :                                        calculate_forces=.FALSE., gapw=gapw)
     420             :             END IF
     421             : 
     422             :          ELSE
     423             : 
     424        8108 :             IF (nspins == 1) THEN
     425        6428 :                IF (.NOT. (lr_triplet)) THEN
     426        6428 :                   CALL pw_scale(v_rspace_new(1), 2.0_dp)
     427        6428 :                   IF (ASSOCIATED(v_xc_tau)) CALL pw_scale(v_xc_tau(1), 2.0_dp)
     428             :                END IF
     429             :                ! add hartree only for SINGLETS
     430             :                !IF (res_etype == tddfpt_singlet) THEN
     431        6428 :                IF (.NOT. lr_triplet) THEN
     432        6428 :                   CALL pw_axpy(v_hartree_rspace, v_rspace_new(1), 2.0_dp)
     433             :                END IF
     434             :             ELSE
     435        1680 :                CALL pw_axpy(v_hartree_rspace, v_rspace_new(ispin), 1.0_dp)
     436             :             END IF
     437             : 
     438        8108 :             IF (lrigpw) THEN
     439          72 :                IF (ASSOCIATED(v_xc_tau)) &
     440           0 :                   CPABORT("metaGGA-functionals not supported with LRI!")
     441             : 
     442          72 :                lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
     443          72 :                CALL get_qs_env(qs_env, nkind=nkind)
     444         216 :                DO ikind = 1, nkind
     445       43008 :                   lri_v_int(ikind)%v_int = 0.0_dp
     446             :                END DO
     447             :                CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
     448          72 :                                                   lri_v_int, .FALSE., "LRI_AUX")
     449         216 :                DO ikind = 1, nkind
     450       85800 :                   CALL para_env%sum(lri_v_int(ikind)%v_int)
     451             :                END DO
     452         144 :                ALLOCATE (k1mat(1))
     453          72 :                k1mat(1)%matrix => kpp1_env%v_ao(ispin)%matrix
     454          72 :                IF (lri_env%exact_1c_terms) THEN
     455             :                   CALL integrate_v_rspace_diagonal(v_rspace_new(ispin), k1mat(1)%matrix, &
     456           0 :                                                    rho_ao(ispin)%matrix, qs_env, .FALSE., "ORB")
     457             :                END IF
     458          72 :                CALL calculate_lri_ks_matrix(lri_env, lri_v_int, k1mat, atomic_kind_set)
     459          72 :                DEALLOCATE (k1mat)
     460             :             ELSE
     461             :                CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
     462             :                                        pmat=rho_ao(ispin), &
     463             :                                        hmat=kpp1_env%v_ao(ispin), &
     464             :                                        qs_env=qs_env, &
     465        8036 :                                        calculate_forces=.FALSE., gapw=gapw)
     466             : 
     467        8036 :                IF (ASSOCIATED(v_xc_tau)) THEN
     468             :                   CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
     469             :                                           pmat=rho_ao(ispin), &
     470             :                                           hmat=kpp1_env%v_ao(ispin), &
     471             :                                           qs_env=qs_env, &
     472             :                                           compute_tau=.TRUE., &
     473         196 :                                           calculate_forces=.FALSE., gapw=gapw)
     474             :                END IF
     475             :             END IF
     476             : 
     477             :          END IF
     478             : 
     479       16096 :          CALL dbcsr_copy(p_env%kpp1(ispin)%matrix, kpp1_env%v_ao(ispin)%matrix)
     480             :       END DO
     481             : 
     482        7620 :       IF (gapw) THEN
     483         792 :          IF (.NOT. ((nspins == 1 .AND. lr_triplet))) THEN
     484             :             CALL Vh_1c_gg_integrals(qs_env, energy_hartree_1c, &
     485             :                                     p_env%hartree_local%ecoul_1c, &
     486             :                                     p_env%local_rho_set, &
     487         792 :                                     para_env, tddft=.TRUE., core_2nd=.TRUE.)
     488             : 
     489             :             CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, &
     490             :                                        calculate_forces=.FALSE., &
     491         792 :                                        local_rho_set=p_env%local_rho_set)
     492             :          END IF
     493             :          ! ***  Add single atom contributions to the KS matrix ***
     494             :          ! remap pointer
     495         792 :          ns = SIZE(p_env%kpp1)
     496         792 :          ksmat(1:ns, 1:1) => p_env%kpp1(1:ns)
     497         792 :          ns = SIZE(rho_ao)
     498         792 :          psmat(1:ns, 1:1) => rho_ao(1:ns)
     499             :          CALL update_ks_atom(qs_env, ksmat, psmat, forces=.FALSE., tddft=.TRUE., &
     500         792 :                              rho_atom_external=p_env%local_rho_set%rho_atom_set)
     501        6828 :       ELSEIF (gapw_xc) THEN
     502         352 :          ns = SIZE(p_env%kpp1)
     503         352 :          ksmat(1:ns, 1:1) => p_env%kpp1(1:ns)
     504         352 :          ns = SIZE(rho_ao)
     505         352 :          psmat(1:ns, 1:1) => rho_ao(1:ns)
     506             :          CALL update_ks_atom(qs_env, ksmat, psmat, forces=.FALSE., tddft=.TRUE., &
     507         352 :                              rho_atom_external=p_env%local_rho_set%rho_atom_set)
     508             :       END IF
     509             : 
     510             :       ! KG embedding, contribution of kinetic energy functional to kernel
     511        7620 :       IF (dft_control%qs_control%do_kg .AND. .NOT. (lr_triplet .OR. gapw .OR. gapw_xc)) THEN
     512          16 :          IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed) THEN
     513             : 
     514          10 :             CALL qs_rho_get(rho1, rho_ao=rho1_ao)
     515          10 :             alpha = 1.0_dp
     516             : 
     517             :             ekin_mol = 0.0_dp
     518          10 :             CALL get_qs_env(qs_env, kg_env=kg_env)
     519             :             CALL kg_ekin_subset(qs_env=qs_env, &
     520             :                                 ks_matrix=p_env%kpp1, &
     521             :                                 ekin_mol=ekin_mol, &
     522             :                                 calc_force=.FALSE., &
     523             :                                 do_kernel=.TRUE., &
     524          10 :                                 pmat_ext=rho1_ao)
     525             :          END IF
     526             :       END IF
     527             : 
     528        7620 :       CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
     529        7620 :       CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
     530       16096 :       DO ispin = 1, nspins
     531       16096 :          CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
     532             :       END DO
     533        7620 :       DEALLOCATE (v_rspace_new)
     534        7620 :       IF (ASSOCIATED(v_xc_tau)) THEN
     535         392 :          DO ispin = 1, nspins
     536         392 :             CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
     537             :          END DO
     538         196 :          DEALLOCATE (v_xc_tau)
     539             :       END IF
     540             : 
     541        7620 :       CALL timestop(handle)
     542             : 
     543        7620 :    END SUBROUTINE apply_op_2_dft
     544             : 
     545             : ! **************************************************************************************************
     546             : !> \brief ...
     547             : !> \param qs_env ...
     548             : !> \param p_env ...
     549             : ! **************************************************************************************************
     550         110 :    SUBROUTINE apply_op_2_xtb(qs_env, p_env)
     551             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     552             :       TYPE(qs_p_env_type)                                :: p_env
     553             : 
     554             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'apply_op_2_xtb'
     555             : 
     556             :       INTEGER                                            :: atom_a, handle, iatom, ikind, is, ispin, &
     557             :                                                             na, natom, natorb, nkind, ns, nsgf, &
     558             :                                                             nspins
     559             :       INTEGER, DIMENSION(25)                             :: lao
     560             :       INTEGER, DIMENSION(5)                              :: occ
     561             :       LOGICAL                                            :: lr_triplet
     562         110 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: mcharge, mcharge1
     563         110 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: aocg, aocg1, charges, charges1
     564         110 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     565         110 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao
     566         110 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_p, matrix_p1, matrix_s
     567             :       TYPE(dft_control_type), POINTER                    :: dft_control
     568             :       TYPE(linres_control_type), POINTER                 :: linres_control
     569             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     570         110 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     571             :       TYPE(pw_env_type), POINTER                         :: pw_env
     572         110 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     573             :       TYPE(qs_kpp1_env_type), POINTER                    :: kpp1_env
     574             :       TYPE(qs_rho_type), POINTER                         :: rho, rho1
     575             :       TYPE(xtb_atom_type), POINTER                       :: xtb_kind
     576             : 
     577         110 :       CALL timeset(routineN, handle)
     578             : 
     579         110 :       CPASSERT(ASSOCIATED(p_env%kpp1_env))
     580         110 :       CPASSERT(ASSOCIATED(p_env%kpp1))
     581         110 :       kpp1_env => p_env%kpp1_env
     582             : 
     583         110 :       rho1 => p_env%rho1
     584         110 :       CPASSERT(ASSOCIATED(rho1))
     585             : 
     586             :       CALL get_qs_env(qs_env=qs_env, &
     587             :                       pw_env=pw_env, &
     588             :                       para_env=para_env, &
     589             :                       rho=rho, &
     590             :                       linres_control=linres_control, &
     591         110 :                       dft_control=dft_control)
     592             : 
     593         110 :       CALL qs_rho_get(rho, rho_ao=rho_ao)
     594             : 
     595         110 :       lr_triplet = linres_control%lr_triplet
     596         110 :       CPASSERT(.NOT. lr_triplet)
     597             : 
     598         110 :       nspins = SIZE(p_env%kpp1)
     599             : 
     600         220 :       DO ispin = 1, nspins
     601         220 :          CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp)
     602             :       END DO
     603             : 
     604         110 :       IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
     605             :          ! Mulliken charges
     606         106 :          CALL get_qs_env(qs_env, particle_set=particle_set, matrix_s_kp=matrix_s)
     607         106 :          natom = SIZE(particle_set)
     608         106 :          CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
     609         106 :          CALL qs_rho_get(rho1, rho_ao_kp=matrix_p1)
     610         530 :          ALLOCATE (mcharge(natom), charges(natom, 5))
     611         318 :          ALLOCATE (mcharge1(natom), charges1(natom, 5))
     612        8066 :          charges = 0.0_dp
     613        8066 :          charges1 = 0.0_dp
     614         106 :          CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
     615         106 :          nkind = SIZE(atomic_kind_set)
     616         106 :          CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
     617         424 :          ALLOCATE (aocg(nsgf, natom))
     618        7536 :          aocg = 0.0_dp
     619         318 :          ALLOCATE (aocg1(nsgf, natom))
     620        7536 :          aocg1 = 0.0_dp
     621         106 :          CALL ao_charges(matrix_p, matrix_s, aocg, para_env)
     622         106 :          CALL ao_charges(matrix_p1, matrix_s, aocg1, para_env)
     623         370 :          DO ikind = 1, nkind
     624         264 :             CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
     625         264 :             CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
     626         264 :             CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
     627        2120 :             DO iatom = 1, na
     628        1486 :                atom_a = atomic_kind_set(ikind)%atom_list(iatom)
     629        8916 :                charges(atom_a, :) = REAL(occ(:), KIND=dp)
     630        5782 :                DO is = 1, natorb
     631        4032 :                   ns = lao(is) + 1
     632        4032 :                   charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
     633        5518 :                   charges1(atom_a, ns) = charges1(atom_a, ns) - aocg1(is, atom_a)
     634             :                END DO
     635             :             END DO
     636             :          END DO
     637         106 :          DEALLOCATE (aocg, aocg1)
     638        1592 :          DO iatom = 1, natom
     639        8916 :             mcharge(iatom) = SUM(charges(iatom, :))
     640        9022 :             mcharge1(iatom) = SUM(charges1(iatom, :))
     641             :          END DO
     642             :          ! Coulomb Kernel
     643         106 :          CALL xtb_coulomb_hessian(qs_env, p_env%kpp1, charges1, mcharge1, mcharge)
     644             :          !
     645         212 :          DEALLOCATE (charges, mcharge, charges1, mcharge1)
     646             :       END IF
     647             : 
     648         110 :       CALL timestop(handle)
     649             : 
     650         220 :    END SUBROUTINE apply_op_2_xtb
     651             : 
     652             : ! **************************************************************************************************
     653             : !> \brief Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
     654             : !> \param qs_env ...
     655             : !> \param p_env ...
     656             : !> \par History
     657             : !>    * 11.2019 adapted from tddfpt_apply_hfx
     658             : ! **************************************************************************************************
     659       16392 :    SUBROUTINE apply_hfx(qs_env, p_env)
     660             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     661             :       TYPE(qs_p_env_type)                                :: p_env
     662             : 
     663             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'apply_hfx'
     664             : 
     665             :       INTEGER                                            :: handle, ispin, nspins
     666             :       LOGICAL                                            :: do_hfx
     667             :       REAL(KIND=dp)                                      :: alpha
     668             :       TYPE(cp_logger_type), POINTER                      :: logger
     669        8196 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: h1_mat, matrix_s, rho1_ao, work
     670             :       TYPE(dft_control_type), POINTER                    :: dft_control
     671             :       TYPE(section_vals_type), POINTER                   :: hfx_section, input
     672             : 
     673        8196 :       CALL timeset(routineN, handle)
     674             : 
     675        8196 :       logger => cp_get_default_logger()
     676             : 
     677             :       CALL get_qs_env(qs_env=qs_env, &
     678             :                       input=input, &
     679             :                       matrix_s=matrix_s, &
     680        8196 :                       dft_control=dft_control)
     681        8196 :       nspins = dft_control%nspins
     682             : 
     683        8196 :       hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
     684        8196 :       CALL section_vals_get(hfx_section, explicit=do_hfx)
     685             : 
     686        8196 :       IF (do_hfx) THEN
     687             : 
     688        3384 :          IF (dft_control%do_admm) THEN
     689        1860 :             IF (dft_control%admm_control%purification_method /= do_admm_purify_none) THEN
     690           0 :                CPABORT("ADMM: Linear Response needs purification_method=none")
     691             :             END IF
     692        1860 :             IF (dft_control%admm_control%scaling_model /= do_admm_exch_scaling_none) THEN
     693           0 :                CPABORT("ADMM: Linear Response needs scaling_model=none")
     694             :             END IF
     695        1860 :             IF (dft_control%admm_control%method /= do_admm_basis_projection) THEN
     696           0 :                CPABORT("ADMM: Linear Response needs admm_method=basis_projection")
     697             :             END IF
     698             :             !
     699        1860 :             rho1_ao => p_env%p1_admm
     700        1860 :             h1_mat => p_env%kpp1_admm
     701             :          ELSE
     702        1524 :             rho1_ao => p_env%p1
     703        1524 :             h1_mat => p_env%kpp1
     704             :          END IF
     705             : 
     706        3384 :          NULLIFY (work)
     707        3384 :          CALL dbcsr_allocate_matrix_set(work, nspins)
     708        7098 :          DO ispin = 1, nspins
     709        3714 :             ALLOCATE (work(ispin)%matrix)
     710        3714 :             CALL dbcsr_create(work(ispin)%matrix, template=h1_mat(ispin)%matrix)
     711        3714 :             CALL dbcsr_copy(work(ispin)%matrix, h1_mat(ispin)%matrix)
     712        7098 :             CALL dbcsr_set(work(ispin)%matrix, 0.0_dp)
     713             :          END DO
     714             : 
     715        3384 :          CALL hfx_matrix(work, rho1_ao, qs_env, hfx_section)
     716             : 
     717        3384 :          alpha = 2.0_dp
     718        3384 :          IF (nspins == 2) alpha = 1.0_dp
     719             : 
     720        7098 :          DO ispin = 1, nspins
     721        7098 :             CALL dbcsr_add(h1_mat(ispin)%matrix, work(ispin)%matrix, 1.0_dp, alpha)
     722             :          END DO
     723             : 
     724        3384 :          CALL dbcsr_deallocate_matrix_set(work)
     725             : 
     726             :       END IF
     727             : 
     728        8196 :       CALL timestop(handle)
     729             : 
     730        8196 :    END SUBROUTINE apply_hfx
     731             : 
     732             : ! **************************************************************************************************
     733             : !> \brief Add the hfx contributions to the Hamiltonian
     734             : !>
     735             : !> \param matrix_ks ...
     736             : !> \param rho_ao ...
     737             : !> \param qs_env ...
     738             : !> \param hfx_sections ...
     739             : !> \param external_x_data ...
     740             : !> \param ex ...
     741             : !> \note
     742             : !>     Simplified version of subroutine hfx_ks_matrix()
     743             : ! **************************************************************************************************
     744        3384 :    SUBROUTINE hfx_matrix(matrix_ks, rho_ao, qs_env, hfx_sections, external_x_data, ex)
     745             :       TYPE(dbcsr_p_type), DIMENSION(:), TARGET           :: matrix_ks, rho_ao
     746             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     747             :       TYPE(section_vals_type), POINTER                   :: hfx_sections
     748             :       TYPE(hfx_type), DIMENSION(:, :), OPTIONAL, TARGET  :: external_x_data
     749             :       REAL(KIND=dp), OPTIONAL                            :: ex
     750             : 
     751             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'hfx_matrix'
     752             : 
     753             :       INTEGER                                            :: handle, irep, ispin, mspin, n_rep_hf, &
     754             :                                                             nspins
     755             :       LOGICAL                                            :: distribute_fock_matrix, &
     756             :                                                             hfx_treat_lsd_in_core, &
     757             :                                                             s_mstruct_changed
     758             :       REAL(KIND=dp)                                      :: eh1, ehfx
     759        3384 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks_kp, rho_ao_kp
     760             :       TYPE(dft_control_type), POINTER                    :: dft_control
     761        3384 :       TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
     762             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     763             : 
     764        3384 :       CALL timeset(routineN, handle)
     765             : 
     766        3384 :       NULLIFY (dft_control, para_env, matrix_ks_kp, rho_ao_kp, x_data)
     767             : 
     768             :       CALL get_qs_env(qs_env=qs_env, &
     769             :                       dft_control=dft_control, &
     770             :                       para_env=para_env, &
     771             :                       s_mstruct_changed=s_mstruct_changed, &
     772        3384 :                       x_data=x_data)
     773             : 
     774        3384 :       IF (PRESENT(external_x_data)) x_data => external_x_data
     775             : 
     776        3384 :       CPASSERT(dft_control%nimages == 1)
     777        3384 :       nspins = dft_control%nspins
     778             : 
     779        3384 :       CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
     780             :       CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
     781        3384 :                                 i_rep_section=1)
     782             : 
     783        3384 :       CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
     784        3384 :       distribute_fock_matrix = .TRUE.
     785             : 
     786        3384 :       mspin = 1
     787        3384 :       IF (hfx_treat_lsd_in_core) mspin = nspins
     788             : 
     789        3384 :       matrix_ks_kp(1:nspins, 1:1) => matrix_ks(1:nspins)
     790        3384 :       rho_ao_kp(1:nspins, 1:1) => rho_ao(1:nspins)
     791             : 
     792        6768 :       DO irep = 1, n_rep_hf
     793        3384 :          ehfx = 0.0_dp
     794             : 
     795        6768 :          IF (x_data(irep, 1)%do_hfx_ri) THEN
     796             :             CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_kp, ehfx, &
     797             :                                   rho_ao=rho_ao_kp, geometry_did_change=s_mstruct_changed, &
     798         170 :                                   nspins=nspins, hf_fraction=x_data(irep, 1)%general_parameter%fraction)
     799             : 
     800             :          ELSE
     801             : 
     802        6428 :             DO ispin = 1, mspin
     803             :                CALL integrate_four_center(qs_env, x_data, matrix_ks_kp, eh1, rho_ao_kp, hfx_sections, para_env, &
     804        3214 :                                           s_mstruct_changed, irep, distribute_fock_matrix, ispin=ispin)
     805        6428 :                ehfx = ehfx + eh1
     806             :             END DO
     807             : 
     808             :          END IF
     809             :       END DO
     810             : 
     811             :       ! Export energy
     812        3384 :       IF (PRESENT(ex)) ex = ehfx
     813             : 
     814        3384 :       CALL timestop(handle)
     815             : 
     816        3384 :    END SUBROUTINE hfx_matrix
     817             : 
     818             : ! **************************************************************************************************
     819             : !> \brief ...
     820             : !> \param qs_env ...
     821             : !> \param p_env ...
     822             : ! **************************************************************************************************
     823       16392 :    SUBROUTINE apply_xc_admm(qs_env, p_env)
     824             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     825             :       TYPE(qs_p_env_type)                                :: p_env
     826             : 
     827             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'apply_xc_admm'
     828             : 
     829             :       CHARACTER(LEN=default_string_length)               :: basis_type
     830             :       INTEGER                                            :: handle, ispin, ns, nspins
     831             :       INTEGER, DIMENSION(2, 3)                           :: bo
     832             :       LOGICAL                                            :: gapw, lsd
     833             :       REAL(KIND=dp)                                      :: alpha
     834             :       TYPE(admm_type), POINTER                           :: admm_env
     835             :       TYPE(dbcsr_p_type)                                 :: xcmat
     836        8196 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     837        8196 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: ksmat, psmat
     838             :       TYPE(dft_control_type), POINTER                    :: dft_control
     839             :       TYPE(linres_control_type), POINTER                 :: linres_control
     840             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     841             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     842        8196 :          POINTER                                         :: sab_aux_fit
     843        8196 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho1_aux_g
     844             :       TYPE(pw_env_type), POINTER                         :: pw_env
     845             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     846       16392 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho1_aux_r, tau_pw, v_xc, v_xc_tau
     847       16392 :       TYPE(rho_atom_type), DIMENSION(:), POINTER         :: rho1_atom_set, rho_atom_set
     848             :       TYPE(section_vals_type), POINTER                   :: xc_fun_section, xc_section
     849             :       TYPE(task_list_type), POINTER                      :: task_list
     850             :       TYPE(xc_rho_cflags_type)                           :: needs
     851             :       TYPE(xc_rho_set_type)                              :: rho1_set
     852             : 
     853        8196 :       CALL timeset(routineN, handle)
     854             : 
     855        8196 :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
     856             : 
     857        8196 :       IF (dft_control%do_admm) THEN
     858        1860 :          IF (qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
     859             :             ! nothing to do
     860             :          ELSE
     861        1068 :             CALL get_qs_env(qs_env=qs_env, linres_control=linres_control)
     862        1068 :             CPASSERT(.NOT. dft_control%qs_control%lrigpw)
     863        1068 :             CPASSERT(.NOT. linres_control%lr_triplet)
     864             : 
     865        1068 :             nspins = dft_control%nspins
     866             : 
     867             :             ! AUX basis contribution
     868        1068 :             CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     869        1068 :             CPASSERT(ASSOCIATED(pw_env))
     870        1068 :             CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     871        1068 :             NULLIFY (tau_pw)
     872             :             ! calculate the xc potential
     873        1068 :             lsd = (nspins == 2)
     874        1068 :             CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s)
     875        1068 :             ALLOCATE (xcmat%matrix)
     876        1068 :             CALL dbcsr_create(xcmat%matrix, template=matrix_s(1)%matrix)
     877             : 
     878        1068 :             CALL get_qs_env(qs_env, admm_env=admm_env)
     879        1068 :             gapw = admm_env%do_gapw
     880             : 
     881        1068 :             CALL qs_rho_get(p_env%rho1_admm, rho_r=rho1_aux_r, rho_g=rho1_aux_g)
     882        1068 :             xc_section => admm_env%xc_section_aux
     883       10680 :             bo = rho1_aux_r(1)%pw_grid%bounds_local
     884             :             ! create the place where to store the argument for the functionals
     885             :             CALL xc_rho_set_create(rho1_set, bo, &
     886             :                                    rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
     887             :                                    drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
     888        1068 :                                    tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
     889             : 
     890        1068 :             xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
     891        1068 :             needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
     892             : 
     893             :             ! calculate the arguments needed by the functionals
     894             :             CALL xc_rho_set_update(rho1_set, rho1_aux_r, rho1_aux_g, tau_pw, needs, &
     895             :                                    section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
     896             :                                    section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
     897        1068 :                                    auxbas_pw_pool)
     898             :             CALL xc_calc_2nd_deriv(v_xc, v_xc_tau, p_env%kpp1_env%deriv_set_admm, p_env%kpp1_env%rho_set_admm, &
     899             :                                    rho1_aux_r, rho1_aux_g, tau_pw, auxbas_pw_pool, gapw=.FALSE., &
     900        1068 :                                    xc_section=xc_section)
     901        1068 :             IF (ASSOCIATED(v_xc_tau)) THEN
     902           0 :                CPABORT("Meta-GGA ADMM functionals not yet supported!")
     903             :             END IF
     904        1068 :             CALL xc_rho_set_release(rho1_set)
     905             : 
     906        1068 :             basis_type = "AUX_FIT"
     907        1068 :             CALL get_qs_env(qs_env, para_env=para_env)
     908        1068 :             CALL get_admm_env(admm_env, task_list_aux_fit=task_list)
     909        1068 :             IF (admm_env%do_gapw) THEN
     910             :                CALL prepare_gapw_den(qs_env, local_rho_set=p_env%local_rho_set_admm, &
     911         160 :                                      do_rho0=.FALSE., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
     912         160 :                rho_atom_set => admm_env%admm_gapw_env%local_rho_set%rho_atom_set
     913         160 :                rho1_atom_set => p_env%local_rho_set_admm%rho_atom_set
     914             :                CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
     915         160 :                                                 kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
     916         160 :                basis_type = "AUX_FIT_SOFT"
     917         160 :                task_list => admm_env%admm_gapw_env%task_list
     918             :             END IF
     919             : 
     920        1068 :             alpha = 1.0_dp
     921        1068 :             IF (nspins == 1) alpha = 2.0_dp
     922             : 
     923        2234 :             DO ispin = 1, nspins
     924        1166 :                CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
     925        1166 :                CALL dbcsr_copy(xcmat%matrix, matrix_s(1)%matrix)
     926        1166 :                CALL dbcsr_set(xcmat%matrix, 0.0_dp)
     927             :                CALL integrate_v_rspace(v_rspace=v_xc(ispin), hmat=xcmat, qs_env=qs_env, &
     928             :                                        calculate_forces=.FALSE., basis_type=basis_type, &
     929        1166 :                                        task_list_external=task_list)
     930        2234 :                CALL dbcsr_add(p_env%kpp1_admm(ispin)%matrix, xcmat%matrix, 1.0_dp, alpha)
     931             :             END DO
     932             : 
     933        1068 :             IF (admm_env%do_gapw) THEN
     934         160 :                CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
     935         160 :                ns = SIZE(p_env%kpp1_admm)
     936         160 :                ksmat(1:ns, 1:1) => p_env%kpp1_admm(1:ns)
     937         160 :                psmat(1:ns, 1:1) => p_env%p1_admm(1:ns)
     938             :                CALL update_ks_atom(qs_env, ksmat, psmat, forces=.FALSE., tddft=.TRUE., &
     939             :                                    rho_atom_external=p_env%local_rho_set_admm%rho_atom_set, &
     940             :                                    kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
     941             :                                    oce_external=admm_env%admm_gapw_env%oce, &
     942         160 :                                    sab_external=sab_aux_fit)
     943             :             END IF
     944             : 
     945        2234 :             DO ispin = 1, nspins
     946        2234 :                CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
     947             :             END DO
     948        1068 :             DEALLOCATE (v_xc)
     949        1068 :             CALL dbcsr_deallocate_matrix(xcmat%matrix)
     950             : 
     951             :          END IF
     952             :       END IF
     953             : 
     954        8196 :       CALL timestop(handle)
     955             : 
     956      180312 :    END SUBROUTINE apply_xc_admm
     957             : 
     958             : END MODULE qs_linres_kernel

Generated by: LCOV version 1.15