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