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