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