LCOV - code coverage report
Current view: top level - src - qs_tddfpt2_fhxc.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 95.1 % 163 155
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
       9              :    USE admm_types,                      ONLY: admm_type
      10              :    USE cp_control_types,                ONLY: dft_control_type,&
      11              :                                               stda_control_type
      12              :    USE cp_dbcsr_api,                    ONLY: &
      13              :         dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_get_info, &
      14              :         dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_symmetric
      15              :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      16              :    USE cp_dbcsr_operations,             ONLY: copy_fm_to_dbcsr,&
      17              :                                               cp_dbcsr_plus_fm_fm_t,&
      18              :                                               cp_dbcsr_sm_fm_multiply
      19              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      20              :                                               cp_fm_get_info,&
      21              :                                               cp_fm_release,&
      22              :                                               cp_fm_type
      23              :    USE input_constants,                 ONLY: do_admm_aux_exch_func_none
      24              :    USE kinds,                           ONLY: default_string_length,&
      25              :                                               dp
      26              :    USE lri_environment_types,           ONLY: lri_kind_type
      27              :    USE message_passing,                 ONLY: mp_para_env_type
      28              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      29              :    USE pw_env_types,                    ONLY: pw_env_get
      30              :    USE pw_methods,                      ONLY: pw_axpy,&
      31              :                                               pw_scale,&
      32              :                                               pw_zero
      33              :    USE pw_pool_types,                   ONLY: pw_pool_type
      34              :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      35              :                                               pw_r3d_rs_type
      36              :    USE qs_environment_types,            ONLY: get_qs_env,&
      37              :                                               qs_environment_type
      38              :    USE qs_gapw_densities,               ONLY: prepare_gapw_den
      39              :    USE qs_integrate_potential,          ONLY: integrate_v_rspace,&
      40              :                                               integrate_v_rspace_one_center
      41              :    USE qs_kernel_types,                 ONLY: full_kernel_env_type
      42              :    USE qs_ks_atom,                      ONLY: update_ks_atom
      43              :    USE qs_rho_atom_types,               ONLY: rho_atom_type
      44              :    USE qs_rho_methods,                  ONLY: qs_rho_update_rho,&
      45              :                                               qs_rho_update_tddfpt
      46              :    USE qs_rho_types,                    ONLY: qs_rho_get
      47              :    USE qs_tddfpt2_densities,            ONLY: tddfpt_construct_aux_fit_density
      48              :    USE qs_tddfpt2_lri_utils,            ONLY: tddfpt2_lri_Amat
      49              :    USE qs_tddfpt2_operators,            ONLY: tddfpt_apply_coulomb,&
      50              :                                               tddfpt_apply_xc,&
      51              :                                               tddfpt_apply_xc_potential
      52              :    USE qs_tddfpt2_stda_types,           ONLY: stda_env_type
      53              :    USE qs_tddfpt2_stda_utils,           ONLY: stda_calculate_kernel
      54              :    USE qs_tddfpt2_subgroups,            ONLY: tddfpt_subgroup_env_type
      55              :    USE qs_tddfpt2_types,                ONLY: tddfpt_work_matrices
      56              :    USE qs_vxc_atom,                     ONLY: calculate_xc_2nd_deriv_atom
      57              :    USE task_list_types,                 ONLY: task_list_type
      58              : #include "./base/base_uses.f90"
      59              : 
      60              :    IMPLICIT NONE
      61              : 
      62              :    PRIVATE
      63              : 
      64              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_fhxc'
      65              : 
      66              :    INTEGER, PARAMETER, PRIVATE          :: maxspins = 2
      67              : 
      68              :    PUBLIC :: fhxc_kernel, stda_kernel
      69              : 
      70              : ! **************************************************************************************************
      71              : 
      72              : CONTAINS
      73              : 
      74              : ! **************************************************************************************************
      75              : !> \brief Compute action matrix-vector products with the FHxc Kernel
      76              : !> \param Aop_evects            action of TDDFPT operator on trial vectors (modified on exit)
      77              : !> \param evects                TDDFPT trial vectors
      78              : !> \param is_rks_triplets       indicates that a triplet excited states calculation using
      79              : !>                              spin-unpolarised molecular orbitals has been requested
      80              : !> \param do_hfx                flag that activates computation of exact-exchange terms
      81              : !> \param do_admm ...
      82              : !> \param qs_env                Quickstep environment
      83              : !> \param kernel_env            kernel environment
      84              : !> \param kernel_env_admm_aux   kernel environment for ADMM correction
      85              : !> \param sub_env               parallel (sub)group environment
      86              : !> \param work_matrices         collection of work matrices (modified on exit)
      87              : !> \param admm_symm             use symmetric definition of ADMM kernel correction
      88              : !> \param admm_xc_correction    use ADMM XC kernel correction
      89              : !> \param do_lrigpw ...
      90              : !> \param tddfpt_mgrid ...
      91              : !> \par History
      92              : !>    * 06.2016 created [Sergey Chulkov]
      93              : !>    * 03.2017 refactored [Sergey Chulkov]
      94              : !>    * 04.2019 refactored [JHU]
      95              : ! **************************************************************************************************
      96         2898 :    SUBROUTINE fhxc_kernel(Aop_evects, evects, is_rks_triplets, &
      97              :                           do_hfx, do_admm, qs_env, kernel_env, kernel_env_admm_aux, &
      98              :                           sub_env, work_matrices, admm_symm, admm_xc_correction, do_lrigpw, &
      99              :                           tddfpt_mgrid)
     100              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT)   :: Aop_evects
     101              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: evects
     102              :       LOGICAL, INTENT(in)                                :: is_rks_triplets, do_hfx, do_admm
     103              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     104              :       TYPE(full_kernel_env_type), POINTER                :: kernel_env, kernel_env_admm_aux
     105              :       TYPE(tddfpt_subgroup_env_type), INTENT(in)         :: sub_env
     106              :       TYPE(tddfpt_work_matrices), INTENT(inout)          :: work_matrices
     107              :       LOGICAL, INTENT(in)                                :: admm_symm, admm_xc_correction, &
     108              :                                                             do_lrigpw, tddfpt_mgrid
     109              : 
     110              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'fhxc_kernel'
     111              : 
     112              :       CHARACTER(LEN=default_string_length)               :: basis_type
     113              :       INTEGER                                            :: handle, ikind, ispin, ivect, nao, &
     114              :                                                             nao_aux, nkind, nspins, nvects
     115         2898 :       INTEGER, DIMENSION(:), POINTER                     :: blk_sizes
     116              :       INTEGER, DIMENSION(maxspins)                       :: nactive
     117              :       LOGICAL                                            :: gapw, gapw_xc
     118              :       TYPE(admm_type), POINTER                           :: admm_env
     119              :       TYPE(cp_fm_type)                                   :: work_aux_orb, work_orb_orb
     120         2898 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: A_xc_munu_sub, rho_ia_ao, &
     121         2898 :                                                             rho_ia_ao_aux_fit
     122              :       TYPE(dbcsr_type), POINTER                          :: dbwork
     123              :       TYPE(dft_control_type), POINTER                    :: dft_control
     124         2898 :       TYPE(lri_kind_type), DIMENSION(:), POINTER         :: lri_v_int
     125              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     126         2898 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_ia_g, rho_ia_g_aux_fit
     127              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     128         2898 :       TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:)    :: V_rspace_sub
     129         2898 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_ia_r, rho_ia_r_aux_fit
     130         2898 :       TYPE(rho_atom_type), DIMENSION(:), POINTER         :: rho1_atom_set, rho_atom_set
     131              :       TYPE(task_list_type), POINTER                      :: task_list
     132              : 
     133         2898 :       CALL timeset(routineN, handle)
     134              : 
     135         2898 :       nspins = SIZE(evects, 1)
     136         2898 :       nvects = SIZE(evects, 2)
     137         2898 :       IF (do_admm) THEN
     138          622 :          CPASSERT(do_hfx)
     139          622 :          CPASSERT(ASSOCIATED(sub_env%admm_A))
     140              :       END IF
     141         2898 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     142              : 
     143         2898 :       gapw = dft_control%qs_control%gapw
     144         2898 :       gapw_xc = dft_control%qs_control%gapw_xc
     145              : 
     146         2898 :       CALL cp_fm_get_info(evects(1, 1), nrow_global=nao)
     147         6318 :       DO ispin = 1, nspins
     148         6318 :          CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive(ispin))
     149              :       END DO
     150              : 
     151              :       CALL qs_rho_get(work_matrices%rho_orb_struct_sub, rho_ao=rho_ia_ao, &
     152         2898 :                       rho_g=rho_ia_g, rho_r=rho_ia_r)
     153         2898 :       IF (do_hfx .AND. do_admm) THEN
     154          622 :          CALL get_qs_env(qs_env, admm_env=admm_env)
     155              :          CALL qs_rho_get(work_matrices%rho_aux_fit_struct_sub, &
     156              :                          rho_ao=rho_ia_ao_aux_fit, rho_g=rho_ia_g_aux_fit, &
     157          622 :                          rho_r=rho_ia_r_aux_fit)
     158              :       END IF
     159              : 
     160         9616 :       DO ivect = 1, nvects
     161         6718 :          IF (ALLOCATED(work_matrices%evects_sub)) THEN
     162           16 :             IF (ASSOCIATED(work_matrices%evects_sub(1, ivect)%matrix_struct)) THEN
     163           16 :                DO ispin = 1, nspins
     164            8 :                   CALL dbcsr_set(rho_ia_ao(ispin)%matrix, 0.0_dp)
     165              :                   CALL cp_dbcsr_plus_fm_fm_t(rho_ia_ao(ispin)%matrix, &
     166              :                                              matrix_v=sub_env%mos_occ(ispin), &
     167              :                                              matrix_g=work_matrices%evects_sub(ispin, ivect), &
     168           16 :                                              ncol=nactive(ispin), symmetry_mode=1)
     169              :                END DO
     170              :             ELSE
     171              :                ! skip trial vectors which are assigned to different parallel groups
     172              :                CYCLE
     173              :             END IF
     174              :          ELSE
     175        14762 :             DO ispin = 1, nspins
     176         8060 :                CALL dbcsr_set(rho_ia_ao(ispin)%matrix, 0.0_dp)
     177              :                CALL cp_dbcsr_plus_fm_fm_t(rho_ia_ao(ispin)%matrix, &
     178              :                                           matrix_v=sub_env%mos_occ(ispin), &
     179              :                                           matrix_g=evects(ispin, ivect), &
     180        14762 :                                           ncol=nactive(ispin), symmetry_mode=1)
     181              :             END DO
     182              :          END IF
     183              : 
     184         6710 :          IF (do_lrigpw) THEN
     185              :             CALL qs_rho_update_tddfpt(work_matrices%rho_orb_struct_sub, qs_env, &
     186              :                                       pw_env_external=sub_env%pw_env, &
     187              :                                       task_list_external=sub_env%task_list_orb, &
     188              :                                       para_env_external=sub_env%para_env, &
     189              :                                       tddfpt_lri_env=kernel_env%lri_env, &
     190          172 :                                       tddfpt_lri_density=kernel_env%lri_density)
     191         6538 :          ELSEIF (dft_control%qs_control%lrigpw .OR. &
     192              :                  dft_control%qs_control%rigpw) THEN
     193              :             CALL qs_rho_update_tddfpt(work_matrices%rho_orb_struct_sub, qs_env, &
     194              :                                       pw_env_external=sub_env%pw_env, &
     195              :                                       task_list_external=sub_env%task_list_orb, &
     196            0 :                                       para_env_external=sub_env%para_env)
     197              :          ELSE
     198         6538 :             IF (gapw) THEN
     199              :                CALL qs_rho_update_rho(work_matrices%rho_orb_struct_sub, qs_env, &
     200              :                                       local_rho_set=work_matrices%local_rho_set, &
     201              :                                       pw_env_external=sub_env%pw_env, &
     202              :                                       task_list_external=sub_env%task_list_orb_soft, &
     203         1106 :                                       para_env_external=sub_env%para_env)
     204              :                CALL prepare_gapw_den(qs_env, work_matrices%local_rho_set, &
     205         1106 :                                      do_rho0=(.NOT. is_rks_triplets), pw_env_sub=sub_env%pw_env)
     206         5432 :             ELSEIF (gapw_xc) THEN
     207              :                CALL qs_rho_update_rho(work_matrices%rho_orb_struct_sub, qs_env, &
     208              :                                       rho_xc_external=work_matrices%rho_xc_struct_sub, &
     209              :                                       local_rho_set=work_matrices%local_rho_set, &
     210              :                                       pw_env_external=sub_env%pw_env, &
     211              :                                       task_list_external=sub_env%task_list_orb, &
     212              :                                       task_list_external_soft=sub_env%task_list_orb_soft, &
     213          218 :                                       para_env_external=sub_env%para_env)
     214              :                CALL prepare_gapw_den(qs_env, work_matrices%local_rho_set, do_rho0=.FALSE., &
     215          218 :                                      pw_env_sub=sub_env%pw_env)
     216              :             ELSE
     217              :                CALL qs_rho_update_rho(work_matrices%rho_orb_struct_sub, qs_env, &
     218              :                                       pw_env_external=sub_env%pw_env, &
     219              :                                       task_list_external=sub_env%task_list_orb, &
     220         5214 :                                       para_env_external=sub_env%para_env)
     221              :             END IF
     222              :          END IF
     223              : 
     224        14778 :          DO ispin = 1, nspins
     225        14778 :             CALL dbcsr_set(work_matrices%A_ia_munu_sub(ispin)%matrix, 0.0_dp)
     226              :          END DO
     227              : 
     228              :          ! electron-hole exchange-correlation interaction
     229        14778 :          DO ispin = 1, nspins
     230        14778 :             CALL pw_zero(work_matrices%A_ia_rspace_sub(ispin))
     231              :          END DO
     232              : 
     233              :          ! C_x d^{2}E_{x}^{DFT}[\rho] / d\rho^2
     234              :          ! + C_{HF} d^{2}E_{x, ADMM}^{DFT}[\rho] / d\rho^2 in case of ADMM calculation
     235         6710 :          IF (gapw_xc) THEN
     236          218 :             IF (kernel_env%do_exck) THEN
     237            0 :                CPABORT("NYA")
     238              :             ELSE
     239              :                CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, &
     240              :                                     rho_ia_struct=work_matrices%rho_xc_struct_sub, &
     241              :                                     is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
     242              :                                     work_v_xc=work_matrices%wpw_rspace_sub, &
     243          218 :                                     work_v_xc_tau=work_matrices%wpw_tau_rspace_sub)
     244              :             END IF
     245          436 :             DO ispin = 1, nspins
     246              :                CALL pw_scale(work_matrices%A_ia_rspace_sub(ispin), &
     247          218 :                              work_matrices%A_ia_rspace_sub(ispin)%pw_grid%dvol)
     248              :                CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), &
     249              :                                        hmat=work_matrices%A_ia_munu_sub(ispin), &
     250              :                                        qs_env=qs_env, calculate_forces=.FALSE., gapw=gapw_xc, &
     251              :                                        pw_env_external=sub_env%pw_env, &
     252          218 :                                        task_list_external=sub_env%task_list_orb_soft)
     253          436 :                CALL pw_zero(work_matrices%A_ia_rspace_sub(ispin))
     254              :             END DO
     255              :          ELSE
     256         6492 :             IF (kernel_env%do_exck) THEN
     257              :                CALL tddfpt_apply_xc_potential(work_matrices%A_ia_rspace_sub, work_matrices%fxc_rspace_sub, &
     258          156 :                                               work_matrices%rho_orb_struct_sub, is_rks_triplets)
     259              :             ELSE
     260              :                CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, &
     261              :                                     rho_ia_struct=work_matrices%rho_orb_struct_sub, &
     262              :                                     is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
     263              :                                     work_v_xc=work_matrices%wpw_rspace_sub, &
     264         6336 :                                     work_v_xc_tau=work_matrices%wpw_tau_rspace_sub)
     265              :             END IF
     266              :          END IF
     267         6710 :          IF (gapw .OR. gapw_xc) THEN
     268         1324 :             rho_atom_set => sub_env%local_rho_set%rho_atom_set
     269         1324 :             rho1_atom_set => work_matrices%local_rho_set%rho_atom_set
     270              :             CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, kernel_env%xc_section, &
     271         1324 :                                              sub_env%para_env, do_tddfpt2=.TRUE., do_triplet=is_rks_triplets)
     272              :          END IF
     273              : 
     274              :          ! ADMM correction
     275         6710 :          IF (do_admm .AND. admm_xc_correction) THEN
     276         1134 :             IF (dft_control%admm_control%aux_exch_func /= do_admm_aux_exch_func_none) THEN
     277              :                CALL tddfpt_construct_aux_fit_density(rho_orb_struct=work_matrices%rho_orb_struct_sub, &
     278              :                                                      rho_aux_fit_struct=work_matrices%rho_aux_fit_struct_sub, &
     279              :                                                      local_rho_set=work_matrices%local_rho_set_admm, &
     280              :                                                      qs_env=qs_env, sub_env=sub_env, &
     281              :                                                      wfm_rho_orb=work_matrices%rho_ao_orb_fm_sub, &
     282              :                                                      wfm_rho_aux_fit=work_matrices%rho_ao_aux_fit_fm_sub, &
     283          806 :                                                      wfm_aux_orb=work_matrices%wfm_aux_orb_sub)
     284              :                ! - C_{HF} d^{2}E_{x, ADMM}^{DFT}[\hat{\rho}] / d\hat{\rho}^2
     285          806 :                IF (admm_symm) THEN
     286          806 :                   CALL dbcsr_get_info(rho_ia_ao_aux_fit(1)%matrix, row_blk_size=blk_sizes)
     287         3224 :                   ALLOCATE (A_xc_munu_sub(nspins))
     288         1612 :                   DO ispin = 1, nspins
     289          806 :                      ALLOCATE (A_xc_munu_sub(ispin)%matrix)
     290              :                      CALL dbcsr_create(matrix=A_xc_munu_sub(ispin)%matrix, name="ADMM_XC", &
     291              :                                        dist=sub_env%dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
     292          806 :                                        row_blk_size=blk_sizes, col_blk_size=blk_sizes)
     293          806 :                      CALL cp_dbcsr_alloc_block_from_nbl(A_xc_munu_sub(ispin)%matrix, sub_env%sab_aux_fit)
     294         1612 :                      CALL dbcsr_set(A_xc_munu_sub(ispin)%matrix, 0.0_dp)
     295              :                   END DO
     296              : 
     297          806 :                   CALL pw_env_get(sub_env%pw_env, auxbas_pw_pool=auxbas_pw_pool)
     298         3224 :                   ALLOCATE (V_rspace_sub(nspins))
     299         1612 :                   DO ispin = 1, nspins
     300          806 :                      CALL auxbas_pw_pool%create_pw(V_rspace_sub(ispin))
     301         1612 :                      CALL pw_zero(V_rspace_sub(ispin))
     302              :                   END DO
     303              : 
     304          806 :                   IF (admm_env%do_gapw) THEN
     305          120 :                      basis_type = "AUX_FIT_SOFT"
     306          120 :                      task_list => sub_env%task_list_aux_fit_soft
     307              :                   ELSE
     308          686 :                      basis_type = "AUX_FIT"
     309          686 :                      task_list => sub_env%task_list_aux_fit
     310              :                   END IF
     311              : 
     312              :                   CALL tddfpt_apply_xc(A_ia_rspace=V_rspace_sub, &
     313              :                                        kernel_env=kernel_env_admm_aux, &
     314              :                                        rho_ia_struct=work_matrices%rho_aux_fit_struct_sub, &
     315              :                                        is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
     316              :                                        work_v_xc=work_matrices%wpw_rspace_sub, &
     317          806 :                                        work_v_xc_tau=work_matrices%wpw_tau_rspace_sub)
     318         1612 :                   DO ispin = 1, nspins
     319          806 :                      CALL pw_scale(V_rspace_sub(ispin), V_rspace_sub(ispin)%pw_grid%dvol)
     320              :                      CALL integrate_v_rspace(v_rspace=V_rspace_sub(ispin), &
     321              :                                              hmat=A_xc_munu_sub(ispin), &
     322              :                                              qs_env=qs_env, calculate_forces=.FALSE., &
     323              :                                              pw_env_external=sub_env%pw_env, &
     324              :                                              basis_type=basis_type, &
     325         1612 :                                              task_list_external=task_list)
     326              :                   END DO
     327          806 :                   IF (admm_env%do_gapw) THEN
     328          120 :                      rho_atom_set => sub_env%local_rho_set_admm%rho_atom_set
     329          120 :                      rho1_atom_set => work_matrices%local_rho_set_admm%rho_atom_set
     330              :                      CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, &
     331              :                                                       kernel_env_admm_aux%xc_section, &
     332              :                                                       sub_env%para_env, do_tddfpt2=.TRUE., &
     333              :                                                       do_triplet=is_rks_triplets, &
     334          120 :                                                       kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
     335              :                      CALL update_ks_atom(qs_env, A_xc_munu_sub, rho_ia_ao_aux_fit, forces=.FALSE., tddft=.TRUE., &
     336              :                                          rho_atom_external=rho1_atom_set, &
     337              :                                          kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
     338              :                                          oce_external=admm_env%admm_gapw_env%oce, &
     339          120 :                                          sab_external=sub_env%sab_aux_fit)
     340              :                   END IF
     341          806 :                   ALLOCATE (dbwork)
     342          806 :                   CALL dbcsr_create(dbwork, template=work_matrices%A_ia_munu_sub(1)%matrix)
     343              :                   CALL cp_fm_create(work_aux_orb, &
     344          806 :                                     matrix_struct=work_matrices%wfm_aux_orb_sub%matrix_struct)
     345              :                   CALL cp_fm_create(work_orb_orb, &
     346          806 :                                     matrix_struct=work_matrices%rho_ao_orb_fm_sub%matrix_struct)
     347          806 :                   CALL cp_fm_get_info(work_aux_orb, nrow_global=nao_aux, ncol_global=nao)
     348         1612 :                   DO ispin = 1, nspins
     349              :                      CALL cp_dbcsr_sm_fm_multiply(A_xc_munu_sub(ispin)%matrix, sub_env%admm_A, &
     350          806 :                                                   work_aux_orb, nao)
     351              :                      CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, sub_env%admm_A, &
     352          806 :                                         work_aux_orb, 0.0_dp, work_orb_orb)
     353          806 :                      CALL dbcsr_copy(dbwork, work_matrices%A_ia_munu_sub(1)%matrix)
     354          806 :                      CALL dbcsr_set(dbwork, 0.0_dp)
     355          806 :                      CALL copy_fm_to_dbcsr(work_orb_orb, dbwork, keep_sparsity=.TRUE.)
     356         1612 :                      CALL dbcsr_add(work_matrices%A_ia_munu_sub(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
     357              :                   END DO
     358          806 :                   CALL dbcsr_release(dbwork)
     359          806 :                   DEALLOCATE (dbwork)
     360         1612 :                   DO ispin = 1, nspins
     361         1612 :                      CALL auxbas_pw_pool%give_back_pw(V_rspace_sub(ispin))
     362              :                   END DO
     363          806 :                   DEALLOCATE (V_rspace_sub)
     364          806 :                   CALL cp_fm_release(work_aux_orb)
     365          806 :                   CALL cp_fm_release(work_orb_orb)
     366         1612 :                   DO ispin = 1, nspins
     367         1612 :                      CALL dbcsr_deallocate_matrix(A_xc_munu_sub(ispin)%matrix)
     368              :                   END DO
     369         1612 :                   DEALLOCATE (A_xc_munu_sub)
     370              :                ELSE
     371              :                   CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, &
     372              :                                        kernel_env=kernel_env_admm_aux, &
     373              :                                        rho_ia_struct=work_matrices%rho_aux_fit_struct_sub, &
     374              :                                        is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
     375              :                                        work_v_xc=work_matrices%wpw_rspace_sub, &
     376            0 :                                        work_v_xc_tau=work_matrices%wpw_tau_rspace_sub)
     377            0 :                   IF (admm_env%do_gapw) THEN
     378            0 :                      CPWARN("GAPW/ADMM needs symmetric ADMM kernel")
     379            0 :                      CPABORT("GAPW/ADMM@TDDFT")
     380              :                   END IF
     381              :                END IF
     382              :             END IF
     383              :          END IF
     384              : 
     385              :          ! electron-hole Coulomb interaction
     386         6710 :          IF (.NOT. is_rks_triplets) THEN
     387              :             ! a sum J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu can be computed by solving
     388              :             ! the Poisson equation for combined density (rho_{ia,alpha} + rho_{ia,beta}) .
     389              :             ! The following action will destroy reciprocal-space grid in spin-unrestricted case.
     390         7544 :             DO ispin = 2, nspins
     391         7544 :                CALL pw_axpy(rho_ia_g(ispin), rho_ia_g(1))
     392              :             END DO
     393              :             CALL tddfpt_apply_coulomb(A_ia_rspace=work_matrices%A_ia_rspace_sub, &
     394              :                                       rho_ia_g=rho_ia_g(1), &
     395              :                                       local_rho_set=work_matrices%local_rho_set, &
     396              :                                       hartree_local=work_matrices%hartree_local, &
     397              :                                       qs_env=qs_env, sub_env=sub_env, gapw=gapw, &
     398              :                                       work_v_gspace=work_matrices%wpw_gspace_sub(1), &
     399              :                                       work_v_rspace=work_matrices%wpw_rspace_sub(1), &
     400         6186 :                                       tddfpt_mgrid=tddfpt_mgrid)
     401              :          END IF
     402              : 
     403              :          ! convert from the plane-wave representation into the Gaussian basis set representation
     404        14778 :          DO ispin = 1, nspins
     405        14778 :             IF (.NOT. do_lrigpw) THEN
     406              :                CALL pw_scale(work_matrices%A_ia_rspace_sub(ispin), &
     407         7896 :                              work_matrices%A_ia_rspace_sub(ispin)%pw_grid%dvol)
     408              : 
     409         7896 :                IF (gapw) THEN
     410              :                   CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), &
     411              :                                           hmat=work_matrices%A_ia_munu_sub(ispin), &
     412              :                                           qs_env=qs_env, calculate_forces=.FALSE., gapw=gapw, &
     413              :                                           pw_env_external=sub_env%pw_env, &
     414         1106 :                                           task_list_external=sub_env%task_list_orb_soft)
     415         6790 :                ELSEIF (gapw_xc) THEN
     416          218 :                   IF (.NOT. is_rks_triplets) THEN
     417              :                      CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), &
     418              :                                              hmat=work_matrices%A_ia_munu_sub(ispin), &
     419              :                                              qs_env=qs_env, calculate_forces=.FALSE., gapw=.FALSE., &
     420          218 :                                              pw_env_external=sub_env%pw_env, task_list_external=sub_env%task_list_orb)
     421              :                   END IF
     422              :                ELSE
     423              :                   CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), &
     424              :                                           hmat=work_matrices%A_ia_munu_sub(ispin), &
     425              :                                           qs_env=qs_env, calculate_forces=.FALSE., gapw=.FALSE., &
     426         6572 :                                           pw_env_external=sub_env%pw_env, task_list_external=sub_env%task_list_orb)
     427              :                END IF
     428              :             ELSE ! for full kernel using lri
     429              :                CALL pw_scale(work_matrices%A_ia_rspace_sub(ispin), &
     430          172 :                              work_matrices%A_ia_rspace_sub(ispin)%pw_grid%dvol)
     431          172 :                lri_v_int => kernel_env%lri_density%lri_coefs(ispin)%lri_kinds
     432          172 :                CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
     433          516 :                DO ikind = 1, nkind
     434       102304 :                   lri_v_int(ikind)%v_int = 0.0_dp
     435              :                END DO
     436              :                CALL integrate_v_rspace_one_center(work_matrices%A_ia_rspace_sub(ispin), &
     437          172 :                                                   qs_env, lri_v_int, .FALSE., "P_LRI_AUX")
     438          516 :                DO ikind = 1, nkind
     439       204092 :                   CALL para_env%sum(lri_v_int(ikind)%v_int)
     440              :                END DO
     441              :             END IF ! for full kernel using lri
     442              :          END DO
     443              : 
     444              :          ! local atom contributions
     445         6710 :          IF (.NOT. do_lrigpw) THEN
     446         6538 :             IF (gapw .OR. gapw_xc) THEN
     447              :                ! rho_ia_ao will not be touched
     448              :                CALL update_ks_atom(qs_env, work_matrices%A_ia_munu_sub, rho_ia_ao, forces=.FALSE., &
     449              :                                    rho_atom_external=work_matrices%local_rho_set%rho_atom_set, &
     450         1324 :                                    tddft=.TRUE.)
     451              :             END IF
     452              :          END IF
     453              : 
     454              :          ! calculate Coulomb contribution to response vector for lrigpw !
     455              :          ! this is restricting lri to Coulomb only at the moment !
     456         6710 :          IF (do_lrigpw .AND. (.NOT. is_rks_triplets)) THEN !
     457          172 :             CALL tddfpt2_lri_Amat(qs_env, sub_env, kernel_env%lri_env, lri_v_int, work_matrices%A_ia_munu_sub)
     458              :          END IF
     459              : 
     460         9608 :          IF (ALLOCATED(work_matrices%evects_sub)) THEN
     461           16 :             DO ispin = 1, nspins
     462              :                CALL cp_dbcsr_sm_fm_multiply(work_matrices%A_ia_munu_sub(ispin)%matrix, &
     463              :                                             sub_env%mos_occ(ispin), &
     464              :                                             work_matrices%Aop_evects_sub(ispin, ivect), &
     465           16 :                                             ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
     466              :             END DO
     467              :          ELSE
     468        14762 :             DO ispin = 1, nspins
     469              :                CALL cp_dbcsr_sm_fm_multiply(work_matrices%A_ia_munu_sub(ispin)%matrix, &
     470              :                                             sub_env%mos_occ(ispin), &
     471              :                                             Aop_evects(ispin, ivect), &
     472        14770 :                                             ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
     473              :             END DO
     474              :          END IF
     475              :       END DO
     476              : 
     477         2898 :       CALL timestop(handle)
     478              : 
     479         5796 :    END SUBROUTINE fhxc_kernel
     480              : 
     481              : ! **************************************************************************************************
     482              : !> \brief Compute action matrix-vector products with the sTDA Kernel
     483              : !> \param Aop_evects            action of TDDFPT operator on trial vectors (modified on exit)
     484              : !> \param evects                TDDFPT trial vectors
     485              : !> \param is_rks_triplets       indicates that a triplet excited states calculation using
     486              : !>                              spin-unpolarised molecular orbitals has been requested
     487              : !> \param qs_env                Quickstep environment
     488              : !> \param stda_control          control parameters for sTDA kernel
     489              : !> \param stda_env ...
     490              : !> \param sub_env               parallel (sub)group environment
     491              : !> \param work_matrices         collection of work matrices (modified on exit)
     492              : !> \par History
     493              : !>    * 04.2019 initial version [JHU]
     494              : ! **************************************************************************************************
     495         2264 :    SUBROUTINE stda_kernel(Aop_evects, evects, is_rks_triplets, &
     496              :                           qs_env, stda_control, stda_env, &
     497              :                           sub_env, work_matrices)
     498              : 
     499              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT)   :: Aop_evects
     500              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: evects
     501              :       LOGICAL, INTENT(in)                                :: is_rks_triplets
     502              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     503              :       TYPE(stda_control_type)                            :: stda_control
     504              :       TYPE(stda_env_type)                                :: stda_env
     505              :       TYPE(tddfpt_subgroup_env_type)                     :: sub_env
     506              :       TYPE(tddfpt_work_matrices), INTENT(inout)          :: work_matrices
     507              : 
     508              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'stda_kernel'
     509              : 
     510              :       INTEGER                                            :: handle, ivect, nvects
     511              : 
     512         2264 :       CALL timeset(routineN, handle)
     513              : 
     514         2264 :       nvects = SIZE(evects, 2)
     515              : 
     516         8726 :       DO ivect = 1, nvects
     517         8726 :          IF (ALLOCATED(work_matrices%evects_sub)) THEN
     518            0 :             IF (ASSOCIATED(work_matrices%evects_sub(1, ivect)%matrix_struct)) THEN
     519              :                CALL stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, work_matrices, &
     520              :                                           is_rks_triplets, work_matrices%evects_sub(:, ivect), &
     521            0 :                                           work_matrices%Aop_evects_sub(:, ivect))
     522              :             ELSE
     523              :                ! skip trial vectors which are assigned to different parallel groups
     524              :                CYCLE
     525              :             END IF
     526              :          ELSE
     527              :             CALL stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, work_matrices, &
     528         6462 :                                        is_rks_triplets, evects(:, ivect), Aop_evects(:, ivect))
     529              :          END IF
     530              :       END DO
     531              : 
     532         2264 :       CALL timestop(handle)
     533              : 
     534         2264 :    END SUBROUTINE stda_kernel
     535              : 
     536              : ! **************************************************************************************************
     537              : 
     538              : END MODULE qs_tddfpt2_fhxc
        

Generated by: LCOV version 2.0-1