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