1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines for RPA with imaginary time 8!> \par History 9!> 10.2015 created [Jan Wilhelm] 10! ************************************************************************************************** 11MODULE rpa_im_time 12 13 USE cell_types, ONLY: cell_type,& 14 get_cell 15 USE cp_control_types, ONLY: dft_control_type 16 USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,& 17 dbcsr_allocate_matrix_set,& 18 dbcsr_deallocate_matrix_set 19 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,& 20 cp_fm_scale 21 USE cp_fm_struct, ONLY: cp_fm_struct_type 22 USE cp_fm_types, ONLY: cp_fm_create,& 23 cp_fm_get_info,& 24 cp_fm_release,& 25 cp_fm_to_fm,& 26 cp_fm_type 27 USE cp_gemm_interface, ONLY: cp_gemm 28 USE cp_para_types, ONLY: cp_para_env_type 29 USE dbcsr_api, ONLY: & 30 dbcsr_add, dbcsr_clear, dbcsr_copy, dbcsr_create, dbcsr_distribution_get, & 31 dbcsr_distribution_type, dbcsr_filter, dbcsr_get_info, dbcsr_init_p, dbcsr_p_type, & 32 dbcsr_release_p, dbcsr_reserve_all_blocks, dbcsr_scalar, dbcsr_scale, dbcsr_set, & 33 dbcsr_type, dbcsr_type_no_symmetry 34 USE dbcsr_tensor_api, ONLY: & 35 dbcsr_t_batched_contract_finalize, dbcsr_t_batched_contract_init, dbcsr_t_contract, & 36 dbcsr_t_copy, dbcsr_t_copy_matrix_to_tensor, dbcsr_t_copy_tensor_to_matrix, & 37 dbcsr_t_create, dbcsr_t_destroy, dbcsr_t_filter, dbcsr_t_get_info, dbcsr_t_nblks_total, & 38 dbcsr_t_nd_mp_comm, dbcsr_t_pgrid_destroy, dbcsr_t_pgrid_type, dbcsr_t_type 39 USE kinds, ONLY: dp,& 40 int_8 41 USE kpoint_types, ONLY: get_kpoint_info,& 42 kpoint_env_type,& 43 kpoint_type 44 USE machine, ONLY: m_walltime 45 USE mathconstants, ONLY: twopi 46 USE message_passing, ONLY: mp_sync 47 USE mp2_types, ONLY: mp2_type 48 USE qs_environment_types, ONLY: get_qs_env,& 49 qs_environment_type 50 USE qs_mo_types, ONLY: get_mo_set,& 51 mo_set_p_type,& 52 mo_set_type 53 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 54 USE qs_tensors, ONLY: get_tensor_occupancy,& 55 tensor_change_pgrid 56 USE qs_tensors_types, ONLY: create_2c_tensor 57#include "./base/base_uses.f90" 58 59 IMPLICIT NONE 60 61 PRIVATE 62 63 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_im_time' 64 65 PUBLIC :: compute_mat_P_omega, & 66 compute_transl_dm, & 67 init_cell_index_rpa, & 68 zero_mat_P_omega 69 70CONTAINS 71 72! ************************************************************************************************** 73!> \brief ... 74!> \param mat_P_omega ... 75!> \param fm_scaled_dm_occ_tau ... 76!> \param fm_scaled_dm_virt_tau ... 77!> \param fm_mo_coeff_occ ... 78!> \param fm_mo_coeff_virt ... 79!> \param fm_mo_coeff_occ_scaled ... 80!> \param fm_mo_coeff_virt_scaled ... 81!> \param mat_P_global ... 82!> \param matrix_s ... 83!> \param ispin ... 84!> \param t_3c_M ... 85!> \param t_3c_O ... 86!> \param starts_array_mc ... 87!> \param ends_array_mc ... 88!> \param starts_array_mc_block ... 89!> \param ends_array_mc_block ... 90!> \param weights_cos_tf_t_to_w ... 91!> \param tj ... 92!> \param tau_tj ... 93!> \param e_fermi ... 94!> \param eps_filter ... 95!> \param alpha ... 96!> \param eps_filter_im_time ... 97!> \param Eigenval ... 98!> \param nmo ... 99!> \param num_integ_points ... 100!> \param cut_memory ... 101!> \param unit_nr ... 102!> \param mp2_env ... 103!> \param para_env ... 104!> \param qs_env ... 105!> \param index_to_cell_3c ... 106!> \param cell_to_index_3c ... 107!> \param has_mat_P_blocks ... 108!> \param do_ri_sos_laplace_mp2 ... 109!> \param dbcsr_time ... 110!> \param dbcsr_nflop ... 111! ************************************************************************************************** 112 SUBROUTINE compute_mat_P_omega(mat_P_omega, fm_scaled_dm_occ_tau, & 113 fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, & 114 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, & 115 mat_P_global, & 116 matrix_s, & 117 ispin, & 118 t_3c_M, t_3c_O, & 119 starts_array_mc, ends_array_mc, & 120 starts_array_mc_block, ends_array_mc_block, & 121 weights_cos_tf_t_to_w, & 122 tj, tau_tj, e_fermi, eps_filter, & 123 alpha, eps_filter_im_time, Eigenval, nmo, & 124 num_integ_points, cut_memory, unit_nr, & 125 mp2_env, para_env, & 126 qs_env, index_to_cell_3c, cell_to_index_3c, & 127 has_mat_P_blocks, do_ri_sos_laplace_mp2, & 128 dbcsr_time, dbcsr_nflop) 129 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega 130 TYPE(cp_fm_type), POINTER :: fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, fm_mo_coeff_occ, & 131 fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled 132 TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_P_global 133 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 134 INTEGER, INTENT(IN) :: ispin 135 TYPE(dbcsr_t_type), INTENT(INOUT) :: t_3c_M 136 TYPE(dbcsr_t_type), DIMENSION(:, :), INTENT(INOUT) :: t_3c_O 137 INTEGER, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc, & 138 starts_array_mc_block, & 139 ends_array_mc_block 140 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 141 INTENT(IN) :: weights_cos_tf_t_to_w 142 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 143 INTENT(IN) :: tj 144 INTEGER, INTENT(IN) :: num_integ_points, nmo 145 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval 146 REAL(KIND=dp), INTENT(IN) :: eps_filter_im_time, alpha, eps_filter, & 147 e_fermi 148 REAL(KIND=dp), DIMENSION(0:num_integ_points), & 149 INTENT(IN) :: tau_tj 150 INTEGER, INTENT(IN) :: cut_memory, unit_nr 151 TYPE(mp2_type), POINTER :: mp2_env 152 TYPE(cp_para_env_type), POINTER :: para_env 153 TYPE(qs_environment_type), POINTER :: qs_env 154 INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell_3c 155 INTEGER, ALLOCATABLE, DIMENSION(:, :, :), & 156 INTENT(IN) :: cell_to_index_3c 157 LOGICAL, DIMENSION(:, :, :, :, :), INTENT(INOUT) :: has_mat_P_blocks 158 LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2 159 REAL(dp), INTENT(INOUT) :: dbcsr_time 160 INTEGER(int_8), INTENT(INOUT) :: dbcsr_nflop 161 162 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_mat_P_omega', & 163 routineP = moduleN//':'//routineN 164 165 INTEGER :: comm_2d, handle, handle2, handle3, i, i_cell, i_cell_R_1, i_cell_R_1_minus_S, & 166 i_cell_R_1_minus_T, i_cell_R_2, i_cell_R_2_minus_S_minus_T, i_cell_S, i_cell_T, i_mem, & 167 iquad, j, j_mem, jquad, num_3c_repl, num_cells_dm, unit_nr_2, unit_nr_dbcsr 168 INTEGER(int_8) :: nze, nze_dm_occ, nze_dm_virt, nze_M_occ, & 169 nze_M_virt, nze_O 170 INTEGER(KIND=int_8) :: flops_1_max_occ, flops_1_max_virt, & 171 flops_1_occ, flops_1_virt, flops_2, & 172 flops_2_max 173 INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_1, dist_2, size_dm, size_P 174 INTEGER, DIMENSION(2) :: pdims_2d 175 INTEGER, DIMENSION(2, 1) :: ibounds_2, jbounds_2 176 INTEGER, DIMENSION(2, 2) :: ibounds_1, jbounds_1 177 INTEGER, DIMENSION(3) :: bounds_3c 178 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm 179 LOGICAL :: do_Gamma_RPA, do_kpoints_cubic_RPA, do_opt_pgrid, first_cycle_im_time, & 180 first_cycle_omega_loop, memory_info, pgrid_1_init_occ, pgrid_1_init_virt, pgrid_2_init, & 181 R_1_minus_S_needed, R_1_minus_T_needed, R_2_minus_S_minus_T_needed 182 REAL(dp) :: occ, occ_dm_occ, occ_dm_virt, occ_M_occ, & 183 occ_M_virt, occ_O, t1_flop 184 REAL(KIND=dp) :: omega, omega_old, t1, t2, tau, weight, & 185 weight_old 186 TYPE(dbcsr_distribution_type) :: dist_P 187 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_occ_global, mat_dm_virt_global 188 TYPE(dbcsr_t_pgrid_type) :: pgrid_1_use_occ, pgrid_1_use_virt, & 189 pgrid_2_use, pgrid_2d 190 TYPE(dbcsr_t_pgrid_type), POINTER :: pgrid_1_opt_occ, pgrid_1_opt_virt, & 191 pgrid_2_opt 192 TYPE(dbcsr_t_type) :: t_3c_M_occ, t_3c_M_occ_tmp, t_3c_M_virt, & 193 t_3c_M_virt_tmp, t_dm, t_dm_tmp, t_P, & 194 t_P_tmp 195 TYPE(dbcsr_t_type), ALLOCATABLE, DIMENSION(:) :: t_dm_occ, t_dm_virt 196 TYPE(dbcsr_t_type), & 197 DIMENSION(SIZE(t_3c_O, 1), SIZE(t_3c_O, 2)) :: t_3c_O_occ, t_3c_O_virt 198 199 NULLIFY (pgrid_1_opt_occ, pgrid_1_opt_virt, pgrid_2_opt) 200 201 CALL timeset(routineN, handle) 202 203 memory_info = mp2_env%ri_rpa_im_time%memory_info 204 IF (memory_info) THEN 205 unit_nr_dbcsr = unit_nr 206 ELSE 207 unit_nr_dbcsr = 0 208 ENDIF 209 210 do_kpoints_cubic_RPA = qs_env%mp2_env%ri_rpa_im_time%do_im_time_kpoints 211 do_Gamma_RPA = .NOT. do_kpoints_cubic_RPA 212 num_3c_repl = MAXVAL(cell_to_index_3c) 213 do_opt_pgrid = qs_env%mp2_env%ri_rpa_im_time%group_size_internal 214 215 first_cycle_im_time = .TRUE. 216 DO i = 1, SIZE(t_3c_O, 1) 217 DO j = 1, SIZE(t_3c_O, 2) 218 CALL dbcsr_t_create(t_3c_O(i, j), t_3c_O_occ(i, j)) 219 CALL dbcsr_t_copy(t_3c_O(i, j), t_3c_O_occ(i, j)) 220 CALL dbcsr_t_create(t_3c_O(i, j), t_3c_O_virt(i, j)) 221 CALL dbcsr_t_copy(t_3c_O(i, j), t_3c_O_virt(i, j)) 222 !CALL dbcsr_t_clear(t_3c_O(i, j)) ! clearing t_3c_O is not safe because it may be used later 223 ENDDO 224 ENDDO 225 226 pgrid_1_init_occ = .FALSE.; pgrid_1_init_virt = .FALSE.; pgrid_2_init = .FALSE. 227 DO jquad = 1, num_integ_points 228 229 CALL mp_sync(para_env%group) 230 t1 = m_walltime() 231 232 flops_1_max_virt = 0; flops_1_max_occ = 0; flops_2_max = 0 233 234 unit_nr_2 = unit_nr_dbcsr 235 IF (pgrid_1_init_occ) THEN 236 DO i = 1, SIZE(t_3c_O, 1) 237 DO j = 1, SIZE(t_3c_O, 2) 238 CALL tensor_change_pgrid(t_3c_O_occ(i, j), pgrid_1_use_occ, & 239 starts_array_mc_block_2=starts_array_mc_block, & 240 ends_array_mc_block_2=ends_array_mc_block, & 241 unit_nr=unit_nr_2) 242 unit_nr_2 = 0 243 ENDDO 244 ENDDO 245 ENDIF 246 247 unit_nr_2 = unit_nr_dbcsr 248 IF (pgrid_1_init_virt) THEN 249 DO i = 1, SIZE(t_3c_O, 1) 250 DO j = 1, SIZE(t_3c_O, 2) 251 CALL tensor_change_pgrid(t_3c_O_virt(i, j), pgrid_1_use_virt, & 252 starts_array_mc_block_2=starts_array_mc_block, & 253 ends_array_mc_block_2=ends_array_mc_block, & 254 unit_nr=unit_nr_2) 255 unit_nr_2 = 0 256 ENDDO 257 ENDDO 258 ENDIF 259 260 IF (pgrid_2_init) THEN 261 CALL tensor_change_pgrid(t_3c_M, pgrid_2_use, & 262 starts_array_mc_block_2=starts_array_mc_block, & 263 ends_array_mc_block_2=ends_array_mc_block, & 264 starts_array_mc_block_3=starts_array_mc_block, & 265 ends_array_mc_block_3=ends_array_mc_block, & 266 nodata=.TRUE., unit_nr=unit_nr_dbcsr) 267 ENDIF 268 269 CALL compute_mat_dm_global(fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, tau_tj, num_integ_points, nmo, & 270 fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, & 271 fm_mo_coeff_virt_scaled, mat_dm_occ_global, mat_dm_virt_global, & 272 matrix_s, ispin, & 273 Eigenval, e_fermi, eps_filter, memory_info, & 274 unit_nr, & 275 jquad, do_kpoints_cubic_RPA, qs_env, & 276 num_cells_dm, index_to_cell_dm) 277 278 ALLOCATE (t_dm_virt(num_cells_dm)) 279 ALLOCATE (t_dm_occ(num_cells_dm)) 280 CALL dbcsr_get_info(mat_P_global%matrix, distribution=dist_P) 281 CALL dbcsr_distribution_get(dist_P, group=comm_2d, nprows=pdims_2d(1), npcols=pdims_2d(2)) 282 283 pgrid_2d = dbcsr_t_nd_mp_comm(comm_2d, [1], [2], pdims_2d=pdims_2d) 284 ALLOCATE (size_P(dbcsr_t_nblks_total(t_3c_M, 1))) 285 CALL dbcsr_t_get_info(t_3c_M, blk_size_1=size_P) 286 287 ALLOCATE (size_dm(dbcsr_t_nblks_total(t_3c_O(1, 1), 3))) 288 CALL dbcsr_t_get_info(t_3c_O(1, 1), blk_size_3=size_dm) 289 CALL create_2c_tensor(t_dm, dist_1, dist_2, pgrid_2d, size_dm, size_dm, name="D (AO | AO)") 290 DEALLOCATE (size_dm) 291 DEALLOCATE (dist_1, dist_2) 292 CALL create_2c_tensor(t_P, dist_1, dist_2, pgrid_2d, size_P, size_P, name="P (RI | RI)") 293 DEALLOCATE (size_P) 294 DEALLOCATE (dist_1, dist_2) 295 CALL dbcsr_t_pgrid_destroy(pgrid_2d) 296 297 DO i_cell = 1, num_cells_dm 298 CALL dbcsr_t_create(t_dm, t_dm_virt(i_cell), name="D virt (AO | AO)") 299 CALL dbcsr_t_create(mat_dm_virt_global(jquad, i_cell)%matrix, t_dm_tmp) 300 CALL dbcsr_t_copy_matrix_to_tensor(mat_dm_virt_global(jquad, i_cell)%matrix, t_dm_tmp) 301 CALL dbcsr_t_copy(t_dm_tmp, t_dm_virt(i_cell), move_data=.TRUE.) 302 CALL dbcsr_clear(mat_dm_virt_global(jquad, i_cell)%matrix) 303 304 CALL dbcsr_t_create(t_dm, t_dm_occ(i_cell), name="D occ (AO | AO)") 305 CALL dbcsr_t_copy_matrix_to_tensor(mat_dm_occ_global(jquad, i_cell)%matrix, t_dm_tmp) 306 CALL dbcsr_t_copy(t_dm_tmp, t_dm_occ(i_cell), move_data=.TRUE.) 307 CALL dbcsr_t_destroy(t_dm_tmp) 308 CALL dbcsr_clear(mat_dm_occ_global(jquad, i_cell)%matrix) 309 ENDDO 310 311 CALL get_tensor_occupancy(t_dm_occ(1), nze_dm_occ, occ_dm_occ) 312 CALL get_tensor_occupancy(t_dm_virt(1), nze_dm_virt, occ_dm_virt) 313 IF (do_Gamma_RPA) CALL get_tensor_occupancy(t_3c_O(1, 1), nze_O, occ_O) 314 315 CALL dbcsr_t_destroy(t_dm) 316 317 CALL dbcsr_t_create(t_3c_O_occ(1, 1), t_3c_M_occ_tmp, name="M (RI AO | AO)") 318 CALL dbcsr_t_create(t_3c_O_virt(1, 1), t_3c_M_virt_tmp, name="M (RI AO | AO)") 319 CALL dbcsr_t_create(t_3c_M, t_3c_M_occ, name="M occ (RI | AO AO)") 320 CALL dbcsr_t_create(t_3c_M, t_3c_M_virt, name="M virt (RI | AO AO)") 321 322 CALL timeset(routineN//"_contract", handle2) 323 324 CALL mp_sync(para_env%group) 325 t1_flop = m_walltime() 326 327 DO i_cell_T = 1, num_cells_dm/2 + 1 328 329 IF (.NOT. ANY(has_mat_P_blocks(i_cell_T, :, :, :, :))) CYCLE 330 331 CALL dbcsr_t_batched_contract_init(t_P) 332 333 IF (do_Gamma_RPA) THEN 334 nze_M_virt = 0 335 nze_M_occ = 0 336 occ_M_virt = 0.0_dp 337 occ_M_occ = 0.0_dp 338 ENDIF 339 340 DO j_mem = 1, cut_memory 341 342 CALL dbcsr_t_get_info(t_3c_O_occ(1, 1), nfull_total=bounds_3c) 343 344 jbounds_1(:, 1) = [1, bounds_3c(1)] 345 jbounds_1(:, 2) = [starts_array_mc(j_mem), ends_array_mc(j_mem)] 346 347 jbounds_2(:, 1) = [starts_array_mc(j_mem), ends_array_mc(j_mem)] 348 349 IF (do_Gamma_RPA) CALL dbcsr_t_batched_contract_init(t_dm_virt(1)) 350 351 DO i_mem = 1, cut_memory 352 353 IF (.NOT. ANY(has_mat_P_blocks(i_cell_T, i_mem, j_mem, :, :))) CYCLE 354 355 ibounds_1(:, 1) = [1, bounds_3c(1)] 356 ibounds_1(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)] 357 358 ibounds_2(:, 1) = [starts_array_mc(i_mem), ends_array_mc(i_mem)] 359 360 IF (unit_nr_dbcsr > 0) WRITE (UNIT=unit_nr_dbcsr, FMT="(T3,A,I3,1X,I3)") & 361 "RPA_LOW_SCALING_INFO| Memory Cut iteration", i_mem, j_mem 362 363 DO i_cell_R_1 = 1, num_3c_repl 364 365 DO i_cell_R_2 = 1, num_3c_repl 366 367 IF (.NOT. has_mat_P_blocks(i_cell_T, i_mem, j_mem, i_cell_R_1, i_cell_R_2)) CYCLE 368 369 CALL get_diff_index_3c(i_cell_R_1, i_cell_T, i_cell_R_1_minus_T, & 370 index_to_cell_3c, cell_to_index_3c, index_to_cell_dm, & 371 R_1_minus_T_needed, do_kpoints_cubic_RPA) 372 DO i_cell_S = 1, num_cells_dm 373 CALL get_diff_index_3c(i_cell_R_1, i_cell_S, i_cell_R_1_minus_S, index_to_cell_3c, & 374 cell_to_index_3c, index_to_cell_dm, R_1_minus_S_needed, & 375 do_kpoints_cubic_RPA) 376 IF (R_1_minus_S_needed) THEN 377 378 CALL timeset(routineN//"_calc_M_occ_t", handle3) 379 380 CALL dbcsr_t_contract(alpha=dbcsr_scalar(1.0_dp), & 381 tensor_1=t_3c_O_occ(i_cell_R_1_minus_S, i_cell_R_2), & 382 tensor_2=t_dm_occ(i_cell_S), & 383 beta=dbcsr_scalar(1.0_dp), & 384 tensor_3=t_3c_M_occ_tmp, & 385 contract_1=[3], notcontract_1=[1, 2], & 386 contract_2=[2], notcontract_2=[1], & 387 map_1=[1, 2], map_2=[3], & 388 bounds_2=jbounds_1, bounds_3=ibounds_2, & 389 pgrid_opt_1=pgrid_1_opt_occ, & 390 filter_eps=eps_filter, unit_nr=unit_nr_dbcsr, & 391 flop=flops_1_occ) 392 CALL timestop(handle3) 393 394 dbcsr_nflop = dbcsr_nflop + flops_1_occ 395 396 IF (do_opt_pgrid) THEN 397 IF (flops_1_occ .GT. flops_1_max_occ) THEN 398 CPASSERT(ASSOCIATED(pgrid_1_opt_occ)) 399 IF (pgrid_1_init_occ) CALL dbcsr_t_pgrid_destroy(pgrid_1_use_occ) 400 pgrid_1_use_occ = pgrid_1_opt_occ 401 DEALLOCATE (pgrid_1_opt_occ) 402 pgrid_1_init_occ = .TRUE. 403 flops_1_max_occ = flops_1_occ 404 ELSEIF (ASSOCIATED(pgrid_1_opt_occ)) THEN 405 CALL dbcsr_t_pgrid_destroy(pgrid_1_opt_occ) 406 DEALLOCATE (pgrid_1_opt_occ) 407 ENDIF 408 ENDIF 409 ENDIF 410 ENDDO 411 412 ! copy matrix to optimal contraction layout - copy is done manually in order 413 ! to better control memory allocations (we can release data of previous 414 ! representation) 415 CALL timeset(routineN//"_copy_M_occ_t", handle3) 416 CALL dbcsr_t_copy(t_3c_M_occ_tmp, t_3c_M_occ, order=[1, 3, 2], move_data=.TRUE.) 417 CALL dbcsr_t_filter(t_3c_M_occ, eps_filter) 418 CALL timestop(handle3) 419 420 IF (do_Gamma_RPA) THEN 421 CALL get_tensor_occupancy(t_3c_M_occ, nze, occ) 422 nze_M_occ = nze_M_occ + nze 423 occ_M_occ = occ_M_occ + occ 424 ENDIF 425 426 DO i_cell_S = 1, num_cells_dm 427 CALL get_diff_diff_index_3c(i_cell_R_2, i_cell_S, i_cell_T, i_cell_R_2_minus_S_minus_T, & 428 index_to_cell_3c, cell_to_index_3c, index_to_cell_dm, & 429 R_2_minus_S_minus_T_needed, do_kpoints_cubic_RPA) 430 431 IF (R_1_minus_T_needed .AND. R_2_minus_S_minus_T_needed) THEN 432 433 CALL timeset(routineN//"_calc_M_virt_t", handle3) 434 CALL dbcsr_t_contract(alpha=dbcsr_scalar(alpha/2.0_dp), & 435 tensor_1=t_3c_O_virt( & 436 i_cell_R_2_minus_S_minus_T, i_cell_R_1_minus_T), & 437 tensor_2=t_dm_virt(i_cell_S), & 438 beta=dbcsr_scalar(1.0_dp), & 439 tensor_3=t_3c_M_virt_tmp, & 440 contract_1=[3], notcontract_1=[1, 2], & 441 contract_2=[2], notcontract_2=[1], & 442 map_1=[1, 2], map_2=[3], & 443 bounds_2=ibounds_1, bounds_3=jbounds_2, & 444 pgrid_opt_1=pgrid_1_opt_virt, & 445 filter_eps=eps_filter, unit_nr=unit_nr_dbcsr, & 446 flop=flops_1_virt) 447 CALL timestop(handle3) 448 449 dbcsr_nflop = dbcsr_nflop + flops_1_virt 450 451 IF (do_opt_pgrid) THEN 452 IF (flops_1_virt .GT. flops_1_max_virt) THEN 453 CPASSERT(ASSOCIATED(pgrid_1_opt_virt)) 454 IF (pgrid_1_init_virt) CALL dbcsr_t_pgrid_destroy(pgrid_1_use_virt) 455 pgrid_1_use_virt = pgrid_1_opt_virt 456 DEALLOCATE (pgrid_1_opt_virt) 457 pgrid_1_init_virt = .TRUE. 458 flops_1_max_virt = flops_1_virt 459 ELSEIF (ASSOCIATED(pgrid_1_opt_virt)) THEN 460 CALL dbcsr_t_pgrid_destroy(pgrid_1_opt_virt) 461 DEALLOCATE (pgrid_1_opt_virt) 462 ENDIF 463 ENDIF 464 ENDIF 465 ENDDO 466 467 CALL timeset(routineN//"_copy_M_virt_t", handle3) 468 CALL dbcsr_t_copy(t_3c_M_virt_tmp, t_3c_M_virt, move_data=.TRUE.) 469 CALL dbcsr_t_filter(t_3c_M_virt, eps_filter) 470 CALL timestop(handle3) 471 472 IF (do_Gamma_RPA) THEN 473 CALL get_tensor_occupancy(t_3c_M_virt, nze, occ) 474 nze_M_virt = nze_M_virt + nze 475 occ_M_virt = occ_M_virt + occ 476 ENDIF 477 478 flops_2 = 0 479 480 CALL timeset(routineN//"_calc_P_t", handle3) 481 482 CALL dbcsr_t_contract(alpha=dbcsr_scalar(1.0_dp), tensor_1=t_3c_M_occ, & 483 tensor_2=t_3c_M_virt, & 484 beta=dbcsr_scalar(1.0_dp), & 485 tensor_3=t_P, & 486 contract_1=[2, 3], notcontract_1=[1], & 487 contract_2=[2, 3], notcontract_2=[1], & 488 map_1=[1], map_2=[2], & 489 pgrid_opt_2=pgrid_2_opt, & 490 filter_eps=eps_filter_im_time/REAL(cut_memory**2, KIND=dp), & 491 flop=flops_2, & 492 move_data=.TRUE., & 493 unit_nr=unit_nr_dbcsr) 494 495 CALL timestop(handle3) 496 497 dbcsr_nflop = dbcsr_nflop + flops_2 498 499 IF (do_opt_pgrid) THEN 500 IF (flops_2 .GT. flops_2_max) THEN 501 CPASSERT(ASSOCIATED(pgrid_2_opt)) 502 IF (pgrid_2_init) CALL dbcsr_t_pgrid_destroy(pgrid_2_use) 503 pgrid_2_use = pgrid_2_opt 504 DEALLOCATE (pgrid_2_opt) 505 pgrid_2_init = .TRUE. 506 flops_2_max = flops_2 507 ELSEIF (ASSOCIATED(pgrid_2_opt)) THEN 508 CALL dbcsr_t_pgrid_destroy(pgrid_2_opt) 509 DEALLOCATE (pgrid_2_opt) 510 ENDIF 511 ENDIF 512 513 first_cycle_im_time = .FALSE. 514 515 IF (jquad == 1 .AND. flops_2 == 0) & 516 has_mat_P_blocks(i_cell_T, i_mem, j_mem, i_cell_R_1, i_cell_R_2) = .FALSE. 517 518 ENDDO 519 ENDDO 520 ENDDO 521 IF (do_Gamma_RPA) CALL dbcsr_t_batched_contract_finalize(t_dm_virt(1)) 522 ENDDO 523 524 CALL dbcsr_t_batched_contract_finalize(t_P, unit_nr=unit_nr_dbcsr) 525 526 CALL dbcsr_t_create(mat_P_global%matrix, t_P_tmp) 527 CALL dbcsr_t_copy(t_P, t_P_tmp, move_data=.TRUE.) 528 CALL dbcsr_t_copy_tensor_to_matrix(t_P_tmp, mat_P_global%matrix) 529 CALL dbcsr_t_destroy(t_P_tmp) 530 531 IF (do_ri_sos_laplace_mp2) THEN 532 ! For RI-SOS-Laplace-MP2 we do not perform a cosine transform, 533 ! but we have to copy P_local to the output matrix 534 535 CALL dbcsr_add(mat_P_omega(jquad, i_cell_T)%matrix, mat_P_global%matrix, 1.0_dp, 1.0_dp) 536 ELSE 537 CALL timeset(routineN//"_Fourier_transform", handle3) 538 539 ! Fourier transform of P(it) to P(iw) 540 first_cycle_omega_loop = .TRUE. 541 542 tau = tau_tj(jquad) 543 544 DO iquad = 1, num_integ_points 545 546 omega = tj(iquad) 547 weight = weights_cos_tf_t_to_w(iquad, jquad) 548 549 IF (first_cycle_omega_loop) THEN 550 ! no multiplication with 2.0 as in Kresses paper (Kaltak, JCTC 10, 2498 (2014), Eq. 12) 551 ! because this factor is already absorbed in the weight w_j 552 CALL dbcsr_scale(mat_P_global%matrix, COS(omega*tau)*weight) 553 ELSE 554 CALL dbcsr_scale(mat_P_global%matrix, COS(omega*tau)/COS(omega_old*tau)*weight/weight_old) 555 END IF 556 557 CALL dbcsr_add(mat_P_omega(iquad, i_cell_T)%matrix, mat_P_global%matrix, 1.0_dp, 1.0_dp) 558 559 first_cycle_omega_loop = .FALSE. 560 561 omega_old = omega 562 weight_old = weight 563 564 END DO 565 566 CALL timestop(handle3) 567 ENDIF 568 569 ENDDO 570 571 CALL timestop(handle2) 572 573 CALL dbcsr_t_destroy(t_P) 574 DO i_cell = 1, num_cells_dm 575 CALL dbcsr_t_destroy(t_dm_virt(i_cell)) 576 CALL dbcsr_t_destroy(t_dm_occ(i_cell)) 577 ENDDO 578 579 CALL dbcsr_t_destroy(t_3c_M_occ_tmp) 580 CALL dbcsr_t_destroy(t_3c_M_virt_tmp) 581 CALL dbcsr_t_destroy(t_3c_M_occ) 582 CALL dbcsr_t_destroy(t_3c_M_virt) 583 DEALLOCATE (t_dm_virt) 584 DEALLOCATE (t_dm_occ) 585 586 CALL mp_sync(para_env%group) 587 t2 = m_walltime() 588 589 dbcsr_time = dbcsr_time + t2 - t1_flop 590 591 IF (unit_nr > 0) THEN 592 WRITE (unit_nr, '(/T3,A,1X,I3)') & 593 'RPA_LOW_SCALING_INFO| Info for time point', jquad 594 WRITE (unit_nr, '(T6,A,T56,F25.6)') & 595 'Time:', t2 - t1 596 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') & 597 'Occupancy of D occ:', REAL(nze_dm_occ, dp), '/', occ_dm_occ*100, '%' 598 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') & 599 'Occupancy of D virt:', REAL(nze_dm_virt, dp), '/', occ_dm_virt*100, '%' 600 IF (do_Gamma_RPA) THEN 601 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') & 602 'Occupancy of 3c ints:', REAL(nze_O, dp), '/', occ_O*100, '%' 603 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') & 604 'Occupancy of M occ:', REAL(nze_M_occ, dp), '/', occ_M_occ*100, '%' 605 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') & 606 'Occupancy of M virt:', REAL(nze_M_virt, dp), '/', occ_M_virt*100, '%' 607 ENDIF 608 WRITE (unit_nr, *) 609 ENDIF 610 611 END DO ! time points 612 613 DO i = 1, SIZE(t_3c_O, 1) 614 DO j = 1, SIZE(t_3c_O, 2) 615 CALL dbcsr_t_destroy(t_3c_O_occ(i, j)) 616 CALL dbcsr_t_destroy(t_3c_O_virt(i, j)) 617 ENDDO 618 ENDDO 619 620 IF (pgrid_1_init_virt) CALL dbcsr_t_pgrid_destroy(pgrid_1_use_virt) 621 IF (pgrid_1_init_occ) CALL dbcsr_t_pgrid_destroy(pgrid_1_use_occ) 622 IF (pgrid_2_init) CALL dbcsr_t_pgrid_destroy(pgrid_2_use) 623 624 CALL clean_up(mat_dm_occ_global, mat_dm_virt_global) 625 626 CALL timestop(handle) 627 628 END SUBROUTINE compute_mat_P_omega 629 630! ************************************************************************************************** 631!> \brief ... 632!> \param mat_P_omega ... 633! ************************************************************************************************** 634 SUBROUTINE zero_mat_P_omega(mat_P_omega) 635 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega 636 637 INTEGER :: i_kp, jquad 638 639 DO jquad = 1, SIZE(mat_P_omega, 1) 640 DO i_kp = 1, SIZE(mat_P_omega, 2) 641 642 CALL dbcsr_set(mat_P_omega(jquad, i_kp)%matrix, 0.0_dp) 643 644 END DO 645 END DO 646 647 END SUBROUTINE zero_mat_P_omega 648 649! ************************************************************************************************** 650!> \brief ... 651!> \param fm_scaled_dm_occ_tau ... 652!> \param fm_scaled_dm_virt_tau ... 653!> \param tau_tj ... 654!> \param num_integ_points ... 655!> \param nmo ... 656!> \param fm_mo_coeff_occ ... 657!> \param fm_mo_coeff_virt ... 658!> \param fm_mo_coeff_occ_scaled ... 659!> \param fm_mo_coeff_virt_scaled ... 660!> \param mat_dm_occ_global ... 661!> \param mat_dm_virt_global ... 662!> \param matrix_s ... 663!> \param ispin ... 664!> \param Eigenval ... 665!> \param e_fermi ... 666!> \param eps_filter ... 667!> \param memory_info ... 668!> \param unit_nr ... 669!> \param jquad ... 670!> \param do_kpoints_cubic_RPA ... 671!> \param qs_env ... 672!> \param num_cells_dm ... 673!> \param index_to_cell_dm ... 674! ************************************************************************************************** 675 SUBROUTINE compute_mat_dm_global(fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, tau_tj, num_integ_points, nmo, & 676 fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, & 677 fm_mo_coeff_virt_scaled, mat_dm_occ_global, mat_dm_virt_global, & 678 matrix_s, ispin, & 679 Eigenval, e_fermi, eps_filter, memory_info, & 680 unit_nr, & 681 jquad, do_kpoints_cubic_RPA, qs_env, & 682 num_cells_dm, index_to_cell_dm) 683 684 TYPE(cp_fm_type), POINTER :: fm_scaled_dm_occ_tau, & 685 fm_scaled_dm_virt_tau 686 INTEGER, INTENT(IN) :: num_integ_points 687 REAL(KIND=dp), DIMENSION(0:num_integ_points), & 688 INTENT(IN) :: tau_tj 689 INTEGER, INTENT(IN) :: nmo 690 TYPE(cp_fm_type), POINTER :: fm_mo_coeff_occ, fm_mo_coeff_virt, & 691 fm_mo_coeff_occ_scaled, & 692 fm_mo_coeff_virt_scaled 693 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_occ_global, mat_dm_virt_global 694 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 695 INTEGER, INTENT(IN) :: ispin 696 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval 697 REAL(KIND=dp), INTENT(IN) :: e_fermi, eps_filter 698 LOGICAL, INTENT(IN) :: memory_info 699 INTEGER, INTENT(IN) :: unit_nr, jquad 700 LOGICAL, INTENT(IN) :: do_kpoints_cubic_RPA 701 TYPE(qs_environment_type), POINTER :: qs_env 702 INTEGER, INTENT(OUT) :: num_cells_dm 703 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm 704 705 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_mat_dm_global', & 706 routineP = moduleN//':'//routineN 707 REAL(KIND=dp), PARAMETER :: stabilize_exp = 70.0_dp 708 709 INTEGER :: handle, i_global, iiB, iquad, jjB, & 710 ncol_local, nrow_local, size_dm_occ, & 711 size_dm_virt 712 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 713 REAL(KIND=dp) :: tau 714 715 CALL timeset(routineN, handle) 716 717 IF (memory_info .AND. unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & 718 "RPA_LOW_SCALING_INFO| Started with time point: ", jquad 719 720 tau = tau_tj(jquad) 721 722 IF (do_kpoints_cubic_RPA) THEN 723 724 CALL compute_transl_dm(mat_dm_occ_global, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, & 725 eps_filter, num_cells_dm, index_to_cell_dm, & 726 remove_occ=.FALSE., remove_virt=.TRUE., first_jquad=1) 727 728 CALL compute_transl_dm(mat_dm_virt_global, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, & 729 eps_filter, num_cells_dm, index_to_cell_dm, & 730 remove_occ=.TRUE., remove_virt=.FALSE., first_jquad=1) 731 732 ELSE 733 734 num_cells_dm = 1 735 736 ! get info of fm_mo_coeff_occ 737 CALL cp_fm_get_info(matrix=fm_mo_coeff_occ, & 738 nrow_local=nrow_local, & 739 ncol_local=ncol_local, & 740 row_indices=row_indices, & 741 col_indices=col_indices) 742 743 ! Multiply the occupied and the virtual MO coefficients with the factor exp((-e_i-e_F)*tau/2). 744 ! Then, we simply get the sum over all occ states and virt. states by a simple matrix-matrix 745 ! multiplication. 746 747 ! first, the occ 748 DO jjB = 1, nrow_local 749 DO iiB = 1, ncol_local 750 i_global = col_indices(iiB) 751 752 ! hard coded: exponential function gets NaN if argument is negative with large absolute value 753 ! use 69, since e^(-69) = 10^(-30) which should be sufficiently small that it does not matter 754 IF (ABS(tau*0.5_dp*(Eigenval(i_global) - e_fermi)) < stabilize_exp) THEN 755 fm_mo_coeff_occ_scaled%local_data(jjB, iiB) = & 756 fm_mo_coeff_occ%local_data(jjB, iiB)*EXP(tau*0.5_dp*(Eigenval(i_global) - e_fermi)) 757 ELSE 758 fm_mo_coeff_occ_scaled%local_data(jjB, iiB) = 0.0_dp 759 END IF 760 761 END DO 762 END DO 763 764 ! get info of fm_mo_coeff_virt 765 CALL cp_fm_get_info(matrix=fm_mo_coeff_virt, & 766 nrow_local=nrow_local, & 767 ncol_local=ncol_local, & 768 row_indices=row_indices, & 769 col_indices=col_indices) 770 771 ! the same for virt 772 DO jjB = 1, nrow_local 773 DO iiB = 1, ncol_local 774 i_global = col_indices(iiB) 775 776 IF (ABS(tau*0.5_dp*(Eigenval(i_global) - e_fermi)) < stabilize_exp) THEN 777 fm_mo_coeff_virt_scaled%local_data(jjB, iiB) = & 778 fm_mo_coeff_virt%local_data(jjB, iiB)*EXP(-tau*0.5_dp*(Eigenval(i_global) - e_fermi)) 779 ELSE 780 fm_mo_coeff_virt_scaled%local_data(jjB, iiB) = 0.0_dp 781 END IF 782 783 END DO 784 END DO 785 786 size_dm_occ = nmo 787 size_dm_virt = nmo 788 789 CALL cp_gemm(transa="N", transb="T", m=size_dm_occ, n=size_dm_occ, k=nmo, alpha=1.0_dp, & 790 matrix_a=fm_mo_coeff_occ_scaled, matrix_b=fm_mo_coeff_occ_scaled, beta=0.0_dp, & 791 matrix_c=fm_scaled_dm_occ_tau) 792 793 CALL cp_gemm(transa="N", transb="T", m=size_dm_virt, n=size_dm_virt, k=nmo, alpha=1.0_dp, & 794 matrix_a=fm_mo_coeff_virt_scaled, matrix_b=fm_mo_coeff_virt_scaled, beta=0.0_dp, & 795 matrix_c=fm_scaled_dm_virt_tau) 796 797 IF (jquad == 1) THEN 798 799 ! transfer fm density matrices to dbcsr matrix 800 NULLIFY (mat_dm_occ_global) 801 CALL dbcsr_allocate_matrix_set(mat_dm_occ_global, num_integ_points, 1) 802 803 DO iquad = 1, num_integ_points 804 805 ALLOCATE (mat_dm_occ_global(iquad, 1)%matrix) 806 CALL dbcsr_create(matrix=mat_dm_occ_global(iquad, 1)%matrix, & 807 template=matrix_s(1)%matrix, & 808 matrix_type=dbcsr_type_no_symmetry) 809 810 END DO 811 812 END IF 813 814 CALL copy_fm_to_dbcsr(fm_scaled_dm_occ_tau, & 815 mat_dm_occ_global(jquad, 1)%matrix, & 816 keep_sparsity=.FALSE.) 817 818 CALL dbcsr_filter(mat_dm_occ_global(jquad, 1)%matrix, eps_filter) 819 820 IF (jquad == 1) THEN 821 822 NULLIFY (mat_dm_virt_global) 823 CALL dbcsr_allocate_matrix_set(mat_dm_virt_global, num_integ_points, 1) 824 825 END IF 826 827 ALLOCATE (mat_dm_virt_global(jquad, 1)%matrix) 828 CALL dbcsr_create(matrix=mat_dm_virt_global(jquad, 1)%matrix, & 829 template=matrix_s(1)%matrix, & 830 matrix_type=dbcsr_type_no_symmetry) 831 CALL copy_fm_to_dbcsr(fm_scaled_dm_virt_tau, & 832 mat_dm_virt_global(jquad, 1)%matrix, & 833 keep_sparsity=.FALSE.) 834 835 CALL dbcsr_filter(mat_dm_virt_global(jquad, 1)%matrix, eps_filter) 836 837 ! release memory 838 IF (jquad > 1) THEN 839 CALL dbcsr_set(mat_dm_occ_global(jquad - 1, 1)%matrix, 0.0_dp) 840 CALL dbcsr_set(mat_dm_virt_global(jquad - 1, 1)%matrix, 0.0_dp) 841 CALL dbcsr_filter(mat_dm_occ_global(jquad - 1, 1)%matrix, 0.0_dp) 842 CALL dbcsr_filter(mat_dm_virt_global(jquad - 1, 1)%matrix, 0.0_dp) 843 END IF 844 845 END IF ! do kpoints 846 847 CALL timestop(handle) 848 END SUBROUTINE compute_mat_dm_global 849 850! ************************************************************************************************** 851!> \brief ... 852!> \param mat_dm_occ_global ... 853!> \param mat_dm_virt_global ... 854! ************************************************************************************************** 855 SUBROUTINE clean_up(mat_dm_occ_global, mat_dm_virt_global) 856 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_occ_global, mat_dm_virt_global 857 858 CALL dbcsr_deallocate_matrix_set(mat_dm_occ_global) 859 CALL dbcsr_deallocate_matrix_set(mat_dm_virt_global) 860 861 END SUBROUTINE clean_up 862 863! ************************************************************************************************** 864!> \brief Calculate kpoint density matrices (rho(k), owned by kpoint groups) 865!> \param kpoint kpoint environment 866!> \param tau ... 867!> \param e_fermi ... 868!> \param remove_occ ... 869!> \param remove_virt ... 870! ************************************************************************************************** 871 SUBROUTINE kpoint_density_matrices_rpa(kpoint, tau, e_fermi, remove_occ, remove_virt) 872 873 TYPE(kpoint_type), POINTER :: kpoint 874 REAL(KIND=dp), INTENT(IN) :: tau, e_fermi 875 LOGICAL, INTENT(IN) :: remove_occ, remove_virt 876 877 CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_density_matrices_rpa', & 878 routineP = moduleN//':'//routineN 879 REAL(KIND=dp), PARAMETER :: stabilize_exp = 70.0_dp 880 881 INTEGER :: handle, i_mo, ikpgr, ispin, kplocal, & 882 nao, nmo, nspin 883 INTEGER, DIMENSION(2) :: kp_range 884 REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, exp_scaling, occupation 885 TYPE(cp_fm_struct_type), POINTER :: matrix_struct 886 TYPE(cp_fm_type), POINTER :: cpmat, fwork, rpmat 887 TYPE(kpoint_env_type), POINTER :: kp 888 TYPE(mo_set_type), POINTER :: mo_set 889 890 CALL timeset(routineN, handle) 891 892 ! only imaginary wavefunctions supported in kpoint cubic scaling RPA 893 CPASSERT(kpoint%use_real_wfn .EQV. .FALSE.) 894 895 ! work matrix 896 mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)%mo_set 897 CALL get_mo_set(mo_set, nao=nao, nmo=nmo) 898 899 ! if this CPASSERT is triggered, please add all virtual MOs to SCF section, 900 ! e.g. ADDED_MOS 1000000 901 CPASSERT(nao == nmo) 902 903 ALLOCATE (exp_scaling(nmo)) 904 905 CALL cp_fm_get_info(mo_set%mo_coeff, matrix_struct=matrix_struct) 906 CALL cp_fm_create(fwork, matrix_struct) 907 908 CALL get_kpoint_info(kpoint, kp_range=kp_range) 909 kplocal = kp_range(2) - kp_range(1) + 1 910 911 DO ikpgr = 1, kplocal 912 kp => kpoint%kp_env(ikpgr)%kpoint_env 913 nspin = SIZE(kp%mos, 2) 914 DO ispin = 1, nspin 915 mo_set => kp%mos(1, ispin)%mo_set 916 CALL get_mo_set(mo_set, eigenvalues=eigenvalues) 917 rpmat => kp%wmat(1, ispin)%matrix 918 cpmat => kp%wmat(2, ispin)%matrix 919 CALL get_mo_set(mo_set, occupation_numbers=occupation) 920 CALL cp_fm_to_fm(mo_set%mo_coeff, fwork) 921 922 IF (remove_virt) THEN 923 CALL cp_fm_column_scale(fwork, occupation) 924 END IF 925 IF (remove_occ) THEN 926 CALL cp_fm_column_scale(fwork, 2.0_dp/REAL(nspin, KIND=dp) - occupation) 927 END IF 928 929 ! proper spin 930 IF (nspin == 1) THEN 931 CALL cp_fm_scale(0.5_dp, fwork) 932 END IF 933 934 DO i_mo = 1, nmo 935 IF (ABS(tau*0.5_dp*(eigenvalues(i_mo) - e_fermi)) < stabilize_exp) THEN 936 exp_scaling(i_mo) = EXP(-ABS(tau*(eigenvalues(i_mo) - e_fermi))) 937 ELSE 938 exp_scaling(i_mo) = 0.0_dp 939 END IF 940 END DO 941 942 CALL cp_fm_column_scale(fwork, exp_scaling) 943 944 ! Re(c)*Re(c) 945 CALL cp_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, rpmat) 946 mo_set => kp%mos(2, ispin)%mo_set 947 ! Im(c)*Re(c) 948 CALL cp_gemm("N", "T", nao, nao, nmo, -1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, cpmat) 949 ! Re(c)*Im(c) 950 CALL cp_gemm("N", "T", nao, nao, nmo, 1.0_dp, fwork, mo_set%mo_coeff, 1.0_dp, cpmat) 951 952 CALL cp_fm_to_fm(mo_set%mo_coeff, fwork) 953 954 IF (remove_virt) THEN 955 CALL cp_fm_column_scale(fwork, occupation) 956 END IF 957 IF (remove_occ) THEN 958 CALL cp_fm_column_scale(fwork, 2.0_dp/REAL(nspin, KIND=dp) - occupation) 959 END IF 960 961 ! proper spin 962 IF (nspin == 1) THEN 963 CALL cp_fm_scale(0.5_dp, fwork) 964 END IF 965 966 DO i_mo = 1, nmo 967 IF (ABS(tau*0.5_dp*(eigenvalues(i_mo) - e_fermi)) < stabilize_exp) THEN 968 exp_scaling(i_mo) = EXP(-ABS(tau*(eigenvalues(i_mo) - e_fermi))) 969 ELSE 970 exp_scaling(i_mo) = 0.0_dp 971 END IF 972 END DO 973 974 CALL cp_fm_column_scale(fwork, exp_scaling) 975 ! Im(c)*Im(c) 976 CALL cp_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 1.0_dp, rpmat) 977 978 END DO 979 980 END DO 981 982 CALL cp_fm_release(fwork) 983 DEALLOCATE (exp_scaling) 984 985 CALL timestop(handle) 986 987 END SUBROUTINE kpoint_density_matrices_rpa 988 989! ************************************************************************************************** 990!> \brief ... 991!> \param mat_dm_global ... 992!> \param qs_env ... 993!> \param ispin ... 994!> \param num_integ_points ... 995!> \param jquad ... 996!> \param e_fermi ... 997!> \param tau ... 998!> \param eps_filter ... 999!> \param num_cells_dm ... 1000!> \param index_to_cell_dm ... 1001!> \param remove_occ ... 1002!> \param remove_virt ... 1003!> \param first_jquad ... 1004! ************************************************************************************************** 1005 SUBROUTINE compute_transl_dm(mat_dm_global, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, & 1006 eps_filter, num_cells_dm, index_to_cell_dm, remove_occ, remove_virt, & 1007 first_jquad) 1008 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_global 1009 TYPE(qs_environment_type), POINTER :: qs_env 1010 INTEGER, INTENT(IN) :: ispin, num_integ_points, jquad 1011 REAL(KIND=dp), INTENT(IN) :: e_fermi, tau, eps_filter 1012 INTEGER, INTENT(OUT) :: num_cells_dm 1013 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm 1014 LOGICAL, INTENT(IN) :: remove_occ, remove_virt 1015 INTEGER, INTENT(IN) :: first_jquad 1016 1017 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_transl_dm', & 1018 routineP = moduleN//':'//routineN 1019 1020 INTEGER :: handle, i_dim, i_img, iquad, jspin, nspin 1021 INTEGER, DIMENSION(3) :: cell_grid_dm 1022 TYPE(cell_type), POINTER :: cell 1023 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_global_work, matrix_s_kp 1024 TYPE(dft_control_type), POINTER :: dft_control 1025 TYPE(kpoint_type), POINTER :: kpoints 1026 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 1027 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 1028 POINTER :: sab_nl 1029 1030 CALL timeset(routineN, handle) 1031 1032 CALL get_qs_env(qs_env, & 1033 matrix_s_kp=matrix_s_kp, & 1034 sab_orb=sab_nl, & 1035 mos=mos, & 1036 dft_control=dft_control, & 1037 cell=cell, & 1038 kpoints=kpoints) 1039 1040 nspin = SIZE(mos) 1041 1042 ! we always use an odd number of image cells 1043 ! CAUTION: also at another point, cell_grid_dm is defined, these definitions have to be identical 1044 DO i_dim = 1, 3 1045 cell_grid_dm(i_dim) = (kpoints%nkp_grid(i_dim)/2)*2 - 1 1046 END DO 1047 1048 num_cells_dm = cell_grid_dm(1)*cell_grid_dm(2)*cell_grid_dm(3) 1049 1050 NULLIFY (mat_dm_global_work) 1051 CALL dbcsr_allocate_matrix_set(mat_dm_global_work, nspin, num_cells_dm) 1052 1053 DO jspin = 1, nspin 1054 1055 DO i_img = 1, num_cells_dm 1056 1057 ALLOCATE (mat_dm_global_work(jspin, i_img)%matrix) 1058 CALL dbcsr_create(matrix=mat_dm_global_work(jspin, i_img)%matrix, & 1059 template=matrix_s_kp(1, 1)%matrix, & 1060 ! matrix_type=dbcsr_type_symmetric) 1061 matrix_type=dbcsr_type_no_symmetry) 1062 1063 CALL dbcsr_reserve_all_blocks(mat_dm_global_work(jspin, i_img)%matrix) 1064 1065 CALL dbcsr_set(mat_dm_global_work(ispin, i_img)%matrix, 0.0_dp) 1066 1067 END DO 1068 1069 END DO 1070 1071 ! density matrices in k-space weighted with EXP(-|e_i-e_F|*t) for occupied orbitals 1072 CALL kpoint_density_matrices_rpa(kpoints, tau, e_fermi, & 1073 remove_occ=remove_occ, remove_virt=remove_virt) 1074 1075 ! overwrite the cell indices in kpoints 1076 CALL init_cell_index_rpa(cell_grid_dm, kpoints%cell_to_index, kpoints%index_to_cell, cell) 1077 1078 ! density matrices in real space, the cell vectors T for transforming are taken from kpoints%index_to_cell 1079 ! (custom made for RPA) and not from sab_nl (which is symmetric and from SCF) 1080 CALL density_matrix_from_kp_to_transl(kpoints, mat_dm_global_work, kpoints%index_to_cell) 1081 1082 ! we need the index to cell for the density matrices later 1083 index_to_cell_dm => kpoints%index_to_cell 1084 1085 ! normally, jquad = 1 to allocate the matrix set, but for GW jquad = 0 is the exchange self-energy 1086 IF (jquad == first_jquad) THEN 1087 1088 NULLIFY (mat_dm_global) 1089! CALL dbcsr_allocate_matrix_set(mat_dm_global, jquad:num_integ_points, num_cells_dm) 1090 ALLOCATE (mat_dm_global(first_jquad:num_integ_points, num_cells_dm)) 1091 1092 DO iquad = first_jquad, num_integ_points 1093 DO i_img = 1, num_cells_dm 1094 NULLIFY (mat_dm_global(iquad, i_img)%matrix) 1095 ALLOCATE (mat_dm_global(iquad, i_img)%matrix) 1096 CALL dbcsr_create(matrix=mat_dm_global(iquad, i_img)%matrix, & 1097 template=matrix_s_kp(1, 1)%matrix, & 1098 matrix_type=dbcsr_type_no_symmetry) 1099 1100 END DO 1101 END DO 1102 1103 END IF 1104 1105 DO i_img = 1, num_cells_dm 1106 1107 ! filter to get rid of the blocks full with zeros on the lower half, otherwise blocks doubled 1108 CALL dbcsr_filter(mat_dm_global_work(ispin, i_img)%matrix, eps_filter) 1109 1110 CALL dbcsr_copy(mat_dm_global(jquad, i_img)%matrix, & 1111 mat_dm_global_work(ispin, i_img)%matrix) 1112 1113 END DO 1114 1115 CALL dbcsr_deallocate_matrix_set(mat_dm_global_work) 1116 1117 CALL timestop(handle) 1118 1119 END SUBROUTINE compute_transl_dm 1120 1121! ************************************************************************************************** 1122!> \brief ... 1123!> \param kpoints ... 1124!> \param mat_dm_global_work ... 1125!> \param index_to_cell ... 1126! ************************************************************************************************** 1127 SUBROUTINE density_matrix_from_kp_to_transl(kpoints, mat_dm_global_work, index_to_cell) 1128 1129 TYPE(kpoint_type), POINTER :: kpoints 1130 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_global_work 1131 INTEGER, DIMENSION(:, :), OPTIONAL, POINTER :: index_to_cell 1132 1133 CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_from_kp_to_transl', & 1134 routineP = moduleN//':'//routineN 1135 1136 INTEGER :: handle, icell, ik, ispin, nkp, nspin, & 1137 xcell, ycell, zcell 1138 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 1139 REAL(KIND=dp) :: arg, coskl, sinkl 1140 REAL(KIND=dp), DIMENSION(:), POINTER :: wkp 1141 REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp 1142 TYPE(cp_fm_type), POINTER :: cpmat, rpmat 1143 TYPE(dbcsr_type), POINTER :: mat_work_im, mat_work_re 1144 TYPE(kpoint_env_type), POINTER :: kp 1145 1146 CALL timeset(routineN, handle) 1147 1148 NULLIFY (cell_to_index, xkp, wkp) 1149 1150 NULLIFY (mat_work_re) 1151 CALL dbcsr_init_p(mat_work_re) 1152 CALL dbcsr_create(matrix=mat_work_re, & 1153 template=mat_dm_global_work(1, 1)%matrix, & 1154 matrix_type=dbcsr_type_no_symmetry) 1155 1156 NULLIFY (mat_work_im) 1157 CALL dbcsr_init_p(mat_work_im) 1158 CALL dbcsr_create(matrix=mat_work_im, & 1159 template=mat_dm_global_work(1, 1)%matrix, & 1160 matrix_type=dbcsr_type_no_symmetry) 1161 1162 CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, wkp=wkp, & 1163 cell_to_index=cell_to_index) 1164 1165 nspin = SIZE(mat_dm_global_work, 1) 1166 1167 CPASSERT(SIZE(mat_dm_global_work, 2) == SIZE(index_to_cell, 2)) 1168 1169 DO ispin = 1, nspin 1170 1171 DO icell = 1, SIZE(mat_dm_global_work, 2) 1172 1173 CALL dbcsr_set(mat_dm_global_work(ispin, icell)%matrix, 0.0_dp) 1174 1175 END DO 1176 1177 END DO 1178 1179 DO ispin = 1, nspin 1180 1181 DO ik = 1, nkp 1182 1183 kp => kpoints%kp_env(ik)%kpoint_env 1184 rpmat => kp%wmat(1, ispin)%matrix 1185 cpmat => kp%wmat(2, ispin)%matrix 1186 1187 CALL copy_fm_to_dbcsr(rpmat, mat_work_re, keep_sparsity=.FALSE.) 1188 CALL copy_fm_to_dbcsr(cpmat, mat_work_im, keep_sparsity=.FALSE.) 1189 1190 DO icell = 1, SIZE(mat_dm_global_work, 2) 1191 1192 xcell = index_to_cell(1, icell) 1193 ycell = index_to_cell(2, icell) 1194 zcell = index_to_cell(3, icell) 1195 1196 arg = REAL(xcell, dp)*xkp(1, ik) + REAL(ycell, dp)*xkp(2, ik) + REAL(zcell, dp)*xkp(3, ik) 1197 coskl = wkp(ik)*COS(twopi*arg) 1198 sinkl = wkp(ik)*SIN(twopi*arg) 1199 1200 CALL dbcsr_add(mat_dm_global_work(ispin, icell)%matrix, mat_work_re, 1.0_dp, coskl) 1201 CALL dbcsr_add(mat_dm_global_work(ispin, icell)%matrix, mat_work_im, 1.0_dp, sinkl) 1202 1203 END DO 1204 1205 END DO 1206 END DO 1207 1208 CALL dbcsr_release_p(mat_work_re) 1209 CALL dbcsr_release_p(mat_work_im) 1210 1211 CALL timestop(handle) 1212 1213 END SUBROUTINE density_matrix_from_kp_to_transl 1214 1215! ************************************************************************************************** 1216!> \brief ... 1217!> \param cell_grid ... 1218!> \param cell_to_index ... 1219!> \param index_to_cell ... 1220!> \param cell ... 1221! ************************************************************************************************** 1222 SUBROUTINE init_cell_index_rpa(cell_grid, cell_to_index, index_to_cell, cell) 1223 INTEGER, DIMENSION(3), INTENT(IN) :: cell_grid 1224 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 1225 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell 1226 TYPE(cell_type), POINTER :: cell 1227 1228 CHARACTER(LEN=*), PARAMETER :: routineN = 'init_cell_index_rpa', & 1229 routineP = moduleN//':'//routineN 1230 1231 INTEGER :: cell_counter, handle, i_cell, & 1232 index_min_dist, num_cells, xcell, & 1233 ycell, zcell 1234 INTEGER, DIMENSION(3) :: itm 1235 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_unsorted 1236 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_unsorted 1237 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: abs_cell_vectors 1238 REAL(KIND=dp), DIMENSION(3) :: cell_vector 1239 REAL(KIND=dp), DIMENSION(3, 3) :: hmat 1240 1241 CALL timeset(routineN, handle) 1242 1243 CALL get_cell(cell=cell, h=hmat) 1244 1245 num_cells = cell_grid(1)*cell_grid(2)*cell_grid(3) 1246 itm(:) = cell_grid(:)/2 1247 1248 ! check that real space super lattice is a (2n+1)x(2m+1)x(2k+1) super lattice with the unit cell 1249 ! in the middle 1250 CPASSERT(cell_grid(1) .NE. itm(1)*2) 1251 CPASSERT(cell_grid(2) .NE. itm(2)*2) 1252 CPASSERT(cell_grid(3) .NE. itm(3)*2) 1253 1254 IF (ASSOCIATED(cell_to_index)) DEALLOCATE (cell_to_index) 1255 IF (ASSOCIATED(index_to_cell)) DEALLOCATE (index_to_cell) 1256 1257 ALLOCATE (cell_to_index_unsorted(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3))) 1258 cell_to_index_unsorted(:, :, :) = 0 1259 1260 ALLOCATE (index_to_cell_unsorted(3, num_cells)) 1261 index_to_cell_unsorted(:, :) = 0 1262 1263 ALLOCATE (cell_to_index(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3))) 1264 cell_to_index(:, :, :) = 0 1265 1266 ALLOCATE (index_to_cell(3, num_cells)) 1267 index_to_cell(:, :) = 0 1268 1269 ALLOCATE (abs_cell_vectors(1:num_cells)) 1270 1271 cell_counter = 0 1272 1273 DO xcell = -itm(1), itm(1) 1274 DO ycell = -itm(2), itm(2) 1275 DO zcell = -itm(3), itm(3) 1276 1277 cell_counter = cell_counter + 1 1278 cell_to_index_unsorted(xcell, ycell, zcell) = cell_counter 1279 1280 index_to_cell_unsorted(1, cell_counter) = xcell 1281 index_to_cell_unsorted(2, cell_counter) = ycell 1282 index_to_cell_unsorted(3, cell_counter) = zcell 1283 1284 cell_vector(1:3) = MATMUL(hmat, REAL(index_to_cell_unsorted(1:3, cell_counter), dp)) 1285 1286 abs_cell_vectors(cell_counter) = SQRT(cell_vector(1)**2 + cell_vector(2)**2 + cell_vector(3)**2) 1287 1288 END DO 1289 END DO 1290 END DO 1291 1292 ! first only do all symmetry non-equivalent cells, we need that because chi^T is computed for 1293 ! cell indices T from index_to_cell(:,1:num_cells/2+1) 1294 DO i_cell = 1, num_cells/2 + 1 1295 1296 index_min_dist = MINLOC(abs_cell_vectors(1:num_cells/2 + 1), DIM=1) 1297 1298 xcell = index_to_cell_unsorted(1, index_min_dist) 1299 ycell = index_to_cell_unsorted(2, index_min_dist) 1300 zcell = index_to_cell_unsorted(3, index_min_dist) 1301 1302 index_to_cell(1, i_cell) = xcell 1303 index_to_cell(2, i_cell) = ycell 1304 index_to_cell(3, i_cell) = zcell 1305 1306 cell_to_index(xcell, ycell, zcell) = i_cell 1307 1308 abs_cell_vectors(index_min_dist) = 10000000000.0_dp 1309 1310 END DO 1311 1312 ! now all the remaining cells 1313 DO i_cell = num_cells/2 + 2, num_cells 1314 1315 index_min_dist = MINLOC(abs_cell_vectors(1:num_cells), DIM=1) 1316 1317 xcell = index_to_cell_unsorted(1, index_min_dist) 1318 ycell = index_to_cell_unsorted(2, index_min_dist) 1319 zcell = index_to_cell_unsorted(3, index_min_dist) 1320 1321 index_to_cell(1, i_cell) = xcell 1322 index_to_cell(2, i_cell) = ycell 1323 index_to_cell(3, i_cell) = zcell 1324 1325 cell_to_index(xcell, ycell, zcell) = i_cell 1326 1327 abs_cell_vectors(index_min_dist) = 10000000000.0_dp 1328 1329 END DO 1330 1331 DEALLOCATE (index_to_cell_unsorted, cell_to_index_unsorted, abs_cell_vectors) 1332 1333 CALL timestop(handle) 1334 1335 END SUBROUTINE init_cell_index_rpa 1336 1337! ************************************************************************************************** 1338!> \brief ... 1339!> \param i_cell_R ... 1340!> \param i_cell_S ... 1341!> \param i_cell_R_minus_S ... 1342!> \param index_to_cell_3c ... 1343!> \param cell_to_index_3c ... 1344!> \param index_to_cell_dm ... 1345!> \param R_minus_S_needed ... 1346!> \param do_kpoints_cubic_RPA ... 1347! ************************************************************************************************** 1348 SUBROUTINE get_diff_index_3c(i_cell_R, i_cell_S, i_cell_R_minus_S, index_to_cell_3c, & 1349 cell_to_index_3c, index_to_cell_dm, R_minus_S_needed, & 1350 do_kpoints_cubic_RPA) 1351 1352 INTEGER, INTENT(IN) :: i_cell_R, i_cell_S 1353 INTEGER, INTENT(OUT) :: i_cell_R_minus_S 1354 INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell_3c 1355 INTEGER, ALLOCATABLE, DIMENSION(:, :, :), & 1356 INTENT(IN) :: cell_to_index_3c 1357 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm 1358 LOGICAL, INTENT(OUT) :: R_minus_S_needed 1359 LOGICAL, INTENT(IN) :: do_kpoints_cubic_RPA 1360 1361 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_diff_index_3c', & 1362 routineP = moduleN//':'//routineN 1363 1364 INTEGER :: handle, x_cell_R, x_cell_R_minus_S, x_cell_S, y_cell_R, y_cell_R_minus_S, & 1365 y_cell_S, z_cell_R, z_cell_R_minus_S, z_cell_S 1366 1367 CALL timeset(routineN, handle) 1368 1369 IF (do_kpoints_cubic_RPA) THEN 1370 1371 x_cell_R = index_to_cell_3c(1, i_cell_R) 1372 y_cell_R = index_to_cell_3c(2, i_cell_R) 1373 z_cell_R = index_to_cell_3c(3, i_cell_R) 1374 1375 x_cell_S = index_to_cell_dm(1, i_cell_S) 1376 y_cell_S = index_to_cell_dm(2, i_cell_S) 1377 z_cell_S = index_to_cell_dm(3, i_cell_S) 1378 1379 x_cell_R_minus_S = x_cell_R - x_cell_S 1380 y_cell_R_minus_S = y_cell_R - y_cell_S 1381 z_cell_R_minus_S = z_cell_R - z_cell_S 1382 1383 IF (x_cell_R_minus_S .GE. LBOUND(cell_to_index_3c, 1) .AND. & 1384 x_cell_R_minus_S .LE. UBOUND(cell_to_index_3c, 1) .AND. & 1385 y_cell_R_minus_S .GE. LBOUND(cell_to_index_3c, 2) .AND. & 1386 y_cell_R_minus_S .LE. UBOUND(cell_to_index_3c, 2) .AND. & 1387 z_cell_R_minus_S .GE. LBOUND(cell_to_index_3c, 3) .AND. & 1388 z_cell_R_minus_S .LE. UBOUND(cell_to_index_3c, 3)) THEN 1389 1390 i_cell_R_minus_S = cell_to_index_3c(x_cell_R_minus_S, y_cell_R_minus_S, z_cell_R_minus_S) 1391 1392 ! 0 means that there is no 3c index with this R-S vector because R-S is too big and the 3c integral is 0 1393 IF (i_cell_R_minus_S == 0) THEN 1394 1395 R_minus_S_needed = .FALSE. 1396 i_cell_R_minus_S = 0 1397 1398 ELSE 1399 1400 R_minus_S_needed = .TRUE. 1401 1402 END IF 1403 1404 ELSE 1405 1406 i_cell_R_minus_S = 0 1407 R_minus_S_needed = .FALSE. 1408 1409 END IF 1410 1411 ELSE ! no k-points 1412 1413 R_minus_S_needed = .TRUE. 1414 i_cell_R_minus_S = 1 1415 1416 END IF 1417 1418 CALL timestop(handle) 1419 1420 END SUBROUTINE get_diff_index_3c 1421 1422! ************************************************************************************************** 1423!> \brief ... 1424!> \param i_cell_R ... 1425!> \param i_cell_S ... 1426!> \param i_cell_T ... 1427!> \param i_cell_R_minus_S_minus_T ... 1428!> \param index_to_cell_3c ... 1429!> \param cell_to_index_3c ... 1430!> \param index_to_cell_dm ... 1431!> \param R_minus_S_minus_T_needed ... 1432!> \param do_kpoints_cubic_RPA ... 1433! ************************************************************************************************** 1434 SUBROUTINE get_diff_diff_index_3c(i_cell_R, i_cell_S, i_cell_T, i_cell_R_minus_S_minus_T, & 1435 index_to_cell_3c, cell_to_index_3c, index_to_cell_dm, & 1436 R_minus_S_minus_T_needed, & 1437 do_kpoints_cubic_RPA) 1438 1439 INTEGER, INTENT(IN) :: i_cell_R, i_cell_S, i_cell_T 1440 INTEGER, INTENT(OUT) :: i_cell_R_minus_S_minus_T 1441 INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell_3c 1442 INTEGER, ALLOCATABLE, DIMENSION(:, :, :), & 1443 INTENT(IN) :: cell_to_index_3c 1444 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm 1445 LOGICAL, INTENT(OUT) :: R_minus_S_minus_T_needed 1446 LOGICAL, INTENT(IN) :: do_kpoints_cubic_RPA 1447 1448 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_diff_diff_index_3c', & 1449 routineP = moduleN//':'//routineN 1450 1451 INTEGER :: handle, x_cell_R, x_cell_R_minus_S_minus_T, x_cell_S, x_cell_T, y_cell_R, & 1452 y_cell_R_minus_S_minus_T, y_cell_S, y_cell_T, z_cell_R, z_cell_R_minus_S_minus_T, & 1453 z_cell_S, z_cell_T 1454 1455 CALL timeset(routineN, handle) 1456 1457 IF (do_kpoints_cubic_RPA) THEN 1458 1459 x_cell_R = index_to_cell_3c(1, i_cell_R) 1460 y_cell_R = index_to_cell_3c(2, i_cell_R) 1461 z_cell_R = index_to_cell_3c(3, i_cell_R) 1462 1463 x_cell_S = index_to_cell_dm(1, i_cell_S) 1464 y_cell_S = index_to_cell_dm(2, i_cell_S) 1465 z_cell_S = index_to_cell_dm(3, i_cell_S) 1466 1467 x_cell_T = index_to_cell_dm(1, i_cell_T) 1468 y_cell_T = index_to_cell_dm(2, i_cell_T) 1469 z_cell_T = index_to_cell_dm(3, i_cell_T) 1470 1471 x_cell_R_minus_S_minus_T = x_cell_R - x_cell_S - x_cell_T 1472 y_cell_R_minus_S_minus_T = y_cell_R - y_cell_S - y_cell_T 1473 z_cell_R_minus_S_minus_T = z_cell_R - z_cell_S - z_cell_T 1474 1475 IF (x_cell_R_minus_S_minus_T .GE. LBOUND(cell_to_index_3c, 1) .AND. & 1476 x_cell_R_minus_S_minus_T .LE. UBOUND(cell_to_index_3c, 1) .AND. & 1477 y_cell_R_minus_S_minus_T .GE. LBOUND(cell_to_index_3c, 2) .AND. & 1478 y_cell_R_minus_S_minus_T .LE. UBOUND(cell_to_index_3c, 2) .AND. & 1479 z_cell_R_minus_S_minus_T .GE. LBOUND(cell_to_index_3c, 3) .AND. & 1480 z_cell_R_minus_S_minus_T .LE. UBOUND(cell_to_index_3c, 3)) THEN 1481 1482 i_cell_R_minus_S_minus_T = cell_to_index_3c(x_cell_R_minus_S_minus_T, & 1483 y_cell_R_minus_S_minus_T, & 1484 z_cell_R_minus_S_minus_T) 1485 1486 ! index 0 means that there are only no 3c matrix elements because R-S-T is too big 1487 IF (i_cell_R_minus_S_minus_T == 0) THEN 1488 1489 R_minus_S_minus_T_needed = .FALSE. 1490 1491 ELSE 1492 1493 R_minus_S_minus_T_needed = .TRUE. 1494 1495 END IF 1496 1497 ELSE 1498 1499 i_cell_R_minus_S_minus_T = 0 1500 R_minus_S_minus_T_needed = .FALSE. 1501 1502 END IF 1503 1504 ! no k-kpoints 1505 ELSE 1506 1507 R_minus_S_minus_T_needed = .TRUE. 1508 i_cell_R_minus_S_minus_T = 1 1509 1510 END IF 1511 1512 CALL timestop(handle) 1513 1514 END SUBROUTINE get_diff_diff_index_3c 1515 1516END MODULE rpa_im_time 1517