LCOV - code coverage report
Current view: top level - src - rpa_rse.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 97.2 % 218 212
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 4 4

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

Generated by: LCOV version 2.0-1