(git:ccc2433)
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-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  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  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  INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_1, dist_2, mc_ranges, size_dm, &
186  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  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  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_occ_global, mat_dm_virt_global
200  TYPE(dbt_pgrid_type) :: pgrid_2d
201  TYPE(dbt_type) :: t_3c_m_occ, t_3c_m_occ_tmp, t_3c_m_virt, &
202  t_3c_m_virt_tmp, t_dm, t_dm_tmp, t_p, &
203  t_p_tmp
204  TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_dm_occ, t_dm_virt
205  TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_o_occ, t_3c_o_virt
206  TYPE(mp_comm_type) :: comm_2d
207 
208  CALL timeset(routinen, handle)
209 
210  memory_info = mp2_env%ri_rpa_im_time%memory_info
211  IF (memory_info) THEN
212  unit_nr_dbcsr = unit_nr
213  ELSE
214  unit_nr_dbcsr = 0
215  END IF
216 
217  do_kpoints_cubic_rpa = qs_env%mp2_env%ri_rpa_im_time%do_im_time_kpoints
218  do_gamma_rpa = .NOT. do_kpoints_cubic_rpa
219  num_3c_repl = maxval(cell_to_index_3c)
220 
221  first_cycle_im_time = .true.
222  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  DO i = 1, SIZE(t_3c_o, 1)
224  DO j = 1, SIZE(t_3c_o, 2)
225  CALL dbt_create(t_3c_o(i, j), t_3c_o_occ(i, j))
226  CALL dbt_create(t_3c_o(i, j), t_3c_o_virt(i, j))
227  END DO
228  END DO
229 
230  CALL dbt_create(t_3c_m, t_3c_m_occ, name="M occ (RI | AO AO)")
231  CALL dbt_create(t_3c_m, t_3c_m_virt, name="M virt (RI | AO AO)")
232 
233  ALLOCATE (mc_ranges(cut_memory + 1))
234  mc_ranges(:cut_memory) = starts_array_mc_block(:)
235  mc_ranges(cut_memory + 1) = ends_array_mc_block(cut_memory) + 1
236 
237  DO jquad = 1, num_integ_points
238 
239  CALL para_env%sync()
240  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  num_cells_dm, index_to_cell_dm, para_env)
249 
250  ALLOCATE (t_dm_virt(num_cells_dm))
251  ALLOCATE (t_dm_occ(num_cells_dm))
252  CALL dbcsr_get_info(mat_p_global%matrix, distribution=dist_p)
253  CALL dbcsr_distribution_get(dist_p, group=comm_2d_handle, nprows=pdims_2d(1), npcols=pdims_2d(2))
254  CALL comm_2d%set_handle(comm_2d_handle)
255 
256  pgrid_2d = dbt_nd_mp_comm(comm_2d, [1], [2], pdims_2d=pdims_2d)
257  ALLOCATE (size_p(dbt_nblks_total(t_3c_m, 1)))
258  CALL dbt_get_info(t_3c_m, blk_size_1=size_p)
259 
260  ALLOCATE (size_dm(dbt_nblks_total(t_3c_o(1, 1), 3)))
261  CALL dbt_get_info(t_3c_o(1, 1), blk_size_3=size_dm)
262  CALL create_2c_tensor(t_dm, dist_1, dist_2, pgrid_2d, size_dm, size_dm, name="D (AO | AO)")
263  DEALLOCATE (size_dm)
264  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  DEALLOCATE (size_p)
267  DEALLOCATE (dist_1, dist_2)
268  CALL dbt_pgrid_destroy(pgrid_2d)
269 
270  DO i_cell = 1, num_cells_dm
271  CALL dbt_create(t_dm, t_dm_virt(i_cell), name="D virt (AO | AO)")
272  CALL dbt_create(mat_dm_virt_global(jquad, i_cell)%matrix, t_dm_tmp)
273  CALL dbt_copy_matrix_to_tensor(mat_dm_virt_global(jquad, i_cell)%matrix, t_dm_tmp)
274  CALL dbt_copy(t_dm_tmp, t_dm_virt(i_cell), move_data=.true.)
275  CALL dbcsr_clear(mat_dm_virt_global(jquad, i_cell)%matrix)
276 
277  CALL dbt_create(t_dm, t_dm_occ(i_cell), name="D occ (AO | AO)")
278  CALL dbt_copy_matrix_to_tensor(mat_dm_occ_global(jquad, i_cell)%matrix, t_dm_tmp)
279  CALL dbt_copy(t_dm_tmp, t_dm_occ(i_cell), move_data=.true.)
280  CALL dbt_destroy(t_dm_tmp)
281  CALL dbcsr_clear(mat_dm_occ_global(jquad, i_cell)%matrix)
282  END DO
283 
284  CALL get_tensor_occupancy(t_dm_occ(1), nze_dm_occ, occ_dm_occ)
285  CALL get_tensor_occupancy(t_dm_virt(1), nze_dm_virt, occ_dm_virt)
286 
287  CALL dbt_destroy(t_dm)
288 
289  CALL dbt_create(t_3c_o_occ(1, 1), t_3c_m_occ_tmp, name="M (RI AO | AO)")
290  CALL dbt_create(t_3c_o_virt(1, 1), t_3c_m_virt_tmp, name="M (RI AO | AO)")
291 
292  CALL timeset(routinen//"_contract", handle2)
293 
294  CALL para_env%sync()
295  t1_flop = m_walltime()
296 
297  DO i = 1, SIZE(t_3c_o_occ, 1)
298  DO j = 1, SIZE(t_3c_o_occ, 2)
299  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  DO i = 1, SIZE(t_3c_o_virt, 1)
303  DO j = 1, SIZE(t_3c_o_virt, 2)
304  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  CALL dbt_batched_contract_init(t_3c_m_occ_tmp, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
308  CALL dbt_batched_contract_init(t_3c_m_virt_tmp, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
309  CALL dbt_batched_contract_init(t_3c_m_occ, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
310  CALL dbt_batched_contract_init(t_3c_m_virt, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
311 
312  DO i_cell_t = 1, num_cells_dm/2 + 1
313 
314  IF (.NOT. any(has_mat_p_blocks(i_cell_t, :, :, :, :))) cycle
315 
316  CALL dbt_batched_contract_init(t_p)
317 
318  IF (do_gamma_rpa) THEN
319  nze_o = 0
320  nze_m_virt = 0
321  nze_m_occ = 0
322  occ_m_virt = 0.0_dp
323  occ_m_occ = 0.0_dp
324  occ_o = 0.0_dp
325  END IF
326 
327  DO j_mem = 1, cut_memory
328 
329  CALL dbt_get_info(t_3c_o_occ(1, 1), nfull_total=bounds_3c)
330 
331  jbounds_1(:, 1) = [1, bounds_3c(1)]
332  jbounds_1(:, 2) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
333 
334  jbounds_2(:, 1) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
335 
336  IF (do_gamma_rpa) CALL dbt_batched_contract_init(t_dm_virt(1))
337 
338  DO i_mem = 1, cut_memory
339 
340  IF (.NOT. any(has_mat_p_blocks(i_cell_t, i_mem, j_mem, :, :))) cycle
341 
342  ibounds_1(:, 1) = [1, bounds_3c(1)]
343  ibounds_1(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
344 
345  ibounds_2(:, 1) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
346 
347  IF (unit_nr_dbcsr > 0) WRITE (unit=unit_nr_dbcsr, fmt="(T3,A,I3,1X,I3)") &
348  "RPA_LOW_SCALING_INFO| Memory Cut iteration", i_mem, j_mem
349 
350  DO i_cell_r_1 = 1, num_3c_repl
351 
352  DO i_cell_r_2 = 1, num_3c_repl
353 
354  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  r_1_minus_t_needed, do_kpoints_cubic_rpa)
359 
360  IF (do_gamma_rpa) CALL dbt_batched_contract_init(t_dm_occ(1))
361  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  do_kpoints_cubic_rpa)
365  IF (r_1_minus_s_needed) THEN
366 
367  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  qs_env%mp2_env%ri_rpa_im_time%eps_compress)
372 
373  IF (do_gamma_rpa .AND. i_mem == 1) THEN
374  CALL get_tensor_occupancy(t_3c_o(1, 1), nze, occ)
375  nze_o = nze_o + nze
376  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  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  flop=flops_1_occ)
393  CALL timestop(handle3)
394 
395  dbcsr_nflop = dbcsr_nflop + flops_1_occ
396 
397  END IF
398  END DO
399 
400  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  CALL timeset(routinen//"_copy_M_occ_t", handle3)
406  CALL dbt_copy(t_3c_m_occ_tmp, t_3c_m_occ, order=[1, 3, 2], move_data=.true.)
407  CALL dbt_filter(t_3c_m_occ, eps_filter)
408  CALL timestop(handle3)
409 
410  IF (do_gamma_rpa) THEN
411  CALL get_tensor_occupancy(t_3c_m_occ, nze, occ)
412  nze_m_occ = nze_m_occ + nze
413  occ_m_occ = occ_m_occ + occ
414  END IF
415 
416  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  r_2_minus_s_minus_t_needed, do_kpoints_cubic_rpa)
420 
421  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  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  t_3c_o_virt(i_cell_r_2_minus_s_minus_t, i_cell_r_1_minus_t), move_data=.true.)
429 
430  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  flop=flops_1_virt)
443  CALL timestop(handle3)
444 
445  dbcsr_nflop = dbcsr_nflop + flops_1_virt
446 
447  END IF
448  END DO
449 
450  CALL timeset(routinen//"_copy_M_virt_t", handle3)
451  CALL dbt_copy(t_3c_m_virt_tmp, t_3c_m_virt, move_data=.true.)
452  CALL dbt_filter(t_3c_m_virt, eps_filter)
453  CALL timestop(handle3)
454 
455  IF (do_gamma_rpa) THEN
456  CALL get_tensor_occupancy(t_3c_m_virt, nze, occ)
457  nze_m_virt = nze_m_virt + nze
458  occ_m_virt = occ_m_virt + occ
459  END IF
460 
461  flops_2 = 0
462 
463  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  unit_nr=unit_nr_dbcsr)
476 
477  CALL timestop(handle3)
478 
479  first_cycle_im_time = .false.
480 
481  IF (jquad == 1 .AND. flops_2 == 0) &
482  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  IF (do_gamma_rpa) CALL dbt_batched_contract_finalize(t_dm_virt(1))
488  END DO
489 
490  CALL dbt_batched_contract_finalize(t_p, unit_nr=unit_nr_dbcsr)
491 
492  CALL dbt_create(mat_p_global%matrix, t_p_tmp)
493  CALL dbt_copy(t_p, t_p_tmp, move_data=.true.)
494  CALL dbt_copy_tensor_to_matrix(t_p_tmp, mat_p_global%matrix)
495  CALL dbt_destroy(t_p_tmp)
496 
497  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  CALL dbcsr_add(mat_p_omega(jquad, i_cell_t)%matrix, mat_p_global%matrix, 1.0_dp, 1.0_dp)
502  ELSE
503  CALL timeset(routinen//"_Fourier_transform", handle3)
504 
505  ! Fourier transform of P(it) to P(iw)
506  first_cycle_omega_loop = .true.
507 
508  tau = tau_tj(jquad)
509 
510  DO iquad = 1, num_integ_points
511 
512  omega = tj(iquad)
513  weight = weights_cos_tf_t_to_w(iquad, jquad)
514 
515  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  CALL dbcsr_scale(mat_p_global%matrix, cos(omega*tau)*weight)
519  ELSE
520  CALL dbcsr_scale(mat_p_global%matrix, cos(omega*tau)/cos(omega_old*tau)*weight/weight_old)
521  END IF
522 
523  CALL dbcsr_add(mat_p_omega(iquad, i_cell_t)%matrix, mat_p_global%matrix, 1.0_dp, 1.0_dp)
524 
525  first_cycle_omega_loop = .false.
526 
527  omega_old = omega
528  weight_old = weight
529 
530  END DO
531 
532  CALL timestop(handle3)
533  END IF
534 
535  END DO
536 
537  CALL timestop(handle2)
538 
539  CALL dbt_batched_contract_finalize(t_3c_m_occ_tmp)
540  CALL dbt_batched_contract_finalize(t_3c_m_virt_tmp)
541  CALL dbt_batched_contract_finalize(t_3c_m_occ)
542  CALL dbt_batched_contract_finalize(t_3c_m_virt)
543 
544  DO i = 1, SIZE(t_3c_o_occ, 1)
545  DO j = 1, SIZE(t_3c_o_occ, 2)
546  CALL dbt_batched_contract_finalize(t_3c_o_occ(i, j))
547  END DO
548  END DO
549 
550  DO i = 1, SIZE(t_3c_o_virt, 1)
551  DO j = 1, SIZE(t_3c_o_virt, 2)
552  CALL dbt_batched_contract_finalize(t_3c_o_virt(i, j))
553  END DO
554  END DO
555 
556  CALL dbt_destroy(t_p)
557  DO i_cell = 1, num_cells_dm
558  CALL dbt_destroy(t_dm_virt(i_cell))
559  CALL dbt_destroy(t_dm_occ(i_cell))
560  END DO
561 
562  CALL dbt_destroy(t_3c_m_occ_tmp)
563  CALL dbt_destroy(t_3c_m_virt_tmp)
564  DEALLOCATE (t_dm_virt)
565  DEALLOCATE (t_dm_occ)
566 
567  CALL para_env%sync()
568  t2 = m_walltime()
569 
570  dbcsr_time = dbcsr_time + t2 - t1_flop
571 
572  IF (unit_nr > 0) THEN
573  WRITE (unit_nr, '(/T3,A,1X,I3)') &
574  'RPA_LOW_SCALING_INFO| Info for time point', jquad
575  WRITE (unit_nr, '(T6,A,T56,F25.1)') &
576  'Execution time (s):', t2 - t1
577  WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
578  '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  'Occupancy of D virt:', real(nze_dm_virt, dp), '/', occ_dm_virt*100, '%'
581  IF (do_gamma_rpa) THEN
582  WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
583  '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  '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  'Occupancy of M virt:', real(nze_m_virt, dp), '/', occ_m_virt*100, '%'
588  END IF
589  WRITE (unit_nr, *)
590  CALL m_flush(unit_nr)
591  END IF
592 
593  END DO ! time points
594 
595  CALL dbt_destroy(t_3c_m_occ)
596  CALL dbt_destroy(t_3c_m_virt)
597 
598  DO i = 1, SIZE(t_3c_o, 1)
599  DO j = 1, SIZE(t_3c_o, 2)
600  CALL dbt_destroy(t_3c_o_occ(i, j))
601  CALL dbt_destroy(t_3c_o_virt(i, j))
602  END DO
603  END DO
604 
605  CALL clean_up(mat_dm_occ_global, mat_dm_virt_global)
606 
607  CALL timestop(handle)
608 
609  END SUBROUTINE compute_mat_p_omega
610 
611 ! **************************************************************************************************
612 !> \brief ...
613 !> \param mat_P_omega ...
614 ! **************************************************************************************************
615  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  DO jquad = 1, SIZE(mat_p_omega, 1)
621  DO i_kp = 1, SIZE(mat_p_omega, 2)
622 
623  CALL dbcsr_set(mat_p_omega(jquad, i_kp)%matrix, 0.0_dp)
624 
625  END DO
626  END DO
627 
628  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  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  matrix_s, ispin, &
662  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  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
696  REAL(kind=dp) :: tau
697 
698  CALL timeset(routinen, handle)
699 
700  IF (memory_info .AND. unit_nr > 0) WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
701  "RPA_LOW_SCALING_INFO| Started with time point: ", jquad
702 
703  tau = tau_tj(jquad)
704 
705  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  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  remove_occ=.true., remove_virt=.false., first_jquad=1)
716 
717  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  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  alloc_dm=(jquad == 1))
728 
729  num_cells_dm = 1
730 
731  ELSE
732 
733  num_cells_dm = 1
734 
735  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  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  DO jjb = 1, nrow_local
750  DO iib = 1, ncol_local
751  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  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  fm_mo_coeff_occ%local_data(jjb, iib)*exp(tau*0.5_dp*(eigenval(i_global) - e_fermi))
758  ELSE
759  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  col_indices=col_indices)
771 
772  ! the same for virt
773  DO jjb = 1, nrow_local
774  DO iib = 1, ncol_local
775  i_global = col_indices(iib)
776 
777  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  fm_mo_coeff_virt%local_data(jjb, iib)*exp(-tau*0.5_dp*(eigenval(i_global) - e_fermi))
780  ELSE
781  fm_mo_coeff_virt_scaled%local_data(jjb, iib) = 0.0_dp
782  END IF
783 
784  END DO
785  END DO
786 
787  CALL para_env%sync()
788 
789  size_dm_occ = nmo
790  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  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  matrix_c=fm_scaled_dm_virt_tau)
799 
800  IF (jquad == 1) THEN
801 
802  ! transfer fm density matrices to dbcsr matrix
803  NULLIFY (mat_dm_occ_global)
804  CALL dbcsr_allocate_matrix_set(mat_dm_occ_global, num_integ_points, 1)
805 
806  DO iquad = 1, num_integ_points
807 
808  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  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  keep_sparsity=.false.)
820 
821  CALL dbcsr_filter(mat_dm_occ_global(jquad, 1)%matrix, eps_filter)
822 
823  IF (jquad == 1) THEN
824 
825  NULLIFY (mat_dm_virt_global)
826  CALL dbcsr_allocate_matrix_set(mat_dm_virt_global, num_integ_points, 1)
827 
828  END IF
829 
830  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  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  keep_sparsity=.false.)
837 
838  CALL dbcsr_filter(mat_dm_virt_global(jquad, 1)%matrix, eps_filter)
839 
840  ! release memory
841  IF (jquad > 1) THEN
842  CALL dbcsr_set(mat_dm_occ_global(jquad - 1, 1)%matrix, 0.0_dp)
843  CALL dbcsr_set(mat_dm_virt_global(jquad - 1, 1)%matrix, 0.0_dp)
844  CALL dbcsr_filter(mat_dm_occ_global(jquad - 1, 1)%matrix, 0.0_dp)
845  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  CALL timestop(handle)
851 
852  END SUBROUTINE compute_mat_dm_global
853 
854 ! **************************************************************************************************
855 !> \brief ...
856 !> \param mat_dm_occ_global ...
857 !> \param mat_dm_virt_global ...
858 ! **************************************************************************************************
859  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  CALL dbcsr_deallocate_matrix_set(mat_dm_occ_global)
863  CALL dbcsr_deallocate_matrix_set(mat_dm_virt_global)
864 
865  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  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  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  CALL timeset(routinen, handle)
895 
896  ! only imaginary wavefunctions supported in kpoint cubic scaling RPA
897  cpassert(kpoint%use_real_wfn .EQV. .false.)
898 
899  ! work matrix
900  mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)
901  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  cpassert(nao == nmo)
906 
907  ALLOCATE (exp_scaling(nmo))
908 
909  CALL cp_fm_get_info(mo_set%mo_coeff, matrix_struct=matrix_struct)
910  CALL cp_fm_create(fwork, matrix_struct)
911 
912  CALL get_kpoint_info(kpoint, kp_range=kp_range)
913  kplocal = kp_range(2) - kp_range(1) + 1
914 
915  DO ikpgr = 1, kplocal
916  kp => kpoint%kp_env(ikpgr)%kpoint_env
917  nspin = SIZE(kp%mos, 2)
918  DO ispin = 1, nspin
919  mo_set => kp%mos(1, ispin)
920  CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
921  rpmat => kp%wmat(1, ispin)
922  cpmat => kp%wmat(2, ispin)
923  CALL get_mo_set(mo_set, occupation_numbers=occupation)
924  CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
925 
926  IF (remove_virt) THEN
927  CALL cp_fm_column_scale(fwork, occupation)
928  END IF
929  IF (remove_occ) THEN
930  CALL cp_fm_column_scale(fwork, 2.0_dp/real(nspin, kind=dp) - occupation)
931  END IF
932 
933  ! proper spin
934  IF (nspin == 1) THEN
935  CALL cp_fm_scale(0.5_dp, fwork)
936  END IF
937 
938  DO i_mo = 1, nmo
939 
940  IF (abs(tau*0.5_dp*(eigenvalues(i_mo) - e_fermi)) < stabilize_exp) THEN
941  exp_scaling(i_mo) = exp(-abs(tau*(eigenvalues(i_mo) - e_fermi)))
942  ELSE
943  exp_scaling(i_mo) = 0.0_dp
944  END IF
945  END DO
946 
947  CALL cp_fm_column_scale(fwork, exp_scaling)
948 
949  ! Re(c)*Re(c)
950  CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, rpmat)
951  mo_set => kp%mos(2, ispin)
952  ! Im(c)*Re(c)
953  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  CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, fwork, mo_set%mo_coeff, 1.0_dp, cpmat)
956 
957  CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
958 
959  IF (remove_virt) THEN
960  CALL cp_fm_column_scale(fwork, occupation)
961  END IF
962  IF (remove_occ) THEN
963  CALL cp_fm_column_scale(fwork, 2.0_dp/real(nspin, kind=dp) - occupation)
964  END IF
965 
966  ! proper spin
967  IF (nspin == 1) THEN
968  CALL cp_fm_scale(0.5_dp, fwork)
969  END IF
970 
971  DO i_mo = 1, nmo
972  IF (abs(tau*0.5_dp*(eigenvalues(i_mo) - e_fermi)) < stabilize_exp) THEN
973  exp_scaling(i_mo) = exp(-abs(tau*(eigenvalues(i_mo) - e_fermi)))
974  ELSE
975  exp_scaling(i_mo) = 0.0_dp
976  END IF
977  END DO
978 
979  CALL cp_fm_column_scale(fwork, exp_scaling)
980  ! Im(c)*Im(c)
981  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  CALL cp_fm_release(fwork)
988  DEALLOCATE (exp_scaling)
989 
990  CALL timestop(handle)
991 
992  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  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  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_global_work, matrix_s_kp
1028  TYPE(kpoint_type), POINTER :: kpoints
1029  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1030 
1031  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  cell=cell)
1038 
1039  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  DO i_dim = 1, 3
1044  cell_grid_dm(i_dim) = (kpoints%nkp_grid(i_dim)/2)*2 - 1
1045  END DO
1046 
1047  num_cells_dm = cell_grid_dm(1)*cell_grid_dm(2)*cell_grid_dm(3)
1048 
1049  NULLIFY (mat_dm_global_work)
1050  CALL dbcsr_allocate_matrix_set(mat_dm_global_work, nspin, num_cells_dm)
1051 
1052  DO jspin = 1, nspin
1053 
1054  DO i_img = 1, num_cells_dm
1055 
1056  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  matrix_type=dbcsr_type_no_symmetry)
1061 
1062  CALL dbcsr_reserve_all_blocks(mat_dm_global_work(jspin, i_img)%matrix)
1063 
1064  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  remove_occ=remove_occ, remove_virt=remove_virt)
1073 
1074  ! overwrite the cell indices in kpoints
1075  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  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  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  IF (jquad == first_jquad) THEN
1086 
1087  NULLIFY (mat_dm_global)
1088  ALLOCATE (mat_dm_global(first_jquad:num_integ_points, num_cells_dm))
1089 
1090  DO iquad = first_jquad, num_integ_points
1091  DO i_img = 1, num_cells_dm
1092  NULLIFY (mat_dm_global(iquad, i_img)%matrix)
1093  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  matrix_type=dbcsr_type_no_symmetry)
1097 
1098  END DO
1099  END DO
1100 
1101  END IF
1102 
1103  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  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  mat_dm_global_work(ispin, i_img)%matrix)
1110 
1111  END DO
1112 
1113  CALL dbcsr_deallocate_matrix_set(mat_dm_global_work)
1114 
1115  CALL timestop(handle)
1116 
1117  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  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  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_global_work, matrix_s_kp
1144  TYPE(kpoint_type), POINTER :: kpoints_g
1145  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1146 
1147  CALL timeset(routinen, handle)
1148 
1149  NULLIFY (matrix_s_kp, mos)
1150 
1151  CALL get_qs_env(qs_env, &
1152  matrix_s_kp=matrix_s_kp, &
1153  mos=mos)
1154 
1155  kpoints_g => qs_env%mp2_env%ri_rpa_im_time%kpoints_G
1156 
1157  nspin = SIZE(mos)
1158 
1159  num_cells_dm = 1
1160 
1161  NULLIFY (mat_dm_global_work)
1162  CALL dbcsr_allocate_matrix_set(mat_dm_global_work, nspin, num_cells_dm)
1163 
1164  ! if necessaray, allocate mat_dm_global
1165  IF (alloc_dm) THEN
1166 
1167  NULLIFY (mat_dm_global)
1168  ALLOCATE (mat_dm_global(1:num_integ_points, num_cells_dm))
1169 
1170  DO iquad = 1, num_integ_points
1171  NULLIFY (mat_dm_global(iquad, 1)%matrix)
1172  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  matrix_type=dbcsr_type_no_symmetry)
1176 
1177  END DO
1178 
1179  END IF
1180 
1181  DO jspin = 1, nspin
1182 
1183  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  matrix_type=dbcsr_type_no_symmetry)
1187 
1188  CALL dbcsr_reserve_all_blocks(mat_dm_global_work(jspin, 1)%matrix)
1189 
1190  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  remove_occ=remove_occ, remove_virt=remove_virt)
1197 
1198  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  mat_dm_global_work(ispin, 1)%matrix)
1202 
1203  CALL dbcsr_deallocate_matrix_set(mat_dm_global_work)
1204 
1205  CALL timestop(handle)
1206 
1207  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  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  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1229  INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
1230  REAL(kind=dp) :: contribution, weight_im, weight_re
1231  REAL(kind=dp), DIMENSION(3, 3) :: hmat
1232  REAL(kind=dp), DIMENSION(:), POINTER :: wkp
1233  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  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1240 
1241  CALL timeset(routinen, handle)
1242 
1243  NULLIFY (xkp, wkp)
1244 
1245  CALL cp_fm_create(fm_mat_work, kpoints_g%kp_env(1)%kpoint_env%wmat(1, 1)%matrix_struct)
1246  CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
1247 
1248  CALL get_kpoint_info(kpoints_g, nkp=nkp, xkp=xkp, wkp=wkp)
1249  index_to_cell => kpoints_g%index_to_cell
1250  num_cells = SIZE(index_to_cell, 2)
1251 
1252  nspin = SIZE(mat_dm_global_work, 1)
1253 
1254  mo_set => kpoints_g%kp_env(1)%kpoint_env%mos(1, 1)
1255  CALL get_mo_set(mo_set, nao=nao)
1256 
1257  ALLOCATE (atom_from_ao_index(nao))
1258 
1259  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  col_indices=col_indices)
1266 
1267  NULLIFY (cell, particle_set)
1268  CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1269  CALL get_cell(cell=cell, h=hmat)
1270 
1271  iatom_old = 0
1272  jatom_old = 0
1273 
1274  DO ispin = 1, nspin
1275 
1276  CALL dbcsr_set(mat_dm_global_work(ispin, 1)%matrix, 0.0_dp)
1277 
1278  DO ik = 1, nkp
1279 
1280  kp => kpoints_g%kp_env(ik)%kpoint_env
1281  rpmat => kp%wmat(1, ispin)
1282  cpmat => kp%wmat(2, ispin)
1283 
1284  DO irow = 1, nrow_local
1285  DO jcol = 1, ncol_local
1286 
1287  iatom = atom_from_ao_index(row_indices(irow))
1288  jatom = atom_from_ao_index(col_indices(jcol))
1289 
1290  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  cell, index_to_cell, hmat, particle_set)
1295 
1296  iatom_old = iatom
1297  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  weight_im*cpmat%local_data(irow, jcol)
1304 
1305  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  CALL copy_fm_to_dbcsr(fm_mat_work, mat_dm_global_work(ispin, 1)%matrix, keep_sparsity=.false.)
1313  CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
1314 
1315  END DO
1316 
1317  CALL cp_fm_release(fm_mat_work)
1318  DEALLOCATE (atom_from_ao_index)
1319 
1320  CALL timestop(handle)
1321 
1322  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  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  REAL(kind=dp), DIMENSION(:), POINTER :: wkp
1342  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  CALL timeset(routinen, handle)
1348 
1349  NULLIFY (xkp, wkp)
1350 
1351  NULLIFY (mat_work_re)
1352  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  matrix_type=dbcsr_type_no_symmetry)
1356 
1357  NULLIFY (mat_work_im)
1358  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  matrix_type=dbcsr_type_no_symmetry)
1362 
1363  CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, wkp=wkp)
1364 
1365  nspin = SIZE(mat_dm_global_work, 1)
1366 
1367  cpassert(SIZE(mat_dm_global_work, 2) == SIZE(index_to_cell, 2))
1368 
1369  DO ispin = 1, nspin
1370 
1371  DO icell = 1, SIZE(mat_dm_global_work, 2)
1372 
1373  CALL dbcsr_set(mat_dm_global_work(ispin, icell)%matrix, 0.0_dp)
1374 
1375  END DO
1376 
1377  END DO
1378 
1379  DO ispin = 1, nspin
1380 
1381  DO ik = 1, nkp
1382 
1383  kp => kpoints%kp_env(ik)%kpoint_env
1384  rpmat => kp%wmat(1, ispin)
1385  cpmat => kp%wmat(2, ispin)
1386 
1387  CALL copy_fm_to_dbcsr(rpmat, mat_work_re, keep_sparsity=.false.)
1388  CALL copy_fm_to_dbcsr(cpmat, mat_work_im, keep_sparsity=.false.)
1389 
1390  DO icell = 1, SIZE(mat_dm_global_work, 2)
1391 
1392  xcell = index_to_cell(1, icell)
1393  ycell = index_to_cell(2, icell)
1394  zcell = index_to_cell(3, icell)
1395 
1396  arg = real(xcell, dp)*xkp(1, ik) + real(ycell, dp)*xkp(2, ik) + real(zcell, dp)*xkp(3, ik)
1397  coskl = wkp(ik)*cos(twopi*arg)
1398  sinkl = wkp(ik)*sin(twopi*arg)
1399 
1400  CALL dbcsr_add(mat_dm_global_work(ispin, icell)%matrix, mat_work_re, 1.0_dp, coskl)
1401  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  CALL dbcsr_release_p(mat_work_re)
1409  CALL dbcsr_release_p(mat_work_im)
1410 
1411  CALL timestop(handle)
1412 
1413  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  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  INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_unsorted
1435  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_unsorted
1436  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  CALL timeset(routinen, handle)
1441 
1442  CALL get_cell(cell=cell, h=hmat)
1443 
1444  num_cells = cell_grid(1)*cell_grid(2)*cell_grid(3)
1445  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  cpassert(cell_grid(1) .NE. itm(1)*2)
1450  cpassert(cell_grid(2) .NE. itm(2)*2)
1451  cpassert(cell_grid(3) .NE. itm(3)*2)
1452 
1453  IF (ASSOCIATED(cell_to_index)) DEALLOCATE (cell_to_index)
1454  IF (ASSOCIATED(index_to_cell)) DEALLOCATE (index_to_cell)
1455 
1456  ALLOCATE (cell_to_index_unsorted(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
1457  cell_to_index_unsorted(:, :, :) = 0
1458 
1459  ALLOCATE (index_to_cell_unsorted(3, num_cells))
1460  index_to_cell_unsorted(:, :) = 0
1461 
1462  ALLOCATE (cell_to_index(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
1463  cell_to_index(:, :, :) = 0
1464 
1465  ALLOCATE (index_to_cell(3, num_cells))
1466  index_to_cell(:, :) = 0
1467 
1468  ALLOCATE (abs_cell_vectors(1:num_cells))
1469 
1470  cell_counter = 0
1471 
1472  DO xcell = -itm(1), itm(1)
1473  DO ycell = -itm(2), itm(2)
1474  DO zcell = -itm(3), itm(3)
1475 
1476  cell_counter = cell_counter + 1
1477  cell_to_index_unsorted(xcell, ycell, zcell) = cell_counter
1478 
1479  index_to_cell_unsorted(1, cell_counter) = xcell
1480  index_to_cell_unsorted(2, cell_counter) = ycell
1481  index_to_cell_unsorted(3, cell_counter) = zcell
1482 
1483  cell_vector(1:3) = matmul(hmat, real(index_to_cell_unsorted(1:3, cell_counter), dp))
1484 
1485  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  DO i_cell = 1, num_cells/2 + 1
1494 
1495  index_min_dist = minloc(abs_cell_vectors(1:num_cells/2 + 1), dim=1)
1496 
1497  xcell = index_to_cell_unsorted(1, index_min_dist)
1498  ycell = index_to_cell_unsorted(2, index_min_dist)
1499  zcell = index_to_cell_unsorted(3, index_min_dist)
1500 
1501  index_to_cell(1, i_cell) = xcell
1502  index_to_cell(2, i_cell) = ycell
1503  index_to_cell(3, i_cell) = zcell
1504 
1505  cell_to_index(xcell, ycell, zcell) = i_cell
1506 
1507  abs_cell_vectors(index_min_dist) = 10000000000.0_dp
1508 
1509  END DO
1510 
1511  ! now all the remaining cells
1512  DO i_cell = num_cells/2 + 2, num_cells
1513 
1514  index_min_dist = minloc(abs_cell_vectors(1:num_cells), dim=1)
1515 
1516  xcell = index_to_cell_unsorted(1, index_min_dist)
1517  ycell = index_to_cell_unsorted(2, index_min_dist)
1518  zcell = index_to_cell_unsorted(3, index_min_dist)
1519 
1520  index_to_cell(1, i_cell) = xcell
1521  index_to_cell(2, i_cell) = ycell
1522  index_to_cell(3, i_cell) = zcell
1523 
1524  cell_to_index(xcell, ycell, zcell) = i_cell
1525 
1526  abs_cell_vectors(index_min_dist) = 10000000000.0_dp
1527 
1528  END DO
1529 
1530  DEALLOCATE (index_to_cell_unsorted, cell_to_index_unsorted, abs_cell_vectors)
1531 
1532  CALL timestop(handle)
1533 
1534  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  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  CALL timeset(routinen, handle)
1566 
1567  IF (do_kpoints_cubic_rpa) THEN
1568 
1569  x_cell_r = index_to_cell_3c(1, i_cell_r)
1570  y_cell_r = index_to_cell_3c(2, i_cell_r)
1571  z_cell_r = index_to_cell_3c(3, i_cell_r)
1572 
1573  x_cell_s = index_to_cell_dm(1, i_cell_s)
1574  y_cell_s = index_to_cell_dm(2, i_cell_s)
1575  z_cell_s = index_to_cell_dm(3, i_cell_s)
1576 
1577  x_cell_r_minus_s = x_cell_r - x_cell_s
1578  y_cell_r_minus_s = y_cell_r - y_cell_s
1579  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  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  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  IF (i_cell_r_minus_s == 0) THEN
1592 
1593  r_minus_s_needed = .false.
1594  i_cell_r_minus_s = 0
1595 
1596  ELSE
1597 
1598  r_minus_s_needed = .true.
1599 
1600  END IF
1601 
1602  ELSE
1603 
1604  i_cell_r_minus_s = 0
1605  r_minus_s_needed = .false.
1606 
1607  END IF
1608 
1609  ELSE ! no k-points
1610 
1611  r_minus_s_needed = .true.
1612  i_cell_r_minus_s = 1
1613 
1614  END IF
1615 
1616  CALL timestop(handle)
1617 
1618  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  SUBROUTINE get_diff_diff_index_3c(i_cell_R, i_cell_S, i_cell_T, i_cell_R_minus_S_minus_T, &
1633  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  CALL timeset(routinen, handle)
1653 
1654  IF (do_kpoints_cubic_rpa) THEN
1655 
1656  x_cell_r = index_to_cell_3c(1, i_cell_r)
1657  y_cell_r = index_to_cell_3c(2, i_cell_r)
1658  z_cell_r = index_to_cell_3c(3, i_cell_r)
1659 
1660  x_cell_s = index_to_cell_dm(1, i_cell_s)
1661  y_cell_s = index_to_cell_dm(2, i_cell_s)
1662  z_cell_s = index_to_cell_dm(3, i_cell_s)
1663 
1664  x_cell_t = index_to_cell_dm(1, i_cell_t)
1665  y_cell_t = index_to_cell_dm(2, i_cell_t)
1666  z_cell_t = index_to_cell_dm(3, i_cell_t)
1667 
1668  x_cell_r_minus_s_minus_t = x_cell_r - x_cell_s - x_cell_t
1669  y_cell_r_minus_s_minus_t = y_cell_r - y_cell_s - y_cell_t
1670  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  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  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  IF (i_cell_r_minus_s_minus_t == 0) THEN
1685 
1686  r_minus_s_minus_t_needed = .false.
1687 
1688  ELSE
1689 
1690  r_minus_s_minus_t_needed = .true.
1691 
1692  END IF
1693 
1694  ELSE
1695 
1696  i_cell_r_minus_s_minus_t = 0
1697  r_minus_s_minus_t_needed = .false.
1698 
1699  END IF
1700 
1701  ! no k-kpoints
1702  ELSE
1703 
1704  r_minus_s_minus_t_needed = .true.
1705  i_cell_r_minus_s_minus_t = 1
1706 
1707  END IF
1708 
1709  CALL timestop(handle)
1710 
1711  END SUBROUTINE get_diff_diff_index_3c
1712 
1713 END MODULE rpa_im_time
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cell_types::get_cell
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
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition: cp_dbcsr_operations.F:118
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition: cp_fm_basic_linalg.F:1882
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition: cp_fm_basic_linalg.F:1680
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
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
Definition: cp_fm_types.F:1016
cp_fm_types::cp_fm_set_all
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
Definition: cp_fm_types.F:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
dbt_api
This is the start of a dbt_api, all publically needed functions are exported here....
Definition: dbt_api.F:17
hfx_types
Types and set/get functions for HFX.
Definition: hfx_types.F:15
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition: kinds.F:54
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
kpoint_types::get_kpoint_info
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.
Definition: kpoint_types.F:333
machine
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition: machine.F:106
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition: machine.F:123
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition: mathconstants.F:112
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
mp2_types
Types needed for MP2 calculations.
Definition: mp2_types.F:14
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
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_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, 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, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
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.
Definition: qs_mo_types.F:397
qs_tensors_types
Utility methods to build 3-center integral tensors of various types.
Definition: qs_tensors_types.F:12
qs_tensors_types::create_2c_tensor
subroutine, public create_2c_tensor(t2c, dist_1, dist_2, pgrid, sizes_1, sizes_2, order, name)
...
Definition: qs_tensors_types.F:382
qs_tensors
Utility methods to build 3-center integral tensors of various types.
Definition: qs_tensors.F:11
qs_tensors::decompress_tensor
subroutine, public decompress_tensor(tensor, blk_indices, compressed, eps)
...
Definition: qs_tensors.F:3788
qs_tensors::get_tensor_occupancy
subroutine, public get_tensor_occupancy(tensor, nze, occ)
...
Definition: qs_tensors.F:3879
rpa_gw_im_time_util
Utility routines for GW with imaginary time.
Definition: rpa_gw_im_time_util.F:13
rpa_gw_im_time_util::compute_weight_re_im
subroutine, public compute_weight_re_im(weight_re, weight_im, num_cells, iatom, jatom, xkp, wkp_W, cell, index_to_cell, hmat, particle_set)
...
Definition: rpa_gw_im_time_util.F:802
rpa_gw_im_time_util::get_atom_index_from_basis_function_index
subroutine, public get_atom_index_from_basis_function_index(qs_env, atom_from_basis_index, basis_size, basis_type, first_bf_from_atom)
...
Definition: rpa_gw_im_time_util.F:747
rpa_im_time
Routines for low-scaling RPA/GW with imaginary time.
Definition: rpa_im_time.F:13
rpa_im_time::compute_mat_p_omega
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)
...
Definition: rpa_im_time.F:139
rpa_im_time::compute_periodic_dm
subroutine, public compute_periodic_dm(mat_dm_global, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, remove_occ, remove_virt, alloc_dm)
...
Definition: rpa_im_time.F:1134
rpa_im_time::compute_mat_dm_global
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)
...
Definition: rpa_im_time.F:665
rpa_im_time::init_cell_index_rpa
subroutine, public init_cell_index_rpa(cell_grid, cell_to_index, index_to_cell, cell)
...
Definition: rpa_im_time.F:1423
rpa_im_time::zero_mat_p_omega
subroutine, public zero_mat_p_omega(mat_P_omega)
...
Definition: rpa_im_time.F:616
rpa_im_time::compute_transl_dm
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)
...
Definition: rpa_im_time.F:1013