(git:73bd903)
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-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Routines for low-scaling RPA/GW with imaginary time
10!> \par History
11!> 10.2015 created [Jan Wilhelm]
12! **************************************************************************************************
14 USE cell_types, ONLY: cell_type,&
18 USE cp_cfm_types, ONLY: cp_cfm_create,&
25 USE cp_dbcsr_api, ONLY: &
28 dbcsr_release_p, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
33 USE cp_fm_types, ONLY: cp_fm_create,&
37 USE dbt_api, ONLY: &
38 dbt_batched_contract_finalize, dbt_batched_contract_init, dbt_contract, dbt_copy, &
39 dbt_copy_matrix_to_tensor, dbt_copy_tensor_to_matrix, dbt_create, dbt_destroy, dbt_filter, &
40 dbt_get_info, dbt_nblks_total, dbt_nd_mp_comm, dbt_pgrid_destroy, dbt_pgrid_type, dbt_type
41 USE hfx_types, ONLY: block_ind_type,&
43 USE kinds, ONLY: dp,&
44 int_8
45 USE kpoint_types, ONLY: get_kpoint_info,&
48 USE machine, ONLY: m_flush,&
50 USE mathconstants, ONLY: twopi
51 USE message_passing, ONLY: mp_comm_type,&
53 USE mp2_types, ONLY: mp2_type
58 USE qs_mo_types, ONLY: get_mo_set,&
60 USE qs_tensors, ONLY: decompress_tensor,&
66#include "./base/base_uses.f90"
67
68 IMPLICIT NONE
69
70 PRIVATE
71
72 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_im_time'
73
74 INTEGER, PARAMETER, PUBLIC :: propagator_sector_occupied = 1
75 INTEGER, PARAMETER, PUBLIC :: propagator_sector_virtual = 2
76 INTEGER, PARAMETER, PUBLIC :: propagator_number_of_sectors = 2
77
78 PUBLIC :: compute_mat_p_omega, &
86
87CONTAINS
88
89! **************************************************************************************************
90!> \brief Creates the sector-, time-, and cell-resolved propagator matrix set.
91!> \param propagator propagator matrix set to create
92!> \param ntime number of time or frequency points
93!> \param matrix_template DBCSR matrix template for every matrix
94!> \param index_to_cell cell coordinates for the third matrix-set index
95! **************************************************************************************************
96 SUBROUTINE create_propagator_matrix_set(propagator, ntime, matrix_template, index_to_cell)
97 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: propagator
98 INTEGER, INTENT(IN) :: ntime
99 TYPE(dbcsr_type), INTENT(IN) :: matrix_template
100 INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell
101
102 INTEGER :: icell, isector, jtime, ncell
103
104 cpassert(ntime > 0)
105 cpassert(SIZE(index_to_cell, 1) == 3)
106 ncell = SIZE(index_to_cell, 2)
107 cpassert(ncell > 0)
108
109 CALL dbcsr_allocate_matrix_set(propagator, propagator_number_of_sectors, ntime, ncell)
110
111 DO icell = 1, ncell
112 DO jtime = 1, ntime
113 DO isector = 1, propagator_number_of_sectors
114 ALLOCATE (propagator(isector, jtime, icell)%matrix)
115 CALL dbcsr_create(matrix=propagator(isector, jtime, icell)%matrix, &
116 template=matrix_template, &
117 matrix_type=dbcsr_type_no_symmetry)
118 END DO
119 END DO
120 END DO
121
122 END SUBROUTINE create_propagator_matrix_set
123
124! **************************************************************************************************
125!> \brief ...
126!> \param mat_P_omega ...
127!> \param cfm_mo_coeff ...
128!> \param homo ...
129!> \param mat_P_global ...
130!> \param matrix_s ...
131!> \param ispin ...
132!> \param t_3c_M ...
133!> \param t_3c_O ...
134!> \param t_3c_O_compressed ...
135!> \param t_3c_O_ind ...
136!> \param starts_array_mc ...
137!> \param ends_array_mc ...
138!> \param starts_array_mc_block ...
139!> \param ends_array_mc_block ...
140!> \param grid ...
141!> \param e_fermi ...
142!> \param eps_filter ...
143!> \param alpha ...
144!> \param eps_filter_im_time ...
145!> \param Eigenval ...
146!> \param nmo ...
147!> \param cut_memory ...
148!> \param unit_nr ...
149!> \param mp2_env ...
150!> \param para_env ...
151!> \param qs_env ...
152!> \param do_kpoints_from_Gamma ...
153!> \param index_to_cell_3c ...
154!> \param cell_to_index_3c ...
155!> \param has_mat_P_blocks ...
156!> \param do_ri_sos_laplace_mp2 ...
157!> \param dbcsr_time ...
158!> \param dbcsr_nflop ...
159! **************************************************************************************************
160 SUBROUTINE compute_mat_p_omega(mat_P_omega, cfm_mo_coeff, homo, &
161 mat_P_global, &
162 matrix_s, &
163 ispin, &
164 t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
165 starts_array_mc, ends_array_mc, &
166 starts_array_mc_block, ends_array_mc_block, &
167 grid, e_fermi, eps_filter, &
168 alpha, eps_filter_im_time, Eigenval, nmo, &
169 cut_memory, unit_nr, &
170 mp2_env, para_env, &
171 qs_env, do_kpoints_from_Gamma, &
172 index_to_cell_3c, cell_to_index_3c, &
173 has_mat_P_blocks, do_ri_sos_laplace_mp2, &
174 dbcsr_time, dbcsr_nflop)
175 TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN) :: mat_p_omega
176 TYPE(cp_cfm_type), INTENT(IN) :: cfm_mo_coeff
177 INTEGER, INTENT(IN) :: homo
178 TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_p_global
179 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
180 INTEGER, INTENT(IN) :: ispin
181 TYPE(dbt_type), INTENT(INOUT) :: t_3c_m
182 TYPE(dbt_type), DIMENSION(:, :), INTENT(INOUT) :: t_3c_o
183 TYPE(hfx_compression_type), DIMENSION(:, :, :), &
184 INTENT(INOUT) :: t_3c_o_compressed
185 TYPE(block_ind_type), DIMENSION(:, :, :), &
186 INTENT(INOUT) :: t_3c_o_ind
187 INTEGER, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc, &
188 starts_array_mc_block, &
189 ends_array_mc_block
190 TYPE(time_frequency_grid_type), INTENT(IN) :: grid
191 REAL(kind=dp), INTENT(IN) :: e_fermi, eps_filter, alpha, &
192 eps_filter_im_time
193 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: eigenval
194 INTEGER, INTENT(IN) :: nmo, cut_memory, unit_nr
195 TYPE(mp2_type) :: mp2_env
196 TYPE(mp_para_env_type), INTENT(IN) :: para_env
197 TYPE(qs_environment_type), POINTER :: qs_env
198 LOGICAL, INTENT(IN) :: do_kpoints_from_gamma
199 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(IN) :: index_to_cell_3c
200 INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
201 INTENT(IN) :: cell_to_index_3c
202 LOGICAL, DIMENSION(:, :, :, :, :), INTENT(INOUT) :: has_mat_p_blocks
203 LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2
204 REAL(dp), INTENT(INOUT) :: dbcsr_time
205 INTEGER(int_8), INTENT(INOUT) :: dbcsr_nflop
206
207 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_mat_P_omega'
208
209 INTEGER :: comm_2d_handle, handle, handle2, handle3, i, i_cell, i_cell_r_1, &
210 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, &
211 i_cell_t, i_mem, iquad, j, j_mem, jquad, num_3c_repl, num_cells_dm, num_integ_points, &
212 unit_nr_dbcsr
213 INTEGER(int_8) :: nze, nze_dm_occ, nze_dm_virt, nze_m_occ, &
214 nze_m_virt, nze_o
215 INTEGER(KIND=int_8) :: flops_1_occ, flops_1_virt, flops_2
216 INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_1, dist_2, mc_ranges, size_dm, &
217 size_p
218 INTEGER, DIMENSION(2) :: pdims_2d
219 INTEGER, DIMENSION(2, 1) :: ibounds_2, jbounds_2
220 INTEGER, DIMENSION(2, 2) :: ibounds_1, jbounds_1
221 INTEGER, DIMENSION(3) :: bounds_3c
222 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm
223 LOGICAL :: do_gamma_rpa, do_kpoints_cubic_rpa, first_cycle_im_time, first_cycle_omega_loop, &
224 memory_info, r_1_minus_s_needed, r_1_minus_t_needed, r_2_minus_s_minus_t_needed
225 REAL(dp) :: occ, occ_dm_occ, occ_dm_virt, occ_m_occ, &
226 occ_m_virt, occ_o, t1_flop
227 REAL(kind=dp) :: omega, omega_old, t1, t2, tau, weight, &
228 weight_old
229 TYPE(dbcsr_distribution_type) :: dist_p
230 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: propagator
231 TYPE(dbt_pgrid_type) :: pgrid_2d
232 TYPE(dbt_type) :: t_3c_m_occ, t_3c_m_occ_tmp, t_3c_m_virt, &
233 t_3c_m_virt_tmp, t_dm, t_dm_tmp, t_p, &
234 t_p_tmp
235 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_dm_occ, t_dm_virt
236 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_o_occ, t_3c_o_virt
237 TYPE(mp_comm_type) :: comm_2d
238
239 CALL timeset(routinen, handle)
240
241 IF (do_ri_sos_laplace_mp2) THEN
242 num_integ_points = SIZE(grid%imaginary_time)
243 ELSE
244 num_integ_points = SIZE(grid%frequency)
245 END IF
246
247 NULLIFY (propagator)
248
249 memory_info = mp2_env%ri_rpa_im_time%memory_info
250 IF (memory_info) THEN
251 unit_nr_dbcsr = unit_nr
252 ELSE
253 unit_nr_dbcsr = 0
254 END IF
255
256 do_kpoints_cubic_rpa = qs_env%mp2_env%ri_rpa_im_time%do_im_time_kpoints
257 do_gamma_rpa = .NOT. do_kpoints_cubic_rpa
258 num_3c_repl = maxval(cell_to_index_3c)
259
260 first_cycle_im_time = .true.
261 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)))
262 DO i = 1, SIZE(t_3c_o, 1)
263 DO j = 1, SIZE(t_3c_o, 2)
264 CALL dbt_create(t_3c_o(i, j), t_3c_o_occ(i, j))
265 CALL dbt_create(t_3c_o(i, j), t_3c_o_virt(i, j))
266 END DO
267 END DO
268
269 CALL dbt_create(t_3c_m, t_3c_m_occ, name="M occ (RI | AO AO)")
270 CALL dbt_create(t_3c_m, t_3c_m_virt, name="M virt (RI | AO AO)")
271
272 ALLOCATE (mc_ranges(cut_memory + 1))
273 mc_ranges(:cut_memory) = starts_array_mc_block(:)
274 mc_ranges(cut_memory + 1) = ends_array_mc_block(cut_memory) + 1
275
276 DO jquad = 1, num_integ_points
277
278 CALL para_env%sync()
279 t1 = m_walltime()
280
281 CALL compute_mat_dm_global(grid, nmo, cfm_mo_coeff, homo, propagator, &
282 matrix_s, ispin, &
283 eigenval, e_fermi, eps_filter, memory_info, unit_nr, &
284 jquad, do_kpoints_cubic_rpa, do_kpoints_from_gamma, qs_env, &
285 num_cells_dm, index_to_cell_dm, para_env)
286
287 ALLOCATE (t_dm_virt(num_cells_dm))
288 ALLOCATE (t_dm_occ(num_cells_dm))
289 CALL dbcsr_get_info(mat_p_global%matrix, distribution=dist_p)
290 CALL dbcsr_distribution_get(dist_p, group=comm_2d_handle, nprows=pdims_2d(1), npcols=pdims_2d(2))
291 CALL comm_2d%set_handle(comm_2d_handle)
292
293 pgrid_2d = dbt_nd_mp_comm(comm_2d, [1], [2], pdims_2d=pdims_2d)
294 ALLOCATE (size_p(dbt_nblks_total(t_3c_m, 1)))
295 CALL dbt_get_info(t_3c_m, blk_size_1=size_p)
296
297 ALLOCATE (size_dm(dbt_nblks_total(t_3c_o(1, 1), 3)))
298 CALL dbt_get_info(t_3c_o(1, 1), blk_size_3=size_dm)
299 CALL create_2c_tensor(t_dm, dist_1, dist_2, pgrid_2d, size_dm, size_dm, name="D (AO | AO)")
300 DEALLOCATE (size_dm)
301 DEALLOCATE (dist_1, dist_2)
302 CALL create_2c_tensor(t_p, dist_1, dist_2, pgrid_2d, size_p, size_p, name="P (RI | RI)")
303 DEALLOCATE (size_p)
304 DEALLOCATE (dist_1, dist_2)
305 CALL dbt_pgrid_destroy(pgrid_2d)
306
307 DO i_cell = 1, num_cells_dm
308 CALL dbt_create(t_dm, t_dm_virt(i_cell), name="D virt (AO | AO)")
309 CALL dbt_create(propagator(propagator_sector_virtual, jquad, i_cell)%matrix, t_dm_tmp)
310 CALL dbt_copy_matrix_to_tensor(propagator(propagator_sector_virtual, jquad, i_cell)%matrix, t_dm_tmp)
311 CALL dbt_copy(t_dm_tmp, t_dm_virt(i_cell), move_data=.true.)
312 CALL dbcsr_clear(propagator(propagator_sector_virtual, jquad, i_cell)%matrix)
313
314 CALL dbt_create(t_dm, t_dm_occ(i_cell), name="D occ (AO | AO)")
315 CALL dbt_copy_matrix_to_tensor(propagator(propagator_sector_occupied, jquad, i_cell)%matrix, t_dm_tmp)
316 CALL dbt_copy(t_dm_tmp, t_dm_occ(i_cell), move_data=.true.)
317 CALL dbt_destroy(t_dm_tmp)
318 CALL dbcsr_clear(propagator(propagator_sector_occupied, jquad, i_cell)%matrix)
319 END DO
320
321 CALL get_tensor_occupancy(t_dm_occ(1), nze_dm_occ, occ_dm_occ)
322 CALL get_tensor_occupancy(t_dm_virt(1), nze_dm_virt, occ_dm_virt)
323
324 CALL dbt_destroy(t_dm)
325
326 CALL dbt_create(t_3c_o_occ(1, 1), t_3c_m_occ_tmp, name="M (RI AO | AO)")
327 CALL dbt_create(t_3c_o_virt(1, 1), t_3c_m_virt_tmp, name="M (RI AO | AO)")
328
329 CALL timeset(routinen//"_contract", handle2)
330
331 CALL para_env%sync()
332 t1_flop = m_walltime()
333
334 DO i = 1, SIZE(t_3c_o_occ, 1)
335 DO j = 1, SIZE(t_3c_o_occ, 2)
336 CALL dbt_batched_contract_init(t_3c_o_occ(i, j), batch_range_2=mc_ranges, batch_range_3=mc_ranges)
337 END DO
338 END DO
339 DO i = 1, SIZE(t_3c_o_virt, 1)
340 DO j = 1, SIZE(t_3c_o_virt, 2)
341 CALL dbt_batched_contract_init(t_3c_o_virt(i, j), batch_range_2=mc_ranges, batch_range_3=mc_ranges)
342 END DO
343 END DO
344 CALL dbt_batched_contract_init(t_3c_m_occ_tmp, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
345 CALL dbt_batched_contract_init(t_3c_m_virt_tmp, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
346 CALL dbt_batched_contract_init(t_3c_m_occ, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
347 CALL dbt_batched_contract_init(t_3c_m_virt, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
348
349 DO i_cell_t = 1, num_cells_dm/2 + 1
350
351 IF (.NOT. any(has_mat_p_blocks(i_cell_t, :, :, :, :))) cycle
352
353 CALL dbt_batched_contract_init(t_p)
354
355 IF (do_gamma_rpa) THEN
356 nze_o = 0
357 nze_m_virt = 0
358 nze_m_occ = 0
359 occ_m_virt = 0.0_dp
360 occ_m_occ = 0.0_dp
361 occ_o = 0.0_dp
362 END IF
363
364 DO j_mem = 1, cut_memory
365
366 CALL dbt_get_info(t_3c_o_occ(1, 1), nfull_total=bounds_3c)
367
368 jbounds_1(:, 1) = [1, bounds_3c(1)]
369 jbounds_1(:, 2) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
370
371 jbounds_2(:, 1) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
372
373 IF (do_gamma_rpa) CALL dbt_batched_contract_init(t_dm_virt(1))
374
375 DO i_mem = 1, cut_memory
376
377 IF (.NOT. any(has_mat_p_blocks(i_cell_t, i_mem, j_mem, :, :))) cycle
378
379 ibounds_1(:, 1) = [1, bounds_3c(1)]
380 ibounds_1(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
381
382 ibounds_2(:, 1) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
383
384 IF (unit_nr_dbcsr > 0) WRITE (unit=unit_nr_dbcsr, fmt="(T3,A,I3,1X,I3)") &
385 "RPA_LOW_SCALING_INFO| Memory Cut iteration", i_mem, j_mem
386
387 DO i_cell_r_1 = 1, num_3c_repl
388
389 DO i_cell_r_2 = 1, num_3c_repl
390
391 IF (.NOT. has_mat_p_blocks(i_cell_t, i_mem, j_mem, i_cell_r_1, i_cell_r_2)) cycle
392
393 CALL get_diff_index_3c(i_cell_r_1, i_cell_t, i_cell_r_1_minus_t, &
394 index_to_cell_3c, cell_to_index_3c, index_to_cell_dm, &
395 r_1_minus_t_needed, do_kpoints_cubic_rpa)
396
397 IF (do_gamma_rpa) CALL dbt_batched_contract_init(t_dm_occ(1))
398 DO i_cell_s = 1, num_cells_dm
399 CALL get_diff_index_3c(i_cell_r_1, i_cell_s, i_cell_r_1_minus_s, index_to_cell_3c, &
400 cell_to_index_3c, index_to_cell_dm, r_1_minus_s_needed, &
401 do_kpoints_cubic_rpa)
402 IF (r_1_minus_s_needed) THEN
403
404 CALL timeset(routinen//"_calc_M_occ_t", handle3)
405 CALL decompress_tensor(t_3c_o(i_cell_r_1_minus_s, i_cell_r_2), &
406 t_3c_o_ind(i_cell_r_1_minus_s, i_cell_r_2, j_mem)%ind, &
407 t_3c_o_compressed(i_cell_r_1_minus_s, i_cell_r_2, j_mem), &
408 qs_env%mp2_env%ri_rpa_im_time%eps_compress)
409
410 IF (do_gamma_rpa .AND. i_mem == 1) THEN
411 CALL get_tensor_occupancy(t_3c_o(1, 1), nze, occ)
412 nze_o = nze_o + nze
413 occ_o = occ_o + occ
414 END IF
415
416 CALL dbt_copy(t_3c_o(i_cell_r_1_minus_s, i_cell_r_2), &
417 t_3c_o_occ(i_cell_r_1_minus_s, i_cell_r_2), move_data=.true.)
418
419 CALL dbt_contract(alpha=1.0_dp, &
420 tensor_1=t_3c_o_occ(i_cell_r_1_minus_s, i_cell_r_2), &
421 tensor_2=t_dm_occ(i_cell_s), &
422 beta=1.0_dp, &
423 tensor_3=t_3c_m_occ_tmp, &
424 contract_1=[3], notcontract_1=[1, 2], &
425 contract_2=[2], notcontract_2=[1], &
426 map_1=[1, 2], map_2=[3], &
427 bounds_2=jbounds_1, bounds_3=ibounds_2, &
428 filter_eps=eps_filter, unit_nr=unit_nr_dbcsr, &
429 flop=flops_1_occ)
430 CALL timestop(handle3)
431
432 dbcsr_nflop = dbcsr_nflop + flops_1_occ
433
434 END IF
435 END DO
436
437 IF (do_gamma_rpa) CALL dbt_batched_contract_finalize(t_dm_occ(1))
438
439 ! copy matrix to optimal contraction layout - copy is done manually in order
440 ! to better control memory allocations (we can release data of previous
441 ! representation)
442 CALL timeset(routinen//"_copy_M_occ_t", handle3)
443 CALL dbt_copy(t_3c_m_occ_tmp, t_3c_m_occ, order=[1, 3, 2], move_data=.true.)
444 CALL dbt_filter(t_3c_m_occ, eps_filter)
445 CALL timestop(handle3)
446
447 IF (do_gamma_rpa) THEN
448 CALL get_tensor_occupancy(t_3c_m_occ, nze, occ)
449 nze_m_occ = nze_m_occ + nze
450 occ_m_occ = occ_m_occ + occ
451 END IF
452
453 DO i_cell_s = 1, num_cells_dm
454 CALL get_diff_diff_index_3c(i_cell_r_2, i_cell_s, i_cell_t, i_cell_r_2_minus_s_minus_t, &
455 index_to_cell_3c, cell_to_index_3c, index_to_cell_dm, &
456 r_2_minus_s_minus_t_needed, do_kpoints_cubic_rpa)
457
458 IF (r_1_minus_t_needed .AND. r_2_minus_s_minus_t_needed) THEN
459 CALL decompress_tensor(t_3c_o(i_cell_r_2_minus_s_minus_t, i_cell_r_1_minus_t), &
460 t_3c_o_ind(i_cell_r_2_minus_s_minus_t, i_cell_r_1_minus_t, i_mem)%ind, &
461 t_3c_o_compressed(i_cell_r_2_minus_s_minus_t, i_cell_r_1_minus_t, i_mem), &
462 qs_env%mp2_env%ri_rpa_im_time%eps_compress)
463
464 CALL dbt_copy(t_3c_o(i_cell_r_2_minus_s_minus_t, i_cell_r_1_minus_t), &
465 t_3c_o_virt(i_cell_r_2_minus_s_minus_t, i_cell_r_1_minus_t), move_data=.true.)
466
467 CALL timeset(routinen//"_calc_M_virt_t", handle3)
468 CALL dbt_contract(alpha=alpha/2.0_dp, &
469 tensor_1=t_3c_o_virt( &
470 i_cell_r_2_minus_s_minus_t, i_cell_r_1_minus_t), &
471 tensor_2=t_dm_virt(i_cell_s), &
472 beta=1.0_dp, &
473 tensor_3=t_3c_m_virt_tmp, &
474 contract_1=[3], notcontract_1=[1, 2], &
475 contract_2=[2], notcontract_2=[1], &
476 map_1=[1, 2], map_2=[3], &
477 bounds_2=ibounds_1, bounds_3=jbounds_2, &
478 filter_eps=eps_filter, unit_nr=unit_nr_dbcsr, &
479 flop=flops_1_virt)
480 CALL timestop(handle3)
481
482 dbcsr_nflop = dbcsr_nflop + flops_1_virt
483
484 END IF
485 END DO
486
487 CALL timeset(routinen//"_copy_M_virt_t", handle3)
488 CALL dbt_copy(t_3c_m_virt_tmp, t_3c_m_virt, move_data=.true.)
489 CALL dbt_filter(t_3c_m_virt, eps_filter)
490 CALL timestop(handle3)
491
492 IF (do_gamma_rpa) THEN
493 CALL get_tensor_occupancy(t_3c_m_virt, nze, occ)
494 nze_m_virt = nze_m_virt + nze
495 occ_m_virt = occ_m_virt + occ
496 END IF
497
498 flops_2 = 0
499
500 CALL timeset(routinen//"_calc_P_t", handle3)
501
502 CALL dbt_contract(alpha=1.0_dp, tensor_1=t_3c_m_occ, &
503 tensor_2=t_3c_m_virt, &
504 beta=1.0_dp, &
505 tensor_3=t_p, &
506 contract_1=[2, 3], notcontract_1=[1], &
507 contract_2=[2, 3], notcontract_2=[1], &
508 map_1=[1], map_2=[2], &
509 filter_eps=eps_filter_im_time/real(cut_memory**2, kind=dp), &
510 flop=flops_2, &
511 move_data=.true., &
512 unit_nr=unit_nr_dbcsr)
513
514 CALL timestop(handle3)
515
516 first_cycle_im_time = .false.
517
518 IF (jquad == 1 .AND. flops_2 == 0) THEN
519 has_mat_p_blocks(i_cell_t, i_mem, j_mem, i_cell_r_1, i_cell_r_2) = .false.
520 END IF
521
522 END DO
523 END DO
524 END DO
525 IF (do_gamma_rpa) CALL dbt_batched_contract_finalize(t_dm_virt(1))
526 END DO
527
528 CALL dbt_batched_contract_finalize(t_p, unit_nr=unit_nr_dbcsr)
529
530 CALL dbt_create(mat_p_global%matrix, t_p_tmp)
531 CALL dbt_copy(t_p, t_p_tmp, move_data=.true.)
532 CALL dbt_copy_tensor_to_matrix(t_p_tmp, mat_p_global%matrix)
533 CALL dbt_destroy(t_p_tmp)
534
535 IF (do_ri_sos_laplace_mp2) THEN
536 ! For RI-SOS-Laplace-MP2 we do not perform a cosine transform,
537 ! but we have to copy P_local to the output matrix
538
539 CALL dbcsr_add(mat_p_omega(jquad, i_cell_t)%matrix, mat_p_global%matrix, 1.0_dp, 1.0_dp)
540 ELSE
541 CALL timeset(routinen//"_Fourier_transform", handle3)
542
543 ! Fourier transform of P(it) to P(iw)
544 first_cycle_omega_loop = .true.
545
546 tau = grid%imaginary_time(jquad)
547
548 DO iquad = 1, num_integ_points
549
550 omega = grid%frequency(iquad)
551 weight = grid%cosine_time_to_frequency_weights(iquad, jquad)
552
553 IF (first_cycle_omega_loop) THEN
554 ! no multiplication with 2.0 as in Kresses paper (Kaltak, JCTC 10, 2498 (2014), Eq. 12)
555 ! because this factor is already absorbed in the weight w_j
556 CALL dbcsr_scale(mat_p_global%matrix, cos(omega*tau)*weight)
557 ELSE
558 CALL dbcsr_scale(mat_p_global%matrix, cos(omega*tau)/cos(omega_old*tau)*weight/weight_old)
559 END IF
560
561 CALL dbcsr_add(mat_p_omega(iquad, i_cell_t)%matrix, mat_p_global%matrix, 1.0_dp, 1.0_dp)
562
563 first_cycle_omega_loop = .false.
564
565 omega_old = omega
566 weight_old = weight
567
568 END DO
569
570 CALL timestop(handle3)
571 END IF
572
573 END DO
574
575 CALL timestop(handle2)
576
577 CALL dbt_batched_contract_finalize(t_3c_m_occ_tmp)
578 CALL dbt_batched_contract_finalize(t_3c_m_virt_tmp)
579 CALL dbt_batched_contract_finalize(t_3c_m_occ)
580 CALL dbt_batched_contract_finalize(t_3c_m_virt)
581
582 DO i = 1, SIZE(t_3c_o_occ, 1)
583 DO j = 1, SIZE(t_3c_o_occ, 2)
584 CALL dbt_batched_contract_finalize(t_3c_o_occ(i, j))
585 END DO
586 END DO
587
588 DO i = 1, SIZE(t_3c_o_virt, 1)
589 DO j = 1, SIZE(t_3c_o_virt, 2)
590 CALL dbt_batched_contract_finalize(t_3c_o_virt(i, j))
591 END DO
592 END DO
593
594 CALL dbt_destroy(t_p)
595 DO i_cell = 1, num_cells_dm
596 CALL dbt_destroy(t_dm_virt(i_cell))
597 CALL dbt_destroy(t_dm_occ(i_cell))
598 END DO
599
600 CALL dbt_destroy(t_3c_m_occ_tmp)
601 CALL dbt_destroy(t_3c_m_virt_tmp)
602 DEALLOCATE (t_dm_virt)
603 DEALLOCATE (t_dm_occ)
604
605 CALL para_env%sync()
606 t2 = m_walltime()
607
608 dbcsr_time = dbcsr_time + t2 - t1_flop
609
610 IF (unit_nr > 0) THEN
611 WRITE (unit_nr, '(/T3,A,1X,I3)') &
612 'RPA_LOW_SCALING_INFO| Info for time point', jquad
613 WRITE (unit_nr, '(T6,A,T56,F25.1)') &
614 'Execution time (s):', t2 - t1
615 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
616 'Occupancy of D occ:', real(nze_dm_occ, dp), '/', occ_dm_occ*100, '%'
617 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
618 'Occupancy of D virt:', real(nze_dm_virt, dp), '/', occ_dm_virt*100, '%'
619 IF (do_gamma_rpa) THEN
620 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
621 'Occupancy of 3c ints:', real(nze_o, dp), '/', occ_o*100, '%'
622 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
623 'Occupancy of M occ:', real(nze_m_occ, dp), '/', occ_m_occ*100, '%'
624 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
625 'Occupancy of M virt:', real(nze_m_virt, dp), '/', occ_m_virt*100, '%'
626 END IF
627 WRITE (unit_nr, *)
628 CALL m_flush(unit_nr)
629 END IF
630
631 END DO ! time points
632
633 CALL dbt_destroy(t_3c_m_occ)
634 CALL dbt_destroy(t_3c_m_virt)
635
636 DO i = 1, SIZE(t_3c_o, 1)
637 DO j = 1, SIZE(t_3c_o, 2)
638 CALL dbt_destroy(t_3c_o_occ(i, j))
639 CALL dbt_destroy(t_3c_o_virt(i, j))
640 END DO
641 END DO
642
643 CALL dbcsr_deallocate_matrix_set(propagator)
644
645 CALL timestop(handle)
646
647 END SUBROUTINE compute_mat_p_omega
648
649! **************************************************************************************************
650!> \brief ...
651!> \param mat_P_omega ...
652! **************************************************************************************************
653 SUBROUTINE zero_mat_p_omega(mat_P_omega)
654 TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN) :: mat_p_omega
655
656 INTEGER :: i_kp, jquad
657
658 DO jquad = 1, SIZE(mat_p_omega, 1)
659 DO i_kp = 1, SIZE(mat_p_omega, 2)
660
661 CALL dbcsr_set(mat_p_omega(jquad, i_kp)%matrix, 0.0_dp)
662
663 END DO
664 END DO
665
666 END SUBROUTINE zero_mat_p_omega
667
668! **************************************************************************************************
669!> \brief ...
670!> \param grid ...
671!> \param nmo ...
672!> \param cfm_mo_coeff ...
673!> \param homo ...
674!> \param propagator ...
675!> \param matrix_s ...
676!> \param ispin ...
677!> \param Eigenval ...
678!> \param e_fermi ...
679!> \param eps_filter ...
680!> \param memory_info ...
681!> \param unit_nr ...
682!> \param jquad ...
683!> \param do_kpoints_cubic_RPA ...
684!> \param do_kpoints_from_Gamma ...
685!> \param qs_env ...
686!> \param num_cells_dm ...
687!> \param index_to_cell_dm ...
688!> \param para_env ...
689! **************************************************************************************************
690 SUBROUTINE compute_mat_dm_global(grid, nmo, cfm_mo_coeff, homo, propagator, &
691 matrix_s, ispin, &
692 Eigenval, e_fermi, eps_filter, memory_info, unit_nr, &
693 jquad, do_kpoints_cubic_RPA, do_kpoints_from_Gamma, qs_env, &
694 num_cells_dm, index_to_cell_dm, para_env)
695
696 TYPE(time_frequency_grid_type), INTENT(IN) :: grid
697 INTEGER, INTENT(IN) :: nmo
698 TYPE(cp_cfm_type), INTENT(IN) :: cfm_mo_coeff
699 INTEGER, INTENT(IN) :: homo
700 TYPE(dbcsr_p_type), DIMENSION(:, :, :), &
701 INTENT(INOUT), POINTER :: propagator
702 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: matrix_s
703 INTEGER, INTENT(IN) :: ispin
704 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: eigenval
705 REAL(kind=dp), INTENT(IN) :: e_fermi, eps_filter
706 LOGICAL, INTENT(IN) :: memory_info
707 INTEGER, INTENT(IN) :: unit_nr, jquad
708 LOGICAL, INTENT(IN) :: do_kpoints_cubic_rpa, &
709 do_kpoints_from_gamma
710 TYPE(qs_environment_type), POINTER :: qs_env
711 INTEGER, INTENT(OUT) :: num_cells_dm
712 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm
713 TYPE(mp_para_env_type), INTENT(IN) :: para_env
714
715 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_mat_dm_global'
716
717 INTEGER :: handle, num_integ_points
718 INTEGER, DIMENSION(3, 1) :: index_to_cell_zero
719 REAL(kind=dp) :: tau
720
721 index_to_cell_zero = 0
722 num_integ_points = SIZE(grid%imaginary_time)
723
724 CALL timeset(routinen, handle)
725
726 IF (memory_info .AND. unit_nr > 0) WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
727 "RPA_LOW_SCALING_INFO| Started with time point: ", jquad
728
729 tau = grid%imaginary_time(jquad)
730
731 IF (do_kpoints_cubic_rpa) THEN
732
733 CALL compute_transl_dm(propagator, qs_env, &
734 ispin, num_integ_points, jquad, e_fermi, tau, &
735 eps_filter, num_cells_dm, index_to_cell_dm, &
737
738 CALL compute_transl_dm(propagator, qs_env, &
739 ispin, num_integ_points, jquad, e_fermi, tau, &
740 eps_filter, num_cells_dm, index_to_cell_dm, &
742
743 ELSE IF (do_kpoints_from_gamma) THEN
744
745 CALL compute_periodic_dm(propagator, qs_env, &
746 ispin, num_integ_points, jquad, e_fermi, tau, &
748
749 CALL compute_periodic_dm(propagator, qs_env, &
750 ispin, num_integ_points, jquad, e_fermi, tau, &
752
753 num_cells_dm = 1
754
755 ELSE
756
757 num_cells_dm = 1
758
759 IF (jquad == 1) THEN
760 CALL create_propagator_matrix_set(propagator, num_integ_points, matrix_s(1)%matrix, &
761 index_to_cell_zero)
762 END IF
763
764 CALL compute_gamma_propagator(propagator, jquad, cfm_mo_coeff, homo, eigenval, nmo, &
765 eps_filter, e_fermi, tau, para_env)
766
767 ! release memory
768 IF (jquad > 1) THEN
769 CALL dbcsr_set(propagator(propagator_sector_occupied, jquad - 1, 1)%matrix, 0.0_dp)
770 CALL dbcsr_set(propagator(propagator_sector_virtual, jquad - 1, 1)%matrix, 0.0_dp)
771 CALL dbcsr_filter(propagator(propagator_sector_occupied, jquad - 1, 1)%matrix, 0.0_dp)
772 CALL dbcsr_filter(propagator(propagator_sector_virtual, jquad - 1, 1)%matrix, 0.0_dp)
773 END IF
774
775 END IF ! do kpoints
776
777 CALL timestop(handle)
778
779 END SUBROUTINE compute_mat_dm_global
780
781! **************************************************************************************************
782!> \brief Builds the Gamma-point occupied and virtual propagators.
783!> \param propagator sector/time/cell propagator matrix set
784!> \param jquad time-point index
785!> \param cfm_mo_coeff complete MO coefficient matrix
786!> \param homo number of occupied orbitals
787!> \param Eigenval orbital energies
788!> \param nmo number of orbitals
789!> \param eps_filter DBCSR filtering threshold
790!> \param e_fermi Fermi energy
791!> \param tau imaginary time
792!> \param para_env parallel environment
793! **************************************************************************************************
794 SUBROUTINE compute_gamma_propagator(propagator, jquad, cfm_mo_coeff, homo, Eigenval, nmo, &
795 eps_filter, e_fermi, tau, para_env)
796 TYPE(dbcsr_p_type), DIMENSION(:, :, :), &
797 INTENT(INOUT), POINTER :: propagator
798 INTEGER, INTENT(IN) :: jquad
799 TYPE(cp_cfm_type), INTENT(IN) :: cfm_mo_coeff
800 INTEGER, INTENT(IN) :: homo
801 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: eigenval
802 INTEGER, INTENT(IN) :: nmo
803 REAL(kind=dp), INTENT(IN) :: eps_filter, e_fermi, tau
804 TYPE(mp_para_env_type), INTENT(IN) :: para_env
805
806 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_gamma_propagator'
807 REAL(kind=dp), PARAMETER :: stabilize_exp = 70.0_dp
808
809 INTEGER :: handle, i_global, iib, jjb, nao, &
810 ncol_local, nrow_local
811 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
812 TYPE(cp_cfm_type) :: cfm_mo_coeff_occ_scaled, &
813 cfm_mo_coeff_virt_scaled, &
814 cfm_scaled_dm_occ_tau, &
815 cfm_scaled_dm_virt_tau
816 TYPE(cp_fm_type) :: fm_scaled_dm_occ_tau, &
817 fm_scaled_dm_virt_tau
818 TYPE(dbcsr_type), POINTER :: mat_occ, mat_virt
819
820 CALL timeset(routinen, handle)
821 CALL para_env%sync()
822
823 mat_occ => propagator(propagator_sector_occupied, jquad, 1)%matrix
824 mat_virt => propagator(propagator_sector_virtual, jquad, 1)%matrix
825
826 CALL cp_cfm_get_info(matrix=cfm_mo_coeff, &
827 nrow_global=nao, &
828 nrow_local=nrow_local, &
829 ncol_local=ncol_local, &
830 row_indices=row_indices, &
831 col_indices=col_indices)
832
833 CALL cp_cfm_create(cfm_scaled_dm_occ_tau, cfm_mo_coeff%matrix_struct, nrow=nao, ncol=nao)
834 CALL cp_cfm_create(cfm_scaled_dm_virt_tau, cfm_mo_coeff%matrix_struct, nrow=nao, ncol=nao)
835 CALL cp_cfm_create(cfm_mo_coeff_occ_scaled, cfm_mo_coeff%matrix_struct)
836 CALL cp_cfm_create(cfm_mo_coeff_virt_scaled, cfm_mo_coeff%matrix_struct)
837
838 DO jjb = 1, nrow_local
839 DO iib = 1, ncol_local
840 i_global = col_indices(iib)
841 IF (i_global <= homo .AND. abs(tau*0.5_dp*(eigenval(i_global) - e_fermi)) < stabilize_exp) THEN
842 cfm_mo_coeff_occ_scaled%local_data(jjb, iib) = &
843 cfm_mo_coeff%local_data(jjb, iib)*exp(tau*0.5_dp*(eigenval(i_global) - e_fermi))
844 ELSE
845 cfm_mo_coeff_occ_scaled%local_data(jjb, iib) = (0.0_dp, 0.0_dp)
846 END IF
847 END DO
848 END DO
849
850 DO jjb = 1, nrow_local
851 DO iib = 1, ncol_local
852 i_global = col_indices(iib)
853 IF (i_global > homo .AND. abs(tau*0.5_dp*(eigenval(i_global) - e_fermi)) < stabilize_exp) THEN
854 cfm_mo_coeff_virt_scaled%local_data(jjb, iib) = &
855 cfm_mo_coeff%local_data(jjb, iib)*exp(-tau*0.5_dp*(eigenval(i_global) - e_fermi))
856 ELSE
857 cfm_mo_coeff_virt_scaled%local_data(jjb, iib) = (0.0_dp, 0.0_dp)
858 END IF
859 END DO
860 END DO
861
862 CALL para_env%sync()
863
864 CALL parallel_gemm(transa="N", transb="C", m=nao, n=nao, k=nmo, alpha=(1.0_dp, 0.0_dp), &
865 matrix_a=cfm_mo_coeff_occ_scaled, matrix_b=cfm_mo_coeff_occ_scaled, beta=(0.0_dp, 0.0_dp), &
866 matrix_c=cfm_scaled_dm_occ_tau)
867 CALL parallel_gemm(transa="N", transb="C", m=nao, n=nao, k=nmo, alpha=(1.0_dp, 0.0_dp), &
868 matrix_a=cfm_mo_coeff_virt_scaled, matrix_b=cfm_mo_coeff_virt_scaled, beta=(0.0_dp, 0.0_dp), &
869 matrix_c=cfm_scaled_dm_virt_tau)
870
871 CALL cp_fm_create(fm_scaled_dm_occ_tau, cfm_scaled_dm_occ_tau%matrix_struct)
872 CALL cp_fm_create(fm_scaled_dm_virt_tau, cfm_scaled_dm_virt_tau%matrix_struct)
873 CALL cp_cfm_to_fm(cfm_scaled_dm_occ_tau, fm_scaled_dm_occ_tau)
874 CALL cp_cfm_to_fm(cfm_scaled_dm_virt_tau, fm_scaled_dm_virt_tau)
875
876 CALL dbcsr_set(mat_occ, 0.0_dp)
877 CALL copy_fm_to_dbcsr(fm_scaled_dm_occ_tau, mat_occ, keep_sparsity=.false.)
878 CALL dbcsr_filter(mat_occ, eps_filter)
879 CALL dbcsr_set(mat_virt, 0.0_dp)
880 CALL copy_fm_to_dbcsr(fm_scaled_dm_virt_tau, mat_virt, keep_sparsity=.false.)
881 CALL dbcsr_filter(mat_virt, eps_filter)
882
883 CALL cp_fm_release(fm_scaled_dm_occ_tau)
884 CALL cp_fm_release(fm_scaled_dm_virt_tau)
885 CALL cp_cfm_release(cfm_scaled_dm_occ_tau)
886 CALL cp_cfm_release(cfm_scaled_dm_virt_tau)
887 CALL cp_cfm_release(cfm_mo_coeff_occ_scaled)
888 CALL cp_cfm_release(cfm_mo_coeff_virt_scaled)
889 CALL timestop(handle)
890
891 END SUBROUTINE compute_gamma_propagator
892
893! **************************************************************************************************
894!> \brief Calculate one kpoint density matrix in the complex AO representation.
895!> \param kpoint kpoint environment
896!> \param ikpgr local kpoint index
897!> \param ispin spin index
898!> \param tau ...
899!> \param e_fermi ...
900!> \param sector ...
901!> \param density complex AO density matrix
902! **************************************************************************************************
903 SUBROUTINE build_kpoint_density_matrix_rpa(kpoint, ikpgr, ispin, tau, e_fermi, sector, density)
904
905 TYPE(kpoint_type), POINTER :: kpoint
906 INTEGER, INTENT(IN) :: ikpgr, ispin
907 REAL(kind=dp), INTENT(IN) :: tau, e_fermi
908 INTEGER, INTENT(IN) :: sector
909 TYPE(cp_cfm_type), INTENT(INOUT) :: density
910
911 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_kpoint_density_matrix_rpa'
912 REAL(kind=dp), PARAMETER :: stabilize_exp = 70.0_dp
913
914 INTEGER :: handle, i_mo, nao, nmo, nspin
915 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues, exp_scaling, occupation
916 TYPE(cp_cfm_type) :: cfm_mo_coeff, cfm_work
917 TYPE(kpoint_env_type), POINTER :: kp
918 TYPE(mo_set_type), POINTER :: mo_set, mo_set_im
919
920 CALL timeset(routinen, handle)
921
922 cpassert(sector == propagator_sector_occupied .OR. sector == propagator_sector_virtual)
923 cpassert(kpoint%use_real_wfn .EQV. .false.)
924
925 kp => kpoint%kp_env(ikpgr)%kpoint_env
926 nspin = SIZE(kp%mos, 2)
927 cpassert(ispin >= 1 .AND. ispin <= nspin)
928 mo_set => kp%mos(1, ispin)
929 mo_set_im => kp%mos(2, ispin)
930 CALL get_mo_set(mo_set, nao=nao, nmo=nmo, eigenvalues=eigenvalues, occupation_numbers=occupation)
931
932 ! if this CPASSERT is triggered, please add all virtual MOs to SCF section,
933 ! e.g. ADDED_MOS 1000000
934 cpassert(nao == nmo)
935
936 ALLOCATE (exp_scaling(nmo))
937 CALL cp_cfm_create(cfm_mo_coeff, mo_set%mo_coeff%matrix_struct)
938 CALL cp_cfm_create(cfm_work, mo_set%mo_coeff%matrix_struct)
939 CALL cp_fm_to_cfm(msourcer=mo_set%mo_coeff, msourcei=mo_set_im%mo_coeff, &
940 mtarget=cfm_mo_coeff)
941 CALL cp_cfm_to_cfm(cfm_mo_coeff, cfm_work)
942
943 IF (sector == propagator_sector_occupied) THEN
944 CALL cp_cfm_column_scale(cfm_work, cmplx(occupation, 0.0_dp, kind=dp))
945 ELSE
946 CALL cp_cfm_column_scale(cfm_work, &
947 cmplx(2.0_dp/real(nspin, kind=dp) - occupation, 0.0_dp, kind=dp))
948 END IF
949
950 IF (nspin == 1) THEN
951 CALL cp_cfm_scale(0.5_dp, cfm_work)
952 END IF
953
954 DO i_mo = 1, nmo
955 IF (abs(tau*0.5_dp*(eigenvalues(i_mo) - e_fermi)) < stabilize_exp) THEN
956 exp_scaling(i_mo) = exp(-abs(tau*(eigenvalues(i_mo) - e_fermi)))
957 ELSE
958 exp_scaling(i_mo) = 0.0_dp
959 END IF
960 END DO
961
962 CALL cp_cfm_column_scale(cfm_work, cmplx(exp_scaling, 0.0_dp, kind=dp))
963 CALL parallel_gemm("N", "C", nao, nao, nmo, (1.0_dp, 0.0_dp), cfm_mo_coeff, cfm_work, &
964 (0.0_dp, 0.0_dp), density)
965
966 CALL cp_cfm_release(cfm_work)
967 CALL cp_cfm_release(cfm_mo_coeff)
968 DEALLOCATE (exp_scaling)
969
970 CALL timestop(handle)
971
972 END SUBROUTINE build_kpoint_density_matrix_rpa
973
974! **************************************************************************************************
975!> \brief ...
976!> \param propagator ...
977!> \param qs_env ...
978!> \param ispin ...
979!> \param num_integ_points ...
980!> \param jquad ...
981!> \param e_fermi ...
982!> \param tau ...
983!> \param eps_filter ...
984!> \param num_cells_dm ...
985!> \param index_to_cell_dm ...
986!> \param sector ...
987! **************************************************************************************************
988 SUBROUTINE compute_transl_dm(propagator, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, &
989 eps_filter, num_cells_dm, index_to_cell_dm, &
990 sector)
991 TYPE(dbcsr_p_type), DIMENSION(:, :, :), &
992 INTENT(INOUT), POINTER :: propagator
993 TYPE(qs_environment_type), POINTER :: qs_env
994 INTEGER, INTENT(IN) :: ispin, num_integ_points, jquad
995 REAL(kind=dp), INTENT(IN) :: e_fermi, tau, eps_filter
996 INTEGER, INTENT(OUT) :: num_cells_dm
997 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm
998 INTEGER, INTENT(IN) :: sector
999
1000 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_transl_dm'
1001
1002 INTEGER :: handle, i_dim, i_img, nao
1003 INTEGER, DIMENSION(3) :: cell_grid_dm
1004 TYPE(cell_type), POINTER :: cell
1005 TYPE(cp_cfm_type) :: cfm_density
1006 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_dm_global_work
1007 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp
1008 TYPE(kpoint_type), POINTER :: kpoints
1009
1010 CALL timeset(routinen, handle)
1011
1012 CALL get_qs_env(qs_env, &
1013 matrix_s_kp=matrix_s_kp, &
1014 kpoints=kpoints, &
1015 cell=cell)
1016
1017 ! we always use an odd number of image cells
1018 ! CAUTION: also at another point, cell_grid_dm is defined, these definitions have to be identical
1019 DO i_dim = 1, 3
1020 cell_grid_dm(i_dim) = (kpoints%nkp_grid(i_dim)/2)*2 - 1
1021 END DO
1022
1023 num_cells_dm = cell_grid_dm(1)*cell_grid_dm(2)*cell_grid_dm(3)
1024
1025 NULLIFY (mat_dm_global_work)
1026 ! Only the requested spin is built.
1027 CALL dbcsr_allocate_matrix_set(mat_dm_global_work, num_cells_dm)
1028
1029 DO i_img = 1, num_cells_dm
1030
1031 ALLOCATE (mat_dm_global_work(i_img)%matrix)
1032 CALL dbcsr_create(matrix=mat_dm_global_work(i_img)%matrix, &
1033 template=matrix_s_kp(1, 1)%matrix, &
1034 ! matrix_type=dbcsr_type_symmetric)
1035 matrix_type=dbcsr_type_no_symmetry)
1036
1037 CALL dbcsr_reserve_all_blocks(mat_dm_global_work(i_img)%matrix)
1038
1039 END DO
1040
1041 ! Reusable complex AO scratch for the k-point density source.
1042 CALL get_mo_set(kpoints%kp_env(1)%kpoint_env%mos(1, 1), nao=nao)
1043 CALL cp_cfm_create(cfm_density, kpoints%kp_env(1)%kpoint_env%mos(1, 1)%mo_coeff%matrix_struct, &
1044 nrow=nao, ncol=nao)
1045
1046 ! overwrite the cell indices in kpoints
1047 CALL init_cell_index_rpa(cell_grid_dm, kpoints%cell_to_index, kpoints%index_to_cell, cell)
1048
1049 ! density matrices in real space, the cell vectors T for transforming are taken from kpoints%index_to_cell
1050 ! (custom made for RPA) and not from sab_nl (which is symmetric and from SCF)
1051 CALL density_matrix_from_kp_to_transl(kpoints, cfm_density, mat_dm_global_work, kpoints%index_to_cell, &
1052 ispin, tau, e_fermi, sector)
1053
1054 CALL cp_cfm_release(cfm_density)
1055
1056 ! we need the index to cell for the density matrices later
1057 index_to_cell_dm => kpoints%index_to_cell
1058
1059 ! The first request owns allocation of the complete time-indexed propagator set.
1060 IF (.NOT. ASSOCIATED(propagator)) THEN
1061 CALL create_propagator_matrix_set(propagator, num_integ_points, matrix_s_kp(1, 1)%matrix, &
1062 kpoints%index_to_cell)
1063 END IF
1064
1065 DO i_img = 1, num_cells_dm
1066
1067 ! filter to get rid of the blocks full with zeros on the lower half, otherwise blocks doubled
1068 CALL dbcsr_filter(mat_dm_global_work(i_img)%matrix, eps_filter)
1069
1070 CALL dbcsr_copy(propagator(sector, jquad, i_img)%matrix, &
1071 mat_dm_global_work(i_img)%matrix)
1072
1073 END DO
1074
1075 CALL dbcsr_deallocate_matrix_set(mat_dm_global_work)
1076
1077 CALL timestop(handle)
1078
1079 END SUBROUTINE compute_transl_dm
1080
1081! **************************************************************************************************
1082!> \brief ...
1083!> \param propagator ...
1084!> \param qs_env ...
1085!> \param ispin ...
1086!> \param num_integ_points ...
1087!> \param jquad ...
1088!> \param e_fermi ...
1089!> \param tau ...
1090!> \param sector ...
1091! **************************************************************************************************
1092 SUBROUTINE compute_periodic_dm(propagator, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, &
1093 sector)
1094 TYPE(dbcsr_p_type), DIMENSION(:, :, :), &
1095 INTENT(INOUT), POINTER :: propagator
1096 TYPE(qs_environment_type), POINTER :: qs_env
1097 INTEGER, INTENT(IN) :: ispin, num_integ_points, jquad
1098 REAL(kind=dp), INTENT(IN) :: e_fermi, tau
1099 INTEGER, INTENT(IN) :: sector
1100
1101 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_periodic_dm'
1102
1103 INTEGER :: handle, nao, num_cells_dm
1104 INTEGER, DIMENSION(3, 1) :: index_to_cell_zero
1105 TYPE(cp_cfm_type) :: cfm_density
1106 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_dm_global_work
1107 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp
1108 TYPE(kpoint_type), POINTER :: kpoints_g
1109
1110 CALL timeset(routinen, handle)
1111
1112 NULLIFY (matrix_s_kp)
1113
1114 CALL get_qs_env(qs_env, &
1115 matrix_s_kp=matrix_s_kp)
1116
1117 kpoints_g => qs_env%mp2_env%ri_rpa_im_time%kpoints_G
1118
1119 num_cells_dm = 1
1120 index_to_cell_zero = 0
1121
1122 NULLIFY (mat_dm_global_work)
1123 ! Only the requested spin is built.
1124 CALL dbcsr_allocate_matrix_set(mat_dm_global_work, num_cells_dm)
1125
1126 ! The first request owns allocation of the complete time-indexed propagator set.
1127 IF (.NOT. ASSOCIATED(propagator)) THEN
1128 CALL create_propagator_matrix_set(propagator, num_integ_points, matrix_s_kp(1, 1)%matrix, &
1129 index_to_cell_zero)
1130 END IF
1131
1132 ALLOCATE (mat_dm_global_work(1)%matrix)
1133 CALL dbcsr_create(matrix=mat_dm_global_work(1)%matrix, &
1134 template=matrix_s_kp(1, 1)%matrix, &
1135 matrix_type=dbcsr_type_no_symmetry)
1136
1137 CALL dbcsr_reserve_all_blocks(mat_dm_global_work(1)%matrix)
1138
1139 CALL dbcsr_set(mat_dm_global_work(1)%matrix, 0.0_dp)
1140
1141 ! Reusable complex AO scratch for the k-point density source.
1142 CALL get_mo_set(kpoints_g%kp_env(1)%kpoint_env%mos(1, 1), nao=nao)
1143 CALL cp_cfm_create(cfm_density, kpoints_g%kp_env(1)%kpoint_env%mos(1, 1)%mo_coeff%matrix_struct, &
1144 nrow=nao, ncol=nao)
1145
1146 CALL density_matrix_from_kp_to_mic(kpoints_g, cfm_density, mat_dm_global_work, qs_env, ispin, tau, e_fermi, sector)
1147
1148 CALL cp_cfm_release(cfm_density)
1149
1150 CALL dbcsr_copy(propagator(sector, jquad, 1)%matrix, &
1151 mat_dm_global_work(1)%matrix)
1152
1153 CALL dbcsr_deallocate_matrix_set(mat_dm_global_work)
1154
1155 CALL timestop(handle)
1156
1157 END SUBROUTINE compute_periodic_dm
1158
1159 ! **************************************************************************************************
1160!> \brief ...
1161!> \param kpoints_G ...
1162!> \param density complex AO density scratch matrix
1163!> \param mat_dm_global_work ...
1164!> \param qs_env ...
1165!> \param ispin ...
1166!> \param tau ...
1167!> \param e_fermi ...
1168!> \param sector ...
1169! **************************************************************************************************
1170 SUBROUTINE density_matrix_from_kp_to_mic(kpoints_G, density, mat_dm_global_work, qs_env, ispin, tau, e_fermi, sector)
1171
1172 TYPE(kpoint_type), POINTER :: kpoints_g
1173 TYPE(cp_cfm_type), INTENT(INOUT) :: density
1174 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_dm_global_work
1175 TYPE(qs_environment_type), POINTER :: qs_env
1176 INTEGER, INTENT(IN) :: ispin
1177 REAL(kind=dp), INTENT(IN) :: tau, e_fermi
1178 INTEGER, INTENT(IN) :: sector
1179
1180 CHARACTER(LEN=*), PARAMETER :: routinen = 'density_matrix_from_kp_to_mic'
1181
1182 COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1183 POINTER :: cfm_data
1184 INTEGER :: handle, iatom, iatom_old, ik, ik_global, &
1185 ik_old, irow, jatom, jatom_old, jcol, &
1186 nao, ncol_local, nrow_local, num_cells
1187 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_ao_index
1188 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1189 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
1190 REAL(kind=dp) :: contribution, weight_im, weight_re
1191 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1192 REAL(kind=dp), DIMENSION(:), POINTER :: wkp
1193 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
1194 TYPE(cell_type), POINTER :: cell
1195 TYPE(cp_fm_type) :: fm_mat_work
1196 TYPE(kpoint_env_type), POINTER :: kp
1197 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1198
1199 CALL timeset(routinen, handle)
1200
1201 NULLIFY (xkp, wkp)
1202
1203 CALL get_kpoint_info(kpoints_g, xkp=xkp, wkp=wkp)
1204 index_to_cell => kpoints_g%index_to_cell
1205 num_cells = SIZE(index_to_cell, 2)
1206
1207 CALL cp_cfm_get_info(density, nrow_global=nao, nrow_local=nrow_local, ncol_local=ncol_local, &
1208 row_indices=row_indices, col_indices=col_indices)
1209 CALL cp_fm_create(fm_mat_work, density%matrix_struct, set_zero=.true.)
1210
1211 ALLOCATE (atom_from_ao_index(nao))
1212
1213 CALL get_atom_index_from_basis_function_index(qs_env, atom_from_ao_index, nao, "ORB")
1214
1215 NULLIFY (cell, particle_set)
1216 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1217 CALL get_cell(cell=cell, h=hmat)
1218
1219 iatom_old = 0
1220 jatom_old = 0
1221 ik_old = 0
1222
1223 CALL dbcsr_set(mat_dm_global_work(1)%matrix, 0.0_dp)
1224
1225 DO ik = 1, SIZE(kpoints_g%kp_env)
1226
1227 kp => kpoints_g%kp_env(ik)%kpoint_env
1228 ik_global = kp%nkpoint
1229 CALL build_kpoint_density_matrix_rpa(kpoints_g, ik, ispin, tau, e_fermi, sector, density)
1230 CALL cp_cfm_get_info(density, local_data=cfm_data)
1231
1232 DO irow = 1, nrow_local
1233 DO jcol = 1, ncol_local
1234
1235 iatom = atom_from_ao_index(row_indices(irow))
1236 jatom = atom_from_ao_index(col_indices(jcol))
1237
1238 IF (ik_global /= ik_old .OR. iatom /= iatom_old .OR. jatom /= jatom_old) THEN
1239
1240 CALL compute_weight_re_im(weight_re, weight_im, &
1241 num_cells, iatom, jatom, xkp(1:3, ik_global), wkp(ik_global), &
1242 cell, index_to_cell, hmat, particle_set)
1243
1244 iatom_old = iatom
1245 jatom_old = jatom
1246 ik_old = ik_global
1247
1248 END IF
1249
1250 contribution = real(cmplx(weight_re, -weight_im, kind=dp)*cfm_data(irow, jcol), kind=dp)
1251
1252 fm_mat_work%local_data(irow, jcol) = fm_mat_work%local_data(irow, jcol) + contribution
1253
1254 END DO
1255 END DO
1256
1257 END DO ! ik
1258
1259 CALL copy_fm_to_dbcsr(fm_mat_work, mat_dm_global_work(1)%matrix, keep_sparsity=.false.)
1260
1261 CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
1262
1263 CALL cp_fm_release(fm_mat_work)
1264 DEALLOCATE (atom_from_ao_index)
1265
1266 CALL timestop(handle)
1267
1268 END SUBROUTINE density_matrix_from_kp_to_mic
1269
1270! **************************************************************************************************
1271!> \brief ...
1272!> \param kpoints ...
1273!> \param density complex AO density scratch matrix
1274!> \param mat_dm_global_work ...
1275!> \param index_to_cell ...
1276!> \param ispin ...
1277!> \param tau ...
1278!> \param e_fermi ...
1279!> \param sector ...
1280! **************************************************************************************************
1281 SUBROUTINE density_matrix_from_kp_to_transl(kpoints, density, mat_dm_global_work, index_to_cell, &
1282 ispin, tau, e_fermi, sector)
1283
1284 TYPE(kpoint_type), POINTER :: kpoints
1285 TYPE(cp_cfm_type), INTENT(INOUT) :: density
1286 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: mat_dm_global_work
1287 INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell
1288 INTEGER, INTENT(IN) :: ispin
1289 REAL(kind=dp), INTENT(IN) :: tau, e_fermi
1290 INTEGER, INTENT(IN) :: sector
1291
1292 CHARACTER(LEN=*), PARAMETER :: routinen = 'density_matrix_from_kp_to_transl'
1293
1294 INTEGER :: handle, icell, ik, ik_global, xcell, &
1295 ycell, zcell
1296 REAL(kind=dp) :: arg, coskl, sinkl
1297 REAL(kind=dp), DIMENSION(:), POINTER :: wkp
1298 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
1299 TYPE(cp_fm_type) :: fm_mat_im, fm_mat_re
1300 TYPE(dbcsr_type), POINTER :: mat_work_im, mat_work_re
1301 TYPE(kpoint_env_type), POINTER :: kp
1302
1303 CALL timeset(routinen, handle)
1304
1305 NULLIFY (xkp, wkp)
1306
1307 NULLIFY (mat_work_re)
1308 CALL dbcsr_init_p(mat_work_re)
1309 CALL dbcsr_create(matrix=mat_work_re, &
1310 template=mat_dm_global_work(1)%matrix, &
1311 matrix_type=dbcsr_type_no_symmetry)
1312
1313 NULLIFY (mat_work_im)
1314 CALL dbcsr_init_p(mat_work_im)
1315 CALL dbcsr_create(matrix=mat_work_im, &
1316 template=mat_dm_global_work(1)%matrix, &
1317 matrix_type=dbcsr_type_no_symmetry)
1318
1319 CALL cp_fm_create(fm_mat_re, density%matrix_struct)
1320 CALL cp_fm_create(fm_mat_im, density%matrix_struct)
1321
1322 CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp)
1323
1324 cpassert(SIZE(mat_dm_global_work) == SIZE(index_to_cell, 2))
1325
1326 DO icell = 1, SIZE(mat_dm_global_work)
1327
1328 CALL dbcsr_set(mat_dm_global_work(icell)%matrix, 0.0_dp)
1329
1330 END DO
1331
1332 DO ik = 1, SIZE(kpoints%kp_env)
1333
1334 kp => kpoints%kp_env(ik)%kpoint_env
1335 ik_global = kp%nkpoint
1336 CALL build_kpoint_density_matrix_rpa(kpoints, ik, ispin, tau, e_fermi, sector, density)
1337 CALL cp_cfm_to_fm(density, fm_mat_re, fm_mat_im)
1338
1339 CALL copy_fm_to_dbcsr(fm_mat_re, mat_work_re, keep_sparsity=.false.)
1340 CALL copy_fm_to_dbcsr(fm_mat_im, mat_work_im, keep_sparsity=.false.)
1341
1342 DO icell = 1, SIZE(mat_dm_global_work)
1343
1344 xcell = index_to_cell(1, icell)
1345 ycell = index_to_cell(2, icell)
1346 zcell = index_to_cell(3, icell)
1347
1348 arg = real(xcell, dp)*xkp(1, ik_global) + real(ycell, dp)*xkp(2, ik_global) + &
1349 REAL(zcell, dp)*xkp(3, ik_global)
1350 coskl = wkp(ik_global)*cos(twopi*arg)
1351 sinkl = wkp(ik_global)*sin(twopi*arg)
1352
1353 CALL dbcsr_add(mat_dm_global_work(icell)%matrix, mat_work_re, 1.0_dp, coskl)
1354 ! The real-space convention is Re(exp(i k R) G(k)).
1355 CALL dbcsr_add(mat_dm_global_work(icell)%matrix, mat_work_im, 1.0_dp, -sinkl)
1356
1357 END DO
1358
1359 END DO
1360
1361 CALL dbcsr_release_p(mat_work_re)
1362 CALL dbcsr_release_p(mat_work_im)
1363 CALL cp_fm_release(fm_mat_re)
1364 CALL cp_fm_release(fm_mat_im)
1365
1366 CALL timestop(handle)
1367
1368 END SUBROUTINE density_matrix_from_kp_to_transl
1369
1370! **************************************************************************************************
1371!> \brief ...
1372!> \param cell_grid ...
1373!> \param cell_to_index ...
1374!> \param index_to_cell ...
1375!> \param cell ...
1376! **************************************************************************************************
1377 SUBROUTINE init_cell_index_rpa(cell_grid, cell_to_index, index_to_cell, cell)
1378 INTEGER, DIMENSION(3), INTENT(IN) :: cell_grid
1379 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1380 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
1381 TYPE(cell_type), INTENT(IN), POINTER :: cell
1382
1383 CHARACTER(LEN=*), PARAMETER :: routinen = 'init_cell_index_rpa'
1384
1385 INTEGER :: cell_counter, handle, i_cell, &
1386 index_min_dist, num_cells, xcell, &
1387 ycell, zcell
1388 INTEGER, DIMENSION(3) :: itm
1389 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_unsorted
1390 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_unsorted
1391 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: abs_cell_vectors
1392 REAL(kind=dp), DIMENSION(3) :: cell_vector
1393 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1394
1395 CALL timeset(routinen, handle)
1396
1397 CALL get_cell(cell=cell, h=hmat)
1398
1399 num_cells = cell_grid(1)*cell_grid(2)*cell_grid(3)
1400 itm(:) = cell_grid(:)/2
1401
1402 ! check that real space super lattice is a (2n+1)x(2m+1)x(2k+1) super lattice with the unit cell
1403 ! in the middle
1404 cpassert(cell_grid(1) /= itm(1)*2)
1405 cpassert(cell_grid(2) /= itm(2)*2)
1406 cpassert(cell_grid(3) /= itm(3)*2)
1407
1408 IF (ASSOCIATED(cell_to_index)) DEALLOCATE (cell_to_index)
1409 IF (ASSOCIATED(index_to_cell)) DEALLOCATE (index_to_cell)
1410
1411 ALLOCATE (cell_to_index_unsorted(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
1412 cell_to_index_unsorted(:, :, :) = 0
1413
1414 ALLOCATE (index_to_cell_unsorted(3, num_cells))
1415 index_to_cell_unsorted(:, :) = 0
1416
1417 ALLOCATE (cell_to_index(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
1418 cell_to_index(:, :, :) = 0
1419
1420 ALLOCATE (index_to_cell(3, num_cells))
1421 index_to_cell(:, :) = 0
1422
1423 ALLOCATE (abs_cell_vectors(1:num_cells))
1424
1425 cell_counter = 0
1426
1427 DO xcell = -itm(1), itm(1)
1428 DO ycell = -itm(2), itm(2)
1429 DO zcell = -itm(3), itm(3)
1430
1431 cell_counter = cell_counter + 1
1432 cell_to_index_unsorted(xcell, ycell, zcell) = cell_counter
1433
1434 index_to_cell_unsorted(1, cell_counter) = xcell
1435 index_to_cell_unsorted(2, cell_counter) = ycell
1436 index_to_cell_unsorted(3, cell_counter) = zcell
1437
1438 cell_vector(1:3) = matmul(hmat, real(index_to_cell_unsorted(1:3, cell_counter), dp))
1439
1440 abs_cell_vectors(cell_counter) = sqrt(cell_vector(1)**2 + cell_vector(2)**2 + cell_vector(3)**2)
1441
1442 END DO
1443 END DO
1444 END DO
1445
1446 ! first only do all symmetry non-equivalent cells, we need that because chi^T is computed for
1447 ! cell indices T from index_to_cell(:,1:num_cells/2+1)
1448 DO i_cell = 1, num_cells/2 + 1
1449
1450 index_min_dist = minloc(abs_cell_vectors(1:num_cells/2 + 1), dim=1)
1451
1452 xcell = index_to_cell_unsorted(1, index_min_dist)
1453 ycell = index_to_cell_unsorted(2, index_min_dist)
1454 zcell = index_to_cell_unsorted(3, index_min_dist)
1455
1456 index_to_cell(1, i_cell) = xcell
1457 index_to_cell(2, i_cell) = ycell
1458 index_to_cell(3, i_cell) = zcell
1459
1460 cell_to_index(xcell, ycell, zcell) = i_cell
1461
1462 abs_cell_vectors(index_min_dist) = 10000000000.0_dp
1463
1464 END DO
1465
1466 ! now all the remaining cells
1467 DO i_cell = num_cells/2 + 2, num_cells
1468
1469 index_min_dist = minloc(abs_cell_vectors(1:num_cells), dim=1)
1470
1471 xcell = index_to_cell_unsorted(1, index_min_dist)
1472 ycell = index_to_cell_unsorted(2, index_min_dist)
1473 zcell = index_to_cell_unsorted(3, index_min_dist)
1474
1475 index_to_cell(1, i_cell) = xcell
1476 index_to_cell(2, i_cell) = ycell
1477 index_to_cell(3, i_cell) = zcell
1478
1479 cell_to_index(xcell, ycell, zcell) = i_cell
1480
1481 abs_cell_vectors(index_min_dist) = 10000000000.0_dp
1482
1483 END DO
1484
1485 DEALLOCATE (index_to_cell_unsorted, cell_to_index_unsorted, abs_cell_vectors)
1486
1487 CALL timestop(handle)
1488
1489 END SUBROUTINE init_cell_index_rpa
1490
1491! **************************************************************************************************
1492!> \brief ...
1493!> \param i_cell_R ...
1494!> \param i_cell_S ...
1495!> \param i_cell_R_minus_S ...
1496!> \param index_to_cell_3c ...
1497!> \param cell_to_index_3c ...
1498!> \param index_to_cell_dm ...
1499!> \param R_minus_S_needed ...
1500!> \param do_kpoints_cubic_RPA ...
1501! **************************************************************************************************
1502 SUBROUTINE get_diff_index_3c(i_cell_R, i_cell_S, i_cell_R_minus_S, index_to_cell_3c, &
1503 cell_to_index_3c, index_to_cell_dm, R_minus_S_needed, &
1504 do_kpoints_cubic_RPA)
1505
1506 INTEGER, INTENT(IN) :: i_cell_r, i_cell_s
1507 INTEGER, INTENT(OUT) :: i_cell_r_minus_s
1508 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(IN) :: index_to_cell_3c
1509 INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
1510 INTENT(IN) :: cell_to_index_3c
1511 INTEGER, DIMENSION(:, :), INTENT(IN), POINTER :: index_to_cell_dm
1512 LOGICAL, INTENT(OUT) :: r_minus_s_needed
1513 LOGICAL, INTENT(IN) :: do_kpoints_cubic_rpa
1514
1515 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_diff_index_3c'
1516
1517 INTEGER :: handle, x_cell_r, x_cell_r_minus_s, x_cell_s, y_cell_r, y_cell_r_minus_s, &
1518 y_cell_s, z_cell_r, z_cell_r_minus_s, z_cell_s
1519
1520 CALL timeset(routinen, handle)
1521
1522 IF (do_kpoints_cubic_rpa) THEN
1523
1524 x_cell_r = index_to_cell_3c(1, i_cell_r)
1525 y_cell_r = index_to_cell_3c(2, i_cell_r)
1526 z_cell_r = index_to_cell_3c(3, i_cell_r)
1527
1528 x_cell_s = index_to_cell_dm(1, i_cell_s)
1529 y_cell_s = index_to_cell_dm(2, i_cell_s)
1530 z_cell_s = index_to_cell_dm(3, i_cell_s)
1531
1532 x_cell_r_minus_s = x_cell_r - x_cell_s
1533 y_cell_r_minus_s = y_cell_r - y_cell_s
1534 z_cell_r_minus_s = z_cell_r - z_cell_s
1535
1536 IF (x_cell_r_minus_s >= lbound(cell_to_index_3c, 1) .AND. &
1537 x_cell_r_minus_s <= ubound(cell_to_index_3c, 1) .AND. &
1538 y_cell_r_minus_s >= lbound(cell_to_index_3c, 2) .AND. &
1539 y_cell_r_minus_s <= ubound(cell_to_index_3c, 2) .AND. &
1540 z_cell_r_minus_s >= lbound(cell_to_index_3c, 3) .AND. &
1541 z_cell_r_minus_s <= ubound(cell_to_index_3c, 3)) THEN
1542
1543 i_cell_r_minus_s = cell_to_index_3c(x_cell_r_minus_s, y_cell_r_minus_s, z_cell_r_minus_s)
1544
1545 ! 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
1546 IF (i_cell_r_minus_s == 0) THEN
1547
1548 r_minus_s_needed = .false.
1549 i_cell_r_minus_s = 0
1550
1551 ELSE
1552
1553 r_minus_s_needed = .true.
1554
1555 END IF
1556
1557 ELSE
1558
1559 i_cell_r_minus_s = 0
1560 r_minus_s_needed = .false.
1561
1562 END IF
1563
1564 ELSE ! no k-points
1565
1566 r_minus_s_needed = .true.
1567 i_cell_r_minus_s = 1
1568
1569 END IF
1570
1571 CALL timestop(handle)
1572
1573 END SUBROUTINE get_diff_index_3c
1574
1575! **************************************************************************************************
1576!> \brief ...
1577!> \param i_cell_R ...
1578!> \param i_cell_S ...
1579!> \param i_cell_T ...
1580!> \param i_cell_R_minus_S_minus_T ...
1581!> \param index_to_cell_3c ...
1582!> \param cell_to_index_3c ...
1583!> \param index_to_cell_dm ...
1584!> \param R_minus_S_minus_T_needed ...
1585!> \param do_kpoints_cubic_RPA ...
1586! **************************************************************************************************
1587 SUBROUTINE get_diff_diff_index_3c(i_cell_R, i_cell_S, i_cell_T, i_cell_R_minus_S_minus_T, &
1588 index_to_cell_3c, cell_to_index_3c, index_to_cell_dm, &
1589 R_minus_S_minus_T_needed, &
1590 do_kpoints_cubic_RPA)
1591
1592 INTEGER, INTENT(IN) :: i_cell_r, i_cell_s, i_cell_t
1593 INTEGER, INTENT(OUT) :: i_cell_r_minus_s_minus_t
1594 INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell_3c
1595 INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
1596 INTENT(IN) :: cell_to_index_3c
1597 INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell_dm
1598 LOGICAL, INTENT(OUT) :: r_minus_s_minus_t_needed
1599 LOGICAL, INTENT(IN) :: do_kpoints_cubic_rpa
1600
1601 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_diff_diff_index_3c'
1602
1603 INTEGER :: handle, x_cell_r, x_cell_r_minus_s_minus_t, x_cell_s, x_cell_t, y_cell_r, &
1604 y_cell_r_minus_s_minus_t, y_cell_s, y_cell_t, z_cell_r, z_cell_r_minus_s_minus_t, &
1605 z_cell_s, z_cell_t
1606
1607 CALL timeset(routinen, handle)
1608
1609 IF (do_kpoints_cubic_rpa) THEN
1610
1611 x_cell_r = index_to_cell_3c(1, i_cell_r)
1612 y_cell_r = index_to_cell_3c(2, i_cell_r)
1613 z_cell_r = index_to_cell_3c(3, i_cell_r)
1614
1615 x_cell_s = index_to_cell_dm(1, i_cell_s)
1616 y_cell_s = index_to_cell_dm(2, i_cell_s)
1617 z_cell_s = index_to_cell_dm(3, i_cell_s)
1618
1619 x_cell_t = index_to_cell_dm(1, i_cell_t)
1620 y_cell_t = index_to_cell_dm(2, i_cell_t)
1621 z_cell_t = index_to_cell_dm(3, i_cell_t)
1622
1623 x_cell_r_minus_s_minus_t = x_cell_r - x_cell_s - x_cell_t
1624 y_cell_r_minus_s_minus_t = y_cell_r - y_cell_s - y_cell_t
1625 z_cell_r_minus_s_minus_t = z_cell_r - z_cell_s - z_cell_t
1626
1627 IF (x_cell_r_minus_s_minus_t >= lbound(cell_to_index_3c, 1) .AND. &
1628 x_cell_r_minus_s_minus_t <= ubound(cell_to_index_3c, 1) .AND. &
1629 y_cell_r_minus_s_minus_t >= lbound(cell_to_index_3c, 2) .AND. &
1630 y_cell_r_minus_s_minus_t <= ubound(cell_to_index_3c, 2) .AND. &
1631 z_cell_r_minus_s_minus_t >= lbound(cell_to_index_3c, 3) .AND. &
1632 z_cell_r_minus_s_minus_t <= ubound(cell_to_index_3c, 3)) THEN
1633
1634 i_cell_r_minus_s_minus_t = cell_to_index_3c(x_cell_r_minus_s_minus_t, &
1635 y_cell_r_minus_s_minus_t, &
1636 z_cell_r_minus_s_minus_t)
1637
1638 ! index 0 means that there are only no 3c matrix elements because R-S-T is too big
1639 IF (i_cell_r_minus_s_minus_t == 0) THEN
1640
1641 r_minus_s_minus_t_needed = .false.
1642
1643 ELSE
1644
1645 r_minus_s_minus_t_needed = .true.
1646
1647 END IF
1648
1649 ELSE
1650
1651 i_cell_r_minus_s_minus_t = 0
1652 r_minus_s_minus_t_needed = .false.
1653
1654 END IF
1655
1656 ! no k-kpoints
1657 ELSE
1658
1659 r_minus_s_minus_t_needed = .true.
1660 i_cell_r_minus_s_minus_t = 1
1661
1662 END IF
1663
1664 CALL timestop(handle)
1665
1666 END SUBROUTINE get_diff_diff_index_3c
1667
1668END 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:233
Basic linear algebra operations for complex full matrices.
subroutine, public cp_cfm_column_scale(matrix_a, scaling)
Scales columns of the full matrix by corresponding factors.
Represents a complex full matrix distributed on many processors.
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
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.
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
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, nrow, ncol, set_zero)
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:16
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, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
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:124
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:141
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, mimic, 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, sab_cneo, 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, matrix_vhxc, 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, xcint_weights, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, do_rixs, tb_tblite)
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 create_propagator_matrix_set(propagator, ntime, matrix_template, index_to_cell)
Creates the sector-, time-, and cell-resolved propagator matrix set.
Definition rpa_im_time.F:97
subroutine, public compute_mat_dm_global(grid, nmo, cfm_mo_coeff, homo, propagator, 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)
...
integer, parameter, public propagator_sector_virtual
Definition rpa_im_time.F:75
subroutine, public compute_periodic_dm(propagator, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, sector)
...
subroutine, public init_cell_index_rpa(cell_grid, cell_to_index, index_to_cell, cell)
...
integer, parameter, public propagator_sector_occupied
Definition rpa_im_time.F:74
subroutine, public compute_transl_dm(propagator, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, eps_filter, num_cells_dm, index_to_cell_dm, sector)
...
subroutine, public compute_gamma_propagator(propagator, jquad, cfm_mo_coeff, homo, eigenval, nmo, eps_filter, e_fermi, tau, para_env)
Builds the Gamma-point occupied and virtual propagators.
subroutine, public compute_mat_p_omega(mat_p_omega, cfm_mo_coeff, homo, 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, grid, e_fermi, eps_filter, alpha, eps_filter_im_time, eigenval, nmo, 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)
...
integer, parameter, public propagator_number_of_sectors
Definition rpa_im_time.F:76
Definition and construction of time/frequency grids for correlation methods.
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
Represent a complex 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