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