LCOV - code coverage report
Current view: top level - src - rpa_gw_im_time_util.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:2fce0f8) Lines: 349 352 99.1 %
Date: 2024-12-21 06:28:57 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 Utility routines for GW with imaginary time
      10             : !> \par History
      11             : !>      06.2019 split from rpa_im_time.F [Frederick Stein]
      12             : ! **************************************************************************************************
      13             : MODULE rpa_gw_im_time_util
      14             : 
      15             :    USE atomic_kind_types,               ONLY: atomic_kind_type
      16             :    USE basis_set_types,                 ONLY: gto_basis_set_p_type
      17             :    USE cell_types,                      ONLY: cell_type,&
      18             :                                               pbc
      19             :    USE cp_dbcsr_api,                    ONLY: &
      20             :         dbcsr_add_on_diag, dbcsr_create, dbcsr_distribution_get, dbcsr_distribution_new, &
      21             :         dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_filter, dbcsr_get_diag, &
      22             :         dbcsr_get_info, dbcsr_get_stored_coordinates, dbcsr_init_p, dbcsr_iterator_blocks_left, &
      23             :         dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
      24             :         dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_release_p, dbcsr_reserve_all_blocks, &
      25             :         dbcsr_set_diag, dbcsr_type, dbcsr_type_no_symmetry
      26             :    USE cp_dbcsr_operations,             ONLY: copy_fm_to_dbcsr,&
      27             :                                               cp_dbcsr_m_by_n_from_row_template
      28             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      29             :                                               cp_fm_release,&
      30             :                                               cp_fm_set_element,&
      31             :                                               cp_fm_to_fm,&
      32             :                                               cp_fm_type
      33             :    USE dbt_api,                         ONLY: &
      34             :         dbt_contract, dbt_copy, dbt_copy_matrix_to_tensor, dbt_create, dbt_default_distvec, &
      35             :         dbt_destroy, dbt_get_info, dbt_pgrid_create, dbt_pgrid_destroy, dbt_pgrid_type, dbt_type
      36             :    USE hfx_types,                       ONLY: alloc_containers,&
      37             :                                               block_ind_type,&
      38             :                                               hfx_compression_type
      39             :    USE kinds,                           ONLY: dp,&
      40             :                                               int_8
      41             :    USE mathconstants,                   ONLY: twopi
      42             :    USE message_passing,                 ONLY: mp_dims_create,&
      43             :                                               mp_para_env_type,&
      44             :                                               mp_request_type
      45             :    USE mp2_types,                       ONLY: integ_mat_buffer_type
      46             :    USE particle_methods,                ONLY: get_particle_set
      47             :    USE particle_types,                  ONLY: particle_type
      48             :    USE qs_environment_types,            ONLY: get_qs_env,&
      49             :                                               qs_environment_type
      50             :    USE qs_integral_utils,               ONLY: basis_set_list_setup
      51             :    USE qs_kind_types,                   ONLY: qs_kind_type
      52             :    USE qs_tensors,                      ONLY: compress_tensor,&
      53             :                                               decompress_tensor,&
      54             :                                               get_tensor_occupancy
      55             :    USE qs_tensors_types,                ONLY: create_2c_tensor,&
      56             :                                               create_3c_tensor,&
      57             :                                               pgf_block_sizes,&
      58             :                                               split_block_sizes
      59             :    USE rpa_communication,               ONLY: communicate_buffer
      60             : #include "./base/base_uses.f90"
      61             : 
      62             :    IMPLICIT NONE
      63             : 
      64             :    PRIVATE
      65             : 
      66             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_gw_im_time_util'
      67             : 
      68             :    PUBLIC :: get_tensor_3c_overl_int_gw, compute_weight_re_im, get_atom_index_from_basis_function_index
      69             : 
      70             : CONTAINS
      71             : 
      72             : ! **************************************************************************************************
      73             : !> \brief ...
      74             : !> \param t_3c_overl_int ...
      75             : !> \param t_3c_O_compressed ...
      76             : !> \param t_3c_O_ind ...
      77             : !> \param t_3c_overl_int_ao_mo ...
      78             : !> \param t_3c_O_mo_compressed ...
      79             : !> \param t_3c_O_mo_ind ...
      80             : !> \param t_3c_overl_int_gw_RI ...
      81             : !> \param t_3c_overl_int_gw_AO ...
      82             : !> \param starts_array_mc ...
      83             : !> \param ends_array_mc ...
      84             : !> \param mo_coeff ...
      85             : !> \param matrix_s ...
      86             : !> \param gw_corr_lev_occ ...
      87             : !> \param gw_corr_lev_virt ...
      88             : !> \param homo ...
      89             : !> \param nmo ...
      90             : !> \param para_env ...
      91             : !> \param do_ic_model ...
      92             : !> \param t_3c_overl_nnP_ic ...
      93             : !> \param t_3c_overl_nnP_ic_reflected ...
      94             : !> \param qs_env ...
      95             : !> \param unit_nr ...
      96             : !> \param do_alpha ...
      97             : ! **************************************************************************************************
      98          56 :    SUBROUTINE get_tensor_3c_overl_int_gw(t_3c_overl_int, t_3c_O_compressed, t_3c_O_ind, &
      99             :                                          t_3c_overl_int_ao_mo, t_3c_O_mo_compressed, t_3c_O_mo_ind, &
     100             :                                          t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, &
     101         112 :                                          starts_array_mc, ends_array_mc, &
     102             :                                          mo_coeff, matrix_s, &
     103             :                                          gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, &
     104             :                                          para_env, &
     105             :                                          do_ic_model, &
     106             :                                          t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, &
     107             :                                          qs_env, unit_nr, do_alpha)
     108             : 
     109             :       TYPE(dbt_type), DIMENSION(:, :)                    :: t_3c_overl_int
     110             :       TYPE(hfx_compression_type), DIMENSION(:, :, :)     :: t_3c_O_compressed
     111             :       TYPE(block_ind_type), DIMENSION(:, :, :)           :: t_3c_O_ind
     112             :       TYPE(dbt_type)                                     :: t_3c_overl_int_ao_mo
     113             :       TYPE(hfx_compression_type)                         :: t_3c_O_mo_compressed
     114             :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: t_3c_O_mo_ind
     115             :       TYPE(dbt_type)                                     :: t_3c_overl_int_gw_RI, &
     116             :                                                             t_3c_overl_int_gw_AO
     117             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: starts_array_mc, ends_array_mc
     118             :       TYPE(cp_fm_type), INTENT(IN)                       :: mo_coeff
     119             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     120             :       INTEGER, INTENT(IN)                                :: gw_corr_lev_occ, gw_corr_lev_virt, homo, &
     121             :                                                             nmo
     122             :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
     123             :       LOGICAL, INTENT(IN)                                :: do_ic_model
     124             :       TYPE(dbt_type)                                     :: t_3c_overl_nnP_ic, &
     125             :                                                             t_3c_overl_nnP_ic_reflected
     126             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     127             :       INTEGER, INTENT(IN)                                :: unit_nr
     128             :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_alpha
     129             : 
     130             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'get_tensor_3c_overl_int_gw'
     131             : 
     132             :       INTEGER :: cut_memory, handle, i_mem, icol_global, imo, irow_global, min_bsize, &
     133             :          min_bsize_mo, nkind, nmo_blk_gw, npcols, nprows, size_MO, unit_nr_prv
     134             :       INTEGER(int_8)                                     :: nze
     135         112 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: dist1, dist2, dist3, sizes_AO, &
     136          56 :                                                             sizes_AO_split, sizes_MO, sizes_MO_1, &
     137         112 :                                                             sizes_RI, sizes_RI_split, tmp
     138             :       INTEGER, DIMENSION(2)                              :: pdims_2d
     139             :       INTEGER, DIMENSION(2, 1)                           :: bounds
     140             :       INTEGER, DIMENSION(2, 3)                           :: ibounds
     141             :       INTEGER, DIMENSION(3)                              :: bounds_3c, pdims
     142         112 :       INTEGER, DIMENSION(:), POINTER                     :: distp_1, distp_2, sizes_MO_blocked, &
     143          56 :                                                             sizes_MO_p1, sizes_MO_p2
     144             :       LOGICAL                                            :: memory_info, my_do_alpha
     145             :       REAL(dp)                                           :: compression_factor, memory_3c, occ
     146          56 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: norm
     147          56 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     148             :       TYPE(cp_fm_type)                                   :: fm_mat_mo_coeff_gw
     149             :       TYPE(dbcsr_distribution_type)                      :: dist, dist_templ
     150             :       TYPE(dbcsr_type)                                   :: mat_mo_coeff_gw_reflected_norm, &
     151             :                                                             mat_norm, mat_norm_diag, mat_work
     152             :       TYPE(dbcsr_type), POINTER                          :: mat_mo_coeff_gw, &
     153             :                                                             mat_mo_coeff_gw_reflected
     154         504 :       TYPE(dbt_pgrid_type)                               :: pgrid_2d, pgrid_AO, pgrid_ic, pgrid_MO
     155         728 :       TYPE(dbt_type)                                     :: mo_coeff_gw_t, mo_coeff_gw_t_tmp, &
     156         392 :                                                             t_3c_overl_int_ao_ao, &
     157         392 :                                                             t_3c_overl_int_mo_ao, &
     158         392 :                                                             t_3c_overl_int_mo_mo
     159          56 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_ao, basis_set_ri_aux
     160          56 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     161          56 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     162             : 
     163          56 :       memory_info = qs_env%mp2_env%ri_rpa_im_time%memory_info
     164          56 :       IF (memory_info) THEN
     165           0 :          unit_nr_prv = unit_nr
     166             :       ELSE
     167          56 :          unit_nr_prv = 0
     168             :       END IF
     169             : 
     170          56 :       my_do_alpha = .FALSE.
     171          56 :       IF (PRESENT(do_alpha)) my_do_alpha = do_alpha
     172             : 
     173          56 :       CALL timeset(routineN, handle)
     174             : 
     175          56 :       CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, particle_set=particle_set, atomic_kind_set=atomic_kind_set)
     176             : 
     177          56 :       CALL cp_fm_create(fm_mat_mo_coeff_gw, mo_coeff%matrix_struct)
     178          56 :       CALL cp_fm_to_fm(mo_coeff, fm_mat_mo_coeff_gw)
     179             : 
     180             :       ! set MO coeffs to zero where
     181        1416 :       DO irow_global = 1, nmo
     182        3142 :          DO icol_global = 1, homo - gw_corr_lev_occ
     183        3142 :             CALL cp_fm_set_element(fm_mat_mo_coeff_gw, irow_global, icol_global, 0.0_dp)
     184             :          END DO
     185       32814 :          DO icol_global = homo + gw_corr_lev_virt + 1, nmo
     186       32758 :             CALL cp_fm_set_element(fm_mat_mo_coeff_gw, irow_global, icol_global, 0.0_dp)
     187             :          END DO
     188             :       END DO
     189             : 
     190          56 :       NULLIFY (mat_mo_coeff_gw)
     191          56 :       CALL dbcsr_init_p(mat_mo_coeff_gw)
     192             : 
     193             :       CALL cp_dbcsr_m_by_n_from_row_template(mat_mo_coeff_gw, template=matrix_s(1)%matrix, n=nmo, &
     194          56 :                                              sym=dbcsr_type_no_symmetry)
     195             : 
     196             :       CALL copy_fm_to_dbcsr(fm_mat_mo_coeff_gw, &
     197             :                             mat_mo_coeff_gw, &
     198          56 :                             keep_sparsity=.FALSE.)
     199             : 
     200             :       ! just remove the blocks which have been set to zero
     201          56 :       CALL dbcsr_filter(mat_mo_coeff_gw, 1.0E-20_dp)
     202             : 
     203          56 :       min_bsize = qs_env%mp2_env%ri_rpa_im_time%min_bsize
     204          56 :       min_bsize_mo = qs_env%mp2_env%ri_rpa_im_time%min_bsize_mo
     205             : 
     206         112 :       CALL split_block_sizes([gw_corr_lev_occ + gw_corr_lev_virt], sizes_MO, min_bsize_mo)
     207         168 :       ALLOCATE (sizes_MO_1(nmo))
     208        1416 :       sizes_MO_1(:) = 1
     209             : 
     210          56 :       nmo_blk_gw = SIZE(sizes_MO)
     211          56 :       CALL move_alloc(sizes_MO, tmp)
     212         168 :       ALLOCATE (sizes_MO(nmo_blk_gw + 2))
     213          56 :       sizes_MO(1) = homo - gw_corr_lev_occ
     214         112 :       sizes_MO(2:SIZE(tmp) + 1) = tmp(:)
     215          56 :       sizes_MO(SIZE(tmp) + 2) = nmo - (homo + gw_corr_lev_virt)
     216             : 
     217         440 :       ALLOCATE (basis_set_ri_aux(nkind), basis_set_ao(nkind))
     218          56 :       CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set)
     219          56 :       CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_RI, basis=basis_set_ri_aux)
     220          56 :       CALL basis_set_list_setup(basis_set_ao, "ORB", qs_kind_set)
     221          56 :       CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_AO, basis=basis_set_ao)
     222             : 
     223             :       CALL pgf_block_sizes(atomic_kind_set, basis_set_ao, min_bsize, sizes_AO_split)
     224             :       CALL pgf_block_sizes(atomic_kind_set, basis_set_ri_aux, min_bsize, sizes_RI_split)
     225             : 
     226          56 :       DEALLOCATE (basis_set_ao, basis_set_ri_aux)
     227             : 
     228          56 :       pdims = 0
     229             :       CALL dbt_pgrid_create(para_env, pdims, pgrid_AO, &
     230         224 :                             tensor_dims=[SIZE(sizes_RI_split), SIZE(sizes_AO_split), SIZE(sizes_AO_split)])
     231             : 
     232          56 :       pdims_2d = 0
     233          56 :       CALL mp_dims_create(para_env%num_pe, pdims_2d)
     234             : 
     235             :       ! we iterate over MO blocks for saving memory during contraction, thus we should not parallelize over MO dimension
     236         224 :       pdims = [pdims_2d(1), pdims_2d(2), 1]
     237             :       CALL dbt_pgrid_create(para_env, pdims, pgrid_MO, &
     238         224 :                             tensor_dims=[SIZE(sizes_RI_split), SIZE(sizes_AO_split), 1])
     239             : 
     240          56 :       pdims_2d = 0
     241             :       CALL dbt_pgrid_create(para_env, pdims_2d, pgrid_2d, &
     242         168 :                             tensor_dims=[SIZE(sizes_AO_split), nmo])
     243             : 
     244             :       CALL create_3c_tensor(t_3c_overl_int_ao_ao, dist1, dist2, dist3, pgrid_AO, &
     245          56 :                             sizes_RI_split, sizes_AO_split, sizes_AO_split, [1, 2], [3], name="(RI AO | AO)")
     246          56 :       DEALLOCATE (dist1, dist2, dist3)
     247             : 
     248          56 :       IF (.NOT. qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
     249             :          CALL create_3c_tensor(t_3c_overl_int_ao_mo, dist1, dist2, dist3, pgrid_AO, &
     250          34 :                                sizes_RI_split, sizes_AO_split, sizes_MO_1, [1, 2], [3], name="(RI AO | MO)")
     251          34 :          DEALLOCATE (dist1, dist2, dist3)
     252             :       END IF
     253             : 
     254             :       CALL create_3c_tensor(t_3c_overl_int_gw_RI, dist1, dist2, dist3, pgrid_MO, &
     255          56 :                             sizes_RI_split, sizes_AO_split, sizes_MO, [1], [2, 3], name="(RI | AO MO)")
     256          56 :       DEALLOCATE (dist1, dist2, dist3)
     257             : 
     258             :       CALL create_3c_tensor(t_3c_overl_int_gw_AO, dist1, dist2, dist3, pgrid_MO, &
     259          56 :                             sizes_AO_split, sizes_RI_split, sizes_MO, [1], [2, 3], name="(AO | RI MO)")
     260          56 :       DEALLOCATE (dist1, dist2, dist3)
     261             : 
     262          56 :       CALL dbt_pgrid_destroy(pgrid_AO)
     263          56 :       CALL dbt_pgrid_destroy(pgrid_MO)
     264             : 
     265             :       CALL create_2c_tensor(mo_coeff_gw_t, dist1, dist2, pgrid_2d, sizes_AO_split, sizes_MO_1, name="(AO|MO)")
     266          56 :       DEALLOCATE (dist1, dist2)
     267          56 :       CALL dbt_pgrid_destroy(pgrid_2d)
     268             : 
     269          56 :       CALL dbt_create(mat_mo_coeff_gw, mo_coeff_gw_t_tmp, name="MO coeffs")
     270          56 :       CALL dbt_copy_matrix_to_tensor(mat_mo_coeff_gw, mo_coeff_gw_t_tmp)
     271             : 
     272          56 :       CALL dbt_copy(mo_coeff_gw_t_tmp, mo_coeff_gw_t)
     273             : 
     274          56 :       bounds(1, 1) = homo - gw_corr_lev_occ + 1
     275          56 :       bounds(2, 1) = homo + gw_corr_lev_virt
     276             : 
     277          56 :       CALL dbt_get_info(t_3c_overl_int_ao_ao, nfull_total=bounds_3c)
     278             : 
     279         168 :       ibounds(:, 1) = [1, bounds_3c(1)]
     280         168 :       ibounds(:, 3) = [1, bounds_3c(3)]
     281             : 
     282          56 :       cut_memory = SIZE(starts_array_mc)
     283             : 
     284          56 :       IF (.NOT. qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
     285          68 :          DO i_mem = 1, cut_memory
     286             :             CALL decompress_tensor(t_3c_overl_int(1, 1), t_3c_O_ind(1, 1, i_mem)%ind, t_3c_O_compressed(1, 1, i_mem), &
     287          34 :                                    qs_env%mp2_env%ri_rpa_im_time%eps_compress)
     288             : 
     289         102 :             ibounds(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
     290             : 
     291          34 :             CALL dbt_copy(t_3c_overl_int(1, 1), t_3c_overl_int_ao_ao, move_data=.TRUE.)
     292             : 
     293             :             CALL dbt_contract(1.0_dp, mo_coeff_gw_t, t_3c_overl_int_ao_ao, 1.0_dp, &
     294             :                               t_3c_overl_int_ao_mo, contract_1=[1], notcontract_1=[2], &
     295             :                               contract_2=[3], notcontract_2=[1, 2], map_1=[3], map_2=[1, 2], &
     296          68 :                               bounds_2=ibounds, move_data=.FALSE., unit_nr=unit_nr_prv)
     297             : 
     298             :          END DO
     299             :       END IF
     300             : 
     301          56 :       CALL cp_fm_release(fm_mat_mo_coeff_gw)
     302             : 
     303          56 :       IF (do_ic_model) THEN
     304           2 :          pdims = 0
     305             :          CALL dbt_pgrid_create(para_env, pdims, pgrid_ic, &
     306           8 :                                tensor_dims=[SIZE(sizes_RI_split), nmo, nmo])
     307             : 
     308             :          CALL create_3c_tensor(t_3c_overl_int_mo_ao, dist1, dist2, dist3, pgrid_ic, &
     309           2 :                                sizes_RI_split, sizes_MO_1, sizes_AO_split, [1, 2], [3], name="(RI MO | AO)")
     310           2 :          DEALLOCATE (dist1, dist2, dist3)
     311             :          CALL create_3c_tensor(t_3c_overl_int_mo_mo, dist1, dist2, dist3, pgrid_ic, &
     312           2 :                                sizes_RI_split, sizes_MO_1, sizes_MO_1, [1, 2], [3], name="(RI MO | MO)")
     313           2 :          DEALLOCATE (dist1, dist2, dist3)
     314           2 :          CALL dbt_create(t_3c_overl_int_mo_mo, t_3c_overl_nnP_ic)
     315             :          CALL create_3c_tensor(t_3c_overl_nnP_ic_reflected, dist1, dist2, dist3, pgrid_ic, &
     316           2 :                                sizes_RI_split, sizes_MO_1, sizes_MO_1, [1], [2, 3], name="(RI | MO MO)")
     317           2 :          DEALLOCATE (dist1, dist2, dist3)
     318             : 
     319           2 :          CALL dbt_pgrid_destroy(pgrid_ic)
     320             : 
     321           2 :          CALL dbt_copy(t_3c_overl_int_ao_mo, t_3c_overl_int_mo_ao, order=[1, 3, 2])
     322             :          CALL dbt_contract(1.0_dp, mo_coeff_gw_t, t_3c_overl_int_mo_ao, 0.0_dp, &
     323             :                            t_3c_overl_int_mo_mo, contract_1=[1], notcontract_1=[2], &
     324             :                            contract_2=[3], notcontract_2=[1, 2], map_1=[3], map_2=[1, 2], &
     325           2 :                            bounds_2=bounds, move_data=.FALSE., unit_nr=unit_nr_prv)
     326           2 :          CALL dbt_copy(t_3c_overl_int_mo_mo, t_3c_overl_nnP_ic)
     327             : 
     328           2 :          NULLIFY (mat_mo_coeff_gw_reflected)
     329           2 :          CALL dbcsr_init_p(mat_mo_coeff_gw_reflected)
     330             : 
     331             :          CALL cp_dbcsr_m_by_n_from_row_template(mat_mo_coeff_gw_reflected, template=matrix_s(1)%matrix, n=nmo, &
     332           2 :                                                 sym=dbcsr_type_no_symmetry)
     333             : 
     334           2 :          CALL reflect_mat_row(mat_mo_coeff_gw_reflected, mat_mo_coeff_gw, para_env, qs_env, unit_nr, do_alpha=my_do_alpha)
     335             : 
     336             :          ! normalize reflected MOs (they are not properly normalized since high angular momentum basis functions
     337             :          ! of the image molecule are not exactly reflected at the image plane (sign problem in p_z function)
     338           2 :          CALL dbcsr_create(matrix=mat_work, template=mat_mo_coeff_gw_reflected, matrix_type=dbcsr_type_no_symmetry)
     339             : 
     340           2 :          CALL dbcsr_get_info(mat_work, distribution=dist_templ, nblkcols_total=size_MO, col_blk_size=sizes_MO_blocked)
     341             : 
     342           2 :          CALL dbcsr_distribution_get(dist_templ, nprows=nprows, npcols=npcols)
     343             : 
     344           8 :          ALLOCATE (distp_1(size_MO), distp_2(size_MO))
     345           2 :          CALL dbt_default_distvec(size_MO, nprows, sizes_MO_blocked, distp_1)
     346           2 :          CALL dbt_default_distvec(size_MO, npcols, sizes_MO_blocked, distp_2)
     347           2 :          CALL dbcsr_distribution_new(dist, template=dist_templ, row_dist=distp_1, col_dist=distp_2, reuse_arrays=.TRUE.)
     348             : 
     349           4 :          ALLOCATE (sizes_MO_p1(size_MO))
     350           4 :          ALLOCATE (sizes_MO_p2(size_MO))
     351          14 :          sizes_MO_p1(:) = sizes_MO_blocked
     352          14 :          sizes_MO_p2(:) = sizes_MO_blocked
     353             :          CALL dbcsr_create(mat_norm, "mo norm", dist, dbcsr_type_no_symmetry, sizes_MO_p1, sizes_MO_p2, &
     354           2 :                            reuse_arrays=.TRUE.)
     355           2 :          CALL dbcsr_distribution_release(dist)
     356             : 
     357           2 :          CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s(1)%matrix, mat_mo_coeff_gw_reflected, 0.0_dp, mat_work)
     358           2 :          CALL dbcsr_multiply("T", "N", 1.0_dp, mat_mo_coeff_gw_reflected, mat_work, 0.0_dp, mat_norm)
     359             : 
     360           2 :          CALL dbcsr_release(mat_work)
     361             : 
     362           6 :          ALLOCATE (norm(nmo))
     363         146 :          norm = 0.0_dp
     364             : 
     365           2 :          CALL dbcsr_get_diag(mat_norm, norm)
     366           2 :          CALL para_env%sum(norm)
     367             : 
     368          14 :          DO imo = bounds(1, 1), bounds(2, 1)
     369          14 :             norm(imo) = 1.0_dp/SQRT(norm(imo))
     370             :          END DO
     371             : 
     372           2 :          CALL dbcsr_create(mat_norm_diag, template=mat_norm)
     373           2 :          CALL dbcsr_release(mat_norm)
     374             : 
     375           2 :          CALL dbcsr_add_on_diag(mat_norm_diag, 1.0_dp)
     376             : 
     377           2 :          CALL dbcsr_set_diag(mat_norm_diag, norm)
     378             : 
     379           2 :          CALL dbcsr_create(mat_mo_coeff_gw_reflected_norm, template=mat_mo_coeff_gw_reflected)
     380           2 :          CALL dbcsr_multiply("N", "N", 1.0_dp, mat_mo_coeff_gw_reflected, mat_norm_diag, 0.0_dp, mat_mo_coeff_gw_reflected_norm)
     381           2 :          CALL dbcsr_release(mat_norm_diag)
     382             : 
     383           2 :          CALL dbcsr_filter(mat_mo_coeff_gw_reflected_norm, 1.0E-20_dp)
     384             : 
     385           2 :          CALL dbt_copy_matrix_to_tensor(mat_mo_coeff_gw_reflected_norm, mo_coeff_gw_t_tmp)
     386           2 :          CALL dbcsr_release(mat_mo_coeff_gw_reflected_norm)
     387           2 :          CALL dbt_copy(mo_coeff_gw_t_tmp, mo_coeff_gw_t)
     388             : 
     389             :          CALL dbt_contract(1.0_dp, mo_coeff_gw_t, t_3c_overl_int_ao_ao, 0.0_dp, &
     390             :                            t_3c_overl_int_ao_mo, contract_1=[1], notcontract_1=[2], &
     391             :                            contract_2=[3], notcontract_2=[1, 2], map_1=[3], map_2=[1, 2], &
     392           2 :                            bounds_2=bounds, move_data=.FALSE., unit_nr=unit_nr_prv)
     393             : 
     394           2 :          CALL dbt_copy(t_3c_overl_int_ao_mo, t_3c_overl_int_mo_ao, order=[1, 3, 2])
     395             :          CALL dbt_contract(1.0_dp, mo_coeff_gw_t, t_3c_overl_int_mo_ao, 0.0_dp, &
     396             :                            t_3c_overl_int_mo_mo, contract_1=[1], notcontract_1=[2], &
     397             :                            contract_2=[3], notcontract_2=[1, 2], map_1=[3], map_2=[1, 2], &
     398           2 :                            bounds_2=bounds, move_data=.FALSE., unit_nr=unit_nr_prv)
     399           2 :          CALL dbt_copy(t_3c_overl_int_mo_mo, t_3c_overl_nnP_ic_reflected)
     400           2 :          CALL dbt_destroy(t_3c_overl_int_mo_ao)
     401           2 :          CALL dbt_destroy(t_3c_overl_int_mo_mo)
     402             : 
     403           4 :          CALL dbcsr_release_p(mat_mo_coeff_gw_reflected)
     404             : 
     405             :       END IF
     406             : 
     407          56 :       IF (.NOT. qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
     408          34 :          CALL alloc_containers(t_3c_O_mo_compressed, 1)
     409          34 :          CALL get_tensor_occupancy(t_3c_overl_int_ao_mo, nze, occ)
     410          34 :          memory_3c = 0.0_dp
     411             : 
     412             :          CALL compress_tensor(t_3c_overl_int_ao_mo, t_3c_O_mo_ind, t_3c_O_mo_compressed, &
     413          34 :                               qs_env%mp2_env%ri_rpa_im_time%eps_compress, memory_3c)
     414             : 
     415          34 :          CALL para_env%sum(memory_3c)
     416          34 :          compression_factor = REAL(nze, dp)*1.0E-06*8.0_dp/memory_3c
     417             : 
     418          34 :          IF (unit_nr > 0) THEN
     419             :             WRITE (UNIT=unit_nr, FMT="((T3,A,T66,F11.2,A4))") &
     420          17 :                "MEMORY_INFO| Memory of MO-contracted tensor (compressed):", memory_3c, ' MiB'
     421             : 
     422             :             WRITE (UNIT=unit_nr, FMT="((T3,A,T60,F21.2))") &
     423          17 :                "MEMORY_INFO| Compression factor:                  ", compression_factor
     424             :          END IF
     425             :       END IF
     426             : 
     427          56 :       CALL dbcsr_release_p(mat_mo_coeff_gw)
     428             : 
     429          56 :       CALL dbt_destroy(t_3c_overl_int_ao_ao)
     430          56 :       CALL dbt_destroy(mo_coeff_gw_t)
     431          56 :       CALL dbt_destroy(mo_coeff_gw_t_tmp)
     432             : 
     433          56 :       CALL timestop(handle)
     434             : 
     435         448 :    END SUBROUTINE
     436             : 
     437             : ! **************************************************************************************************
     438             : !> \brief reflect from V = (A,B|B,A) to V_reflected = (B,A|A,B) where A belongs to the block of the molecule
     439             : !>        and B to the off diagonal block between molecule and image of the molecule
     440             : !> \param mat_reflected ...
     441             : !> \param mat_orig ...
     442             : !> \param para_env ...
     443             : !> \param qs_env ...
     444             : !> \param unit_nr ...
     445             : !> \param do_alpha ...
     446             : ! **************************************************************************************************
     447           2 :    SUBROUTINE reflect_mat_row(mat_reflected, mat_orig, para_env, qs_env, unit_nr, do_alpha)
     448             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: mat_reflected
     449             :       TYPE(dbcsr_type), INTENT(IN)                       :: mat_orig
     450             :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
     451             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     452             :       INTEGER, INTENT(IN)                                :: unit_nr
     453             :       LOGICAL, INTENT(IN)                                :: do_alpha
     454             : 
     455             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'reflect_mat_row'
     456             : 
     457             :       INTEGER :: block, block_size, col, col_rec, col_size, handle, i_atom, i_block, imepos, &
     458             :          j_atom, natom, nblkcols_total, nblkrows_total, offset, row, row_rec, row_reflected, &
     459             :          row_size
     460           2 :       INTEGER, ALLOCATABLE, DIMENSION(:) :: block_counter, entry_counter, image_atom, &
     461           2 :          num_blocks_rec, num_blocks_send, num_entries_rec, num_entries_send, sizes_rec, sizes_send
     462           2 :       INTEGER, DIMENSION(:), POINTER                     :: col_blk_sizes, row_blk_sizes
     463             :       LOGICAL                                            :: found_image_atom
     464             :       REAL(KIND=dp)                                      :: avg_z_dist, delta, eps_dist2, &
     465             :                                                             min_z_dist, ra(3), rb(3), sum_z, &
     466             :                                                             z_reflection
     467           2 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: data_block
     468             :       TYPE(cell_type), POINTER                           :: cell
     469             :       TYPE(dbcsr_iterator_type)                          :: iter
     470             :       TYPE(integ_mat_buffer_type), ALLOCATABLE, &
     471           2 :          DIMENSION(:)                                    :: buffer_rec, buffer_send
     472           2 :       TYPE(mp_request_type), DIMENSION(:, :), POINTER    :: req_array
     473           2 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     474             : 
     475           2 :       CALL timeset(routineN, handle)
     476             : 
     477           2 :       CALL dbcsr_reserve_all_blocks(mat_reflected)
     478             : 
     479             :       CALL get_qs_env(qs_env, cell=cell, &
     480           2 :                       particle_set=particle_set)
     481             : 
     482             :       ! first check, whether we have an image molecule
     483             :       CALL dbcsr_get_info(mat_orig, &
     484             :                           nblkrows_total=nblkrows_total, &
     485             :                           nblkcols_total=nblkcols_total, &
     486             :                           row_blk_size=row_blk_sizes, &
     487           2 :                           col_blk_size=col_blk_sizes)
     488             : 
     489           2 :       natom = SIZE(particle_set)
     490           2 :       CPASSERT(natom == nblkrows_total)
     491             : 
     492           2 :       eps_dist2 = qs_env%mp2_env%ri_g0w0%eps_dist
     493           2 :       eps_dist2 = eps_dist2*eps_dist2
     494             : 
     495           2 :       sum_z = 0.0_dp
     496             : 
     497          18 :       DO i_atom = 1, natom
     498             : 
     499          16 :          ra(:) = pbc(particle_set(i_atom)%r, cell)
     500             : 
     501          18 :          sum_z = sum_z + ra(3)
     502             : 
     503             :       END DO
     504             : 
     505           2 :       z_reflection = sum_z/REAL(natom, KIND=dp)
     506             : 
     507           2 :       sum_z = 0.0_dp
     508             : 
     509          18 :       DO i_atom = 1, natom
     510             : 
     511          16 :          ra(:) = pbc(particle_set(i_atom)%r, cell)
     512             : 
     513          18 :          sum_z = sum_z + ABS(ra(3) - z_reflection)
     514             : 
     515             :       END DO
     516             : 
     517           2 :       avg_z_dist = sum_z/REAL(natom, KIND=dp)
     518             : 
     519           2 :       min_z_dist = avg_z_dist
     520             : 
     521          18 :       DO i_atom = 1, natom
     522             : 
     523          16 :          ra(:) = pbc(particle_set(i_atom)%r, cell)
     524             : 
     525          18 :          IF (ABS(ra(3) - z_reflection) < min_z_dist) THEN
     526           2 :             min_z_dist = ABS(ra(3) - z_reflection)
     527             :          END IF
     528             : 
     529             :       END DO
     530             : 
     531           2 :       IF (unit_nr > 0 .AND. do_alpha) THEN
     532           1 :          WRITE (unit_nr, '(T3,A,T70,F9.2,A2)') 'IC_MODEL| Average distance of the molecule to the image plane:', &
     533           2 :             avg_z_dist*0.529_dp, ' A'
     534           1 :          WRITE (unit_nr, '(T3,A,T70,F9.2,A2)') 'IC_MODEL| Minimum distance of the molecule to the image plane:', &
     535           2 :             min_z_dist*0.529_dp, ' A'
     536             :       END IF
     537             : 
     538           6 :       ALLOCATE (image_atom(nblkrows_total))
     539          18 :       image_atom = 0
     540             : 
     541          18 :       DO i_atom = 1, natom
     542             : 
     543          16 :          found_image_atom = .FALSE.
     544             : 
     545          16 :          ra(:) = pbc(particle_set(i_atom)%r, cell)
     546             : 
     547         144 :          DO j_atom = 1, natom
     548             : 
     549         128 :             rb(:) = pbc(particle_set(j_atom)%r, cell)
     550             : 
     551         128 :             delta = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) + rb(3) - 2.0_dp*z_reflection)**2
     552             : 
     553             :             ! SQRT(delta) < eps_dist
     554         144 :             IF (delta < eps_dist2) THEN
     555             :                ! this CPASSERT ensures that there is at most one image atom for each atom
     556          16 :                CPASSERT(.NOT. found_image_atom)
     557          16 :                image_atom(i_atom) = j_atom
     558          16 :                found_image_atom = .TRUE.
     559             :                ! check whether we have the same basis at the image atom
     560             :                ! if this is wrong, check whether you have the same basis sets for the molecule and the image
     561          16 :                CPASSERT(row_blk_sizes(i_atom) == row_blk_sizes(j_atom))
     562             :             END IF
     563             : 
     564             :          END DO
     565             : 
     566             :          ! this CPASSERT ensures that there is at least one image atom for each atom
     567          18 :          CPASSERT(found_image_atom)
     568             : 
     569             :       END DO
     570             : 
     571          10 :       ALLOCATE (buffer_rec(0:para_env%num_pe - 1))
     572           8 :       ALLOCATE (buffer_send(0:para_env%num_pe - 1))
     573             : 
     574           6 :       ALLOCATE (num_entries_rec(0:para_env%num_pe - 1))
     575           4 :       ALLOCATE (num_blocks_rec(0:para_env%num_pe - 1))
     576           4 :       ALLOCATE (num_entries_send(0:para_env%num_pe - 1))
     577           4 :       ALLOCATE (num_blocks_send(0:para_env%num_pe - 1))
     578           6 :       num_entries_rec = 0
     579           6 :       num_blocks_rec = 0
     580           6 :       num_entries_send = 0
     581           6 :       num_blocks_send = 0
     582             : 
     583           2 :       CALL dbcsr_iterator_start(iter, mat_orig)
     584          10 :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     585             : 
     586             :          CALL dbcsr_iterator_next_block(iter, row, col, data_block, &
     587           8 :                                         row_size=row_size, col_size=col_size)
     588             : 
     589           8 :          row_reflected = image_atom(row)
     590             : 
     591           8 :          CALL dbcsr_get_stored_coordinates(mat_reflected, row_reflected, col, imepos)
     592             : 
     593           8 :          num_entries_send(imepos) = num_entries_send(imepos) + row_size*col_size
     594           8 :          num_blocks_send(imepos) = num_blocks_send(imepos) + 1
     595             : 
     596             :       END DO
     597             : 
     598           2 :       CALL dbcsr_iterator_stop(iter)
     599             : 
     600           2 :       IF (para_env%num_pe > 1) THEN
     601             : 
     602           6 :          ALLOCATE (sizes_rec(0:2*para_env%num_pe - 1))
     603           4 :          ALLOCATE (sizes_send(0:2*para_env%num_pe - 1))
     604             : 
     605           6 :          DO imepos = 0, para_env%num_pe - 1
     606             : 
     607           4 :             sizes_send(2*imepos) = num_entries_send(imepos)
     608           6 :             sizes_send(2*imepos + 1) = num_blocks_send(imepos)
     609             : 
     610             :          END DO
     611             : 
     612           2 :          CALL para_env%alltoall(sizes_send, sizes_rec, 2)
     613             : 
     614           6 :          DO imepos = 0, para_env%num_pe - 1
     615           4 :             num_entries_rec(imepos) = sizes_rec(2*imepos)
     616           6 :             num_blocks_rec(imepos) = sizes_rec(2*imepos + 1)
     617             :          END DO
     618             : 
     619           2 :          DEALLOCATE (sizes_rec, sizes_send)
     620             : 
     621             :       ELSE
     622             : 
     623           0 :          num_entries_rec(0) = num_entries_send(0)
     624           0 :          num_blocks_rec(0) = num_blocks_send(0)
     625             : 
     626             :       END IF
     627             : 
     628             :       ! allocate data message and corresponding indices
     629           6 :       DO imepos = 0, para_env%num_pe - 1
     630             : 
     631          10 :          ALLOCATE (buffer_rec(imepos)%msg(num_entries_rec(imepos)))
     632        2308 :          buffer_rec(imepos)%msg = 0.0_dp
     633             : 
     634          10 :          ALLOCATE (buffer_send(imepos)%msg(num_entries_send(imepos)))
     635        2308 :          buffer_send(imepos)%msg = 0.0_dp
     636             : 
     637          10 :          ALLOCATE (buffer_rec(imepos)%indx(num_blocks_rec(imepos), 3))
     638          40 :          buffer_rec(imepos)%indx = 0
     639             : 
     640          10 :          ALLOCATE (buffer_send(imepos)%indx(num_blocks_send(imepos), 3))
     641          42 :          buffer_send(imepos)%indx = 0
     642             : 
     643             :       END DO
     644             : 
     645           4 :       ALLOCATE (block_counter(0:para_env%num_pe - 1))
     646           6 :       block_counter(:) = 0
     647             : 
     648           4 :       ALLOCATE (entry_counter(0:para_env%num_pe - 1))
     649           6 :       entry_counter(:) = 0
     650             : 
     651           2 :       CALL dbcsr_iterator_start(iter, mat_orig)
     652          10 :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     653             : 
     654             :          CALL dbcsr_iterator_next_block(iter, row, col, data_block, &
     655           8 :                                         row_size=row_size, col_size=col_size)
     656             : 
     657           8 :          row_reflected = image_atom(row)
     658             : 
     659           8 :          CALL dbcsr_get_stored_coordinates(mat_reflected, row_reflected, col, imepos)
     660             : 
     661           8 :          block_size = row_size*col_size
     662             : 
     663           8 :          offset = entry_counter(imepos)
     664             : 
     665             :          buffer_send(imepos)%msg(offset + 1:offset + block_size) = &
     666        2320 :             RESHAPE(data_block(1:row_size, 1:col_size), (/block_size/))
     667             : 
     668           8 :          block = block_counter(imepos) + 1
     669             : 
     670           8 :          buffer_send(imepos)%indx(block, 1) = row_reflected
     671           8 :          buffer_send(imepos)%indx(block, 2) = col
     672           8 :          buffer_send(imepos)%indx(block, 3) = offset
     673             : 
     674           8 :          entry_counter(imepos) = entry_counter(imepos) + block_size
     675             : 
     676           8 :          block_counter(imepos) = block_counter(imepos) + 1
     677             : 
     678             :       END DO
     679             : 
     680           2 :       CALL dbcsr_iterator_stop(iter)
     681             : 
     682          30 :       ALLOCATE (req_array(1:para_env%num_pe, 4))
     683             : 
     684           2 :       CALL communicate_buffer(para_env, num_entries_rec, num_entries_send, buffer_rec, buffer_send, req_array)
     685             : 
     686           2 :       DEALLOCATE (req_array)
     687             : 
     688             :       ! fill the reflected matrix
     689           6 :       DO imepos = 0, para_env%num_pe - 1
     690             : 
     691          14 :          DO i_block = 1, num_blocks_rec(imepos)
     692             : 
     693           8 :             row_rec = buffer_rec(imepos)%indx(i_block, 1)
     694           8 :             col_rec = buffer_rec(imepos)%indx(i_block, 2)
     695             : 
     696           8 :             CALL dbcsr_iterator_start(iter, mat_reflected)
     697         104 :             DO WHILE (dbcsr_iterator_blocks_left(iter))
     698             : 
     699             :                CALL dbcsr_iterator_next_block(iter, row, col, data_block, &
     700          96 :                                               row_size=row_size, col_size=col_size)
     701             : 
     702         104 :                IF (row_rec == row .AND. col_rec == col) THEN
     703             : 
     704           8 :                   offset = buffer_rec(imepos)%indx(i_block, 3)
     705             : 
     706             :                   data_block(:, :) = RESHAPE(buffer_rec(imepos)%msg(offset + 1:offset + row_size*col_size), &
     707          24 :                                              (/row_size, col_size/))
     708             : 
     709             :                END IF
     710             : 
     711             :             END DO
     712             : 
     713          20 :             CALL dbcsr_iterator_stop(iter)
     714             : 
     715             :          END DO
     716             : 
     717             :       END DO
     718             : 
     719           6 :       DO imepos = 0, para_env%num_pe - 1
     720           4 :          DEALLOCATE (buffer_rec(imepos)%msg)
     721           4 :          DEALLOCATE (buffer_rec(imepos)%indx)
     722           4 :          DEALLOCATE (buffer_send(imepos)%msg)
     723           6 :          DEALLOCATE (buffer_send(imepos)%indx)
     724             :       END DO
     725             : 
     726          10 :       DEALLOCATE (buffer_rec, buffer_send)
     727           2 :       DEALLOCATE (block_counter, entry_counter)
     728           2 :       DEALLOCATE (num_entries_rec)
     729           2 :       DEALLOCATE (num_blocks_rec)
     730           2 :       DEALLOCATE (num_entries_send)
     731           2 :       DEALLOCATE (num_blocks_send)
     732             : 
     733           2 :       CALL timestop(handle)
     734             : 
     735          10 :    END SUBROUTINE
     736             : 
     737             : ! **************************************************************************************************
     738             : !> \brief ...
     739             : !> \param qs_env ...
     740             : !> \param atom_from_basis_index ...
     741             : !> \param basis_size ...
     742             : !> \param basis_type ...
     743             : !> \param first_bf_from_atom ...
     744             : ! **************************************************************************************************
     745        9696 :    SUBROUTINE get_atom_index_from_basis_function_index(qs_env, atom_from_basis_index, basis_size, &
     746             :                                                        basis_type, first_bf_from_atom)
     747             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     748             :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_from_basis_index
     749             :       INTEGER                                            :: basis_size
     750             :       CHARACTER(LEN=*)                                   :: basis_type
     751             :       INTEGER, ALLOCATABLE, DIMENSION(:), OPTIONAL       :: first_bf_from_atom
     752             : 
     753             :       INTEGER                                            :: iatom, LLL, natom, nkind
     754             :       INTEGER, DIMENSION(:), POINTER                     :: row_blk_end, row_blk_start
     755        9696 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set
     756        9696 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     757        9696 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     758             : 
     759        9696 :       NULLIFY (qs_kind_set, particle_set)
     760             :       CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, natom=natom, nkind=nkind, &
     761        9696 :                       particle_set=particle_set)
     762             : 
     763       29088 :       ALLOCATE (row_blk_start(natom))
     764       29088 :       ALLOCATE (row_blk_end(natom))
     765       46832 :       ALLOCATE (basis_set(nkind))
     766        9696 :       CALL basis_set_list_setup(basis_set, basis_type, qs_kind_set)
     767             :       CALL get_particle_set(particle_set, qs_kind_set, first_sgf=row_blk_start, last_sgf=row_blk_end, &
     768        9696 :                             basis=basis_set)
     769      254616 :       DO LLL = 1, basis_size
     770      950264 :          DO iatom = 1, natom
     771      940568 :             IF (LLL >= row_blk_start(iatom) .AND. LLL <= row_blk_end(iatom)) THEN
     772      244920 :                atom_from_basis_index(LLL) = iatom
     773             :             END IF
     774             :          END DO
     775             :       END DO
     776             : 
     777        9696 :       IF (PRESENT(first_bf_from_atom)) first_bf_from_atom(1:natom) = row_blk_start(:)
     778             : 
     779        9696 :       DEALLOCATE (basis_set)
     780        9696 :       DEALLOCATE (row_blk_start)
     781        9696 :       DEALLOCATE (row_blk_end)
     782             : 
     783        9696 :    END SUBROUTINE get_atom_index_from_basis_function_index
     784             : 
     785             : ! **************************************************************************************************
     786             : !> \brief ...
     787             : !> \param weight_re ...
     788             : !> \param weight_im ...
     789             : !> \param num_cells ...
     790             : !> \param iatom ...
     791             : !> \param jatom ...
     792             : !> \param xkp ...
     793             : !> \param wkp_W ...
     794             : !> \param cell ...
     795             : !> \param index_to_cell ...
     796             : !> \param hmat ...
     797             : !> \param particle_set ...
     798             : ! **************************************************************************************************
     799      315352 :    SUBROUTINE compute_weight_re_im(weight_re, weight_im, &
     800             :                                    num_cells, iatom, jatom, xkp, wkp_W, &
     801             :                                    cell, index_to_cell, hmat, particle_set)
     802             : 
     803             :       REAL(KIND=dp)                                      :: weight_re, weight_im
     804             :       INTEGER                                            :: num_cells, iatom, jatom
     805             :       REAL(KIND=dp), DIMENSION(3)                        :: xkp
     806             :       REAL(KIND=dp)                                      :: wkp_W
     807             :       TYPE(cell_type), POINTER                           :: cell
     808             :       INTEGER, DIMENSION(:, :), POINTER                  :: index_to_cell
     809             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat
     810             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     811             : 
     812             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_weight_re_im'
     813             : 
     814             :       INTEGER                                            :: handle, icell, n_equidistant_cells, &
     815             :                                                             xcell, ycell, zcell
     816             :       REAL(KIND=dp)                                      :: abs_rab_cell, abs_rab_cell_min, arg
     817             :       REAL(KIND=dp), DIMENSION(3)                        :: cell_vector, rab_cell_i
     818             : 
     819      315352 :       CALL timeset(routineN, handle)
     820             : 
     821      315352 :       weight_re = 0.0_dp
     822      315352 :       weight_im = 0.0_dp
     823             : 
     824      315352 :       abs_rab_cell_min = 1.0E10_dp
     825             : 
     826      315352 :       n_equidistant_cells = 0
     827             : 
     828     3091512 :       DO icell = 1, num_cells
     829             : 
     830     2776160 :          xcell = index_to_cell(1, icell)
     831     2776160 :          ycell = index_to_cell(2, icell)
     832     2776160 :          zcell = index_to_cell(3, icell)
     833             : 
     834    44418560 :          cell_vector(1:3) = MATMUL(hmat, REAL((/xcell, ycell, zcell/), dp))
     835             : 
     836             :          rab_cell_i(1:3) = pbc(particle_set(iatom)%r(1:3), cell) - &
     837    11104640 :                            (pbc(particle_set(jatom)%r(1:3), cell) + cell_vector(1:3))
     838             : 
     839     2776160 :          abs_rab_cell = SQRT(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2)
     840             : 
     841     3091512 :          IF (abs_rab_cell < abs_rab_cell_min) THEN
     842      478126 :             abs_rab_cell_min = abs_rab_cell
     843             :          END IF
     844             : 
     845             :       END DO
     846             : 
     847     3091512 :       DO icell = 1, num_cells
     848             : 
     849     2776160 :          xcell = index_to_cell(1, icell)
     850     2776160 :          ycell = index_to_cell(2, icell)
     851     2776160 :          zcell = index_to_cell(3, icell)
     852             : 
     853    44418560 :          cell_vector(1:3) = MATMUL(hmat, REAL((/xcell, ycell, zcell/), dp))
     854             : 
     855             :          rab_cell_i(1:3) = pbc(particle_set(iatom)%r(1:3), cell) - &
     856    11104640 :                            (pbc(particle_set(jatom)%r(1:3), cell) + cell_vector(1:3))
     857             : 
     858     2776160 :          abs_rab_cell = SQRT(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2)
     859             : 
     860     3091512 :          IF (abs_rab_cell < abs_rab_cell_min + 0.1_dp) THEN
     861             : 
     862      315352 :             arg = REAL(xcell, dp)*xkp(1) + REAL(ycell, dp)*xkp(2) + REAL(zcell, dp)*xkp(3)
     863             : 
     864      315352 :             weight_re = weight_re + wkp_W*COS(twopi*arg)
     865      315352 :             weight_im = weight_im + wkp_W*SIN(twopi*arg)
     866             : 
     867      315352 :             n_equidistant_cells = n_equidistant_cells + 1
     868             : 
     869             :          END IF
     870             : 
     871             :       END DO
     872             : 
     873      315352 :       weight_re = weight_re/REAL(n_equidistant_cells, KIND=dp)
     874      315352 :       weight_im = weight_im/REAL(n_equidistant_cells, KIND=dp)
     875             : 
     876      315352 :       CALL timestop(handle)
     877             : 
     878      315352 :    END SUBROUTINE compute_weight_re_im
     879             : 
     880             : END MODULE rpa_gw_im_time_util

Generated by: LCOV version 1.15