(git:ed6f26b)
Loading...
Searching...
No Matches
rpa_im_time.F
Go to the documentation of this file.
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! **************************************************************************************************
14 USE cell_types, ONLY: cell_type,&
16 USE cp_dbcsr_api, ONLY: &
19 dbcsr_release_p, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
27 USE cp_fm_types, ONLY: cp_fm_create,&
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,&
39 USE kinds, ONLY: dp,&
40 int_8
41 USE kpoint_types, ONLY: get_kpoint_info,&
44 USE machine, ONLY: m_flush,&
46 USE mathconstants, ONLY: twopi
47 USE message_passing, ONLY: mp_comm_type,&
49 USE mp2_types, ONLY: mp2_type
54 USE qs_mo_types, ONLY: get_mo_set,&
56 USE qs_tensors, ONLY: decompress_tensor,&
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, &
75
76CONTAINS
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 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 t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
128 starts_array_mc, ends_array_mc, &
129 starts_array_mc_block, ends_array_mc_block, &
130 weights_cos_tf_t_to_w, &
131 tj, tau_tj, e_fermi, eps_filter, &
132 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_1, dist_2, mc_ranges, size_dm, &
185 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 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 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_occ_global, mat_dm_virt_global
199 TYPE(dbt_pgrid_type) :: pgrid_2d
200 TYPE(dbt_type) :: t_3c_m_occ, t_3c_m_occ_tmp, t_3c_m_virt, &
201 t_3c_m_virt_tmp, t_dm, t_dm_tmp, t_p, &
202 t_p_tmp
203 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_dm_occ, t_dm_virt
204 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_o_occ, t_3c_o_virt
205 TYPE(mp_comm_type) :: comm_2d
206
207 CALL timeset(routinen, handle)
208
209 memory_info = mp2_env%ri_rpa_im_time%memory_info
210 IF (memory_info) THEN
211 unit_nr_dbcsr = unit_nr
212 ELSE
213 unit_nr_dbcsr = 0
214 END IF
215
216 do_kpoints_cubic_rpa = qs_env%mp2_env%ri_rpa_im_time%do_im_time_kpoints
217 do_gamma_rpa = .NOT. do_kpoints_cubic_rpa
218 num_3c_repl = maxval(cell_to_index_3c)
219
220 first_cycle_im_time = .true.
221 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 DO i = 1, SIZE(t_3c_o, 1)
223 DO j = 1, SIZE(t_3c_o, 2)
224 CALL dbt_create(t_3c_o(i, j), t_3c_o_occ(i, j))
225 CALL dbt_create(t_3c_o(i, j), t_3c_o_virt(i, j))
226 END DO
227 END DO
228
229 CALL dbt_create(t_3c_m, t_3c_m_occ, name="M occ (RI | AO AO)")
230 CALL dbt_create(t_3c_m, t_3c_m_virt, name="M virt (RI | AO AO)")
231
232 ALLOCATE (mc_ranges(cut_memory + 1))
233 mc_ranges(:cut_memory) = starts_array_mc_block(:)
234 mc_ranges(cut_memory + 1) = ends_array_mc_block(cut_memory) + 1
235
236 DO jquad = 1, num_integ_points
237
238 CALL para_env%sync()
239 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 num_cells_dm, index_to_cell_dm, para_env)
248
249 ALLOCATE (t_dm_virt(num_cells_dm))
250 ALLOCATE (t_dm_occ(num_cells_dm))
251 CALL dbcsr_get_info(mat_p_global%matrix, distribution=dist_p)
252 CALL dbcsr_distribution_get(dist_p, group=comm_2d_handle, nprows=pdims_2d(1), npcols=pdims_2d(2))
253 CALL comm_2d%set_handle(comm_2d_handle)
254
255 pgrid_2d = dbt_nd_mp_comm(comm_2d, [1], [2], pdims_2d=pdims_2d)
256 ALLOCATE (size_p(dbt_nblks_total(t_3c_m, 1)))
257 CALL dbt_get_info(t_3c_m, blk_size_1=size_p)
258
259 ALLOCATE (size_dm(dbt_nblks_total(t_3c_o(1, 1), 3)))
260 CALL dbt_get_info(t_3c_o(1, 1), blk_size_3=size_dm)
261 CALL create_2c_tensor(t_dm, dist_1, dist_2, pgrid_2d, size_dm, size_dm, name="D (AO | AO)")
262 DEALLOCATE (size_dm)
263 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 DEALLOCATE (size_p)
266 DEALLOCATE (dist_1, dist_2)
267 CALL dbt_pgrid_destroy(pgrid_2d)
268
269 DO i_cell = 1, num_cells_dm
270 CALL dbt_create(t_dm, t_dm_virt(i_cell), name="D virt (AO | AO)")
271 CALL dbt_create(mat_dm_virt_global(jquad, i_cell)%matrix, t_dm_tmp)
272 CALL dbt_copy_matrix_to_tensor(mat_dm_virt_global(jquad, i_cell)%matrix, t_dm_tmp)
273 CALL dbt_copy(t_dm_tmp, t_dm_virt(i_cell), move_data=.true.)
274 CALL dbcsr_clear(mat_dm_virt_global(jquad, i_cell)%matrix)
275
276 CALL dbt_create(t_dm, t_dm_occ(i_cell), name="D occ (AO | AO)")
277 CALL dbt_copy_matrix_to_tensor(mat_dm_occ_global(jquad, i_cell)%matrix, t_dm_tmp)
278 CALL dbt_copy(t_dm_tmp, t_dm_occ(i_cell), move_data=.true.)
279 CALL dbt_destroy(t_dm_tmp)
280 CALL dbcsr_clear(mat_dm_occ_global(jquad, i_cell)%matrix)
281 END DO
282
283 CALL get_tensor_occupancy(t_dm_occ(1), nze_dm_occ, occ_dm_occ)
284 CALL get_tensor_occupancy(t_dm_virt(1), nze_dm_virt, occ_dm_virt)
285
286 CALL dbt_destroy(t_dm)
287
288 CALL dbt_create(t_3c_o_occ(1, 1), t_3c_m_occ_tmp, name="M (RI AO | AO)")
289 CALL dbt_create(t_3c_o_virt(1, 1), t_3c_m_virt_tmp, name="M (RI AO | AO)")
290
291 CALL timeset(routinen//"_contract", handle2)
292
293 CALL para_env%sync()
294 t1_flop = m_walltime()
295
296 DO i = 1, SIZE(t_3c_o_occ, 1)
297 DO j = 1, SIZE(t_3c_o_occ, 2)
298 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 DO i = 1, SIZE(t_3c_o_virt, 1)
302 DO j = 1, SIZE(t_3c_o_virt, 2)
303 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 CALL dbt_batched_contract_init(t_3c_m_occ_tmp, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
307 CALL dbt_batched_contract_init(t_3c_m_virt_tmp, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
308 CALL dbt_batched_contract_init(t_3c_m_occ, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
309 CALL dbt_batched_contract_init(t_3c_m_virt, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
310
311 DO i_cell_t = 1, num_cells_dm/2 + 1
312
313 IF (.NOT. any(has_mat_p_blocks(i_cell_t, :, :, :, :))) cycle
314
315 CALL dbt_batched_contract_init(t_p)
316
317 IF (do_gamma_rpa) THEN
318 nze_o = 0
319 nze_m_virt = 0
320 nze_m_occ = 0
321 occ_m_virt = 0.0_dp
322 occ_m_occ = 0.0_dp
323 occ_o = 0.0_dp
324 END IF
325
326 DO j_mem = 1, cut_memory
327
328 CALL dbt_get_info(t_3c_o_occ(1, 1), nfull_total=bounds_3c)
329
330 jbounds_1(:, 1) = [1, bounds_3c(1)]
331 jbounds_1(:, 2) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
332
333 jbounds_2(:, 1) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
334
335 IF (do_gamma_rpa) CALL dbt_batched_contract_init(t_dm_virt(1))
336
337 DO i_mem = 1, cut_memory
338
339 IF (.NOT. any(has_mat_p_blocks(i_cell_t, i_mem, j_mem, :, :))) cycle
340
341 ibounds_1(:, 1) = [1, bounds_3c(1)]
342 ibounds_1(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
343
344 ibounds_2(:, 1) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
345
346 IF (unit_nr_dbcsr > 0) WRITE (unit=unit_nr_dbcsr, fmt="(T3,A,I3,1X,I3)") &
347 "RPA_LOW_SCALING_INFO| Memory Cut iteration", i_mem, j_mem
348
349 DO i_cell_r_1 = 1, num_3c_repl
350
351 DO i_cell_r_2 = 1, num_3c_repl
352
353 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 r_1_minus_t_needed, do_kpoints_cubic_rpa)
358
359 IF (do_gamma_rpa) CALL dbt_batched_contract_init(t_dm_occ(1))
360 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 do_kpoints_cubic_rpa)
364 IF (r_1_minus_s_needed) THEN
365
366 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 qs_env%mp2_env%ri_rpa_im_time%eps_compress)
371
372 IF (do_gamma_rpa .AND. i_mem == 1) THEN
373 CALL get_tensor_occupancy(t_3c_o(1, 1), nze, occ)
374 nze_o = nze_o + nze
375 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 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 flop=flops_1_occ)
392 CALL timestop(handle3)
393
394 dbcsr_nflop = dbcsr_nflop + flops_1_occ
395
396 END IF
397 END DO
398
399 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 CALL timeset(routinen//"_copy_M_occ_t", handle3)
405 CALL dbt_copy(t_3c_m_occ_tmp, t_3c_m_occ, order=[1, 3, 2], move_data=.true.)
406 CALL dbt_filter(t_3c_m_occ, eps_filter)
407 CALL timestop(handle3)
408
409 IF (do_gamma_rpa) THEN
410 CALL get_tensor_occupancy(t_3c_m_occ, nze, occ)
411 nze_m_occ = nze_m_occ + nze
412 occ_m_occ = occ_m_occ + occ
413 END IF
414
415 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 r_2_minus_s_minus_t_needed, do_kpoints_cubic_rpa)
419
420 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 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 t_3c_o_virt(i_cell_r_2_minus_s_minus_t, i_cell_r_1_minus_t), move_data=.true.)
428
429 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 flop=flops_1_virt)
442 CALL timestop(handle3)
443
444 dbcsr_nflop = dbcsr_nflop + flops_1_virt
445
446 END IF
447 END DO
448
449 CALL timeset(routinen//"_copy_M_virt_t", handle3)
450 CALL dbt_copy(t_3c_m_virt_tmp, t_3c_m_virt, move_data=.true.)
451 CALL dbt_filter(t_3c_m_virt, eps_filter)
452 CALL timestop(handle3)
453
454 IF (do_gamma_rpa) THEN
455 CALL get_tensor_occupancy(t_3c_m_virt, nze, occ)
456 nze_m_virt = nze_m_virt + nze
457 occ_m_virt = occ_m_virt + occ
458 END IF
459
460 flops_2 = 0
461
462 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 unit_nr=unit_nr_dbcsr)
475
476 CALL timestop(handle3)
477
478 first_cycle_im_time = .false.
479
480 IF (jquad == 1 .AND. flops_2 == 0) &
481 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 IF (do_gamma_rpa) CALL dbt_batched_contract_finalize(t_dm_virt(1))
487 END DO
488
489 CALL dbt_batched_contract_finalize(t_p, unit_nr=unit_nr_dbcsr)
490
491 CALL dbt_create(mat_p_global%matrix, t_p_tmp)
492 CALL dbt_copy(t_p, t_p_tmp, move_data=.true.)
493 CALL dbt_copy_tensor_to_matrix(t_p_tmp, mat_p_global%matrix)
494 CALL dbt_destroy(t_p_tmp)
495
496 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 CALL dbcsr_add(mat_p_omega(jquad, i_cell_t)%matrix, mat_p_global%matrix, 1.0_dp, 1.0_dp)
501 ELSE
502 CALL timeset(routinen//"_Fourier_transform", handle3)
503
504 ! Fourier transform of P(it) to P(iw)
505 first_cycle_omega_loop = .true.
506
507 tau = tau_tj(jquad)
508
509 DO iquad = 1, num_integ_points
510
511 omega = tj(iquad)
512 weight = weights_cos_tf_t_to_w(iquad, jquad)
513
514 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 CALL dbcsr_scale(mat_p_global%matrix, cos(omega*tau)*weight)
518 ELSE
519 CALL dbcsr_scale(mat_p_global%matrix, cos(omega*tau)/cos(omega_old*tau)*weight/weight_old)
520 END IF
521
522 CALL dbcsr_add(mat_p_omega(iquad, i_cell_t)%matrix, mat_p_global%matrix, 1.0_dp, 1.0_dp)
523
524 first_cycle_omega_loop = .false.
525
526 omega_old = omega
527 weight_old = weight
528
529 END DO
530
531 CALL timestop(handle3)
532 END IF
533
534 END DO
535
536 CALL timestop(handle2)
537
538 CALL dbt_batched_contract_finalize(t_3c_m_occ_tmp)
539 CALL dbt_batched_contract_finalize(t_3c_m_virt_tmp)
540 CALL dbt_batched_contract_finalize(t_3c_m_occ)
541 CALL dbt_batched_contract_finalize(t_3c_m_virt)
542
543 DO i = 1, SIZE(t_3c_o_occ, 1)
544 DO j = 1, SIZE(t_3c_o_occ, 2)
545 CALL dbt_batched_contract_finalize(t_3c_o_occ(i, j))
546 END DO
547 END DO
548
549 DO i = 1, SIZE(t_3c_o_virt, 1)
550 DO j = 1, SIZE(t_3c_o_virt, 2)
551 CALL dbt_batched_contract_finalize(t_3c_o_virt(i, j))
552 END DO
553 END DO
554
555 CALL dbt_destroy(t_p)
556 DO i_cell = 1, num_cells_dm
557 CALL dbt_destroy(t_dm_virt(i_cell))
558 CALL dbt_destroy(t_dm_occ(i_cell))
559 END DO
560
561 CALL dbt_destroy(t_3c_m_occ_tmp)
562 CALL dbt_destroy(t_3c_m_virt_tmp)
563 DEALLOCATE (t_dm_virt)
564 DEALLOCATE (t_dm_occ)
565
566 CALL para_env%sync()
567 t2 = m_walltime()
568
569 dbcsr_time = dbcsr_time + t2 - t1_flop
570
571 IF (unit_nr > 0) THEN
572 WRITE (unit_nr, '(/T3,A,1X,I3)') &
573 'RPA_LOW_SCALING_INFO| Info for time point', jquad
574 WRITE (unit_nr, '(T6,A,T56,F25.1)') &
575 'Execution time (s):', t2 - t1
576 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
577 '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 'Occupancy of D virt:', real(nze_dm_virt, dp), '/', occ_dm_virt*100, '%'
580 IF (do_gamma_rpa) THEN
581 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
582 '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 '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 'Occupancy of M virt:', real(nze_m_virt, dp), '/', occ_m_virt*100, '%'
587 END IF
588 WRITE (unit_nr, *)
589 CALL m_flush(unit_nr)
590 END IF
591
592 END DO ! time points
593
594 CALL dbt_destroy(t_3c_m_occ)
595 CALL dbt_destroy(t_3c_m_virt)
596
597 DO i = 1, SIZE(t_3c_o, 1)
598 DO j = 1, SIZE(t_3c_o, 2)
599 CALL dbt_destroy(t_3c_o_occ(i, j))
600 CALL dbt_destroy(t_3c_o_virt(i, j))
601 END DO
602 END DO
603
604 CALL clean_up(mat_dm_occ_global, mat_dm_virt_global)
605
606 CALL timestop(handle)
607
608 END SUBROUTINE compute_mat_p_omega
609
610! **************************************************************************************************
611!> \brief ...
612!> \param mat_P_omega ...
613! **************************************************************************************************
614 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 DO jquad = 1, SIZE(mat_p_omega, 1)
620 DO i_kp = 1, SIZE(mat_p_omega, 2)
621
622 CALL dbcsr_set(mat_p_omega(jquad, i_kp)%matrix, 0.0_dp)
623
624 END DO
625 END DO
626
627 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 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 matrix_s, ispin, &
661 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 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
693 REAL(kind=dp) :: tau
694
695 CALL timeset(routinen, handle)
696
697 IF (memory_info .AND. unit_nr > 0) WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
698 "RPA_LOW_SCALING_INFO| Started with time point: ", jquad
699
700 tau = tau_tj(jquad)
701
702 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 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 remove_occ=.true., remove_virt=.false., first_jquad=1)
713
714 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 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 alloc_dm=(jquad == 1))
725
726 num_cells_dm = 1
727
728 ELSE
729
730 num_cells_dm = 1
731
732 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 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 DO jjb = 1, nrow_local
747 DO iib = 1, ncol_local
748 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 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 fm_mo_coeff_occ%local_data(jjb, iib)*exp(tau*0.5_dp*(eigenval(i_global) - e_fermi))
755 ELSE
756 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 col_indices=col_indices)
768
769 ! the same for virt
770 DO jjb = 1, nrow_local
771 DO iib = 1, ncol_local
772 i_global = col_indices(iib)
773
774 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 fm_mo_coeff_virt%local_data(jjb, iib)*exp(-tau*0.5_dp*(eigenval(i_global) - e_fermi))
777 ELSE
778 fm_mo_coeff_virt_scaled%local_data(jjb, iib) = 0.0_dp
779 END IF
780
781 END DO
782 END DO
783
784 CALL para_env%sync()
785
786 size_dm_occ = nmo
787 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 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 matrix_c=fm_scaled_dm_virt_tau)
796
797 IF (jquad == 1) THEN
798
799 ! transfer fm density matrices to dbcsr matrix
800 NULLIFY (mat_dm_occ_global)
801 CALL dbcsr_allocate_matrix_set(mat_dm_occ_global, num_integ_points, 1)
802
803 DO iquad = 1, num_integ_points
804
805 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 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 keep_sparsity=.false.)
817
818 CALL dbcsr_filter(mat_dm_occ_global(jquad, 1)%matrix, eps_filter)
819
820 IF (jquad == 1) THEN
821
822 NULLIFY (mat_dm_virt_global)
823 CALL dbcsr_allocate_matrix_set(mat_dm_virt_global, num_integ_points, 1)
824
825 END IF
826
827 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 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 keep_sparsity=.false.)
834
835 CALL dbcsr_filter(mat_dm_virt_global(jquad, 1)%matrix, eps_filter)
836
837 ! release memory
838 IF (jquad > 1) THEN
839 CALL dbcsr_set(mat_dm_occ_global(jquad - 1, 1)%matrix, 0.0_dp)
840 CALL dbcsr_set(mat_dm_virt_global(jquad - 1, 1)%matrix, 0.0_dp)
841 CALL dbcsr_filter(mat_dm_occ_global(jquad - 1, 1)%matrix, 0.0_dp)
842 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 CALL timestop(handle)
848
849 END SUBROUTINE compute_mat_dm_global
850
851! **************************************************************************************************
852!> \brief ...
853!> \param mat_dm_occ_global ...
854!> \param mat_dm_virt_global ...
855! **************************************************************************************************
856 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 CALL dbcsr_deallocate_matrix_set(mat_dm_occ_global)
860 CALL dbcsr_deallocate_matrix_set(mat_dm_virt_global)
861
862 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 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 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 CALL timeset(routinen, handle)
892
893 ! only imaginary wavefunctions supported in kpoint cubic scaling RPA
894 cpassert(kpoint%use_real_wfn .EQV. .false.)
895
896 ! work matrix
897 mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)
898 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 cpassert(nao == nmo)
903
904 ALLOCATE (exp_scaling(nmo))
905
906 CALL cp_fm_get_info(mo_set%mo_coeff, matrix_struct=matrix_struct)
907 CALL cp_fm_create(fwork, matrix_struct)
908
909 CALL get_kpoint_info(kpoint, kp_range=kp_range)
910 kplocal = kp_range(2) - kp_range(1) + 1
911
912 DO ikpgr = 1, kplocal
913 kp => kpoint%kp_env(ikpgr)%kpoint_env
914 nspin = SIZE(kp%mos, 2)
915 DO ispin = 1, nspin
916 mo_set => kp%mos(1, ispin)
917 CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
918 rpmat => kp%wmat(1, ispin)
919 cpmat => kp%wmat(2, ispin)
920 CALL get_mo_set(mo_set, occupation_numbers=occupation)
921 CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
922
923 IF (remove_virt) THEN
924 CALL cp_fm_column_scale(fwork, occupation)
925 END IF
926 IF (remove_occ) THEN
927 CALL cp_fm_column_scale(fwork, 2.0_dp/real(nspin, kind=dp) - occupation)
928 END IF
929
930 ! proper spin
931 IF (nspin == 1) THEN
932 CALL cp_fm_scale(0.5_dp, fwork)
933 END IF
934
935 DO i_mo = 1, nmo
936
937 IF (abs(tau*0.5_dp*(eigenvalues(i_mo) - e_fermi)) < stabilize_exp) THEN
938 exp_scaling(i_mo) = exp(-abs(tau*(eigenvalues(i_mo) - e_fermi)))
939 ELSE
940 exp_scaling(i_mo) = 0.0_dp
941 END IF
942 END DO
943
944 CALL cp_fm_column_scale(fwork, exp_scaling)
945
946 ! Re(c)*Re(c)
947 CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, rpmat)
948 mo_set => kp%mos(2, ispin)
949 ! Im(c)*Re(c)
950 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 CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, fwork, mo_set%mo_coeff, 1.0_dp, cpmat)
953
954 CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
955
956 IF (remove_virt) THEN
957 CALL cp_fm_column_scale(fwork, occupation)
958 END IF
959 IF (remove_occ) THEN
960 CALL cp_fm_column_scale(fwork, 2.0_dp/real(nspin, kind=dp) - occupation)
961 END IF
962
963 ! proper spin
964 IF (nspin == 1) THEN
965 CALL cp_fm_scale(0.5_dp, fwork)
966 END IF
967
968 DO i_mo = 1, nmo
969 IF (abs(tau*0.5_dp*(eigenvalues(i_mo) - e_fermi)) < stabilize_exp) THEN
970 exp_scaling(i_mo) = exp(-abs(tau*(eigenvalues(i_mo) - e_fermi)))
971 ELSE
972 exp_scaling(i_mo) = 0.0_dp
973 END IF
974 END DO
975
976 CALL cp_fm_column_scale(fwork, exp_scaling)
977 ! Im(c)*Im(c)
978 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 CALL cp_fm_release(fwork)
985 DEALLOCATE (exp_scaling)
986
987 CALL timestop(handle)
988
989 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 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 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_global_work, matrix_s_kp
1025 TYPE(kpoint_type), POINTER :: kpoints
1026 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1027
1028 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 cell=cell)
1035
1036 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 DO i_dim = 1, 3
1041 cell_grid_dm(i_dim) = (kpoints%nkp_grid(i_dim)/2)*2 - 1
1042 END DO
1043
1044 num_cells_dm = cell_grid_dm(1)*cell_grid_dm(2)*cell_grid_dm(3)
1045
1046 NULLIFY (mat_dm_global_work)
1047 CALL dbcsr_allocate_matrix_set(mat_dm_global_work, nspin, num_cells_dm)
1048
1049 DO jspin = 1, nspin
1050
1051 DO i_img = 1, num_cells_dm
1052
1053 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 matrix_type=dbcsr_type_no_symmetry)
1058
1059 CALL dbcsr_reserve_all_blocks(mat_dm_global_work(jspin, i_img)%matrix)
1060
1061 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 remove_occ=remove_occ, remove_virt=remove_virt)
1070
1071 ! overwrite the cell indices in kpoints
1072 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 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 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 IF (jquad == first_jquad) THEN
1083
1084 NULLIFY (mat_dm_global)
1085 ALLOCATE (mat_dm_global(first_jquad:num_integ_points, num_cells_dm))
1086
1087 DO iquad = first_jquad, num_integ_points
1088 DO i_img = 1, num_cells_dm
1089 NULLIFY (mat_dm_global(iquad, i_img)%matrix)
1090 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 matrix_type=dbcsr_type_no_symmetry)
1094
1095 END DO
1096 END DO
1097
1098 END IF
1099
1100 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 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 mat_dm_global_work(ispin, i_img)%matrix)
1107
1108 END DO
1109
1110 CALL dbcsr_deallocate_matrix_set(mat_dm_global_work)
1111
1112 CALL timestop(handle)
1113
1114 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 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 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_global_work, matrix_s_kp
1141 TYPE(kpoint_type), POINTER :: kpoints_g
1142 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1143
1144 CALL timeset(routinen, handle)
1145
1146 NULLIFY (matrix_s_kp, mos)
1147
1148 CALL get_qs_env(qs_env, &
1149 matrix_s_kp=matrix_s_kp, &
1150 mos=mos)
1151
1152 kpoints_g => qs_env%mp2_env%ri_rpa_im_time%kpoints_G
1153
1154 nspin = SIZE(mos)
1155
1156 num_cells_dm = 1
1157
1158 NULLIFY (mat_dm_global_work)
1159 CALL dbcsr_allocate_matrix_set(mat_dm_global_work, nspin, num_cells_dm)
1160
1161 ! if necessaray, allocate mat_dm_global
1162 IF (alloc_dm) THEN
1163
1164 NULLIFY (mat_dm_global)
1165 ALLOCATE (mat_dm_global(1:num_integ_points, num_cells_dm))
1166
1167 DO iquad = 1, num_integ_points
1168 NULLIFY (mat_dm_global(iquad, 1)%matrix)
1169 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 matrix_type=dbcsr_type_no_symmetry)
1173
1174 END DO
1175
1176 END IF
1177
1178 DO jspin = 1, nspin
1179
1180 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 matrix_type=dbcsr_type_no_symmetry)
1184
1185 CALL dbcsr_reserve_all_blocks(mat_dm_global_work(jspin, 1)%matrix)
1186
1187 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 remove_occ=remove_occ, remove_virt=remove_virt)
1194
1195 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 mat_dm_global_work(ispin, 1)%matrix)
1199
1200 CALL dbcsr_deallocate_matrix_set(mat_dm_global_work)
1201
1202 CALL timestop(handle)
1203
1204 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 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 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1226 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
1227 REAL(kind=dp) :: contribution, weight_im, weight_re
1228 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1229 REAL(kind=dp), DIMENSION(:), POINTER :: wkp
1230 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1237
1238 CALL timeset(routinen, handle)
1239
1240 NULLIFY (xkp, wkp)
1241
1242 CALL cp_fm_create(fm_mat_work, kpoints_g%kp_env(1)%kpoint_env%wmat(1, 1)%matrix_struct)
1243 CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
1244
1245 CALL get_kpoint_info(kpoints_g, nkp=nkp, xkp=xkp, wkp=wkp)
1246 index_to_cell => kpoints_g%index_to_cell
1247 num_cells = SIZE(index_to_cell, 2)
1248
1249 nspin = SIZE(mat_dm_global_work, 1)
1250
1251 mo_set => kpoints_g%kp_env(1)%kpoint_env%mos(1, 1)
1252 CALL get_mo_set(mo_set, nao=nao)
1253
1254 ALLOCATE (atom_from_ao_index(nao))
1255
1256 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 col_indices=col_indices)
1263
1264 NULLIFY (cell, particle_set)
1265 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1266 CALL get_cell(cell=cell, h=hmat)
1267
1268 iatom_old = 0
1269 jatom_old = 0
1270
1271 DO ispin = 1, nspin
1272
1273 CALL dbcsr_set(mat_dm_global_work(ispin, 1)%matrix, 0.0_dp)
1274
1275 DO ik = 1, nkp
1276
1277 kp => kpoints_g%kp_env(ik)%kpoint_env
1278 rpmat => kp%wmat(1, ispin)
1279 cpmat => kp%wmat(2, ispin)
1280
1281 DO irow = 1, nrow_local
1282 DO jcol = 1, ncol_local
1283
1284 iatom = atom_from_ao_index(row_indices(irow))
1285 jatom = atom_from_ao_index(col_indices(jcol))
1286
1287 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 cell, index_to_cell, hmat, particle_set)
1292
1293 iatom_old = iatom
1294 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 weight_im*cpmat%local_data(irow, jcol)
1301
1302 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 CALL copy_fm_to_dbcsr(fm_mat_work, mat_dm_global_work(ispin, 1)%matrix, keep_sparsity=.false.)
1310 CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
1311
1312 END DO
1313
1314 CALL cp_fm_release(fm_mat_work)
1315 DEALLOCATE (atom_from_ao_index)
1316
1317 CALL timestop(handle)
1318
1319 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: wkp
1339 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 CALL timeset(routinen, handle)
1345
1346 NULLIFY (xkp, wkp)
1347
1348 NULLIFY (mat_work_re)
1349 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 matrix_type=dbcsr_type_no_symmetry)
1353
1354 NULLIFY (mat_work_im)
1355 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 matrix_type=dbcsr_type_no_symmetry)
1359
1360 CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, wkp=wkp)
1361
1362 nspin = SIZE(mat_dm_global_work, 1)
1363
1364 cpassert(SIZE(mat_dm_global_work, 2) == SIZE(index_to_cell, 2))
1365
1366 DO ispin = 1, nspin
1367
1368 DO icell = 1, SIZE(mat_dm_global_work, 2)
1369
1370 CALL dbcsr_set(mat_dm_global_work(ispin, icell)%matrix, 0.0_dp)
1371
1372 END DO
1373
1374 END DO
1375
1376 DO ispin = 1, nspin
1377
1378 DO ik = 1, nkp
1379
1380 kp => kpoints%kp_env(ik)%kpoint_env
1381 rpmat => kp%wmat(1, ispin)
1382 cpmat => kp%wmat(2, ispin)
1383
1384 CALL copy_fm_to_dbcsr(rpmat, mat_work_re, keep_sparsity=.false.)
1385 CALL copy_fm_to_dbcsr(cpmat, mat_work_im, keep_sparsity=.false.)
1386
1387 DO icell = 1, SIZE(mat_dm_global_work, 2)
1388
1389 xcell = index_to_cell(1, icell)
1390 ycell = index_to_cell(2, icell)
1391 zcell = index_to_cell(3, icell)
1392
1393 arg = real(xcell, dp)*xkp(1, ik) + real(ycell, dp)*xkp(2, ik) + real(zcell, dp)*xkp(3, ik)
1394 coskl = wkp(ik)*cos(twopi*arg)
1395 sinkl = wkp(ik)*sin(twopi*arg)
1396
1397 CALL dbcsr_add(mat_dm_global_work(ispin, icell)%matrix, mat_work_re, 1.0_dp, coskl)
1398 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 CALL dbcsr_release_p(mat_work_re)
1406 CALL dbcsr_release_p(mat_work_im)
1407
1408 CALL timestop(handle)
1409
1410 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 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 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_unsorted
1432 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_unsorted
1433 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 CALL timeset(routinen, handle)
1438
1439 CALL get_cell(cell=cell, h=hmat)
1440
1441 num_cells = cell_grid(1)*cell_grid(2)*cell_grid(3)
1442 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 cpassert(cell_grid(1) .NE. itm(1)*2)
1447 cpassert(cell_grid(2) .NE. itm(2)*2)
1448 cpassert(cell_grid(3) .NE. itm(3)*2)
1449
1450 IF (ASSOCIATED(cell_to_index)) DEALLOCATE (cell_to_index)
1451 IF (ASSOCIATED(index_to_cell)) DEALLOCATE (index_to_cell)
1452
1453 ALLOCATE (cell_to_index_unsorted(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
1454 cell_to_index_unsorted(:, :, :) = 0
1455
1456 ALLOCATE (index_to_cell_unsorted(3, num_cells))
1457 index_to_cell_unsorted(:, :) = 0
1458
1459 ALLOCATE (cell_to_index(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
1460 cell_to_index(:, :, :) = 0
1461
1462 ALLOCATE (index_to_cell(3, num_cells))
1463 index_to_cell(:, :) = 0
1464
1465 ALLOCATE (abs_cell_vectors(1:num_cells))
1466
1467 cell_counter = 0
1468
1469 DO xcell = -itm(1), itm(1)
1470 DO ycell = -itm(2), itm(2)
1471 DO zcell = -itm(3), itm(3)
1472
1473 cell_counter = cell_counter + 1
1474 cell_to_index_unsorted(xcell, ycell, zcell) = cell_counter
1475
1476 index_to_cell_unsorted(1, cell_counter) = xcell
1477 index_to_cell_unsorted(2, cell_counter) = ycell
1478 index_to_cell_unsorted(3, cell_counter) = zcell
1479
1480 cell_vector(1:3) = matmul(hmat, real(index_to_cell_unsorted(1:3, cell_counter), dp))
1481
1482 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 DO i_cell = 1, num_cells/2 + 1
1491
1492 index_min_dist = minloc(abs_cell_vectors(1:num_cells/2 + 1), dim=1)
1493
1494 xcell = index_to_cell_unsorted(1, index_min_dist)
1495 ycell = index_to_cell_unsorted(2, index_min_dist)
1496 zcell = index_to_cell_unsorted(3, index_min_dist)
1497
1498 index_to_cell(1, i_cell) = xcell
1499 index_to_cell(2, i_cell) = ycell
1500 index_to_cell(3, i_cell) = zcell
1501
1502 cell_to_index(xcell, ycell, zcell) = i_cell
1503
1504 abs_cell_vectors(index_min_dist) = 10000000000.0_dp
1505
1506 END DO
1507
1508 ! now all the remaining cells
1509 DO i_cell = num_cells/2 + 2, num_cells
1510
1511 index_min_dist = minloc(abs_cell_vectors(1:num_cells), dim=1)
1512
1513 xcell = index_to_cell_unsorted(1, index_min_dist)
1514 ycell = index_to_cell_unsorted(2, index_min_dist)
1515 zcell = index_to_cell_unsorted(3, index_min_dist)
1516
1517 index_to_cell(1, i_cell) = xcell
1518 index_to_cell(2, i_cell) = ycell
1519 index_to_cell(3, i_cell) = zcell
1520
1521 cell_to_index(xcell, ycell, zcell) = i_cell
1522
1523 abs_cell_vectors(index_min_dist) = 10000000000.0_dp
1524
1525 END DO
1526
1527 DEALLOCATE (index_to_cell_unsorted, cell_to_index_unsorted, abs_cell_vectors)
1528
1529 CALL timestop(handle)
1530
1531 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 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 CALL timeset(routinen, handle)
1563
1564 IF (do_kpoints_cubic_rpa) THEN
1565
1566 x_cell_r = index_to_cell_3c(1, i_cell_r)
1567 y_cell_r = index_to_cell_3c(2, i_cell_r)
1568 z_cell_r = index_to_cell_3c(3, i_cell_r)
1569
1570 x_cell_s = index_to_cell_dm(1, i_cell_s)
1571 y_cell_s = index_to_cell_dm(2, i_cell_s)
1572 z_cell_s = index_to_cell_dm(3, i_cell_s)
1573
1574 x_cell_r_minus_s = x_cell_r - x_cell_s
1575 y_cell_r_minus_s = y_cell_r - y_cell_s
1576 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 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 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 IF (i_cell_r_minus_s == 0) THEN
1589
1590 r_minus_s_needed = .false.
1591 i_cell_r_minus_s = 0
1592
1593 ELSE
1594
1595 r_minus_s_needed = .true.
1596
1597 END IF
1598
1599 ELSE
1600
1601 i_cell_r_minus_s = 0
1602 r_minus_s_needed = .false.
1603
1604 END IF
1605
1606 ELSE ! no k-points
1607
1608 r_minus_s_needed = .true.
1609 i_cell_r_minus_s = 1
1610
1611 END IF
1612
1613 CALL timestop(handle)
1614
1615 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 SUBROUTINE get_diff_diff_index_3c(i_cell_R, i_cell_S, i_cell_T, i_cell_R_minus_S_minus_T, &
1630 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 CALL timeset(routinen, handle)
1650
1651 IF (do_kpoints_cubic_rpa) THEN
1652
1653 x_cell_r = index_to_cell_3c(1, i_cell_r)
1654 y_cell_r = index_to_cell_3c(2, i_cell_r)
1655 z_cell_r = index_to_cell_3c(3, i_cell_r)
1656
1657 x_cell_s = index_to_cell_dm(1, i_cell_s)
1658 y_cell_s = index_to_cell_dm(2, i_cell_s)
1659 z_cell_s = index_to_cell_dm(3, i_cell_s)
1660
1661 x_cell_t = index_to_cell_dm(1, i_cell_t)
1662 y_cell_t = index_to_cell_dm(2, i_cell_t)
1663 z_cell_t = index_to_cell_dm(3, i_cell_t)
1664
1665 x_cell_r_minus_s_minus_t = x_cell_r - x_cell_s - x_cell_t
1666 y_cell_r_minus_s_minus_t = y_cell_r - y_cell_s - y_cell_t
1667 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 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 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 IF (i_cell_r_minus_s_minus_t == 0) THEN
1682
1683 r_minus_s_minus_t_needed = .false.
1684
1685 ELSE
1686
1687 r_minus_s_minus_t_needed = .true.
1688
1689 END IF
1690
1691 ELSE
1692
1693 i_cell_r_minus_s_minus_t = 0
1694 r_minus_s_minus_t_needed = .false.
1695
1696 END IF
1697
1698 ! no k-kpoints
1699 ELSE
1700
1701 r_minus_s_minus_t_needed = .true.
1702 i_cell_r_minus_s_minus_t = 1
1703
1704 END IF
1705
1706 CALL timestop(handle)
1707
1708 END SUBROUTINE get_diff_diff_index_3c
1709
1710END MODULE rpa_im_time
Handles all functions related to the CELL.
Definition cell_types.F:15
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:195
subroutine, public dbcsr_release_p(matrix)
...
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_init_p(matrix)
...
subroutine, public dbcsr_filter(matrix, eps)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_clear(matrix)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
subroutine, public dbcsr_reserve_all_blocks(matrix)
Reserves all blocks.
DBCSR operations in CP2K.
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
basic linear algebra operations for full matrices
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
represent the structure of a full matrix
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
This is the start of a dbt_api, all publically needed functions are exported here....
Definition dbt_api.F:17
Types and set/get functions for HFX.
Definition hfx_types.F:15
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public int_8
Definition kinds.F:54
integer, parameter, public dp
Definition kinds.F:34
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:130
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:147
Definition of mathematical constants and functions.
real(kind=dp), parameter, public twopi
Interface to the message passing library MPI.
Types needed for MP2 calculations.
Definition mp2_types.F:14
basic linear algebra operations for full matrixes
Define the data structure for the particle information.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Utility methods to build 3-center integral tensors of various types.
subroutine, public create_2c_tensor(t2c, dist_1, dist_2, pgrid, sizes_1, sizes_2, order, name)
...
Utility methods to build 3-center integral tensors of various types.
Definition qs_tensors.F:11
subroutine, public decompress_tensor(tensor, blk_indices, compressed, eps)
...
subroutine, public get_tensor_occupancy(tensor, nze, occ)
...
Utility routines for GW with imaginary time.
subroutine, public compute_weight_re_im(weight_re, weight_im, num_cells, iatom, jatom, xkp, wkp_w, cell, index_to_cell, hmat, particle_set)
...
subroutine, public get_atom_index_from_basis_function_index(qs_env, atom_from_basis_index, basis_size, basis_type, first_bf_from_atom)
...
Routines for low-scaling RPA/GW with imaginary time.
Definition rpa_im_time.F:13
subroutine, public zero_mat_p_omega(mat_p_omega)
...
subroutine, public compute_mat_dm_global(fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, tau_tj, num_integ_points, nmo, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm_occ_global, mat_dm_virt_global, matrix_s, ispin, eigenval, e_fermi, eps_filter, memory_info, unit_nr, jquad, do_kpoints_cubic_rpa, do_kpoints_from_gamma, qs_env, num_cells_dm, index_to_cell_dm, para_env)
...
subroutine, public compute_mat_p_omega(mat_p_omega, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_p_global, matrix_s, ispin, t_3c_m, t_3c_o, t_3c_o_compressed, t_3c_o_ind, starts_array_mc, ends_array_mc, starts_array_mc_block, ends_array_mc_block, weights_cos_tf_t_to_w, tj, tau_tj, e_fermi, eps_filter, alpha, eps_filter_im_time, eigenval, nmo, num_integ_points, cut_memory, unit_nr, mp2_env, para_env, qs_env, do_kpoints_from_gamma, index_to_cell_3c, cell_to_index_3c, has_mat_p_blocks, do_ri_sos_laplace_mp2, dbcsr_time, dbcsr_nflop)
...
subroutine, public compute_periodic_dm(mat_dm_global, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, remove_occ, remove_virt, alloc_dm)
...
subroutine, public init_cell_index_rpa(cell_grid, cell_to_index, index_to_cell, cell)
...
subroutine, public compute_transl_dm(mat_dm_global, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, eps_filter, num_cells_dm, index_to_cell_dm, remove_occ, remove_virt, first_jquad)
...
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
keeps the information about the structure of a full matrix
represent a full matrix
Keeps information about a specific k-point.
Contains information about kpoints.
stores all the informations relevant to an mpi environment