LCOV - code coverage report
Current view: top level - src - rpa_gw_sigma_x.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 188 407 46.2 %
Date: 2024-03-28 07:31:50 Functions: 1 4 25.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Routines to calculate EXX within GW
      10             : !> \par History
      11             : !>      07.2020 separated from mp2.F [F. Stein, code by Jan Wilhelm]
      12             : !> \author Jan Wilhelm, Frederick Stein
      13             : ! **************************************************************************************************
      14             : MODULE rpa_gw_sigma_x
      15             :    USE admm_methods,                    ONLY: admm_mo_merge_ks_matrix
      16             :    USE admm_types,                      ONLY: admm_type,&
      17             :                                               get_admm_env
      18             :    USE cp_cfm_basic_linalg,             ONLY: cp_cfm_scale_and_add_fm
      19             :    USE cp_cfm_types,                    ONLY: cp_cfm_create,&
      20             :                                               cp_cfm_get_info,&
      21             :                                               cp_cfm_release,&
      22             :                                               cp_cfm_type
      23             :    USE cp_control_types,                ONLY: dft_control_type
      24             :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      25             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      26             :                                               copy_fm_to_dbcsr,&
      27             :                                               dbcsr_allocate_matrix_set,&
      28             :                                               dbcsr_deallocate_matrix_set
      29             :    USE cp_files,                        ONLY: close_file,&
      30             :                                               open_file
      31             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_type
      32             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      33             :                                               cp_fm_get_info,&
      34             :                                               cp_fm_release,&
      35             :                                               cp_fm_type
      36             :    USE dbcsr_api,                       ONLY: &
      37             :         dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_get_diag, dbcsr_multiply, &
      38             :         dbcsr_p_type, dbcsr_release, dbcsr_release_p, dbcsr_set, dbcsr_type, &
      39             :         dbcsr_type_antisymmetric, dbcsr_type_symmetric
      40             :    USE hfx_energy_potential,            ONLY: integrate_four_center
      41             :    USE hfx_exx,                         ONLY: calc_exx_admm_xc_contributions,&
      42             :                                               exx_post_hfx,&
      43             :                                               exx_pre_hfx
      44             :    USE hfx_ri,                          ONLY: hfx_ri_update_ks
      45             :    USE input_constants,                 ONLY: do_admm_basis_projection,&
      46             :                                               do_admm_purify_none,&
      47             :                                               gw_print_exx,&
      48             :                                               gw_read_exx,&
      49             :                                               xc_none
      50             :    USE input_section_types,             ONLY: section_vals_get,&
      51             :                                               section_vals_get_subs_vals,&
      52             :                                               section_vals_type,&
      53             :                                               section_vals_val_get,&
      54             :                                               section_vals_val_set
      55             :    USE kinds,                           ONLY: dp
      56             :    USE kpoint_methods,                  ONLY: rskp_transform
      57             :    USE kpoint_types,                    ONLY: get_kpoint_info,&
      58             :                                               kpoint_env_type,&
      59             :                                               kpoint_type
      60             :    USE machine,                         ONLY: m_walltime
      61             :    USE mathconstants,                   ONLY: gaussi,&
      62             :                                               z_one,&
      63             :                                               z_zero
      64             :    USE message_passing,                 ONLY: mp_para_env_type
      65             :    USE mp2_integrals,                   ONLY: compute_kpoints
      66             :    USE mp2_ri_2c,                       ONLY: setup_trunc_coulomb_pot_for_exchange_self_energy
      67             :    USE mp2_types,                       ONLY: mp2_type
      68             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      69             :    USE physcon,                         ONLY: evolt
      70             :    USE qs_energy_types,                 ONLY: qs_energy_type
      71             :    USE qs_environment_types,            ONLY: get_qs_env,&
      72             :                                               qs_environment_type
      73             :    USE qs_ks_methods,                   ONLY: qs_ks_build_kohn_sham_matrix
      74             :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      75             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      76             :                                               mo_set_type
      77             :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
      78             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      79             :                                               qs_rho_type
      80             :    USE rpa_gw,                          ONLY: compute_minus_vxc_kpoints,&
      81             :                                               trafo_to_mo_and_kpoints
      82             :    USE rpa_gw_kpoints_util,             ONLY: get_bandstruc_and_k_dependent_MOs
      83             : 
      84             : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
      85             : 
      86             : #include "./base/base_uses.f90"
      87             : 
      88             :    IMPLICIT NONE
      89             : 
      90             :    PRIVATE
      91             : 
      92             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_gw_sigma_x'
      93             : 
      94             :    PUBLIC :: compute_vec_Sigma_x_minus_vxc_gw
      95             : 
      96             : CONTAINS
      97             : 
      98             : ! **************************************************************************************************
      99             : !> \brief ...
     100             : !> \param qs_env ...
     101             : !> \param mp2_env ...
     102             : !> \param mos_mp2 ...
     103             : !> \param energy_ex ...
     104             : !> \param energy_xc_admm ...
     105             : !> \param t3 ...
     106             : !> \param unit_nr ...
     107             : !> \par History
     108             : !>      04.2015 created
     109             : !> \author Jan Wilhelm
     110             : ! **************************************************************************************************
     111          72 :    SUBROUTINE compute_vec_Sigma_x_minus_vxc_gw(qs_env, mp2_env, mos_mp2, energy_ex, energy_xc_admm, t3, unit_nr)
     112             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     113             :       TYPE(mp2_type)                                     :: mp2_env
     114             :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mos_mp2
     115             :       REAL(KIND=dp), INTENT(OUT)                         :: energy_ex, energy_xc_admm(2), t3
     116             :       INTEGER, INTENT(IN)                                :: unit_nr
     117             : 
     118             :       CHARACTER(len=*), PARAMETER :: routineN = 'compute_vec_Sigma_x_minus_vxc_gw'
     119             : 
     120             :       CHARACTER(4)                                       :: occ_virt
     121             :       CHARACTER(LEN=40)                                  :: line
     122             :       INTEGER :: dimen, gw_corr_lev_occ, gw_corr_lev_tot, gw_corr_lev_virt, handle, homo, i_img, &
     123             :          ikp, irep, ispin, iunit, myfun, myfun_aux, myfun_prim, n_level_gw, n_level_gw_ref, &
     124             :          n_rep_hf, nkp, nkp_Sigma, nmo, nspins, print_exx
     125             :       LOGICAL :: calc_ints, charge_constrain_tmp, do_admm_rpa, do_hfx, do_kpoints_cubic_RPA, &
     126             :          do_kpoints_from_Gamma, do_ri_Sigma_x, really_read_line
     127             :       REAL(KIND=dp) :: E_GAP_GW, E_HOMO_GW, E_LUMO_GW, eh1, ehfx, eigval_dft, eigval_hf_at_dft, &
     128             :          energy_exc, energy_exc1, energy_exc1_aux_fit, energy_exc_aux_fit, energy_total, &
     129             :          exx_minus_vxc, hfx_fraction, min_direct_HF_at_DFT_gap, t1, t2, tmp
     130          72 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: Eigenval_kp_HF_at_DFT, vec_Sigma_x
     131          72 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: Eigenval_kp, vec_Sigma_x_minus_vxc_gw, &
     132          72 :                                                             vec_Sigma_x_minus_vxc_gw_im
     133             :       TYPE(admm_type), POINTER                           :: admm_env
     134             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     135          72 :       TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:)      :: mat_exchange_for_kp_from_gamma
     136          72 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_ks_aux_fit, &
     137          72 :                                                             matrix_ks_aux_fit_hfx, rho_ao, &
     138          72 :                                                             rho_ao_aux_fit
     139          72 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_2d, matrix_ks_kp_im, &
     140          72 :          matrix_ks_kp_re, matrix_ks_transl, matrix_sigma_x_minus_vxc, matrix_sigma_x_minus_vxc_im, &
     141          72 :          rho_ao_2d
     142             :       TYPE(dbcsr_type)                                   :: matrix_tmp, matrix_tmp_2, mo_coeff_b
     143             :       TYPE(dft_control_type), POINTER                    :: dft_control
     144             :       TYPE(kpoint_type), POINTER                         :: kpoints, kpoints_Sigma
     145             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     146             :       TYPE(qs_energy_type), POINTER                      :: energy
     147             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     148             :       TYPE(qs_rho_type), POINTER                         :: rho, rho_aux_fit
     149             :       TYPE(section_vals_type), POINTER                   :: hfx_sections, input, xc_section, &
     150             :                                                             xc_section_admm_aux, &
     151             :                                                             xc_section_admm_prim
     152             : 
     153          72 :       NULLIFY (admm_env, matrix_ks, matrix_ks_aux_fit, rho_ao, matrix_sigma_x_minus_vxc, input, &
     154          72 :                xc_section, xc_section_admm_aux, xc_section_admm_prim, hfx_sections, rho, &
     155          72 :                dft_control, para_env, ks_env, mo_coeff, matrix_sigma_x_minus_vxc_im, matrix_ks_aux_fit_hfx, &
     156          72 :                rho_aux_fit, rho_ao_aux_fit)
     157             : 
     158          72 :       CALL timeset(routineN, handle)
     159             : 
     160          72 :       t1 = m_walltime()
     161             : 
     162          72 :       do_admm_rpa = mp2_env%ri_rpa%do_admm
     163          72 :       do_ri_Sigma_x = mp2_env%ri_g0w0%do_ri_Sigma_x
     164          72 :       do_kpoints_cubic_RPA = qs_env%mp2_env%ri_rpa_im_time%do_im_time_kpoints
     165          72 :       do_kpoints_from_Gamma = qs_env%mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma
     166          72 :       print_exx = mp2_env%ri_g0w0%print_exx
     167             : 
     168          72 :       IF (do_kpoints_cubic_RPA) THEN
     169           0 :          CPASSERT(do_ri_Sigma_x)
     170             :       END IF
     171             : 
     172          72 :       IF (do_kpoints_cubic_RPA) THEN
     173             : 
     174             :          CALL get_qs_env(qs_env, &
     175             :                          admm_env=admm_env, &
     176             :                          matrix_ks_kp=matrix_ks_transl, &
     177             :                          rho=rho, &
     178             :                          input=input, &
     179             :                          dft_control=dft_control, &
     180             :                          para_env=para_env, &
     181             :                          kpoints=kpoints, &
     182             :                          ks_env=ks_env, &
     183           0 :                          energy=energy)
     184           0 :          nkp = kpoints%nkp
     185             : 
     186             :       ELSE
     187             : 
     188             :          CALL get_qs_env(qs_env, &
     189             :                          admm_env=admm_env, &
     190             :                          matrix_ks=matrix_ks, &
     191             :                          rho=rho, &
     192             :                          input=input, &
     193             :                          dft_control=dft_control, &
     194             :                          para_env=para_env, &
     195             :                          ks_env=ks_env, &
     196          72 :                          energy=energy)
     197          72 :          nkp = 1
     198          72 :          CALL qs_rho_get(rho, rho_ao=rho_ao)
     199             : 
     200          72 :          IF (do_admm_rpa) THEN
     201             :             CALL get_admm_env(admm_env, matrix_ks_aux_fit=matrix_ks_aux_fit, rho_aux_fit=rho_aux_fit, &
     202           8 :                               matrix_ks_aux_fit_hfx=matrix_ks_aux_fit_hfx)
     203           8 :             CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux_fit)
     204             : 
     205             :             ! RPA/GW with ADMM for EXX or the exchange self-energy only implemented for
     206             :             ! ADMM_PURIFICATION_METHOD  NONE
     207             :             ! METHOD                    BASIS_PROJECTION
     208             :             ! in the admm section
     209           8 :             CPASSERT(admm_env%purification_method == do_admm_purify_none)
     210           8 :             CPASSERT(dft_control%admm_control%method == do_admm_basis_projection)
     211             :          END IF
     212             :       END IF
     213             : 
     214          72 :       nspins = dft_control%nspins
     215             : 
     216             :       ! safe ks matrix for later: we will transform matrix_ks
     217             :       ! to T-cell index and then to k-points for band structure calculation
     218          72 :       IF (do_kpoints_from_Gamma) THEN
     219             :          ! not yet there: open shell
     220          76 :          ALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_ks(nspins))
     221          40 :          DO ispin = 1, nspins
     222          22 :             NULLIFY (qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix)
     223          22 :             ALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix)
     224             :             CALL dbcsr_create(qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix, &
     225          22 :                               template=matrix_ks(ispin)%matrix)
     226             :             CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, &
     227          40 :                                     qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix)
     228             : 
     229             :          END DO
     230             :       END IF
     231             : 
     232          72 :       IF (do_kpoints_cubic_RPA) THEN
     233             : 
     234           0 :          CALL allocate_matrix_ks_kp(matrix_ks_transl, matrix_ks_kp_re, matrix_ks_kp_im, kpoints)
     235           0 :          CALL transform_matrix_ks_to_kp(matrix_ks_transl, matrix_ks_kp_re, matrix_ks_kp_im, kpoints)
     236             : 
     237           0 :          DO ispin = 1, nspins
     238           0 :          DO i_img = 1, SIZE(matrix_ks_transl, 2)
     239           0 :             CALL dbcsr_set(matrix_ks_transl(ispin, i_img)%matrix, 0.0_dp)
     240             :          END DO
     241             :          END DO
     242             : 
     243             :       END IF
     244             : 
     245             :       ! initialize matrix_sigma_x_minus_vxc
     246          72 :       NULLIFY (matrix_sigma_x_minus_vxc)
     247          72 :       CALL dbcsr_allocate_matrix_set(matrix_sigma_x_minus_vxc, nspins, nkp)
     248          72 :       IF (do_kpoints_cubic_RPA) THEN
     249           0 :          NULLIFY (matrix_sigma_x_minus_vxc_im)
     250           0 :          CALL dbcsr_allocate_matrix_set(matrix_sigma_x_minus_vxc_im, nspins, nkp)
     251             :       END IF
     252             : 
     253         158 :       DO ispin = 1, nspins
     254         244 :          DO ikp = 1, nkp
     255             : 
     256         172 :             IF (do_kpoints_cubic_RPA) THEN
     257             : 
     258           0 :                ALLOCATE (matrix_sigma_x_minus_vxc(ispin, ikp)%matrix)
     259             :                CALL dbcsr_create(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, &
     260             :                                  template=matrix_ks_kp_re(1, 1)%matrix, &
     261           0 :                                  matrix_type=dbcsr_type_symmetric)
     262             : 
     263           0 :                CALL dbcsr_copy(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, matrix_ks_kp_re(ispin, ikp)%matrix)
     264           0 :                CALL dbcsr_set(matrix_ks_kp_re(ispin, ikp)%matrix, 0.0_dp)
     265             : 
     266           0 :                ALLOCATE (matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix)
     267             :                CALL dbcsr_create(matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix, &
     268             :                                  template=matrix_ks_kp_im(1, 1)%matrix, &
     269           0 :                                  matrix_type=dbcsr_type_antisymmetric)
     270             : 
     271           0 :                CALL dbcsr_copy(matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix, matrix_ks_kp_im(ispin, ikp)%matrix)
     272           0 :                CALL dbcsr_set(matrix_ks_kp_im(ispin, ikp)%matrix, 0.0_dp)
     273             : 
     274             :             ELSE
     275             : 
     276          86 :                ALLOCATE (matrix_sigma_x_minus_vxc(ispin, ikp)%matrix)
     277             :                CALL dbcsr_create(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, &
     278          86 :                                  template=matrix_ks(1)%matrix)
     279             : 
     280          86 :                CALL dbcsr_copy(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, matrix_ks(ispin)%matrix)
     281          86 :                CALL dbcsr_set(matrix_ks(ispin)%matrix, 0.0_dp)
     282             : 
     283             :             END IF
     284             : 
     285             :          END DO
     286             :       END DO
     287             : 
     288             :       ! set DFT functional to none and hfx_fraction to zero
     289          72 :       hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
     290          72 :       CALL section_vals_get(hfx_sections, explicit=do_hfx)
     291             : 
     292          72 :       IF (do_hfx) THEN
     293          16 :          hfx_fraction = qs_env%x_data(1, 1)%general_parameter%fraction
     294          48 :          qs_env%x_data(:, :)%general_parameter%fraction = 0.0_dp
     295             :       END IF
     296          72 :       xc_section => section_vals_get_subs_vals(input, "DFT%XC")
     297             :       CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
     298          72 :                                 i_val=myfun)
     299             :       CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
     300          72 :                                 i_val=xc_none)
     301             : 
     302             :       ! in ADMM, also set the XC functional for ADMM correction to none
     303             :       ! do not do this if we do ADMM for Sigma_x
     304          72 :       IF (dft_control%do_admm) THEN
     305             :          xc_section_admm_aux => section_vals_get_subs_vals(admm_env%xc_section_aux, &
     306           8 :                                                            "XC_FUNCTIONAL")
     307             :          CALL section_vals_val_get(xc_section_admm_aux, "_SECTION_PARAMETERS_", &
     308           8 :                                    i_val=myfun_aux)
     309             :          CALL section_vals_val_set(xc_section_admm_aux, "_SECTION_PARAMETERS_", &
     310           8 :                                    i_val=xc_none)
     311             : 
     312             :          ! the same for the primary basis
     313             :          xc_section_admm_prim => section_vals_get_subs_vals(admm_env%xc_section_primary, &
     314           8 :                                                             "XC_FUNCTIONAL")
     315             :          CALL section_vals_val_get(xc_section_admm_prim, "_SECTION_PARAMETERS_", &
     316           8 :                                    i_val=myfun_prim)
     317             :          CALL section_vals_val_set(xc_section_admm_prim, "_SECTION_PARAMETERS_", &
     318           8 :                                    i_val=xc_none)
     319             : 
     320             :          ! for ADMMQ/S, set the charge_constrain to false (otherwise wrong results)
     321           8 :          charge_constrain_tmp = .FALSE.
     322           8 :          IF (admm_env%charge_constrain) THEN
     323           0 :             admm_env%charge_constrain = .FALSE.
     324           0 :             charge_constrain_tmp = .TRUE.
     325             :          END IF
     326             : 
     327             :       END IF
     328             : 
     329             :       ! if we do ADMM for Sigma_x, we write the ADMM correction into matrix_ks_aux_fit
     330             :       ! and therefore we should set it to zero
     331          72 :       IF (do_admm_rpa) THEN
     332          18 :          DO ispin = 1, nspins
     333          18 :             CALL dbcsr_set(matrix_ks_aux_fit(ispin)%matrix, 0.0_dp)
     334             :          END DO
     335             :       END IF
     336             : 
     337          72 :       IF (.NOT. mp2_env%ri_g0w0%update_xc_energy) THEN
     338          46 :          energy_total = energy%total
     339          46 :          energy_exc = energy%exc
     340          46 :          energy_exc1 = energy%exc1
     341          46 :          energy_exc_aux_fit = energy%ex
     342          46 :          energy_exc1_aux_fit = energy%exc_aux_fit
     343          46 :          energy_ex = energy%exc1_aux_fit
     344             :       END IF
     345             : 
     346             :       ! Remove the Exchange-correlation energy contributions from the total energy
     347             :       energy%total = energy%total - (energy%exc + energy%exc1 + energy%ex + &
     348          72 :                                      energy%exc_aux_fit + energy%exc1_aux_fit)
     349             : 
     350             :       ! calculate KS-matrix without XC and without HF
     351             :       CALL qs_ks_build_kohn_sham_matrix(qs_env=qs_env, calculate_forces=.FALSE., &
     352          72 :                                         just_energy=.FALSE.)
     353             : 
     354          72 :       IF (.NOT. mp2_env%ri_g0w0%update_xc_energy) THEN
     355          46 :          energy%exc = energy_exc
     356          46 :          energy%exc1 = energy_exc1
     357          46 :          energy%exc_aux_fit = energy_ex
     358          46 :          energy%exc1_aux_fit = energy_exc_aux_fit
     359          46 :          energy%ex = energy_exc1_aux_fit
     360          46 :          energy%total = energy_total
     361             :       END IF
     362             : 
     363             :       ! set the DFT functional and HF fraction back
     364             :       CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
     365          72 :                                 i_val=myfun)
     366          72 :       IF (do_hfx) THEN
     367          48 :          qs_env%x_data(:, :)%general_parameter%fraction = hfx_fraction
     368             :       END IF
     369             : 
     370          72 :       IF (dft_control%do_admm) THEN
     371             :          xc_section_admm_aux => section_vals_get_subs_vals(admm_env%xc_section_aux, &
     372           8 :                                                            "XC_FUNCTIONAL")
     373             :          xc_section_admm_prim => section_vals_get_subs_vals(admm_env%xc_section_primary, &
     374           8 :                                                             "XC_FUNCTIONAL")
     375             : 
     376             :          CALL section_vals_val_set(xc_section_admm_aux, "_SECTION_PARAMETERS_", &
     377           8 :                                    i_val=myfun_aux)
     378             :          CALL section_vals_val_set(xc_section_admm_prim, "_SECTION_PARAMETERS_", &
     379           8 :                                    i_val=myfun_prim)
     380           8 :          IF (charge_constrain_tmp) THEN
     381           0 :             admm_env%charge_constrain = .TRUE.
     382             :          END IF
     383             :       END IF
     384             : 
     385          72 :       IF (do_kpoints_cubic_RPA) THEN
     386           0 :          CALL transform_matrix_ks_to_kp(matrix_ks_transl, matrix_ks_kp_re, matrix_ks_kp_im, kpoints)
     387             :       END IF
     388             : 
     389             :       ! remove the single-particle part (kin. En + Hartree pot) and change the sign
     390         158 :       DO ispin = 1, nspins
     391         158 :          IF (do_kpoints_cubic_RPA) THEN
     392           0 :             DO ikp = 1, nkp
     393           0 :                CALL dbcsr_add(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, matrix_ks_kp_re(ispin, ikp)%matrix, -1.0_dp, 1.0_dp)
     394           0 :                CALL dbcsr_add(matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix, matrix_ks_kp_im(ispin, ikp)%matrix, -1.0_dp, 1.0_dp)
     395             :             END DO
     396             :          ELSE
     397          86 :             CALL dbcsr_add(matrix_sigma_x_minus_vxc(ispin, 1)%matrix, matrix_ks(ispin)%matrix, -1.0_dp, 1.0_dp)
     398             :          END IF
     399             :       END DO
     400             : 
     401          72 :       IF (do_kpoints_cubic_RPA) THEN
     402             : 
     403             :          CALL transform_sigma_x_minus_vxc_to_MO_basis(kpoints, matrix_sigma_x_minus_vxc, &
     404             :                                                       matrix_sigma_x_minus_vxc_im, &
     405             :                                                       vec_Sigma_x_minus_vxc_gw, &
     406             :                                                       vec_Sigma_x_minus_vxc_gw_im, &
     407           0 :                                                       para_env, nmo, mp2_env)
     408             : 
     409             :       ELSE
     410             : 
     411         158 :          DO ispin = 1, nspins
     412          86 :             CALL dbcsr_set(matrix_ks(ispin)%matrix, 0.0_dp)
     413         158 :             IF (do_admm_rpa) THEN
     414          10 :                CALL dbcsr_set(matrix_ks_aux_fit(ispin)%matrix, 0.0_dp)
     415             :             END IF
     416             :          END DO
     417             : 
     418          72 :          hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
     419             : 
     420          72 :          CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
     421             : 
     422             :          ! in most cases, we calculate the exchange self-energy here. But if we do only RI for
     423             :          ! the exchange self-energy, we do not calculate exchange here
     424          72 :          ehfx = 0.0_dp
     425          72 :          IF (.NOT. do_ri_Sigma_x) THEN
     426             : 
     427          46 :             CALL exx_pre_hfx(hfx_sections, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx)
     428          46 :             calc_ints = .NOT. qs_env%mp2_env%ri_rpa%reuse_hfx
     429             : 
     430             :             ! add here HFX (=Sigma_exchange) to matrix_sigma_x_minus_vxc
     431          92 :             DO irep = 1, n_rep_hf
     432          46 :                IF (do_admm_rpa) THEN
     433           8 :                   matrix_ks_2d(1:nspins, 1:1) => matrix_ks_aux_fit(1:nspins)
     434           8 :                   rho_ao_2d(1:nspins, 1:1) => rho_ao_aux_fit(1:nspins)
     435             :                ELSE
     436          38 :                   matrix_ks_2d(1:nspins, 1:1) => matrix_ks(1:nspins)
     437          38 :                   rho_ao_2d(1:nspins, 1:1) => rho_ao(1:nspins)
     438             :                END IF
     439             : 
     440          92 :                IF (qs_env%mp2_env%ri_rpa%x_data(irep, 1)%do_hfx_ri) THEN
     441             :                   CALL hfx_ri_update_ks(qs_env, qs_env%mp2_env%ri_rpa%x_data(irep, 1)%ri_data, matrix_ks_2d, ehfx, &
     442             :                                         rho_ao=rho_ao_2d, geometry_did_change=calc_ints, nspins=nspins, &
     443           0 :                                         hf_fraction=qs_env%mp2_env%ri_rpa%x_data(irep, 1)%general_parameter%fraction)
     444             : 
     445           0 :                   IF (do_admm_rpa) THEN
     446             :                      !for ADMMS, we need the exchange matrix k(d) for both spins
     447           0 :                      DO ispin = 1, nspins
     448             :                         CALL dbcsr_copy(matrix_ks_aux_fit_hfx(ispin)%matrix, matrix_ks_2d(ispin, 1)%matrix, &
     449           0 :                                         name="HF exch. part of matrix_ks_aux_fit for ADMMS")
     450             :                      END DO
     451             :                   END IF
     452             :                ELSE
     453             :                   CALL integrate_four_center(qs_env, qs_env%mp2_env%ri_rpa%x_data, matrix_ks_2d, eh1, &
     454             :                                              rho_ao_2d, hfx_sections, &
     455             :                                              para_env, calc_ints, irep, .TRUE., &
     456          46 :                                              ispin=1)
     457          46 :                   ehfx = ehfx + eh1
     458             :                END IF
     459             :             END DO
     460             : 
     461             :             !ADMM XC correction
     462          46 :             IF (do_admm_rpa) THEN
     463             :                CALL calc_exx_admm_xc_contributions(qs_env=qs_env, &
     464             :                                                    matrix_prim=matrix_ks, &
     465             :                                                    matrix_aux=matrix_ks_aux_fit, &
     466             :                                                    x_data=qs_env%mp2_env%ri_rpa%x_data, &
     467             :                                                    exc=energy_xc_admm(1), &
     468             :                                                    exc_aux_fit=energy_xc_admm(2), &
     469             :                                                    calc_forces=.FALSE., &
     470           8 :                                                    use_virial=.FALSE.)
     471             :             END IF
     472             : 
     473          46 :             IF (do_kpoints_from_Gamma .AND. print_exx == gw_print_exx) THEN
     474             :                ! JW not yet there: open shell
     475           0 :                ALLOCATE (mat_exchange_for_kp_from_gamma(1))
     476             : 
     477           0 :                DO ispin = 1, 1
     478           0 :                   NULLIFY (mat_exchange_for_kp_from_gamma(ispin)%matrix)
     479           0 :                   ALLOCATE (mat_exchange_for_kp_from_gamma(ispin)%matrix)
     480           0 :                   CALL dbcsr_create(mat_exchange_for_kp_from_gamma(ispin)%matrix, template=matrix_ks(ispin)%matrix)
     481           0 :                   CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, mat_exchange_for_kp_from_gamma(ispin)%matrix)
     482             :                END DO
     483             : 
     484             :             END IF
     485             : 
     486          46 :             CALL exx_post_hfx(qs_env, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx)
     487             :          END IF
     488             : 
     489          72 :          energy_ex = ehfx
     490             : 
     491             :          ! transform Fock-Matrix (calculated in integrate_four_center, written in matrix_ks_aux_fit in case
     492             :          ! of ADMM) from ADMM basis to primary basis
     493          72 :          IF (do_admm_rpa) THEN
     494           8 :             CALL admm_mo_merge_ks_matrix(qs_env)
     495             :          END IF
     496             : 
     497         158 :          DO ispin = 1, nspins
     498         158 :             CALL dbcsr_add(matrix_sigma_x_minus_vxc(ispin, 1)%matrix, matrix_ks(ispin)%matrix, 1.0_dp, 1.0_dp)
     499             :          END DO
     500             : 
     501             :          ! safe matrix_sigma_x_minus_vxc for later: for example, we will transform matrix_sigma_x_minus_vxc
     502             :          ! to T-cell index and then to k-points for band structure calculation
     503          72 :          IF (do_kpoints_from_Gamma) THEN
     504             :             ! not yet there: open shell
     505          76 :             ALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(nspins))
     506          40 :             DO ispin = 1, nspins
     507          22 :                NULLIFY (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix)
     508          22 :                ALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix)
     509             :                CALL dbcsr_create(qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix, &
     510          22 :                                  template=matrix_ks(ispin)%matrix)
     511             : 
     512             :                CALL dbcsr_desymmetrize(matrix_sigma_x_minus_vxc(ispin, 1)%matrix, &
     513          40 :                                        qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix)
     514             : 
     515             :             END DO
     516             :          END IF
     517             : 
     518          72 :          CALL dbcsr_desymmetrize(matrix_ks(1)%matrix, mo_coeff_b)
     519          72 :          CALL dbcsr_set(mo_coeff_b, 0.0_dp)
     520             : 
     521             :          ! Transform matrix_sigma_x_minus_vxc to MO basis
     522         158 :          DO ispin = 1, nspins
     523             : 
     524             :             CALL get_mo_set(mo_set=mos_mp2(ispin), &
     525             :                             mo_coeff=mo_coeff, &
     526             :                             nmo=nmo, &
     527             :                             homo=homo, &
     528          86 :                             nao=dimen)
     529             : 
     530          86 :             IF (ispin == 1) THEN
     531             : 
     532         360 :                ALLOCATE (vec_Sigma_x_minus_vxc_gw(nmo, nspins, nkp))
     533        2276 :                vec_Sigma_x_minus_vxc_gw = 0.0_dp
     534             :             END IF
     535             : 
     536          86 :             CALL dbcsr_set(mo_coeff_b, 0.0_dp)
     537          86 :             CALL copy_fm_to_dbcsr(mo_coeff, mo_coeff_b, keep_sparsity=.FALSE.)
     538             : 
     539             :             ! initialize matrix_tmp and matrix_tmp2
     540          86 :             IF (ispin == 1) THEN
     541          72 :                CALL dbcsr_create(matrix_tmp, template=mo_coeff_b)
     542          72 :                CALL dbcsr_copy(matrix_tmp, mo_coeff_b)
     543          72 :                CALL dbcsr_set(matrix_tmp, 0.0_dp)
     544             : 
     545          72 :                CALL dbcsr_create(matrix_tmp_2, template=mo_coeff_b)
     546          72 :                CALL dbcsr_copy(matrix_tmp_2, mo_coeff_b)
     547          72 :                CALL dbcsr_set(matrix_tmp_2, 0.0_dp)
     548             :             END IF
     549             : 
     550          86 :             gw_corr_lev_occ = mp2_env%ri_g0w0%corr_mos_occ
     551          86 :             gw_corr_lev_virt = mp2_env%ri_g0w0%corr_mos_virt
     552             :             ! if requested number of occ/virt levels for correction either exceed the number of
     553             :             ! occ/virt levels or the requested number is negative, default to correct all
     554             :             ! occ/virt level energies
     555          86 :             IF (gw_corr_lev_occ > homo .OR. gw_corr_lev_occ < 0) gw_corr_lev_occ = homo
     556          86 :             IF (gw_corr_lev_virt > dimen - homo .OR. gw_corr_lev_virt < 0) gw_corr_lev_virt = dimen - homo
     557          86 :             IF (ispin == 1) THEN
     558          72 :                mp2_env%ri_g0w0%corr_mos_occ = gw_corr_lev_occ
     559          72 :                mp2_env%ri_g0w0%corr_mos_virt = gw_corr_lev_virt
     560          14 :             ELSE IF (ispin == 2) THEN
     561             :                ! ensure that the total number of corrected MOs is the same for alpha and beta, important
     562             :                ! for parallelization
     563          14 :                IF (mp2_env%ri_g0w0%corr_mos_occ + mp2_env%ri_g0w0%corr_mos_virt /= &
     564             :                    gw_corr_lev_occ + gw_corr_lev_virt) THEN
     565          10 :                   gw_corr_lev_virt = mp2_env%ri_g0w0%corr_mos_occ + mp2_env%ri_g0w0%corr_mos_virt - gw_corr_lev_occ
     566             :                END IF
     567          14 :                mp2_env%ri_g0w0%corr_mos_occ_beta = gw_corr_lev_occ
     568          14 :                mp2_env%ri_g0w0%corr_mos_virt_beta = gw_corr_lev_virt
     569             : 
     570             :             END IF
     571             : 
     572             :             CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_sigma_x_minus_vxc(ispin, 1)%matrix, &
     573             :                                 mo_coeff_b, 0.0_dp, matrix_tmp, first_column=homo + 1 - gw_corr_lev_occ, &
     574          86 :                                 last_column=homo + gw_corr_lev_virt)
     575             : 
     576             :             CALL dbcsr_multiply('T', 'N', 1.0_dp, mo_coeff_b, &
     577             :                                 matrix_tmp, 0.0_dp, matrix_tmp_2, first_row=homo + 1 - gw_corr_lev_occ, &
     578          86 :                                 last_row=homo + gw_corr_lev_virt)
     579             : 
     580          86 :             CALL dbcsr_get_diag(matrix_tmp_2, vec_Sigma_x_minus_vxc_gw(:, ispin, 1))
     581             : 
     582          86 :             CALL dbcsr_set(matrix_tmp, 0.0_dp)
     583         244 :             CALL dbcsr_set(matrix_tmp_2, 0.0_dp)
     584             : 
     585             :          END DO
     586             : 
     587          72 :          CALL para_env%sum(vec_Sigma_x_minus_vxc_gw)
     588             : 
     589             :       END IF
     590             : 
     591          72 :       CALL dbcsr_release(mo_coeff_b)
     592          72 :       CALL dbcsr_release(matrix_tmp)
     593          72 :       CALL dbcsr_release(matrix_tmp_2)
     594          72 :       IF (do_kpoints_cubic_RPA) THEN
     595           0 :          CALL dbcsr_deallocate_matrix_set(matrix_ks_kp_re)
     596           0 :          CALL dbcsr_deallocate_matrix_set(matrix_ks_kp_im)
     597             :       END IF
     598             : 
     599         158 :       DO ispin = 1, nspins
     600         244 :          DO ikp = 1, nkp
     601          86 :             CALL dbcsr_release_p(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix)
     602         172 :             IF (do_kpoints_cubic_RPA) THEN
     603           0 :                CALL dbcsr_release_p(matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix)
     604             :             END IF
     605             :          END DO
     606             :       END DO
     607             : 
     608         360 :       ALLOCATE (mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(nmo, nspins, nkp))
     609             : 
     610          72 :       IF (print_exx == gw_print_exx) THEN
     611             : 
     612           0 :          IF (do_kpoints_from_Gamma) THEN
     613             : 
     614           0 :             gw_corr_lev_tot = gw_corr_lev_occ + gw_corr_lev_virt
     615             : 
     616             :             CALL get_qs_env(qs_env=qs_env, &
     617           0 :                             kpoints=kpoints)
     618             : 
     619           0 :             CALL setup_trunc_coulomb_pot_for_exchange_self_energy(qs_env)
     620             : 
     621           0 :             CALL compute_kpoints(qs_env, kpoints, unit_nr)
     622             : 
     623           0 :             ALLOCATE (Eigenval_kp(nmo, 1, nspins))
     624             : 
     625           0 :             CALL get_bandstruc_and_k_dependent_MOs(qs_env, Eigenval_kp)
     626             : 
     627           0 :             CALL compute_minus_vxc_kpoints(qs_env)
     628             : 
     629           0 :             nkp_Sigma = SIZE(Eigenval_kp, 2)
     630             : 
     631           0 :             ALLOCATE (vec_Sigma_x(nmo, nkp_Sigma))
     632           0 :             vec_Sigma_x(:, :) = 0.0_dp
     633             : 
     634             :             CALL trafo_to_mo_and_kpoints(qs_env, &
     635             :                                          mat_exchange_for_kp_from_gamma(1)%matrix, &
     636             :                                          vec_Sigma_x(homo - gw_corr_lev_occ + 1:homo + gw_corr_lev_virt, :), &
     637           0 :                                          homo, gw_corr_lev_occ, gw_corr_lev_virt, 1)
     638             : 
     639           0 :             CALL dbcsr_release(mat_exchange_for_kp_from_gamma(1)%matrix)
     640           0 :             DEALLOCATE (mat_exchange_for_kp_from_gamma(1)%matrix)
     641           0 :             DEALLOCATE (mat_exchange_for_kp_from_gamma)
     642             : 
     643           0 :             DEALLOCATE (vec_Sigma_x_minus_vxc_gw)
     644             : 
     645           0 :             ALLOCATE (vec_Sigma_x_minus_vxc_gw(nmo, nspins, nkp_Sigma))
     646             : 
     647             :             vec_Sigma_x_minus_vxc_gw(:, 1, :) = vec_Sigma_x(:, :) + &
     648           0 :                                                 qs_env%mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, :)
     649             : 
     650           0 :             kpoints_Sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma
     651             : 
     652             :          ELSE
     653             : 
     654           0 :             nkp_Sigma = 1
     655             : 
     656             :          END IF
     657             : 
     658           0 :          IF (unit_nr > 0) THEN
     659             : 
     660           0 :             ALLOCATE (Eigenval_kp_HF_at_DFT(nmo, nkp_Sigma))
     661           0 :             Eigenval_kp_HF_at_DFT(:, :) = Eigenval_kp(:, :, 1) + vec_Sigma_x_minus_vxc_gw(:, 1, :)
     662             : 
     663           0 :             min_direct_HF_at_DFT_gap = 100.0_dp
     664             : 
     665           0 :             WRITE (unit_nr, '(T3,A)') ''
     666           0 :             WRITE (unit_nr, '(T3,A)') 'Exchange energies'
     667           0 :             WRITE (unit_nr, '(T3,A)') '-----------------'
     668           0 :             WRITE (unit_nr, '(T3,A)') ''
     669           0 :             WRITE (unit_nr, '(T6,2A)') 'MO                        e_n^DFT          Sigma_x-vxc           e_n^HF@DFT'
     670           0 :             DO ikp = 1, nkp_Sigma
     671           0 :                IF (nkp_Sigma > 1) THEN
     672           0 :                   WRITE (unit_nr, '(T3,A)') ''
     673           0 :                   WRITE (unit_nr, '(T3,A7,I3,A3,I3,A8,3F7.3,A12,3F7.3)') 'Kpoint ', ikp, '  /', nkp_Sigma, &
     674           0 :                      '   xkp =', kpoints_Sigma%xkp(1, ikp), kpoints_Sigma%xkp(2, ikp), &
     675           0 :                      kpoints_Sigma%xkp(3, ikp), '  and  xkp =', -kpoints_Sigma%xkp(1, ikp), &
     676           0 :                      -kpoints_Sigma%xkp(2, ikp), -kpoints_Sigma%xkp(3, ikp)
     677           0 :                   WRITE (unit_nr, '(T3,A)') ''
     678             :                END IF
     679           0 :                DO n_level_gw = 1, gw_corr_lev_occ + gw_corr_lev_virt
     680             : 
     681           0 :                   n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
     682           0 :                   IF (n_level_gw <= gw_corr_lev_occ) THEN
     683           0 :                      occ_virt = 'occ'
     684             :                   ELSE
     685           0 :                      occ_virt = 'vir'
     686             :                   END IF
     687             : 
     688           0 :                   eigval_dft = Eigenval_kp(n_level_gw_ref, ikp, 1)*evolt
     689           0 :                   exx_minus_vxc = REAL(vec_Sigma_x_minus_vxc_gw(n_level_gw_ref, 1, ikp)*evolt, kind=dp)
     690           0 :                   eigval_hf_at_dft = Eigenval_kp_HF_at_DFT(n_level_gw_ref, ikp)*evolt
     691             : 
     692             :                   WRITE (unit_nr, '(T4,I4,3A,3F21.3,3F21.3,3F21.3)') &
     693           0 :                      n_level_gw_ref, ' ( ', occ_virt, ')  ', eigval_dft, exx_minus_vxc, eigval_hf_at_dft
     694             : 
     695             :                END DO
     696           0 :                E_HOMO_GW = MAXVAL(Eigenval_kp_HF_at_DFT(homo - gw_corr_lev_occ + 1:homo, ikp))
     697           0 :                E_LUMO_GW = MINVAL(Eigenval_kp_HF_at_DFT(homo + 1:homo + gw_corr_lev_virt, ikp))
     698           0 :                E_GAP_GW = E_LUMO_GW - E_HOMO_GW
     699           0 :                IF (E_GAP_GW < min_direct_HF_at_DFT_gap) min_direct_HF_at_DFT_gap = E_GAP_GW
     700           0 :                WRITE (unit_nr, '(T3,A)') ''
     701           0 :                WRITE (unit_nr, '(T3,A,F53.2)') 'HF@DFT HOMO-LUMO gap (eV)', E_GAP_GW*evolt
     702           0 :                WRITE (unit_nr, '(T3,A)') ''
     703             :             END DO
     704             : 
     705           0 :             WRITE (unit_nr, '(T3,A)') ''
     706           0 :             WRITE (unit_nr, '(T3,A)') ''
     707           0 :             WRITE (unit_nr, '(T3,A,F63.3)') 'HF@DFT direct bandgap (eV)', min_direct_HF_at_DFT_gap*evolt
     708             : 
     709           0 :             WRITE (unit_nr, '(T3,A)') ''
     710           0 :             WRITE (unit_nr, '(T3,A)') 'End of exchange energies'
     711           0 :             WRITE (unit_nr, '(T3,A)') '------------------------'
     712           0 :             WRITE (unit_nr, '(T3,A)') ''
     713             : 
     714           0 :             CPABORT('Stop after printing exchange energies.')
     715             : 
     716             :          ELSE
     717           0 :             CALL para_env%sync()
     718             :          END IF
     719             : 
     720             :       END IF
     721             : 
     722          72 :       IF (print_exx == gw_read_exx) THEN
     723             : 
     724           0 :          CALL open_file(unit_number=iunit, file_name="exx.out")
     725             : 
     726           0 :          really_read_line = .FALSE.
     727             : 
     728             :          DO WHILE (.TRUE.)
     729             : 
     730           0 :             READ (iunit, '(A)') line
     731             : 
     732           0 :             IF (line == "  End of exchange energies              ") EXIT
     733             : 
     734           0 :             IF (really_read_line) THEN
     735             : 
     736           0 :                READ (line(1:7), *) n_level_gw_ref
     737           0 :                READ (line(17:40), *) tmp
     738             : 
     739           0 :                DO ikp = 1, SIZE(vec_Sigma_x_minus_vxc_gw, 3)
     740           0 :                   vec_Sigma_x_minus_vxc_gw(n_level_gw_ref, 1, ikp) = tmp/evolt
     741             :                END DO
     742             : 
     743             :             END IF
     744             : 
     745           0 :             IF (line == "     MO                    Sigma_x-vxc  ") really_read_line = .TRUE.
     746             : 
     747             :          END DO
     748             : 
     749           0 :          CALL close_file(iunit)
     750             : 
     751             :       END IF
     752             : 
     753             :       ! store vec_Sigma_x_minus_vxc_gw in the mp2_environment
     754        2276 :       mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, :, :) = vec_Sigma_x_minus_vxc_gw(:, :, :)
     755             : 
     756             :       ! clean up
     757          72 :       DEALLOCATE (matrix_sigma_x_minus_vxc, vec_Sigma_x_minus_vxc_gw)
     758          72 :       IF (do_kpoints_cubic_RPA) THEN
     759           0 :          DEALLOCATE (matrix_sigma_x_minus_vxc_im)
     760             :       END IF
     761             : 
     762          72 :       t2 = m_walltime()
     763             : 
     764          72 :       t3 = t2 - t1
     765             : 
     766          72 :       CALL timestop(handle)
     767             : 
     768         288 :    END SUBROUTINE compute_vec_Sigma_x_minus_vxc_gw
     769             : 
     770             : ! **************************************************************************************************
     771             : !> \brief ...
     772             : !> \param kpoints ...
     773             : !> \param matrix_sigma_x_minus_vxc ...
     774             : !> \param matrix_sigma_x_minus_vxc_im ...
     775             : !> \param vec_Sigma_x_minus_vxc_gw ...
     776             : !> \param vec_Sigma_x_minus_vxc_gw_im ...
     777             : !> \param para_env ...
     778             : !> \param nmo ...
     779             : !> \param mp2_env ...
     780             : ! **************************************************************************************************
     781           0 :    SUBROUTINE transform_sigma_x_minus_vxc_to_MO_basis(kpoints, matrix_sigma_x_minus_vxc, &
     782             :                                                       matrix_sigma_x_minus_vxc_im, vec_Sigma_x_minus_vxc_gw, &
     783             :                                                       vec_Sigma_x_minus_vxc_gw_im, para_env, nmo, mp2_env)
     784             : 
     785             :       TYPE(kpoint_type), POINTER                         :: kpoints
     786             :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_sigma_x_minus_vxc, &
     787             :                                                             matrix_sigma_x_minus_vxc_im
     788             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: vec_Sigma_x_minus_vxc_gw, &
     789             :                                                             vec_Sigma_x_minus_vxc_gw_im
     790             :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
     791             :       INTEGER                                            :: nmo
     792             :       TYPE(mp2_type)                                     :: mp2_env
     793             : 
     794             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_sigma_x_minus_vxc_to_MO_basis'
     795             : 
     796             :       INTEGER :: dimen, gw_corr_lev_occ, gw_corr_lev_virt, handle, homo, i_global, iiB, ikp, &
     797             :          ispin, j_global, jjB, ncol_local, nkp, nrow_local, nspins
     798             :       INTEGER, DIMENSION(2)                              :: kp_range
     799           0 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     800             :       REAL(KIND=dp)                                      :: imval, reval
     801             :       TYPE(cp_cfm_type)                                  :: cfm_mos, cfm_sigma_x_minus_vxc, &
     802             :                                                             cfm_sigma_x_minus_vxc_mo_basis, cfm_tmp
     803             :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct
     804             :       TYPE(cp_fm_type)                                   :: fwork_im, fwork_re
     805             :       TYPE(kpoint_env_type), POINTER                     :: kp
     806             :       TYPE(mo_set_type), POINTER                         :: mo_set, mo_set_im, mo_set_re
     807             : 
     808           0 :       CALL timeset(routineN, handle)
     809             : 
     810           0 :       mo_set => kpoints%kp_env(1)%kpoint_env%mos(1, 1)
     811           0 :       CALL get_mo_set(mo_set, nmo=nmo)
     812             : 
     813           0 :       nspins = SIZE(matrix_sigma_x_minus_vxc, 1)
     814           0 :       CALL get_kpoint_info(kpoints, nkp=nkp, kp_range=kp_range)
     815             : 
     816             :       ! if this CPASSERT is wrong, please make sure that the kpoint group size PARALLEL_GROUP_SIZE
     817             :       ! in the kpoint environment &DFT &KPOINTS is -1
     818           0 :       CPASSERT(kp_range(1) == 1 .AND. kp_range(2) == nkp)
     819             : 
     820           0 :       ALLOCATE (vec_Sigma_x_minus_vxc_gw(nmo, nspins, nkp))
     821           0 :       vec_Sigma_x_minus_vxc_gw = 0.0_dp
     822             : 
     823           0 :       ALLOCATE (vec_Sigma_x_minus_vxc_gw_im(nmo, nspins, nkp))
     824           0 :       vec_Sigma_x_minus_vxc_gw_im = 0.0_dp
     825             : 
     826           0 :       CALL cp_fm_get_info(mo_set%mo_coeff, matrix_struct=matrix_struct)
     827           0 :       CALL cp_fm_create(fwork_re, matrix_struct)
     828           0 :       CALL cp_fm_create(fwork_im, matrix_struct)
     829           0 :       CALL cp_cfm_create(cfm_mos, matrix_struct)
     830           0 :       CALL cp_cfm_create(cfm_sigma_x_minus_vxc, matrix_struct)
     831           0 :       CALL cp_cfm_create(cfm_sigma_x_minus_vxc_mo_basis, matrix_struct)
     832           0 :       CALL cp_cfm_create(cfm_tmp, matrix_struct)
     833             : 
     834             :       CALL cp_cfm_get_info(matrix=cfm_sigma_x_minus_vxc_mo_basis, &
     835             :                            nrow_local=nrow_local, &
     836             :                            ncol_local=ncol_local, &
     837             :                            row_indices=row_indices, &
     838           0 :                            col_indices=col_indices)
     839             : 
     840             :       ! Transform matrix_sigma_x_minus_vxc to MO basis
     841           0 :       DO ikp = 1, nkp
     842             : 
     843           0 :          kp => kpoints%kp_env(ikp)%kpoint_env
     844             : 
     845           0 :          DO ispin = 1, nspins
     846             : 
     847             :             ! v_xc_n to fm matrix
     848           0 :             CALL copy_dbcsr_to_fm(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, fwork_re)
     849           0 :             CALL copy_dbcsr_to_fm(matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix, fwork_im)
     850             : 
     851           0 :             CALL cp_cfm_scale_and_add_fm(z_zero, cfm_sigma_x_minus_vxc, z_one, fwork_re)
     852           0 :             CALL cp_cfm_scale_and_add_fm(z_one, cfm_sigma_x_minus_vxc, gaussi, fwork_im)
     853             : 
     854             :             ! get real part (1) and imag. part (2) of the mo coeffs
     855           0 :             mo_set_re => kp%mos(1, ispin)
     856           0 :             mo_set_im => kp%mos(2, ispin)
     857             : 
     858           0 :             CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mos, z_one, mo_set_re%mo_coeff)
     859           0 :             CALL cp_cfm_scale_and_add_fm(z_one, cfm_mos, gaussi, mo_set_im%mo_coeff)
     860             : 
     861             :             ! tmp = V(k)*C(k)
     862             :             CALL parallel_gemm('N', 'N', nmo, nmo, nmo, z_one, cfm_sigma_x_minus_vxc, &
     863           0 :                                cfm_mos, z_zero, cfm_tmp)
     864             : 
     865             :             ! V_n(k) = C^H(k)*tmp
     866             :             CALL parallel_gemm('C', 'N', nmo, nmo, nmo, z_one, cfm_mos, cfm_tmp, &
     867           0 :                                z_zero, cfm_sigma_x_minus_vxc_mo_basis)
     868             : 
     869           0 :             DO jjB = 1, ncol_local
     870             : 
     871           0 :                j_global = col_indices(jjB)
     872             : 
     873           0 :                DO iiB = 1, nrow_local
     874             : 
     875           0 :                   i_global = row_indices(iiB)
     876             : 
     877           0 :                   IF (j_global == i_global .AND. i_global <= nmo) THEN
     878             : 
     879           0 :                      reval = REAL(cfm_sigma_x_minus_vxc_mo_basis%local_data(iiB, jjB), kind=dp)
     880           0 :                      imval = AIMAG(cfm_sigma_x_minus_vxc_mo_basis%local_data(iiB, jjB))
     881             : 
     882           0 :                      vec_Sigma_x_minus_vxc_gw(i_global, ispin, ikp) = reval
     883           0 :                      vec_Sigma_x_minus_vxc_gw_im(i_global, ispin, ikp) = imval
     884             : 
     885             :                   END IF
     886             : 
     887             :                END DO
     888             : 
     889             :             END DO
     890             : 
     891             :          END DO
     892             : 
     893             :       END DO
     894             : 
     895           0 :       CALL para_env%sum(vec_Sigma_x_minus_vxc_gw)
     896           0 :       CALL para_env%sum(vec_Sigma_x_minus_vxc_gw_im)
     897             : 
     898             :       ! also adjust in the case of kpoints too big gw_corr_lev_occ and gw_corr_lev_virt
     899           0 :       DO ispin = 1, nspins
     900             :          CALL get_mo_set(mo_set=kpoints%kp_env(1)%kpoint_env%mos(ispin, 1), &
     901           0 :                          homo=homo, nao=dimen)
     902           0 :          gw_corr_lev_occ = mp2_env%ri_g0w0%corr_mos_occ
     903           0 :          gw_corr_lev_virt = mp2_env%ri_g0w0%corr_mos_virt
     904             :          ! if corrected occ/virt levels exceed the number of occ/virt levels,
     905             :          ! correct all occ/virt level energies
     906           0 :          IF (gw_corr_lev_occ > homo) gw_corr_lev_occ = homo
     907           0 :          IF (gw_corr_lev_virt > dimen - homo) gw_corr_lev_virt = dimen - homo
     908           0 :          IF (ispin == 1) THEN
     909           0 :             mp2_env%ri_g0w0%corr_mos_occ = gw_corr_lev_occ
     910           0 :             mp2_env%ri_g0w0%corr_mos_virt = gw_corr_lev_virt
     911           0 :          ELSE IF (ispin == 2) THEN
     912             :             ! ensure that the total number of corrected MOs is the same for alpha and beta, important
     913             :             ! for parallelization
     914           0 :             IF (mp2_env%ri_g0w0%corr_mos_occ + mp2_env%ri_g0w0%corr_mos_virt /= &
     915             :                 gw_corr_lev_occ + gw_corr_lev_virt) THEN
     916           0 :                gw_corr_lev_virt = mp2_env%ri_g0w0%corr_mos_occ + mp2_env%ri_g0w0%corr_mos_virt - gw_corr_lev_occ
     917             :             END IF
     918           0 :             mp2_env%ri_g0w0%corr_mos_occ_beta = gw_corr_lev_occ
     919           0 :             mp2_env%ri_g0w0%corr_mos_virt_beta = gw_corr_lev_virt
     920             :          END IF
     921             :       END DO
     922             : 
     923           0 :       CALL cp_fm_release(fwork_re)
     924           0 :       CALL cp_fm_release(fwork_im)
     925           0 :       CALL cp_cfm_release(cfm_mos)
     926           0 :       CALL cp_cfm_release(cfm_sigma_x_minus_vxc)
     927           0 :       CALL cp_cfm_release(cfm_sigma_x_minus_vxc_mo_basis)
     928           0 :       CALL cp_cfm_release(cfm_tmp)
     929             : 
     930           0 :       CALL timestop(handle)
     931             : 
     932           0 :    END SUBROUTINE
     933             : 
     934             : ! **************************************************************************************************
     935             : !> \brief ...
     936             : !> \param matrix_ks_transl ...
     937             : !> \param matrix_ks_kp_re ...
     938             : !> \param matrix_ks_kp_im ...
     939             : !> \param kpoints ...
     940             : ! **************************************************************************************************
     941           0 :    SUBROUTINE transform_matrix_ks_to_kp(matrix_ks_transl, matrix_ks_kp_re, matrix_ks_kp_im, kpoints)
     942             : 
     943             :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks_transl, matrix_ks_kp_re, &
     944             :                                                             matrix_ks_kp_im
     945             :       TYPE(kpoint_type), POINTER                         :: kpoints
     946             : 
     947             :       CHARACTER(len=*), PARAMETER :: routineN = 'transform_matrix_ks_to_kp'
     948             : 
     949             :       INTEGER                                            :: handle, ikp, ispin, nkp, nspin
     950           0 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
     951           0 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: xkp
     952             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     953           0 :          POINTER                                         :: sab_nl
     954             : 
     955           0 :       CALL timeset(routineN, handle)
     956             : 
     957           0 :       NULLIFY (sab_nl)
     958           0 :       CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, sab_nl=sab_nl, cell_to_index=cell_to_index)
     959             : 
     960           0 :       CPASSERT(ASSOCIATED(sab_nl))
     961             : 
     962           0 :       nspin = SIZE(matrix_ks_transl, 1)
     963             : 
     964           0 :       DO ikp = 1, nkp
     965           0 :          DO ispin = 1, nspin
     966             : 
     967           0 :             CALL dbcsr_set(matrix_ks_kp_re(ispin, ikp)%matrix, 0.0_dp)
     968           0 :             CALL dbcsr_set(matrix_ks_kp_im(ispin, ikp)%matrix, 0.0_dp)
     969             :             CALL rskp_transform(rmatrix=matrix_ks_kp_re(ispin, ikp)%matrix, &
     970             :                                 cmatrix=matrix_ks_kp_im(ispin, ikp)%matrix, &
     971             :                                 rsmat=matrix_ks_transl, ispin=ispin, &
     972           0 :                                 xkp=xkp(1:3, ikp), cell_to_index=cell_to_index, sab_nl=sab_nl)
     973             : 
     974             :          END DO
     975             :       END DO
     976             : 
     977           0 :       CALL timestop(handle)
     978             : 
     979           0 :    END SUBROUTINE
     980             : 
     981             : ! **************************************************************************************************
     982             : !> \brief ...
     983             : !> \param matrix_ks_transl ...
     984             : !> \param matrix_ks_kp_re ...
     985             : !> \param matrix_ks_kp_im ...
     986             : !> \param kpoints ...
     987             : ! **************************************************************************************************
     988           0 :    SUBROUTINE allocate_matrix_ks_kp(matrix_ks_transl, matrix_ks_kp_re, matrix_ks_kp_im, kpoints)
     989             : 
     990             :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks_transl, matrix_ks_kp_re, &
     991             :                                                             matrix_ks_kp_im
     992             :       TYPE(kpoint_type), POINTER                         :: kpoints
     993             : 
     994             :       CHARACTER(len=*), PARAMETER :: routineN = 'allocate_matrix_ks_kp'
     995             : 
     996             :       INTEGER                                            :: handle, ikp, ispin, nkp, nspin
     997           0 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
     998           0 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: xkp
     999             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1000           0 :          POINTER                                         :: sab_nl
    1001             : 
    1002           0 :       CALL timeset(routineN, handle)
    1003             : 
    1004           0 :       NULLIFY (sab_nl)
    1005           0 :       CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, sab_nl=sab_nl, cell_to_index=cell_to_index)
    1006             : 
    1007           0 :       CPASSERT(ASSOCIATED(sab_nl))
    1008             : 
    1009           0 :       nspin = SIZE(matrix_ks_transl, 1)
    1010             : 
    1011           0 :       NULLIFY (matrix_ks_kp_re, matrix_ks_kp_im)
    1012           0 :       CALL dbcsr_allocate_matrix_set(matrix_ks_kp_re, nspin, nkp)
    1013           0 :       CALL dbcsr_allocate_matrix_set(matrix_ks_kp_im, nspin, nkp)
    1014             : 
    1015           0 :       DO ikp = 1, nkp
    1016           0 :       DO ispin = 1, nspin
    1017             : 
    1018           0 :          ALLOCATE (matrix_ks_kp_re(ispin, ikp)%matrix)
    1019           0 :          ALLOCATE (matrix_ks_kp_im(ispin, ikp)%matrix)
    1020             : 
    1021             :          CALL dbcsr_create(matrix_ks_kp_re(ispin, ikp)%matrix, &
    1022             :                            template=matrix_ks_transl(1, 1)%matrix, &
    1023           0 :                            matrix_type=dbcsr_type_symmetric)
    1024             :          CALL dbcsr_create(matrix_ks_kp_im(ispin, ikp)%matrix, &
    1025             :                            template=matrix_ks_transl(1, 1)%matrix, &
    1026           0 :                            matrix_type=dbcsr_type_antisymmetric)
    1027             : 
    1028           0 :          CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_kp_re(ispin, ikp)%matrix, sab_nl)
    1029           0 :          CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_kp_im(ispin, ikp)%matrix, sab_nl)
    1030             : 
    1031           0 :          CALL dbcsr_set(matrix_ks_kp_re(ispin, ikp)%matrix, 0.0_dp)
    1032           0 :          CALL dbcsr_set(matrix_ks_kp_im(ispin, ikp)%matrix, 0.0_dp)
    1033             : 
    1034             :       END DO
    1035             :       END DO
    1036             : 
    1037           0 :       CALL timestop(handle)
    1038             : 
    1039           0 :    END SUBROUTINE
    1040             : 
    1041             : END MODULE rpa_gw_sigma_x
    1042             : 

Generated by: LCOV version 1.15