LCOV - code coverage report
Current view: top level - src - rpa_rse.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 210 216 97.2 %
Date: 2024-03-28 07:31:50 Functions: 4 4 100.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 compute singles correction to RPA (RSE)
      10             : !> \par History
      11             : !>      08.2019 created [Vladimir Rybkin]
      12             : !> \author Vladimir Rybkin
      13             : ! **************************************************************************************************
      14             : MODULE rpa_rse
      15             : 
      16             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      17             :    USE cp_control_types,                ONLY: dft_control_type
      18             :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      19             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      20             :                                               copy_fm_to_dbcsr,&
      21             :                                               dbcsr_allocate_matrix_set
      22             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_scale_and_add
      23             :    USE cp_fm_diag,                      ONLY: choose_eigv_solver
      24             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      25             :                                               cp_fm_struct_release,&
      26             :                                               cp_fm_struct_type
      27             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      28             :                                               cp_fm_get_diag,&
      29             :                                               cp_fm_get_info,&
      30             :                                               cp_fm_release,&
      31             :                                               cp_fm_set_all,&
      32             :                                               cp_fm_to_fm_submat,&
      33             :                                               cp_fm_type
      34             :    USE dbcsr_api,                       ONLY: dbcsr_copy,&
      35             :                                               dbcsr_create,&
      36             :                                               dbcsr_init_p,&
      37             :                                               dbcsr_p_type,&
      38             :                                               dbcsr_release,&
      39             :                                               dbcsr_scale,&
      40             :                                               dbcsr_set,&
      41             :                                               dbcsr_type_symmetric
      42             :    USE hfx_energy_potential,            ONLY: integrate_four_center
      43             :    USE hfx_exx,                         ONLY: exx_post_hfx,&
      44             :                                               exx_pre_hfx
      45             :    USE hfx_ri,                          ONLY: hfx_ri_update_ks
      46             :    USE input_section_types,             ONLY: section_vals_get,&
      47             :                                               section_vals_get_subs_vals,&
      48             :                                               section_vals_type,&
      49             :                                               section_vals_val_get
      50             :    USE kinds,                           ONLY: dp
      51             :    USE message_passing,                 ONLY: mp_para_env_type
      52             :    USE mp2_types,                       ONLY: mp2_type
      53             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      54             :    USE pw_types,                        ONLY: pw_r3d_rs_type
      55             :    USE qs_energy_types,                 ONLY: qs_energy_type
      56             :    USE qs_environment_types,            ONLY: get_qs_env,&
      57             :                                               qs_environment_type
      58             :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      59             :    USE qs_ks_utils,                     ONLY: compute_matrix_vxc
      60             :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
      61             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      62             :                                               qs_rho_type
      63             :    USE qs_vxc,                          ONLY: qs_vxc_create
      64             : 
      65             : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
      66             : 
      67             : #include "./base/base_uses.f90"
      68             : 
      69             :    IMPLICIT NONE
      70             : 
      71             :    PRIVATE
      72             : 
      73             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_rse'
      74             : 
      75             :    PUBLIC :: rse_energy
      76             : 
      77             : CONTAINS
      78             : 
      79             : ! **************************************************************************************************
      80             : !> \brief Single excitations energy corrections for RPA
      81             : !> \param qs_env ...
      82             : !> \param mp2_env ...
      83             : !> \param para_env ...
      84             : !> \param dft_control ...
      85             : !> \param mo_coeff ...
      86             : !> \param nmo ...
      87             : !> \param homo ...
      88             : !> \param Eigenval ...
      89             : !> \author Vladimir Rybkin, 08/2019
      90             : ! **************************************************************************************************
      91           4 :    SUBROUTINE rse_energy(qs_env, mp2_env, para_env, dft_control, &
      92           4 :                          mo_coeff, nmo, homo, Eigenval)
      93             :       TYPE(qs_environment_type), INTENT(IN), POINTER     :: qs_env
      94             :       TYPE(mp2_type), INTENT(INOUT)                      :: mp2_env
      95             :       TYPE(mp_para_env_type), INTENT(IN), POINTER        :: para_env
      96             :       TYPE(dft_control_type), INTENT(IN), POINTER        :: dft_control
      97             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: mo_coeff
      98             :       INTEGER, INTENT(IN)                                :: nmo
      99             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: homo
     100             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: Eigenval
     101             : 
     102             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'rse_energy'
     103             : 
     104             :       INTEGER                                            :: dimen, handle, i_global, iiB, ispin, &
     105             :                                                             j_global, jjB, n_rep_hf, ncol_local, &
     106             :                                                             nrow_local, nspins
     107           4 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     108             :       LOGICAL                                            :: do_hfx, hfx_treat_lsd_in_core
     109             :       REAL(KIND=dp)                                      :: coeff, corr, rse_corr
     110           4 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: diag_diff
     111             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     112             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     113             :       TYPE(cp_fm_type)                                   :: fm_ao, fm_ao_mo
     114           4 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_P_mu_nu, fm_X_mo, fm_XC_mo
     115           4 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mat_mu_nu, matrix_s, rho_ao
     116             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     117           4 :          POINTER                                         :: sab_orb
     118             :       TYPE(qs_energy_type), POINTER                      :: energy
     119             :       TYPE(qs_rho_type), POINTER                         :: rho
     120             :       TYPE(section_vals_type), POINTER                   :: hfx_sections, input
     121             : 
     122           4 :       CALL timeset(routineN, handle)
     123             : 
     124           4 :       nspins = dft_control%nspins
     125             : 
     126             :       ! Pick the diagonal terms
     127             :       CALL cp_fm_get_info(matrix=mo_coeff(1), &
     128             :                           nrow_local=nrow_local, &
     129             :                           ncol_local=ncol_local, &
     130             :                           row_indices=row_indices, &
     131           4 :                           col_indices=col_indices)
     132             : 
     133             :       ! start collecting stuff
     134           4 :       dimen = nmo
     135           4 :       NULLIFY (input, matrix_s, blacs_env, rho, energy, sab_orb)
     136             :       CALL get_qs_env(qs_env, &
     137             :                       input=input, &
     138             :                       matrix_s=matrix_s, &
     139             :                       blacs_env=blacs_env, &
     140             :                       rho=rho, &
     141             :                       energy=energy, &
     142           4 :                       sab_orb=sab_orb)
     143             : 
     144           4 :       CALL qs_rho_get(rho, rho_ao=rho_ao)
     145             : 
     146             :       ! hfx section
     147           4 :       NULLIFY (hfx_sections)
     148           4 :       hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
     149           4 :       CALL section_vals_get(hfx_sections, explicit=do_hfx, n_repetition=n_rep_hf)
     150           4 :       IF (do_hfx) THEN
     151             :          CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
     152           4 :                                    i_rep_section=1)
     153             :       END IF
     154             : 
     155             :       ! create work array
     156           4 :       NULLIFY (mat_mu_nu)
     157           4 :       CALL dbcsr_allocate_matrix_set(mat_mu_nu, nspins)
     158          10 :       DO ispin = 1, nspins
     159           6 :          ALLOCATE (mat_mu_nu(ispin)%matrix)
     160             :          CALL dbcsr_create(matrix=mat_mu_nu(ispin)%matrix, template=matrix_s(1)%matrix, name="T_mu_nu", &
     161           6 :                            matrix_type=dbcsr_type_symmetric, nze=0)
     162           6 :          CALL cp_dbcsr_alloc_block_from_nbl(mat_mu_nu(ispin)%matrix, sab_orb)
     163          10 :          CALL dbcsr_set(mat_mu_nu(ispin)%matrix, 0.0_dp)
     164             :       END DO
     165             : 
     166             :       ! Dense (full) matrices
     167          18 :       ALLOCATE (fm_P_mu_nu(nspins))
     168           4 :       NULLIFY (fm_struct_tmp)
     169             :       CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
     170           4 :                                nrow_global=dimen, ncol_global=dimen)
     171          10 :       DO ispin = 1, nspins
     172           6 :          CALL cp_fm_create(fm_P_mu_nu(ispin), fm_struct_tmp, name="P_mu_nu")
     173          10 :          CALL cp_fm_set_all(fm_P_mu_nu(ispin), 0.0_dp)
     174             :       END DO
     175           4 :       CALL cp_fm_struct_release(fm_struct_tmp)
     176             : 
     177           4 :       NULLIFY (fm_struct_tmp)
     178             :       CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
     179           4 :                                nrow_global=dimen, ncol_global=dimen)
     180          32 :       ALLOCATE (fm_X_mo(nspins), fm_XC_mo(nspins))
     181          10 :       DO ispin = 1, nspins
     182           6 :          CALL cp_fm_create(fm_X_mo(ispin), fm_struct_tmp, name="f_X_mo")
     183           6 :          CALL cp_fm_create(fm_XC_mo(ispin), fm_struct_tmp, name="f_XC_mo")
     184           6 :          CALL cp_fm_set_all(fm_X_mo(ispin), 0.0_dp)
     185          10 :          CALL cp_fm_set_all(fm_XC_mo(ispin), 0.0_dp)
     186             :       END DO
     187           4 :       CALL cp_fm_struct_release(fm_struct_tmp)
     188             : 
     189             :       CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
     190           4 :                                nrow_global=dimen, ncol_global=dimen)
     191           4 :       CALL cp_fm_create(fm_ao, fm_struct_tmp, name="f_ao")
     192           4 :       CALL cp_fm_struct_release(fm_struct_tmp)
     193           4 :       CALL cp_fm_set_all(fm_ao, 0.0_dp)
     194             :       CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
     195           4 :                                nrow_global=dimen, ncol_global=dimen)
     196           4 :       CALL cp_fm_create(fm_ao_mo, fm_struct_tmp, name="f_ao_mo")
     197           4 :       CALL cp_fm_struct_release(fm_struct_tmp)
     198           4 :       CALL cp_fm_set_all(fm_ao_mo, 0.0_dp)
     199             : 
     200             :       !
     201             :       !     Ready with preparations, do the real staff
     202             :       !
     203             : 
     204             :       ! Obtain density matrix like quantity
     205             : 
     206           4 :       coeff = 1.0_dp
     207           4 :       IF (nspins == 1) coeff = 2.0_dp
     208          10 :       DO ispin = 1, nspins
     209             :          CALL parallel_gemm(transa='N', transb='T', m=dimen, n=dimen, k=homo(ispin), alpha=coeff, &
     210             :                             matrix_a=mo_coeff(ispin), matrix_b=mo_coeff(ispin), &
     211          10 :                             beta=0.0_dp, matrix_c=fm_P_mu_nu(ispin))
     212             :       END DO
     213             : 
     214             :       ! Calculate exact exchange contribution
     215             :       CALL exchange_contribution(qs_env, para_env, dimen, mo_coeff, &
     216             :                                  hfx_sections, n_rep_hf, &
     217             :                                  rho, mat_mu_nu, fm_P_mu_nu, &
     218             :                                  fm_ao, fm_X_mo, fm_ao_mo, &
     219           4 :                                  mp2_env%ri_grad%free_hfx_buffer)
     220             : 
     221             :       ! Calculate DFT exchange-correlation contribution
     222           4 :       CALL xc_contribution(qs_env, fm_ao, fm_ao_mo, fm_XC_mo, mo_coeff, dimen)
     223             : 
     224          12 :       ALLOCATE (diag_diff(dimen))
     225           4 :       rse_corr = 0.0_dp
     226             : 
     227          10 :       DO ispin = 1, nspins
     228             :          ! Compute the correction matrix: it is stored in fm_X_mo
     229           6 :          CALL cp_fm_scale_and_add(1.0_dp, fm_X_mo(ispin), -1.0_dp, fm_XC_mo(ispin))
     230             : 
     231             :          ! Pick the diagonal terms
     232           6 :          CALL cp_fm_get_diag(fm_X_mo(ispin), diag_diff)
     233             : 
     234             :          ! Compute the correction
     235             :          CALL cp_fm_get_info(matrix=fm_X_mo(ispin), &
     236             :                              nrow_local=nrow_local, &
     237             :                              ncol_local=ncol_local, &
     238             :                              row_indices=row_indices, &
     239           6 :                              col_indices=col_indices)
     240             : 
     241           6 :          corr = 0.0_dp
     242             : 
     243             : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
     244             : !$OMP             REDUCTION(+: corr) &
     245           6 : !$OMP             SHARED(ncol_local,nrow_local,col_indices,row_indices,diag_diff,eigenval,fm_X_mo,homo,ispin)
     246             :          DO jjB = 1, ncol_local
     247             :             j_global = col_indices(jjB)
     248             :             DO iiB = 1, nrow_local
     249             :                i_global = row_indices(iiB)
     250             :                IF ((i_global .LE. homo(ispin)) .AND. (j_global .GT. homo(ispin))) THEN
     251             :                   corr = corr + fm_X_mo(ispin)%local_data(iib, jjb)**2.0_dp/ &
     252             :                          (eigenval(i_global, ispin) - eigenval(j_global, ispin) - diag_diff(i_global) + diag_diff(j_global))
     253             :                END IF
     254             :             END DO
     255             :          END DO
     256             : !$OMP END PARALLEL DO
     257             : 
     258          10 :          rse_corr = rse_corr + corr
     259             :       END DO
     260             : 
     261           4 :       CALL para_env%sum(rse_corr)
     262             : 
     263           4 :       IF (nspins == 1) rse_corr = rse_corr*2.0_dp
     264             : 
     265           4 :       mp2_env%ri_rpa%rse_corr_diag = rse_corr
     266             : 
     267           4 :       CALL non_diag_rse(fm_X_mo, eigenval, dimen, homo, para_env, blacs_env, rse_corr)
     268             : 
     269           4 :       IF (nspins == 1) rse_corr = rse_corr*2.0_dp
     270             : 
     271           4 :       mp2_env%ri_rpa%rse_corr = rse_corr
     272             : 
     273             :       ! Release staff
     274           4 :       DEALLOCATE (diag_diff)
     275           4 :       CALL cp_fm_release(fm_ao)
     276           4 :       CALL cp_fm_release(fm_ao_mo)
     277           4 :       CALL cp_fm_release(fm_P_mu_nu)
     278           4 :       CALL cp_fm_release(fm_X_mo)
     279           4 :       CALL cp_fm_release(fm_XC_mo)
     280          10 :       DO ispin = 1, nspins
     281           6 :          CALL dbcsr_release(mat_mu_nu(ispin)%matrix)
     282          10 :          DEALLOCATE (mat_mu_nu(ispin)%matrix)
     283             :       END DO
     284           4 :       DEALLOCATE (mat_mu_nu)
     285             : 
     286           4 :       CALL timestop(handle)
     287             : 
     288          12 :    END SUBROUTINE rse_energy
     289             : 
     290             : ! **************************************************************************************************
     291             : !> \brief HF exchange occupied-virtual matrix
     292             : !> \param qs_env ...
     293             : !> \param para_env ...
     294             : !> \param dimen ...
     295             : !> \param mo_coeff ...
     296             : !> \param hfx_sections ...
     297             : !> \param n_rep_hf ...
     298             : !> \param rho_work ...
     299             : !> \param mat_mu_nu ...
     300             : !> \param fm_P_mu_nu ...
     301             : !> \param fm_X_ao ...
     302             : !> \param fm_X_mo ...
     303             : !> \param fm_X_ao_mo ...
     304             : !> \param recalc_hfx_integrals ...
     305             : ! **************************************************************************************************
     306          12 :    SUBROUTINE exchange_contribution(qs_env, para_env, dimen, mo_coeff, &
     307             :                                     hfx_sections, n_rep_hf, &
     308           4 :                                     rho_work, mat_mu_nu, fm_P_mu_nu, &
     309           4 :                                     fm_X_ao, fm_X_mo, fm_X_ao_mo, &
     310             :                                     recalc_hfx_integrals)
     311             :       TYPE(qs_environment_type), INTENT(IN), POINTER     :: qs_env
     312             :       TYPE(mp_para_env_type), INTENT(IN), POINTER        :: para_env
     313             :       INTEGER, INTENT(IN)                                :: dimen
     314             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: mo_coeff
     315             :       TYPE(section_vals_type), INTENT(IN), POINTER       :: hfx_sections
     316             :       INTEGER, INTENT(IN)                                :: n_rep_hf
     317             :       TYPE(qs_rho_type), INTENT(IN), POINTER             :: rho_work
     318             :       TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
     319             :          POINTER                                         :: mat_mu_nu
     320             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: fm_P_mu_nu
     321             :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_X_ao
     322             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: fm_X_mo
     323             :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_X_ao_mo
     324             :       LOGICAL, INTENT(IN)                                :: recalc_hfx_integrals
     325             : 
     326             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'exchange_contribution'
     327             : 
     328             :       INTEGER                                            :: handle, irep, is, ns
     329             :       LOGICAL                                            :: my_recalc_hfx_integrals
     330             :       REAL(KIND=dp)                                      :: ehfx
     331           4 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: P_mu_nu, rho_work_ao
     332           4 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_2d, rho_ao_2d
     333             : 
     334           4 :       CALL timeset(routineN, handle)
     335             : 
     336           4 :       my_recalc_hfx_integrals = recalc_hfx_integrals
     337             : 
     338           4 :       CALL qs_rho_get(rho_work, rho_ao=rho_work_ao)
     339           4 :       ns = SIZE(rho_work_ao)
     340           4 :       NULLIFY (P_mu_nu)
     341           4 :       CALL dbcsr_allocate_matrix_set(P_mu_nu, ns)
     342          10 :       DO is = 1, ns
     343           6 :          CALL dbcsr_init_p(P_mu_nu(is)%matrix)
     344           6 :          CALL dbcsr_create(P_mu_nu(is)%matrix, template=rho_work_ao(1)%matrix)
     345           6 :          CALL dbcsr_copy(P_mu_nu(is)%matrix, rho_work_ao(1)%matrix)
     346          10 :          CALL dbcsr_set(P_mu_nu(is)%matrix, 0.0_dp)
     347             :       END DO
     348             : 
     349           4 :       CALL exx_pre_hfx(hfx_sections, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx)
     350          10 :       DO is = 1, ns
     351           6 :          CALL copy_fm_to_dbcsr(fm_P_mu_nu(is), P_mu_nu(1)%matrix, keep_sparsity=.TRUE.)
     352             : 
     353           6 :          CALL dbcsr_set(mat_mu_nu(1)%matrix, 0.0_dp)
     354             : 
     355           6 :          IF (qs_env%mp2_env%ri_rpa%x_data(1, 1)%do_hfx_ri) THEN
     356             : 
     357           0 :             DO irep = 1, n_rep_hf
     358           0 :                rho_ao_2d(1:ns, 1:1) => P_mu_nu(1:ns)
     359           0 :                mat_2d(1:ns, 1:1) => mat_mu_nu(1:ns)
     360             :                CALL hfx_ri_update_ks(qs_env, qs_env%mp2_env%ri_rpa%x_data(irep, 1)%ri_data, mat_2d, ehfx, &
     361             :                                      rho_ao=rho_ao_2d, geometry_did_change=my_recalc_hfx_integrals, nspins=1, &
     362           0 :                                      hf_fraction=qs_env%mp2_env%ri_rpa%x_data(irep, 1)%general_parameter%fraction)
     363             : 
     364           0 :                IF (ns == 2) CALL dbcsr_scale(mat_mu_nu(1)%matrix, 2.0_dp)
     365           0 :                my_recalc_hfx_integrals = .FALSE.
     366             :             END DO
     367             : 
     368             :          ELSE
     369             : 
     370          12 :             DO irep = 1, n_rep_hf
     371           6 :                rho_ao_2d(1:ns, 1:1) => P_mu_nu(1:ns)
     372           6 :                mat_2d(1:ns, 1:1) => mat_mu_nu(1:ns)
     373             :                CALL integrate_four_center(qs_env, qs_env%mp2_env%ri_rpa%x_data, mat_2d, ehfx, rho_ao_2d, hfx_sections, &
     374             :                                           para_env, my_recalc_hfx_integrals, irep, .TRUE., &
     375           6 :                                           ispin=1)
     376             : 
     377          12 :                my_recalc_hfx_integrals = .FALSE.
     378             :             END DO
     379             :          END IF
     380             : 
     381             :          ! copy back to fm
     382           6 :          CALL cp_fm_set_all(fm_X_ao, 0.0_dp)
     383           6 :          CALL copy_dbcsr_to_fm(matrix=mat_mu_nu(1)%matrix, fm=fm_X_ao)
     384           6 :          CALL cp_fm_set_all(fm_X_mo(is), 0.0_dp)
     385             : 
     386             :          ! First index
     387             :          CALL parallel_gemm('T', 'N', dimen, dimen, dimen, 1.0_dp, &
     388           6 :                             mo_coeff(is), fm_X_ao, 0.0_dp, fm_X_ao_mo)
     389             : 
     390             :          ! Second index
     391             :          CALL parallel_gemm('N', 'N', dimen, dimen, dimen, 1.0_dp, &
     392          10 :                             fm_X_ao_mo, mo_coeff(is), 1.0_dp, fm_X_mo(is))
     393             : 
     394             :       END DO
     395           4 :       CALL exx_post_hfx(qs_env, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx)
     396             : 
     397             :       ! Release dbcsr objects
     398          10 :       DO is = 1, SIZE(P_mu_nu)
     399           6 :          CALL dbcsr_release(P_mu_nu(is)%matrix)
     400          10 :          DEALLOCATE (P_mu_nu(is)%matrix)
     401             :       END DO
     402           4 :       DEALLOCATE (P_mu_nu)
     403             : 
     404           4 :       CALL timestop(handle)
     405             : 
     406           4 :    END SUBROUTINE exchange_contribution
     407             : 
     408             : ! **************************************************************************************************
     409             : !> \brief Exchange-correlation occupied-virtual matrix
     410             : !> \param qs_env ...
     411             : !> \param fm_XC_ao ...
     412             : !> \param fm_XC_ao_mo ...
     413             : !> \param fm_XC_mo ...
     414             : !> \param mo_coeff ...
     415             : !> \param dimen ...
     416             : ! **************************************************************************************************
     417           4 :    SUBROUTINE xc_contribution(qs_env, fm_XC_ao, fm_XC_ao_mo, fm_XC_mo, mo_coeff, dimen)
     418             :       TYPE(qs_environment_type), INTENT(IN), POINTER     :: qs_env
     419             :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_XC_ao, fm_XC_ao_mo
     420             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: fm_XC_mo, mo_coeff
     421             :       INTEGER, INTENT(IN)                                :: dimen
     422             : 
     423             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'xc_contribution'
     424             : 
     425             :       INTEGER                                            :: handle, i
     426             :       REAL(KIND=dp)                                      :: exc
     427           4 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_vxc
     428           4 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: tau_rspace, v_rspace
     429             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     430             :       TYPE(qs_rho_type), POINTER                         :: rho
     431             :       TYPE(section_vals_type), POINTER                   :: input, xc_section
     432             : 
     433           4 :       CALL timeset(routineN, handle)
     434             : 
     435           4 :       NULLIFY (matrix_vxc, v_rspace, tau_rspace, input, xc_section, ks_env, &
     436           4 :                rho)
     437           4 :       CALL get_qs_env(qs_env, matrix_vxc=matrix_vxc, input=input, ks_env=ks_env, rho=rho)
     438           4 :       xc_section => section_vals_get_subs_vals(input, "DFT%XC")
     439             : 
     440             :       ! Compute XC matrix in AO basis
     441             :       CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=xc_section, &
     442           4 :                          vxc_rho=v_rspace, vxc_tau=tau_rspace, exc=exc)
     443             : 
     444           4 :       IF (ASSOCIATED(v_rspace)) THEN
     445           4 :          CALL compute_matrix_vxc(qs_env=qs_env, v_rspace=v_rspace, matrix_vxc=matrix_vxc)
     446             : 
     447          10 :          DO i = 1, SIZE(v_rspace)
     448          10 :             CALL v_rspace(i)%release()
     449             :          END DO
     450           4 :          DEALLOCATE (v_rspace)
     451             : 
     452          10 :          DO i = 1, SIZE(matrix_vxc)
     453           6 :             CALL cp_fm_set_all(fm_XC_ao, 0.0_dp)
     454           6 :             CALL copy_dbcsr_to_fm(matrix=matrix_vxc(i)%matrix, fm=fm_XC_ao)
     455           6 :             CALL cp_fm_set_all(fm_XC_mo(i), 0.0_dp)
     456             : 
     457             :             ! First index
     458             :             CALL parallel_gemm('T', 'N', dimen, dimen, dimen, 1.0_dp, &
     459           6 :                                mo_coeff(i), fm_XC_ao, 0.0_dp, fm_XC_ao_mo)
     460             : 
     461             :             ! Second index
     462             :             CALL parallel_gemm('N', 'N', dimen, dimen, dimen, 1.0_dp, &
     463          10 :                                fm_XC_ao_mo, mo_coeff(i), 1.0_dp, fm_XC_mo(i))
     464             : 
     465             :          END DO
     466             : 
     467          10 :          DO i = 1, SIZE(matrix_vxc)
     468           6 :             CALL dbcsr_release(matrix_vxc(i)%matrix)
     469          10 :             DEALLOCATE (matrix_vxc(i)%matrix)
     470             :          END DO
     471           4 :          DEALLOCATE (matrix_vxc)
     472             :       END IF
     473             : 
     474           4 :       CALL timestop(handle)
     475             : 
     476           4 :    END SUBROUTINE xc_contribution
     477             : 
     478             : ! **************************************************************************************************
     479             : !> \brief ...
     480             : !> \param fm_F_mo ...
     481             : !> \param eigenval ...
     482             : !> \param dimen ...
     483             : !> \param homo ...
     484             : !> \param para_env ...
     485             : !> \param blacs_env ...
     486             : !> \param rse_corr ...
     487             : ! **************************************************************************************************
     488           4 :    SUBROUTINE non_diag_rse(fm_F_mo, eigenval, dimen, homo, para_env, &
     489             :                            blacs_env, rse_corr)
     490             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: fm_F_mo
     491             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: Eigenval
     492             :       INTEGER, INTENT(IN)                                :: dimen
     493             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: homo
     494             :       TYPE(mp_para_env_type), INTENT(IN), POINTER        :: para_env
     495             :       TYPE(cp_blacs_env_type), INTENT(IN), POINTER       :: blacs_env
     496             :       REAL(KIND=dp), INTENT(OUT)                         :: rse_corr
     497             : 
     498             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'non_diag_rse'
     499             : 
     500             :       INTEGER                                            :: handle, i_global, iiB, ispin, j_global, &
     501             :                                                             jjB, ncol_local, nrow_local, nspins, &
     502             :                                                             virtual
     503           4 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     504             :       REAL(KIND=dp)                                      :: corr
     505           4 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eig_o, eig_semi_can, eig_v
     506             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     507             :       TYPE(cp_fm_type)                                   :: fm_F_oo, fm_F_ov, fm_F_vv, fm_O, fm_tmp, &
     508             :                                                             fm_U
     509             : 
     510           4 :       CALL timeset(routineN, handle)
     511             : 
     512           4 :       nspins = SIZE(fm_f_mo)
     513             : 
     514          10 :       DO ispin = 1, nspins
     515             :          ! Add eigenvalues on the diagonal
     516             :          CALL cp_fm_get_info(matrix=fm_F_mo(ispin), &
     517             :                              nrow_local=nrow_local, &
     518             :                              ncol_local=ncol_local, &
     519             :                              row_indices=row_indices, &
     520           6 :                              col_indices=col_indices)
     521             : 
     522             : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
     523          10 : !$OMP             SHARED(ncol_local,nrow_local,col_indices,row_indices,fm_F_mo,eigenval,ispin)
     524             :          DO jjB = 1, ncol_local
     525             :             j_global = col_indices(jjB)
     526             :             DO iiB = 1, nrow_local
     527             :                i_global = row_indices(iiB)
     528             :                IF (i_global .EQ. j_global) fm_F_mo(ispin)%local_data(iib, jjb) = &
     529             :                   fm_F_mo(ispin)%local_data(iib, jjb) + eigenval(i_global, ispin)
     530             :             END DO
     531             :          END DO
     532             : !$OMP END PARALLEL DO
     533             :       END DO
     534             : 
     535           4 :       rse_corr = 0.0_dp
     536             : 
     537          10 :       DO ispin = 1, nspins
     538             :          ! Create the occupied-occupied and virtual-virtual blocks, eigenvectors
     539           6 :          NULLIFY (fm_struct_tmp)
     540             :          CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
     541           6 :                                   nrow_global=homo(ispin), ncol_global=homo(ispin))
     542           6 :          CALL cp_fm_create(fm_F_oo, fm_struct_tmp, name="F_oo")
     543           6 :          CALL cp_fm_create(fm_O, fm_struct_tmp, name="O")
     544           6 :          CALL cp_fm_set_all(fm_F_oo, 0.0_dp)
     545           6 :          CALL cp_fm_set_all(fm_O, 0.0_dp)
     546           6 :          CALL cp_fm_struct_release(fm_struct_tmp)
     547             : 
     548             :          CALL cp_fm_to_fm_submat(msource=fm_F_mo(ispin), mtarget=fm_F_oo, &
     549             :                                  nrow=homo(ispin), ncol=homo(ispin), &
     550             :                                  s_firstrow=1, s_firstcol=1, &
     551           6 :                                  t_firstrow=1, t_firstcol=1)
     552           6 :          virtual = dimen - homo(ispin)
     553           6 :          NULLIFY (fm_struct_tmp)
     554             :          CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
     555           6 :                                   nrow_global=virtual, ncol_global=virtual)
     556           6 :          CALL cp_fm_create(fm_F_vv, fm_struct_tmp, name="F_vv")
     557           6 :          CALL cp_fm_create(fm_U, fm_struct_tmp, name="U")
     558           6 :          CALL cp_fm_set_all(fm_F_vv, 0.0_dp)
     559           6 :          CALL cp_fm_set_all(fm_U, 0.0_dp)
     560           6 :          CALL cp_fm_struct_release(fm_struct_tmp)
     561             : 
     562             :          CALL cp_fm_to_fm_submat(msource=fm_F_mo(ispin), mtarget=fm_F_vv, &
     563             :                                  nrow=virtual, ncol=virtual, &
     564             :                                  s_firstrow=homo(ispin) + 1, s_firstcol=homo(ispin) + 1, &
     565           6 :                                  t_firstrow=1, t_firstcol=1)
     566             : 
     567             :          ! Diagonalize occupied-occupied and virtual-virtual matrices
     568          18 :          ALLOCATE (eig_o(homo(ispin)))
     569          18 :          ALLOCATE (eig_v(virtual))
     570          58 :          eig_v = 0.0_dp
     571          14 :          eig_o = 0.0_dp
     572           6 :          CALL choose_eigv_solver(fm_F_oo, fm_O, eig_o)
     573           6 :          CALL choose_eigv_solver(fm_F_vv, fm_U, eig_v)
     574             : 
     575             :          ! Collect the eigenvalues to one array
     576          18 :          ALLOCATE (eig_semi_can(dimen))
     577          66 :          eig_semi_can = 0.0_dp
     578          14 :          eig_semi_can(1:homo(ispin)) = eig_o(:)
     579          58 :          eig_semi_can(homo(ispin) + 1:dimen) = eig_v(:)
     580             : 
     581             :          ! Create occupied-virtual block
     582           6 :          NULLIFY (fm_struct_tmp)
     583             :          CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
     584           6 :                                   nrow_global=homo(ispin), ncol_global=virtual)
     585           6 :          CALL cp_fm_create(fm_F_ov, fm_struct_tmp, name="F_ov")
     586           6 :          CALL cp_fm_create(fm_tmp, fm_struct_tmp, name="tmp")
     587           6 :          CALL cp_fm_set_all(fm_F_ov, 0.0_dp)
     588           6 :          CALL cp_fm_set_all(fm_tmp, 0.0_dp)
     589           6 :          CALL cp_fm_struct_release(fm_struct_tmp)
     590             : 
     591             :          CALL cp_fm_to_fm_submat(msource=fm_F_mo(ispin), mtarget=fm_F_ov, &
     592             :                                  nrow=homo(ispin), ncol=virtual, &
     593             :                                  s_firstrow=1, s_firstcol=homo(ispin) + 1, &
     594           6 :                                  t_firstrow=1, t_firstcol=1)
     595             : 
     596             :          CALL parallel_gemm(transa='T', transb='N', m=homo(ispin), n=virtual, k=homo(ispin), alpha=1.0_dp, &
     597           6 :                             matrix_a=fm_O, matrix_b=fm_F_ov, beta=0.0_dp, matrix_c=fm_tmp)
     598             : 
     599             :          CALL parallel_gemm(transa='N', transb='N', m=homo(ispin), n=virtual, k=virtual, alpha=1.0_dp, &
     600           6 :                             matrix_a=fm_tmp, matrix_b=fm_U, beta=0.0_dp, matrix_c=fm_F_ov)
     601             : 
     602             :          ! Compute the correction
     603             :          CALL cp_fm_get_info(matrix=fm_F_ov, &
     604             :                              nrow_local=nrow_local, &
     605             :                              ncol_local=ncol_local, &
     606             :                              row_indices=row_indices, &
     607           6 :                              col_indices=col_indices)
     608           6 :          corr = 0.0_dp
     609             : !$OMP    PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
     610             : !$OMP             REDUCTION(+:corr) &
     611           6 : !$OMP                SHARED(ncol_local,nrow_local,col_indices,row_indices,fm_F_ov,eig_semi_can,homo,ispin)
     612             :          DO jjB = 1, ncol_local
     613             :             j_global = col_indices(jjB)
     614             :             DO iiB = 1, nrow_local
     615             :                i_global = row_indices(iiB)
     616             :                corr = corr + fm_F_ov%local_data(iib, jjb)**2.0_dp/ &
     617             :                       (eig_semi_can(i_global) - eig_semi_can(j_global + homo(ispin)))
     618             :             END DO
     619             :          END DO
     620             : !$OMP    END PARALLEL DO
     621             : 
     622           6 :          rse_corr = rse_corr + corr
     623             : 
     624             :          ! Release
     625           6 :          DEALLOCATE (eig_semi_can)
     626           6 :          DEALLOCATE (eig_o)
     627           6 :          DEALLOCATE (eig_v)
     628             : 
     629           6 :          CALL cp_fm_release(fm_F_ov)
     630           6 :          CALL cp_fm_release(fm_F_oo)
     631           6 :          CALL cp_fm_release(fm_F_vv)
     632           6 :          CALL cp_fm_release(fm_U)
     633           6 :          CALL cp_fm_release(fm_O)
     634          16 :          CALL cp_fm_release(fm_tmp)
     635             : 
     636             :       END DO
     637             : 
     638           4 :       CALL para_env%sum(rse_corr)
     639             : 
     640           4 :       CALL timestop(handle)
     641             : 
     642           8 :    END SUBROUTINE non_diag_rse
     643             : 
     644             : END MODULE rpa_rse

Generated by: LCOV version 1.15