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