Line data Source code
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 : ! **************************************************************************************************
13 : MODULE rpa_im_time
14 : USE cell_types, ONLY: cell_type,&
15 : get_cell
16 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,&
17 : cp_cfm_scale
18 : USE cp_cfm_types, ONLY: cp_cfm_create,&
19 : cp_cfm_get_info,&
20 : cp_cfm_release,&
21 : cp_cfm_to_cfm,&
22 : cp_cfm_to_fm,&
23 : cp_cfm_type
24 : USE cp_dbcsr_api, ONLY: &
25 : dbcsr_add, dbcsr_clear, dbcsr_copy, dbcsr_create, dbcsr_distribution_get, &
26 : dbcsr_distribution_type, dbcsr_filter, dbcsr_get_info, dbcsr_init_p, dbcsr_p_type, &
27 : dbcsr_release_p, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
28 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_all_blocks
29 : USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,&
30 : dbcsr_allocate_matrix_set,&
31 : dbcsr_deallocate_matrix_set
32 : USE cp_fm_types, ONLY: cp_fm_create,&
33 : cp_fm_get_info,&
34 : cp_fm_release,&
35 : cp_fm_set_all,&
36 : cp_fm_type
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,&
42 : hfx_compression_type
43 : USE kinds, ONLY: dp,&
44 : int_8
45 : USE kpoint_types, ONLY: get_kpoint_info,&
46 : kpoint_env_type,&
47 : kpoint_type
48 : USE machine, ONLY: m_flush,&
49 : m_walltime
50 : USE mathconstants, ONLY: twopi
51 : USE message_passing, ONLY: mp_comm_type,&
52 : mp_para_env_type
53 : USE mp2_types, ONLY: mp2_type
54 : USE parallel_gemm_api, ONLY: parallel_gemm
55 : USE particle_types, ONLY: particle_type
56 : USE qs_environment_types, ONLY: get_qs_env,&
57 : qs_environment_type
58 : USE qs_mo_types, ONLY: get_mo_set,&
59 : mo_set_type
60 : USE qs_tensors, ONLY: decompress_tensor,&
61 : get_tensor_occupancy
62 : USE qs_tensors_types, ONLY: create_2c_tensor
63 : USE rpa_gw_im_time_util, ONLY: build_ao_propagator,&
64 : compute_weight_re_im,&
65 : get_atom_index_from_basis_function_index
66 : USE time_frequency_grids, ONLY: time_frequency_grid_type
67 : #include "./base/base_uses.f90"
68 :
69 : IMPLICIT NONE
70 :
71 : PRIVATE
72 :
73 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_im_time'
74 :
75 : INTEGER, PARAMETER, PUBLIC :: propagator_sector_occupied = 1
76 : INTEGER, PARAMETER, PUBLIC :: propagator_sector_virtual = 2
77 : INTEGER, PARAMETER, PUBLIC :: propagator_number_of_sectors = 2
78 :
79 : PUBLIC :: compute_mat_P_omega, &
80 : compute_transl_dm, &
81 : init_cell_index_rpa, &
82 : compute_periodic_dm, &
83 : compute_mat_dm_global, &
84 : compute_gamma_propagator, &
85 : create_propagator_matrix_set
86 :
87 : CONTAINS
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 328 : 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 328 : CPASSERT(ntime > 0)
105 328 : CPASSERT(SIZE(index_to_cell, 1) == 3)
106 328 : ncell = SIZE(index_to_cell, 2)
107 328 : CPASSERT(ncell > 0)
108 :
109 328 : CALL dbcsr_allocate_matrix_set(propagator, propagator_number_of_sectors, ntime, ncell)
110 :
111 668 : DO icell = 1, ncell
112 2562 : DO jtime = 1, ntime
113 6022 : DO isector = 1, propagator_number_of_sectors
114 3788 : ALLOCATE (propagator(isector, jtime, icell)%matrix)
115 : CALL dbcsr_create(matrix=propagator(isector, jtime, icell)%matrix, &
116 : template=matrix_template, &
117 5682 : matrix_type=dbcsr_type_no_symmetry)
118 : END DO
119 : END DO
120 : END DO
121 :
122 328 : END SUBROUTINE create_propagator_matrix_set
123 :
124 : ! **************************************************************************************************
125 : !> \brief ...
126 : !> \param mat_P_omega ...
127 : !> \param fm_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 188 : SUBROUTINE compute_mat_P_omega(mat_P_omega, fm_mo_coeff, homo, &
161 : mat_P_global, &
162 : matrix_s, &
163 : ispin, &
164 376 : t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
165 188 : starts_array_mc, ends_array_mc, &
166 188 : starts_array_mc_block, ends_array_mc_block, &
167 : grid, e_fermi, eps_filter, &
168 188 : 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 188 : 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_fm_type), INTENT(IN) :: fm_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 376 : INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_1, dist_2, mc_ranges, size_dm, &
217 188 : 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 188 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm
223 : LOGICAL :: do_Gamma_RPA, do_kpoints_cubic_RPA, first_cycle_im_time, memory_info, &
224 : 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) :: t1, t2, weight
228 : TYPE(dbcsr_distribution_type) :: dist_P
229 188 : TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: propagator
230 564 : TYPE(dbt_pgrid_type) :: pgrid_2d
231 3572 : TYPE(dbt_type) :: t_3c_M_occ, t_3c_M_occ_tmp, t_3c_M_virt, &
232 5264 : t_3c_M_virt_tmp, t_dm, t_dm_tmp, t_P, &
233 1316 : t_P_tmp
234 376 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_dm_occ, t_dm_virt
235 188 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_O_occ, t_3c_O_virt
236 : TYPE(mp_comm_type) :: comm_2d
237 :
238 188 : CALL timeset(routineN, handle)
239 :
240 1670 : DO jquad = 1, SIZE(mat_P_omega, 1)
241 5888 : DO i_cell_T = 1, SIZE(mat_P_omega, 2)
242 5700 : CALL dbcsr_set(mat_P_omega(jquad, i_cell_T)%matrix, 0.0_dp)
243 : END DO
244 : END DO
245 :
246 188 : IF (do_ri_sos_laplace_mp2) THEN
247 48 : num_integ_points = SIZE(grid%imaginary_time)
248 : ELSE
249 140 : num_integ_points = SIZE(grid%frequency)
250 : END IF
251 :
252 188 : NULLIFY (propagator)
253 :
254 188 : memory_info = mp2_env%ri_rpa_im_time%memory_info
255 188 : IF (memory_info) THEN
256 0 : unit_nr_dbcsr = unit_nr
257 : ELSE
258 188 : unit_nr_dbcsr = 0
259 : END IF
260 :
261 188 : do_kpoints_cubic_RPA = qs_env%mp2_env%ri_rpa_im_time%do_im_time_kpoints
262 188 : do_Gamma_RPA = .NOT. do_kpoints_cubic_RPA
263 824 : num_3c_repl = MAXVAL(cell_to_index_3c)
264 :
265 188 : first_cycle_im_time = .TRUE.
266 4472 : 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)))
267 400 : DO i = 1, SIZE(t_3c_O, 1)
268 732 : DO j = 1, SIZE(t_3c_O, 2)
269 332 : CALL dbt_create(t_3c_O(i, j), t_3c_O_occ(i, j))
270 544 : CALL dbt_create(t_3c_O(i, j), t_3c_O_virt(i, j))
271 : END DO
272 : END DO
273 :
274 188 : CALL dbt_create(t_3c_M, t_3c_M_occ, name="M occ (RI | AO AO)")
275 188 : CALL dbt_create(t_3c_M, t_3c_M_virt, name="M virt (RI | AO AO)")
276 :
277 564 : ALLOCATE (mc_ranges(cut_memory + 1))
278 540 : mc_ranges(:cut_memory) = starts_array_mc_block(:)
279 188 : mc_ranges(cut_memory + 1) = ends_array_mc_block(cut_memory) + 1
280 :
281 1670 : DO jquad = 1, num_integ_points
282 :
283 1482 : CALL para_env%sync()
284 1482 : t1 = m_walltime()
285 :
286 : CALL compute_mat_dm_global(grid, nmo, fm_mo_coeff, homo, propagator, &
287 : matrix_s, ispin, &
288 : Eigenval, e_fermi, eps_filter, memory_info, unit_nr, &
289 : jquad, do_kpoints_cubic_RPA, do_kpoints_from_Gamma, qs_env, &
290 1482 : num_cells_dm, index_to_cell_dm, para_env)
291 :
292 14892 : ALLOCATE (t_dm_virt(num_cells_dm))
293 13410 : ALLOCATE (t_dm_occ(num_cells_dm))
294 1482 : CALL dbcsr_get_info(mat_P_global%matrix, distribution=dist_P)
295 1482 : CALL dbcsr_distribution_get(dist_P, group=comm_2d_handle, nprows=pdims_2d(1), npcols=pdims_2d(2))
296 1482 : CALL comm_2d%set_handle(comm_2d_handle)
297 :
298 1482 : pgrid_2d = dbt_nd_mp_comm(comm_2d, [1], [2], pdims_2d=pdims_2d)
299 4446 : ALLOCATE (size_P(dbt_nblks_total(t_3c_M, 1)))
300 1482 : CALL dbt_get_info(t_3c_M, blk_size_1=size_P)
301 :
302 4446 : ALLOCATE (size_dm(dbt_nblks_total(t_3c_O(1, 1), 3)))
303 1482 : CALL dbt_get_info(t_3c_O(1, 1), blk_size_3=size_dm)
304 1482 : CALL create_2c_tensor(t_dm, dist_1, dist_2, pgrid_2d, size_dm, size_dm, name="D (AO | AO)")
305 1482 : DEALLOCATE (size_dm)
306 1482 : DEALLOCATE (dist_1, dist_2)
307 1482 : CALL create_2c_tensor(t_P, dist_1, dist_2, pgrid_2d, size_P, size_P, name="P (RI | RI)")
308 1482 : DEALLOCATE (size_P)
309 1482 : DEALLOCATE (dist_1, dist_2)
310 1482 : CALL dbt_pgrid_destroy(pgrid_2d)
311 :
312 3036 : DO i_cell = 1, num_cells_dm
313 1554 : CALL dbt_create(t_dm, t_dm_virt(i_cell), name="D virt (AO | AO)")
314 1554 : CALL dbt_create(propagator(propagator_sector_virtual, jquad, i_cell)%matrix, t_dm_tmp)
315 1554 : CALL dbt_copy_matrix_to_tensor(propagator(propagator_sector_virtual, jquad, i_cell)%matrix, t_dm_tmp)
316 1554 : CALL dbt_copy(t_dm_tmp, t_dm_virt(i_cell), move_data=.TRUE.)
317 1554 : CALL dbcsr_clear(propagator(propagator_sector_virtual, jquad, i_cell)%matrix)
318 :
319 1554 : CALL dbt_create(t_dm, t_dm_occ(i_cell), name="D occ (AO | AO)")
320 1554 : CALL dbt_copy_matrix_to_tensor(propagator(propagator_sector_occupied, jquad, i_cell)%matrix, t_dm_tmp)
321 1554 : CALL dbt_copy(t_dm_tmp, t_dm_occ(i_cell), move_data=.TRUE.)
322 1554 : CALL dbt_destroy(t_dm_tmp)
323 3036 : CALL dbcsr_clear(propagator(propagator_sector_occupied, jquad, i_cell)%matrix)
324 : END DO
325 :
326 1482 : CALL get_tensor_occupancy(t_dm_occ(1), nze_dm_occ, occ_dm_occ)
327 1482 : CALL get_tensor_occupancy(t_dm_virt(1), nze_dm_virt, occ_dm_virt)
328 :
329 1482 : CALL dbt_destroy(t_dm)
330 :
331 1482 : CALL dbt_create(t_3c_O_occ(1, 1), t_3c_M_occ_tmp, name="M (RI AO | AO)")
332 1482 : CALL dbt_create(t_3c_O_virt(1, 1), t_3c_M_virt_tmp, name="M (RI AO | AO)")
333 :
334 1482 : CALL timeset(routineN//"_contract", handle2)
335 :
336 1482 : CALL para_env%sync()
337 1482 : t1_flop = m_walltime()
338 :
339 3108 : DO i = 1, SIZE(t_3c_O_occ, 1)
340 5454 : DO j = 1, SIZE(t_3c_O_occ, 2)
341 3972 : CALL dbt_batched_contract_init(t_3c_O_occ(i, j), batch_range_2=mc_ranges, batch_range_3=mc_ranges)
342 : END DO
343 : END DO
344 3108 : DO i = 1, SIZE(t_3c_O_virt, 1)
345 5454 : DO j = 1, SIZE(t_3c_O_virt, 2)
346 3972 : CALL dbt_batched_contract_init(t_3c_O_virt(i, j), batch_range_2=mc_ranges, batch_range_3=mc_ranges)
347 : END DO
348 : END DO
349 1482 : CALL dbt_batched_contract_init(t_3c_M_occ_tmp, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
350 1482 : CALL dbt_batched_contract_init(t_3c_M_virt_tmp, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
351 1482 : CALL dbt_batched_contract_init(t_3c_M_occ, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
352 1482 : CALL dbt_batched_contract_init(t_3c_M_virt, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
353 :
354 3000 : DO i_cell_T = 1, num_cells_dm/2 + 1
355 :
356 1518 : IF (.NOT. ANY(has_mat_P_blocks(i_cell_T, :, :, :, :))) CYCLE
357 :
358 1518 : CALL dbt_batched_contract_init(t_P)
359 :
360 1518 : IF (do_Gamma_RPA) THEN
361 1446 : nze_O = 0
362 1446 : nze_M_virt = 0
363 1446 : nze_M_occ = 0
364 1446 : occ_M_virt = 0.0_dp
365 1446 : occ_M_occ = 0.0_dp
366 1446 : occ_O = 0.0_dp
367 : END IF
368 :
369 3798 : DO j_mem = 1, cut_memory
370 :
371 2280 : CALL dbt_get_info(t_3c_O_occ(1, 1), nfull_total=bounds_3c)
372 :
373 6840 : jbounds_1(:, 1) = [1, bounds_3c(1)]
374 6840 : jbounds_1(:, 2) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
375 :
376 6840 : jbounds_2(:, 1) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
377 :
378 2280 : IF (do_Gamma_RPA) CALL dbt_batched_contract_init(t_dm_virt(1))
379 :
380 6924 : DO i_mem = 1, cut_memory
381 :
382 4804 : IF (.NOT. ANY(has_mat_P_blocks(i_cell_T, i_mem, j_mem, :, :))) CYCLE
383 :
384 13692 : ibounds_1(:, 1) = [1, bounds_3c(1)]
385 13692 : ibounds_1(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
386 :
387 13692 : ibounds_2(:, 1) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
388 :
389 4564 : IF (unit_nr_dbcsr > 0) WRITE (UNIT=unit_nr_dbcsr, FMT="(T3,A,I3,1X,I3)") &
390 0 : "RPA_LOW_SCALING_INFO| Memory Cut iteration", i_mem, j_mem
391 :
392 11984 : DO i_cell_R_1 = 1, num_3c_repl
393 :
394 17804 : DO i_cell_R_2 = 1, num_3c_repl
395 :
396 8020 : IF (.NOT. has_mat_P_blocks(i_cell_T, i_mem, j_mem, i_cell_R_1, i_cell_R_2)) CYCLE
397 :
398 : CALL get_diff_index_3c(i_cell_R_1, i_cell_T, i_cell_R_1_minus_T, &
399 : index_to_cell_3c, cell_to_index_3c, index_to_cell_dm, &
400 5680 : R_1_minus_T_needed, do_kpoints_cubic_RPA)
401 :
402 5680 : IF (do_Gamma_RPA) CALL dbt_batched_contract_init(t_dm_occ(1))
403 13880 : DO i_cell_S = 1, num_cells_dm
404 : CALL get_diff_index_3c(i_cell_R_1, i_cell_S, i_cell_R_1_minus_S, index_to_cell_3c, &
405 : cell_to_index_3c, index_to_cell_dm, R_1_minus_S_needed, &
406 8200 : do_kpoints_cubic_RPA)
407 13880 : IF (R_1_minus_S_needed) THEN
408 :
409 7960 : CALL timeset(routineN//"_calc_M_occ_t", handle3)
410 : CALL decompress_tensor(t_3c_O(i_cell_R_1_minus_S, i_cell_R_2), &
411 : t_3c_O_ind(i_cell_R_1_minus_S, i_cell_R_2, j_mem)%ind, &
412 : t_3c_O_compressed(i_cell_R_1_minus_S, i_cell_R_2, j_mem), &
413 7960 : qs_env%mp2_env%ri_rpa_im_time%eps_compress)
414 :
415 7960 : IF (do_Gamma_RPA .AND. i_mem == 1) THEN
416 2164 : CALL get_tensor_occupancy(t_3c_O(1, 1), nze, occ)
417 2164 : nze_O = nze_O + nze
418 2164 : occ_O = occ_O + occ
419 : END IF
420 :
421 : CALL dbt_copy(t_3c_O(i_cell_R_1_minus_S, i_cell_R_2), &
422 7960 : t_3c_O_occ(i_cell_R_1_minus_S, i_cell_R_2), move_data=.TRUE.)
423 :
424 : CALL dbt_contract(alpha=1.0_dp, &
425 : tensor_1=t_3c_O_occ(i_cell_R_1_minus_S, i_cell_R_2), &
426 : tensor_2=t_dm_occ(i_cell_S), &
427 : beta=1.0_dp, &
428 : tensor_3=t_3c_M_occ_tmp, &
429 : contract_1=[3], notcontract_1=[1, 2], &
430 : contract_2=[2], notcontract_2=[1], &
431 : map_1=[1, 2], map_2=[3], &
432 : bounds_2=jbounds_1, bounds_3=ibounds_2, &
433 : filter_eps=eps_filter, unit_nr=unit_nr_dbcsr, &
434 7960 : flop=flops_1_occ)
435 7960 : CALL timestop(handle3)
436 :
437 7960 : dbcsr_nflop = dbcsr_nflop + flops_1_occ
438 :
439 : END IF
440 : END DO
441 :
442 5680 : IF (do_Gamma_RPA) CALL dbt_batched_contract_finalize(t_dm_occ(1))
443 :
444 : ! copy matrix to optimal contraction layout - copy is done manually in order
445 : ! to better control memory allocations (we can release data of previous
446 : ! representation)
447 5680 : CALL timeset(routineN//"_copy_M_occ_t", handle3)
448 5680 : CALL dbt_copy(t_3c_M_occ_tmp, t_3c_M_occ, order=[1, 3, 2], move_data=.TRUE.)
449 5680 : CALL dbt_filter(t_3c_M_occ, eps_filter)
450 5680 : CALL timestop(handle3)
451 :
452 5680 : IF (do_Gamma_RPA) THEN
453 4420 : CALL get_tensor_occupancy(t_3c_M_occ, nze, occ)
454 4420 : nze_M_occ = nze_M_occ + nze
455 4420 : occ_M_occ = occ_M_occ + occ
456 : END IF
457 :
458 13880 : DO i_cell_S = 1, num_cells_dm
459 : CALL get_diff_diff_index_3c(i_cell_R_2, i_cell_S, i_cell_T, i_cell_R_2_minus_S_minus_T, &
460 : index_to_cell_3c, cell_to_index_3c, index_to_cell_dm, &
461 8200 : R_2_minus_S_minus_T_needed, do_kpoints_cubic_RPA)
462 :
463 13880 : IF (R_1_minus_T_needed .AND. R_2_minus_S_minus_T_needed) THEN
464 : CALL decompress_tensor(t_3c_O(i_cell_R_2_minus_S_minus_T, i_cell_R_1_minus_T), &
465 : t_3c_O_ind(i_cell_R_2_minus_S_minus_T, i_cell_R_1_minus_T, i_mem)%ind, &
466 : t_3c_O_compressed(i_cell_R_2_minus_S_minus_T, i_cell_R_1_minus_T, i_mem), &
467 7756 : qs_env%mp2_env%ri_rpa_im_time%eps_compress)
468 :
469 : CALL dbt_copy(t_3c_O(i_cell_R_2_minus_S_minus_T, i_cell_R_1_minus_T), &
470 7756 : t_3c_O_virt(i_cell_R_2_minus_S_minus_T, i_cell_R_1_minus_T), move_data=.TRUE.)
471 :
472 7756 : CALL timeset(routineN//"_calc_M_virt_t", handle3)
473 : CALL dbt_contract(alpha=alpha/2.0_dp, &
474 : tensor_1=t_3c_O_virt( &
475 : i_cell_R_2_minus_S_minus_T, i_cell_R_1_minus_T), &
476 : tensor_2=t_dm_virt(i_cell_S), &
477 : beta=1.0_dp, &
478 : tensor_3=t_3c_M_virt_tmp, &
479 : contract_1=[3], notcontract_1=[1, 2], &
480 : contract_2=[2], notcontract_2=[1], &
481 : map_1=[1, 2], map_2=[3], &
482 : bounds_2=ibounds_1, bounds_3=jbounds_2, &
483 : filter_eps=eps_filter, unit_nr=unit_nr_dbcsr, &
484 7756 : flop=flops_1_virt)
485 7756 : CALL timestop(handle3)
486 :
487 7756 : dbcsr_nflop = dbcsr_nflop + flops_1_virt
488 :
489 : END IF
490 : END DO
491 :
492 5680 : CALL timeset(routineN//"_copy_M_virt_t", handle3)
493 5680 : CALL dbt_copy(t_3c_M_virt_tmp, t_3c_M_virt, move_data=.TRUE.)
494 5680 : CALL dbt_filter(t_3c_M_virt, eps_filter)
495 5680 : CALL timestop(handle3)
496 :
497 5680 : IF (do_Gamma_RPA) THEN
498 4420 : CALL get_tensor_occupancy(t_3c_M_virt, nze, occ)
499 4420 : nze_M_virt = nze_M_virt + nze
500 4420 : occ_M_virt = occ_M_virt + occ
501 : END IF
502 :
503 : flops_2 = 0
504 :
505 5680 : CALL timeset(routineN//"_calc_P_t", handle3)
506 :
507 : CALL dbt_contract(alpha=1.0_dp, tensor_1=t_3c_M_occ, &
508 : tensor_2=t_3c_M_virt, &
509 : beta=1.0_dp, &
510 : tensor_3=t_P, &
511 : contract_1=[2, 3], notcontract_1=[1], &
512 : contract_2=[2, 3], notcontract_2=[1], &
513 : map_1=[1], map_2=[2], &
514 : filter_eps=eps_filter_im_time/REAL(cut_memory**2, KIND=dp), &
515 : flop=flops_2, &
516 : move_data=.TRUE., &
517 5680 : unit_nr=unit_nr_dbcsr)
518 :
519 5680 : CALL timestop(handle3)
520 :
521 5680 : first_cycle_im_time = .FALSE.
522 :
523 33540 : IF (jquad == 1 .AND. flops_2 == 0) THEN
524 484 : has_mat_P_blocks(i_cell_T, i_mem, j_mem, i_cell_R_1, i_cell_R_2) = .FALSE.
525 : END IF
526 :
527 : END DO
528 : END DO
529 : END DO
530 3798 : IF (do_Gamma_RPA) CALL dbt_batched_contract_finalize(t_dm_virt(1))
531 : END DO
532 :
533 1518 : CALL dbt_batched_contract_finalize(t_P, unit_nr=unit_nr_dbcsr)
534 :
535 1518 : CALL dbt_create(mat_P_global%matrix, t_P_tmp)
536 1518 : CALL dbt_copy(t_P, t_P_tmp, move_data=.TRUE.)
537 1518 : CALL dbt_copy_tensor_to_matrix(t_P_tmp, mat_P_global%matrix)
538 1518 : CALL dbt_destroy(t_P_tmp)
539 :
540 3000 : IF (do_ri_sos_laplace_mp2) THEN
541 : ! For RI-SOS-Laplace-MP2 we do not perform a cosine transform,
542 : ! but we have to copy P_local to the output matrix
543 :
544 188 : CALL dbcsr_add(mat_P_omega(jquad, i_cell_T)%matrix, mat_P_global%matrix, 1.0_dp, 1.0_dp)
545 : ELSE
546 1330 : CALL timeset(routineN//"_Fourier_transform", handle3)
547 :
548 : ! Fourier transform of P(it) to P(iw)
549 24560 : DO iquad = 1, num_integ_points
550 :
551 23230 : weight = grid%cosine_time_to_frequency(iquad, jquad)
552 :
553 : ! no multiplication with 2.0 as in Kresses paper (Kaltak, JCTC 10, 2498 (2014), Eq. 12)
554 : ! because this factor is already absorbed in the weight w_j
555 : CALL dbcsr_add(mat_P_omega(iquad, i_cell_T)%matrix, mat_P_global%matrix, &
556 24560 : 1.0_dp, weight)
557 :
558 : END DO
559 :
560 1330 : CALL timestop(handle3)
561 : END IF
562 :
563 : END DO
564 :
565 1482 : CALL timestop(handle2)
566 :
567 1482 : CALL dbt_batched_contract_finalize(t_3c_M_occ_tmp)
568 1482 : CALL dbt_batched_contract_finalize(t_3c_M_virt_tmp)
569 1482 : CALL dbt_batched_contract_finalize(t_3c_M_occ)
570 1482 : CALL dbt_batched_contract_finalize(t_3c_M_virt)
571 :
572 3108 : DO i = 1, SIZE(t_3c_O_occ, 1)
573 5454 : DO j = 1, SIZE(t_3c_O_occ, 2)
574 3972 : CALL dbt_batched_contract_finalize(t_3c_O_occ(i, j))
575 : END DO
576 : END DO
577 :
578 3108 : DO i = 1, SIZE(t_3c_O_virt, 1)
579 5454 : DO j = 1, SIZE(t_3c_O_virt, 2)
580 3972 : CALL dbt_batched_contract_finalize(t_3c_O_virt(i, j))
581 : END DO
582 : END DO
583 :
584 1482 : CALL dbt_destroy(t_P)
585 3036 : DO i_cell = 1, num_cells_dm
586 1554 : CALL dbt_destroy(t_dm_virt(i_cell))
587 3036 : CALL dbt_destroy(t_dm_occ(i_cell))
588 : END DO
589 :
590 1482 : CALL dbt_destroy(t_3c_M_occ_tmp)
591 1482 : CALL dbt_destroy(t_3c_M_virt_tmp)
592 3036 : DEALLOCATE (t_dm_virt)
593 3036 : DEALLOCATE (t_dm_occ)
594 :
595 1482 : CALL para_env%sync()
596 1482 : t2 = m_walltime()
597 :
598 1482 : dbcsr_time = dbcsr_time + t2 - t1_flop
599 :
600 6116 : IF (unit_nr > 0) THEN
601 : WRITE (unit_nr, '(/T3,A,1X,I3)') &
602 741 : 'RPA_LOW_SCALING_INFO| Info for time point', jquad
603 : WRITE (unit_nr, '(T6,A,T56,F25.1)') &
604 741 : 'Execution time (s):', t2 - t1
605 : WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
606 741 : 'Occupancy of D occ:', REAL(nze_dm_occ, dp), '/', occ_dm_occ*100, '%'
607 : WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
608 741 : 'Occupancy of D virt:', REAL(nze_dm_virt, dp), '/', occ_dm_virt*100, '%'
609 741 : IF (do_Gamma_RPA) THEN
610 : WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
611 723 : 'Occupancy of 3c ints:', REAL(nze_O, dp), '/', occ_O*100, '%'
612 : WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
613 723 : 'Occupancy of M occ:', REAL(nze_M_occ, dp), '/', occ_M_occ*100, '%'
614 : WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
615 723 : 'Occupancy of M virt:', REAL(nze_M_virt, dp), '/', occ_M_virt*100, '%'
616 : END IF
617 741 : WRITE (unit_nr, *)
618 741 : CALL m_flush(unit_nr)
619 : END IF
620 :
621 : END DO ! time points
622 :
623 188 : CALL dbt_destroy(t_3c_M_occ)
624 188 : CALL dbt_destroy(t_3c_M_virt)
625 :
626 400 : DO i = 1, SIZE(t_3c_O, 1)
627 732 : DO j = 1, SIZE(t_3c_O, 2)
628 332 : CALL dbt_destroy(t_3c_O_occ(i, j))
629 544 : CALL dbt_destroy(t_3c_O_virt(i, j))
630 : END DO
631 : END DO
632 :
633 188 : CALL dbcsr_deallocate_matrix_set(propagator)
634 :
635 188 : CALL timestop(handle)
636 :
637 1416 : END SUBROUTINE compute_mat_P_omega
638 :
639 : ! **************************************************************************************************
640 : !> \brief ...
641 : !> \param grid ...
642 : !> \param nmo ...
643 : !> \param fm_mo_coeff ...
644 : !> \param homo ...
645 : !> \param propagator ...
646 : !> \param matrix_s ...
647 : !> \param ispin ...
648 : !> \param Eigenval ...
649 : !> \param e_fermi ...
650 : !> \param eps_filter ...
651 : !> \param memory_info ...
652 : !> \param unit_nr ...
653 : !> \param jquad ...
654 : !> \param do_kpoints_cubic_RPA ...
655 : !> \param do_kpoints_from_Gamma ...
656 : !> \param qs_env ...
657 : !> \param num_cells_dm ...
658 : !> \param index_to_cell_dm ...
659 : !> \param para_env ...
660 : ! **************************************************************************************************
661 3304 : SUBROUTINE compute_mat_dm_global(grid, nmo, fm_mo_coeff, homo, propagator, &
662 1652 : matrix_s, ispin, &
663 1652 : Eigenval, e_fermi, eps_filter, memory_info, unit_nr, &
664 : jquad, do_kpoints_cubic_RPA, do_kpoints_from_Gamma, qs_env, &
665 : num_cells_dm, index_to_cell_dm, para_env)
666 :
667 : TYPE(time_frequency_grid_type), INTENT(IN) :: grid
668 : INTEGER, INTENT(IN) :: nmo
669 : TYPE(cp_fm_type), INTENT(IN) :: fm_mo_coeff
670 : INTEGER, INTENT(IN) :: homo
671 : TYPE(dbcsr_p_type), DIMENSION(:, :, :), &
672 : INTENT(INOUT), POINTER :: propagator
673 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: matrix_s
674 : INTEGER, INTENT(IN) :: ispin
675 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval
676 : REAL(KIND=dp), INTENT(IN) :: e_fermi, eps_filter
677 : LOGICAL, INTENT(IN) :: memory_info
678 : INTEGER, INTENT(IN) :: unit_nr, jquad
679 : LOGICAL, INTENT(IN) :: do_kpoints_cubic_RPA, &
680 : do_kpoints_from_Gamma
681 : TYPE(qs_environment_type), POINTER :: qs_env
682 : INTEGER, INTENT(OUT) :: num_cells_dm
683 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm
684 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
685 :
686 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_mat_dm_global'
687 :
688 : INTEGER :: handle, num_integ_points
689 : INTEGER, DIMENSION(3, 1) :: index_to_cell_zero
690 : REAL(KIND=dp) :: tau
691 :
692 1652 : index_to_cell_zero = 0
693 1652 : num_integ_points = SIZE(grid%imaginary_time)
694 :
695 1652 : CALL timeset(routineN, handle)
696 :
697 1652 : IF (memory_info .AND. unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
698 0 : "RPA_LOW_SCALING_INFO| Started with time point: ", jquad
699 :
700 1652 : tau = grid%imaginary_time(jquad)
701 :
702 1652 : IF (do_kpoints_cubic_RPA) THEN
703 :
704 : CALL compute_transl_dm(propagator, qs_env, &
705 : ispin, num_integ_points, jquad, e_fermi, tau, &
706 : eps_filter, num_cells_dm, index_to_cell_dm, &
707 36 : sector=propagator_sector_occupied)
708 :
709 : CALL compute_transl_dm(propagator, qs_env, &
710 : ispin, num_integ_points, jquad, e_fermi, tau, &
711 : eps_filter, num_cells_dm, index_to_cell_dm, &
712 36 : sector=propagator_sector_virtual)
713 :
714 1616 : ELSE IF (do_kpoints_from_Gamma) THEN
715 :
716 : CALL compute_periodic_dm(propagator, qs_env, &
717 : ispin, num_integ_points, jquad, e_fermi, tau, &
718 108 : sector=propagator_sector_occupied)
719 :
720 : CALL compute_periodic_dm(propagator, qs_env, &
721 : ispin, num_integ_points, jquad, e_fermi, tau, &
722 108 : sector=propagator_sector_virtual)
723 :
724 108 : num_cells_dm = 1
725 :
726 : ELSE
727 :
728 1508 : num_cells_dm = 1
729 :
730 1508 : IF (jquad == 1) THEN
731 : CALL create_propagator_matrix_set(propagator, num_integ_points, matrix_s(1)%matrix, &
732 226 : index_to_cell_zero)
733 : END IF
734 :
735 : CALL compute_gamma_propagator(propagator, jquad, fm_mo_coeff, homo, Eigenval, nmo, &
736 1508 : eps_filter, e_fermi, tau, para_env)
737 :
738 : ! release memory
739 1508 : IF (jquad > 1) THEN
740 1282 : CALL dbcsr_set(propagator(propagator_sector_occupied, jquad - 1, 1)%matrix, 0.0_dp)
741 1282 : CALL dbcsr_set(propagator(propagator_sector_virtual, jquad - 1, 1)%matrix, 0.0_dp)
742 1282 : CALL dbcsr_filter(propagator(propagator_sector_occupied, jquad - 1, 1)%matrix, 0.0_dp)
743 1282 : CALL dbcsr_filter(propagator(propagator_sector_virtual, jquad - 1, 1)%matrix, 0.0_dp)
744 : END IF
745 :
746 : END IF ! do kpoints
747 :
748 1652 : CALL timestop(handle)
749 :
750 1652 : END SUBROUTINE compute_mat_dm_global
751 :
752 : ! **************************************************************************************************
753 : !> \brief Builds the Gamma-point occupied and virtual propagators.
754 : !> \param propagator sector/time/cell propagator matrix set
755 : !> \param jquad time-point index
756 : !> \param fm_mo_coeff real Gamma-point MO coefficient matrix
757 : !> \param homo number of occupied orbitals
758 : !> \param Eigenval orbital energies
759 : !> \param nmo number of orbitals
760 : !> \param eps_filter DBCSR filtering threshold
761 : !> \param e_fermi Fermi energy
762 : !> \param tau imaginary time
763 : !> \param para_env parallel environment
764 : ! **************************************************************************************************
765 2456 : SUBROUTINE compute_gamma_propagator(propagator, jquad, fm_mo_coeff, homo, Eigenval, nmo, &
766 : eps_filter, e_fermi, tau, para_env)
767 : TYPE(dbcsr_p_type), DIMENSION(:, :, :), &
768 : INTENT(INOUT), POINTER :: propagator
769 : INTEGER, INTENT(IN) :: jquad
770 : TYPE(cp_fm_type), INTENT(IN) :: fm_mo_coeff
771 : INTEGER, INTENT(IN) :: homo
772 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval
773 : INTEGER, INTENT(IN) :: nmo
774 : REAL(KIND=dp), INTENT(IN) :: eps_filter, e_fermi, tau
775 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
776 :
777 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_gamma_propagator'
778 : REAL(KIND=dp), PARAMETER :: stabilize_exp = 70.0_dp
779 :
780 : INTEGER :: first_mo, handle, i_mo, isector, &
781 : last_mo, nao
782 4912 : REAL(KIND=dp) :: column_factors(nmo), tau_E
783 : TYPE(cp_fm_type) :: fm_mo_coeff_scaled, fm_scaled_dm_tau
784 : TYPE(dbcsr_type), POINTER :: mat
785 :
786 2456 : CALL timeset(routineN, handle)
787 2456 : CALL para_env%sync()
788 :
789 2456 : CALL cp_fm_get_info(fm_mo_coeff, nrow_global=nao)
790 2456 : CALL cp_fm_create(fm_scaled_dm_tau, fm_mo_coeff%matrix_struct, nrow=nao, ncol=nao)
791 2456 : CALL cp_fm_create(fm_mo_coeff_scaled, fm_mo_coeff%matrix_struct)
792 :
793 64556 : DO i_mo = 1, nmo
794 62100 : tau_E = tau*0.5_dp*(Eigenval(i_mo) - e_fermi)
795 64556 : IF (ABS(tau_E) < stabilize_exp) THEN
796 50888 : IF (i_mo > homo) tau_E = -tau_E
797 50888 : column_factors(i_mo) = EXP(tau_E)
798 : ELSE
799 11212 : column_factors(i_mo) = 0.0_dp
800 : END IF
801 : END DO
802 :
803 2456 : CALL para_env%sync()
804 :
805 7368 : DO isector = propagator_sector_occupied, propagator_sector_virtual
806 4912 : first_mo = 1
807 4912 : last_mo = homo
808 4912 : IF (isector == propagator_sector_virtual) THEN
809 2456 : first_mo = homo + 1
810 2456 : last_mo = nmo
811 : END IF
812 : CALL build_ao_propagator(fm_mo_coeff, column_factors(first_mo:last_mo), first_mo, &
813 4912 : fm_mo_coeff_scaled, fm_scaled_dm_tau)
814 :
815 : ! Each sector is copied before the dense coefficient and AO workspaces are reused.
816 4912 : mat => propagator(isector, jquad, 1)%matrix
817 4912 : CALL dbcsr_set(mat, 0.0_dp)
818 4912 : CALL copy_fm_to_dbcsr(fm_scaled_dm_tau, mat, keep_sparsity=.FALSE.)
819 7368 : CALL dbcsr_filter(mat, eps_filter)
820 : END DO
821 :
822 2456 : CALL cp_fm_release(fm_scaled_dm_tau)
823 2456 : CALL cp_fm_release(fm_mo_coeff_scaled)
824 2456 : CALL timestop(handle)
825 :
826 2456 : END SUBROUTINE compute_gamma_propagator
827 :
828 : ! **************************************************************************************************
829 : !> \brief Calculate one kpoint density matrix in the complex AO representation.
830 : !> \param kpoint kpoint environment
831 : !> \param ikpgr local kpoint index
832 : !> \param ispin spin index
833 : !> \param tau ...
834 : !> \param e_fermi ...
835 : !> \param sector ...
836 : !> \param density complex AO density matrix
837 : ! **************************************************************************************************
838 604 : SUBROUTINE build_kpoint_density_matrix_rpa(kpoint, ikpgr, ispin, tau, e_fermi, sector, density)
839 :
840 : TYPE(kpoint_type), POINTER :: kpoint
841 : INTEGER, INTENT(IN) :: ikpgr, ispin
842 : REAL(KIND=dp), INTENT(IN) :: tau, e_fermi
843 : INTEGER, INTENT(IN) :: sector
844 : TYPE(cp_cfm_type), INTENT(INOUT) :: density
845 :
846 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_kpoint_density_matrix_rpa'
847 : REAL(KIND=dp), PARAMETER :: stabilize_exp = 70.0_dp
848 :
849 : INTEGER :: handle, i_mo, nao, nmo, nspin
850 604 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, exp_scaling, occupation
851 : TYPE(cp_cfm_type) :: cfm_work
852 : TYPE(cp_cfm_type), POINTER :: cfm_mo_coeff
853 : TYPE(kpoint_env_type), POINTER :: kp
854 : TYPE(mo_set_type), POINTER :: mo_set
855 :
856 604 : CALL timeset(routineN, handle)
857 :
858 604 : CPASSERT(sector == propagator_sector_occupied .OR. sector == propagator_sector_virtual)
859 604 : CPASSERT(kpoint%use_real_wfn .EQV. .FALSE.)
860 :
861 604 : kp => kpoint%kp_env(ikpgr)%kpoint_env
862 604 : nspin = SIZE(kp%mos)
863 604 : CPASSERT(ispin >= 1 .AND. ispin <= nspin)
864 604 : mo_set => kp%mos(ispin)
865 : CALL get_mo_set(mo_set, cmo_coeff=cfm_mo_coeff, nao=nao, nmo=nmo, &
866 604 : eigenvalues=eigenvalues, occupation_numbers=occupation)
867 :
868 : ! if this CPASSERT is triggered, please add all virtual MOs to SCF section,
869 : ! e.g. ADDED_MOS 1000000
870 604 : CPASSERT(nao == nmo)
871 :
872 1812 : ALLOCATE (exp_scaling(nmo))
873 604 : CALL cp_cfm_create(cfm_work, cfm_mo_coeff%matrix_struct)
874 604 : CALL cp_cfm_to_cfm(cfm_mo_coeff, cfm_work)
875 :
876 604 : IF (sector == propagator_sector_occupied) THEN
877 6468 : CALL cp_cfm_column_scale(cfm_work, CMPLX(occupation, 0.0_dp, KIND=dp))
878 : ELSE
879 : CALL cp_cfm_column_scale(cfm_work, &
880 6468 : CMPLX(2.0_dp/REAL(nspin, KIND=dp) - occupation, 0.0_dp, KIND=dp))
881 : END IF
882 :
883 604 : IF (nspin == 1) THEN
884 500 : CALL cp_cfm_scale(0.5_dp, cfm_work)
885 : END IF
886 :
887 12936 : DO i_mo = 1, nmo
888 12936 : IF (ABS(tau*0.5_dp*(eigenvalues(i_mo) - e_fermi)) < stabilize_exp) THEN
889 12332 : exp_scaling(i_mo) = EXP(-ABS(tau*(eigenvalues(i_mo) - e_fermi)))
890 : ELSE
891 0 : exp_scaling(i_mo) = 0.0_dp
892 : END IF
893 : END DO
894 :
895 12936 : CALL cp_cfm_column_scale(cfm_work, CMPLX(exp_scaling, 0.0_dp, KIND=dp))
896 : CALL parallel_gemm("N", "C", nao, nao, nmo, (1.0_dp, 0.0_dp), cfm_mo_coeff, cfm_work, &
897 604 : (0.0_dp, 0.0_dp), density)
898 :
899 604 : CALL cp_cfm_release(cfm_work)
900 604 : DEALLOCATE (exp_scaling)
901 :
902 604 : CALL timestop(handle)
903 :
904 1812 : END SUBROUTINE build_kpoint_density_matrix_rpa
905 :
906 : ! **************************************************************************************************
907 : !> \brief ...
908 : !> \param propagator ...
909 : !> \param qs_env ...
910 : !> \param ispin ...
911 : !> \param num_integ_points ...
912 : !> \param jquad ...
913 : !> \param e_fermi ...
914 : !> \param tau ...
915 : !> \param eps_filter ...
916 : !> \param num_cells_dm ...
917 : !> \param index_to_cell_dm ...
918 : !> \param sector ...
919 : ! **************************************************************************************************
920 216 : SUBROUTINE compute_transl_dm(propagator, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, &
921 : eps_filter, num_cells_dm, index_to_cell_dm, &
922 : sector)
923 : TYPE(dbcsr_p_type), DIMENSION(:, :, :), &
924 : INTENT(INOUT), POINTER :: propagator
925 : TYPE(qs_environment_type), POINTER :: qs_env
926 : INTEGER, INTENT(IN) :: ispin, num_integ_points, jquad
927 : REAL(KIND=dp), INTENT(IN) :: e_fermi, tau, eps_filter
928 : INTEGER, INTENT(OUT) :: num_cells_dm
929 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm
930 : INTEGER, INTENT(IN) :: sector
931 :
932 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_transl_dm'
933 :
934 : INTEGER :: handle, i_dim, i_img, nao
935 : INTEGER, DIMENSION(3) :: cell_grid_dm
936 : TYPE(cell_type), POINTER :: cell
937 : TYPE(cp_cfm_type) :: cfm_density
938 72 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_dm_global_work
939 72 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp
940 : TYPE(kpoint_type), POINTER :: kpoints
941 :
942 72 : CALL timeset(routineN, handle)
943 :
944 : CALL get_qs_env(qs_env, &
945 : matrix_s_kp=matrix_s_kp, &
946 : kpoints=kpoints, &
947 72 : cell=cell)
948 :
949 : ! we always use an odd number of image cells
950 : ! CAUTION: also at another point, cell_grid_dm is defined, these definitions have to be identical
951 288 : DO i_dim = 1, 3
952 288 : cell_grid_dm(i_dim) = (kpoints%nkp_grid(i_dim)/2)*2 - 1
953 : END DO
954 :
955 72 : num_cells_dm = cell_grid_dm(1)*cell_grid_dm(2)*cell_grid_dm(3)
956 :
957 72 : NULLIFY (mat_dm_global_work)
958 : ! Only the requested spin is built.
959 72 : CALL dbcsr_allocate_matrix_set(mat_dm_global_work, num_cells_dm)
960 :
961 288 : DO i_img = 1, num_cells_dm
962 :
963 216 : ALLOCATE (mat_dm_global_work(i_img)%matrix)
964 : CALL dbcsr_create(matrix=mat_dm_global_work(i_img)%matrix, &
965 : template=matrix_s_kp(1, 1)%matrix, &
966 : ! matrix_type=dbcsr_type_symmetric)
967 216 : matrix_type=dbcsr_type_no_symmetry)
968 :
969 288 : CALL dbcsr_reserve_all_blocks(mat_dm_global_work(i_img)%matrix)
970 :
971 : END DO
972 :
973 : ! Reusable complex AO scratch for the k-point density source.
974 72 : CALL get_mo_set(kpoints%kp_env(1)%kpoint_env%mos(1), nao=nao)
975 : CALL cp_cfm_create(cfm_density, kpoints%kp_env(1)%kpoint_env%mos(1)%cmo_coeff%matrix_struct, &
976 72 : nrow=nao, ncol=nao)
977 :
978 : ! overwrite the cell indices in kpoints
979 72 : CALL init_cell_index_rpa(cell_grid_dm, kpoints%cell_to_index, kpoints%index_to_cell, cell)
980 :
981 : ! density matrices in real space, the cell vectors T for transforming are taken from kpoints%index_to_cell
982 : ! (custom made for RPA) and not from sab_nl (which is symmetric and from SCF)
983 : CALL density_matrix_from_kp_to_transl(kpoints, cfm_density, mat_dm_global_work, kpoints%index_to_cell, &
984 72 : ispin, tau, e_fermi, sector)
985 :
986 72 : CALL cp_cfm_release(cfm_density)
987 :
988 : ! we need the index to cell for the density matrices later
989 72 : index_to_cell_dm => kpoints%index_to_cell
990 :
991 : ! The first request owns allocation of the complete time-indexed propagator set.
992 72 : IF (.NOT. ASSOCIATED(propagator)) THEN
993 : CALL create_propagator_matrix_set(propagator, num_integ_points, matrix_s_kp(1, 1)%matrix, &
994 6 : kpoints%index_to_cell)
995 : END IF
996 :
997 288 : DO i_img = 1, num_cells_dm
998 :
999 : ! filter to get rid of the blocks full with zeros on the lower half, otherwise blocks doubled
1000 216 : CALL dbcsr_filter(mat_dm_global_work(i_img)%matrix, eps_filter)
1001 :
1002 : CALL dbcsr_copy(propagator(sector, jquad, i_img)%matrix, &
1003 288 : mat_dm_global_work(i_img)%matrix)
1004 :
1005 : END DO
1006 :
1007 72 : CALL dbcsr_deallocate_matrix_set(mat_dm_global_work)
1008 :
1009 72 : CALL timestop(handle)
1010 :
1011 72 : END SUBROUTINE compute_transl_dm
1012 :
1013 : ! **************************************************************************************************
1014 : !> \brief ...
1015 : !> \param propagator ...
1016 : !> \param qs_env ...
1017 : !> \param ispin ...
1018 : !> \param num_integ_points ...
1019 : !> \param jquad ...
1020 : !> \param e_fermi ...
1021 : !> \param tau ...
1022 : !> \param sector ...
1023 : ! **************************************************************************************************
1024 1380 : SUBROUTINE compute_periodic_dm(propagator, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, &
1025 : sector)
1026 : TYPE(dbcsr_p_type), DIMENSION(:, :, :), &
1027 : INTENT(INOUT), POINTER :: propagator
1028 : TYPE(qs_environment_type), POINTER :: qs_env
1029 : INTEGER, INTENT(IN) :: ispin, num_integ_points, jquad
1030 : REAL(KIND=dp), INTENT(IN) :: e_fermi, tau
1031 : INTEGER, INTENT(IN) :: sector
1032 :
1033 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_periodic_dm'
1034 :
1035 : INTEGER :: handle, nao, num_cells_dm
1036 : INTEGER, DIMENSION(3, 1) :: index_to_cell_zero
1037 : TYPE(cp_cfm_type) :: cfm_density
1038 460 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_dm_global_work
1039 460 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp
1040 : TYPE(kpoint_type), POINTER :: kpoints_G
1041 :
1042 460 : CALL timeset(routineN, handle)
1043 :
1044 460 : NULLIFY (matrix_s_kp)
1045 :
1046 : CALL get_qs_env(qs_env, &
1047 460 : matrix_s_kp=matrix_s_kp)
1048 :
1049 460 : kpoints_G => qs_env%mp2_env%ri_rpa_im_time%kpoints_G
1050 :
1051 460 : num_cells_dm = 1
1052 460 : index_to_cell_zero = 0
1053 :
1054 460 : NULLIFY (mat_dm_global_work)
1055 : ! Only the requested spin is built.
1056 460 : CALL dbcsr_allocate_matrix_set(mat_dm_global_work, num_cells_dm)
1057 :
1058 : ! The first request owns allocation of the complete time-indexed propagator set.
1059 460 : IF (.NOT. ASSOCIATED(propagator)) THEN
1060 : CALL create_propagator_matrix_set(propagator, num_integ_points, matrix_s_kp(1, 1)%matrix, &
1061 34 : index_to_cell_zero)
1062 : END IF
1063 :
1064 460 : ALLOCATE (mat_dm_global_work(1)%matrix)
1065 : CALL dbcsr_create(matrix=mat_dm_global_work(1)%matrix, &
1066 : template=matrix_s_kp(1, 1)%matrix, &
1067 460 : matrix_type=dbcsr_type_no_symmetry)
1068 :
1069 460 : CALL dbcsr_reserve_all_blocks(mat_dm_global_work(1)%matrix)
1070 :
1071 460 : CALL dbcsr_set(mat_dm_global_work(1)%matrix, 0.0_dp)
1072 :
1073 : ! Reusable complex AO scratch for the k-point density source.
1074 460 : CALL get_mo_set(kpoints_G%kp_env(1)%kpoint_env%mos(1), nao=nao)
1075 : CALL cp_cfm_create(cfm_density, kpoints_G%kp_env(1)%kpoint_env%mos(1)%cmo_coeff%matrix_struct, &
1076 460 : nrow=nao, ncol=nao)
1077 :
1078 460 : CALL density_matrix_from_kp_to_mic(kpoints_G, cfm_density, mat_dm_global_work, qs_env, ispin, tau, e_fermi, sector)
1079 :
1080 460 : CALL cp_cfm_release(cfm_density)
1081 :
1082 : CALL dbcsr_copy(propagator(sector, jquad, 1)%matrix, &
1083 460 : mat_dm_global_work(1)%matrix)
1084 :
1085 460 : CALL dbcsr_deallocate_matrix_set(mat_dm_global_work)
1086 :
1087 460 : CALL timestop(handle)
1088 :
1089 460 : END SUBROUTINE compute_periodic_dm
1090 :
1091 : ! **************************************************************************************************
1092 : !> \brief ...
1093 : !> \param kpoints_G ...
1094 : !> \param density complex AO density scratch matrix
1095 : !> \param mat_dm_global_work ...
1096 : !> \param qs_env ...
1097 : !> \param ispin ...
1098 : !> \param tau ...
1099 : !> \param e_fermi ...
1100 : !> \param sector ...
1101 : ! **************************************************************************************************
1102 460 : SUBROUTINE density_matrix_from_kp_to_mic(kpoints_G, density, mat_dm_global_work, qs_env, ispin, tau, e_fermi, sector)
1103 :
1104 : TYPE(kpoint_type), POINTER :: kpoints_G
1105 : TYPE(cp_cfm_type), INTENT(INOUT) :: density
1106 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_dm_global_work
1107 : TYPE(qs_environment_type), POINTER :: qs_env
1108 : INTEGER, INTENT(IN) :: ispin
1109 : REAL(KIND=dp), INTENT(IN) :: tau, e_fermi
1110 : INTEGER, INTENT(IN) :: sector
1111 :
1112 : CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_from_kp_to_mic'
1113 :
1114 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1115 460 : POINTER :: cfm_data
1116 : INTEGER :: handle, iatom, iatom_old, ik, ik_global, &
1117 : ik_old, irow, jatom, jatom_old, jcol, &
1118 : nao, ncol_local, nrow_local, num_cells
1119 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_ao_index
1120 460 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1121 460 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
1122 : REAL(KIND=dp) :: contribution, weight_im, weight_re
1123 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1124 460 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
1125 460 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1126 : TYPE(cell_type), POINTER :: cell
1127 : TYPE(cp_fm_type) :: fm_mat_work
1128 : TYPE(kpoint_env_type), POINTER :: kp
1129 460 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1130 :
1131 460 : CALL timeset(routineN, handle)
1132 :
1133 460 : NULLIFY (xkp, wkp)
1134 :
1135 460 : CALL get_kpoint_info(kpoints_G, xkp=xkp, wkp=wkp)
1136 460 : index_to_cell => kpoints_G%index_to_cell
1137 460 : num_cells = SIZE(index_to_cell, 2)
1138 :
1139 : CALL cp_cfm_get_info(density, nrow_global=nao, nrow_local=nrow_local, ncol_local=ncol_local, &
1140 460 : row_indices=row_indices, col_indices=col_indices)
1141 460 : CALL cp_fm_create(fm_mat_work, density%matrix_struct, set_zero=.TRUE.)
1142 :
1143 1380 : ALLOCATE (atom_from_ao_index(nao))
1144 :
1145 460 : CALL get_atom_index_from_basis_function_index(qs_env, atom_from_ao_index, nao, "ORB")
1146 :
1147 460 : NULLIFY (cell, particle_set)
1148 460 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1149 460 : CALL get_cell(cell=cell, h=hmat)
1150 :
1151 460 : iatom_old = 0
1152 460 : jatom_old = 0
1153 460 : ik_old = 0
1154 :
1155 460 : CALL dbcsr_set(mat_dm_global_work(1)%matrix, 0.0_dp)
1156 :
1157 920 : DO ik = 1, SIZE(kpoints_G%kp_env)
1158 :
1159 460 : kp => kpoints_G%kp_env(ik)%kpoint_env
1160 460 : ik_global = kp%nkpoint
1161 460 : CALL build_kpoint_density_matrix_rpa(kpoints_G, ik, ispin, tau, e_fermi, sector, density)
1162 460 : CALL cp_cfm_get_info(density, local_data=cfm_data)
1163 :
1164 5430 : DO irow = 1, nrow_local
1165 99184 : DO jcol = 1, ncol_local
1166 :
1167 94214 : iatom = atom_from_ao_index(row_indices(irow))
1168 94214 : jatom = atom_from_ao_index(col_indices(jcol))
1169 :
1170 94214 : IF (ik_global /= ik_old .OR. iatom /= iatom_old .OR. jatom /= jatom_old) THEN
1171 :
1172 : CALL compute_weight_re_im(weight_re, weight_im, &
1173 : num_cells, iatom, jatom, xkp(1:3, ik_global), wkp(ik_global), &
1174 13296 : cell, index_to_cell, hmat, particle_set)
1175 :
1176 13296 : iatom_old = iatom
1177 13296 : jatom_old = jatom
1178 13296 : ik_old = ik_global
1179 :
1180 : END IF
1181 :
1182 94214 : contribution = REAL(CMPLX(weight_re, -weight_im, KIND=dp)*cfm_data(irow, jcol), KIND=dp)
1183 :
1184 98724 : fm_mat_work%local_data(irow, jcol) = fm_mat_work%local_data(irow, jcol) + contribution
1185 :
1186 : END DO
1187 : END DO
1188 :
1189 : END DO ! ik
1190 :
1191 460 : CALL copy_fm_to_dbcsr(fm_mat_work, mat_dm_global_work(1)%matrix, keep_sparsity=.FALSE.)
1192 :
1193 460 : CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
1194 :
1195 460 : CALL cp_fm_release(fm_mat_work)
1196 460 : DEALLOCATE (atom_from_ao_index)
1197 :
1198 460 : CALL timestop(handle)
1199 :
1200 920 : END SUBROUTINE density_matrix_from_kp_to_mic
1201 :
1202 : ! **************************************************************************************************
1203 : !> \brief ...
1204 : !> \param kpoints ...
1205 : !> \param density complex AO density scratch matrix
1206 : !> \param mat_dm_global_work ...
1207 : !> \param index_to_cell ...
1208 : !> \param ispin ...
1209 : !> \param tau ...
1210 : !> \param e_fermi ...
1211 : !> \param sector ...
1212 : ! **************************************************************************************************
1213 72 : SUBROUTINE density_matrix_from_kp_to_transl(kpoints, density, mat_dm_global_work, index_to_cell, &
1214 : ispin, tau, e_fermi, sector)
1215 :
1216 : TYPE(kpoint_type), POINTER :: kpoints
1217 : TYPE(cp_cfm_type), INTENT(INOUT) :: density
1218 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: mat_dm_global_work
1219 : INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell
1220 : INTEGER, INTENT(IN) :: ispin
1221 : REAL(KIND=dp), INTENT(IN) :: tau, e_fermi
1222 : INTEGER, INTENT(IN) :: sector
1223 :
1224 : CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_from_kp_to_transl'
1225 :
1226 : INTEGER :: handle, icell, ik, ik_global, xcell, &
1227 : ycell, zcell
1228 : REAL(KIND=dp) :: arg, coskl, sinkl
1229 72 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
1230 72 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1231 : TYPE(cp_fm_type) :: fm_mat_im, fm_mat_re
1232 : TYPE(dbcsr_type), POINTER :: mat_work_im, mat_work_re
1233 : TYPE(kpoint_env_type), POINTER :: kp
1234 :
1235 72 : CALL timeset(routineN, handle)
1236 :
1237 72 : NULLIFY (xkp, wkp)
1238 :
1239 72 : NULLIFY (mat_work_re)
1240 72 : CALL dbcsr_init_p(mat_work_re)
1241 : CALL dbcsr_create(matrix=mat_work_re, &
1242 : template=mat_dm_global_work(1)%matrix, &
1243 72 : matrix_type=dbcsr_type_no_symmetry)
1244 :
1245 72 : NULLIFY (mat_work_im)
1246 72 : CALL dbcsr_init_p(mat_work_im)
1247 : CALL dbcsr_create(matrix=mat_work_im, &
1248 : template=mat_dm_global_work(1)%matrix, &
1249 72 : matrix_type=dbcsr_type_no_symmetry)
1250 :
1251 72 : CALL cp_fm_create(fm_mat_re, density%matrix_struct)
1252 72 : CALL cp_fm_create(fm_mat_im, density%matrix_struct)
1253 :
1254 72 : CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp)
1255 :
1256 72 : CPASSERT(SIZE(mat_dm_global_work) == SIZE(index_to_cell, 2))
1257 :
1258 288 : DO icell = 1, SIZE(mat_dm_global_work)
1259 :
1260 288 : CALL dbcsr_set(mat_dm_global_work(icell)%matrix, 0.0_dp)
1261 :
1262 : END DO
1263 :
1264 216 : DO ik = 1, SIZE(kpoints%kp_env)
1265 :
1266 144 : kp => kpoints%kp_env(ik)%kpoint_env
1267 144 : ik_global = kp%nkpoint
1268 144 : CALL build_kpoint_density_matrix_rpa(kpoints, ik, ispin, tau, e_fermi, sector, density)
1269 144 : CALL cp_cfm_to_fm(density, fm_mat_re, fm_mat_im)
1270 :
1271 144 : CALL copy_fm_to_dbcsr(fm_mat_re, mat_work_re, keep_sparsity=.FALSE.)
1272 144 : CALL copy_fm_to_dbcsr(fm_mat_im, mat_work_im, keep_sparsity=.FALSE.)
1273 :
1274 648 : DO icell = 1, SIZE(mat_dm_global_work)
1275 :
1276 432 : xcell = index_to_cell(1, icell)
1277 432 : ycell = index_to_cell(2, icell)
1278 432 : zcell = index_to_cell(3, icell)
1279 :
1280 : arg = REAL(xcell, dp)*xkp(1, ik_global) + REAL(ycell, dp)*xkp(2, ik_global) + &
1281 432 : REAL(zcell, dp)*xkp(3, ik_global)
1282 432 : coskl = wkp(ik_global)*COS(twopi*arg)
1283 432 : sinkl = wkp(ik_global)*SIN(twopi*arg)
1284 :
1285 432 : CALL dbcsr_add(mat_dm_global_work(icell)%matrix, mat_work_re, 1.0_dp, coskl)
1286 : ! The real-space convention is Re(exp(i k R) G(k)).
1287 576 : CALL dbcsr_add(mat_dm_global_work(icell)%matrix, mat_work_im, 1.0_dp, -sinkl)
1288 :
1289 : END DO
1290 :
1291 : END DO
1292 :
1293 72 : CALL dbcsr_release_p(mat_work_re)
1294 72 : CALL dbcsr_release_p(mat_work_im)
1295 72 : CALL cp_fm_release(fm_mat_re)
1296 72 : CALL cp_fm_release(fm_mat_im)
1297 :
1298 72 : CALL timestop(handle)
1299 :
1300 72 : END SUBROUTINE density_matrix_from_kp_to_transl
1301 :
1302 : ! **************************************************************************************************
1303 : !> \brief ...
1304 : !> \param cell_grid ...
1305 : !> \param cell_to_index ...
1306 : !> \param index_to_cell ...
1307 : !> \param cell ...
1308 : ! **************************************************************************************************
1309 330 : SUBROUTINE init_cell_index_rpa(cell_grid, cell_to_index, index_to_cell, cell)
1310 : INTEGER, DIMENSION(3), INTENT(IN) :: cell_grid
1311 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1312 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
1313 : TYPE(cell_type), INTENT(IN), POINTER :: cell
1314 :
1315 : CHARACTER(LEN=*), PARAMETER :: routineN = 'init_cell_index_rpa'
1316 :
1317 : INTEGER :: cell_counter, handle, i_cell, &
1318 : index_min_dist, num_cells, xcell, &
1319 : ycell, zcell
1320 : INTEGER, DIMENSION(3) :: itm
1321 330 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_unsorted
1322 330 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_unsorted
1323 330 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: abs_cell_vectors
1324 : REAL(KIND=dp), DIMENSION(3) :: cell_vector
1325 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1326 :
1327 330 : CALL timeset(routineN, handle)
1328 :
1329 330 : CALL get_cell(cell=cell, h=hmat)
1330 :
1331 330 : num_cells = cell_grid(1)*cell_grid(2)*cell_grid(3)
1332 1320 : itm(:) = cell_grid(:)/2
1333 :
1334 : ! check that real space super lattice is a (2n+1)x(2m+1)x(2k+1) super lattice with the unit cell
1335 : ! in the middle
1336 330 : CPASSERT(cell_grid(1) /= itm(1)*2)
1337 330 : CPASSERT(cell_grid(2) /= itm(2)*2)
1338 330 : CPASSERT(cell_grid(3) /= itm(3)*2)
1339 :
1340 330 : IF (ASSOCIATED(cell_to_index)) DEALLOCATE (cell_to_index)
1341 330 : IF (ASSOCIATED(index_to_cell)) DEALLOCATE (index_to_cell)
1342 :
1343 1650 : ALLOCATE (cell_to_index_unsorted(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
1344 6036 : cell_to_index_unsorted(:, :, :) = 0
1345 :
1346 990 : ALLOCATE (index_to_cell_unsorted(3, num_cells))
1347 10482 : index_to_cell_unsorted(:, :) = 0
1348 :
1349 1320 : ALLOCATE (cell_to_index(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
1350 6036 : cell_to_index(:, :, :) = 0
1351 :
1352 660 : ALLOCATE (index_to_cell(3, num_cells))
1353 10482 : index_to_cell(:, :) = 0
1354 :
1355 990 : ALLOCATE (abs_cell_vectors(1:num_cells))
1356 :
1357 780 : cell_counter = 0
1358 :
1359 780 : DO xcell = -itm(1), itm(1)
1360 1626 : DO ycell = -itm(2), itm(2)
1361 3834 : DO zcell = -itm(3), itm(3)
1362 :
1363 2538 : cell_counter = cell_counter + 1
1364 2538 : cell_to_index_unsorted(xcell, ycell, zcell) = cell_counter
1365 :
1366 2538 : index_to_cell_unsorted(1, cell_counter) = xcell
1367 2538 : index_to_cell_unsorted(2, cell_counter) = ycell
1368 2538 : index_to_cell_unsorted(3, cell_counter) = zcell
1369 :
1370 40608 : cell_vector(1:3) = MATMUL(hmat, REAL(index_to_cell_unsorted(1:3, cell_counter), dp))
1371 :
1372 3384 : abs_cell_vectors(cell_counter) = SQRT(cell_vector(1)**2 + cell_vector(2)**2 + cell_vector(3)**2)
1373 :
1374 : END DO
1375 : END DO
1376 : END DO
1377 :
1378 : ! first only do all symmetry non-equivalent cells, we need that because chi^T is computed for
1379 : ! cell indices T from index_to_cell(:,1:num_cells/2+1)
1380 1764 : DO i_cell = 1, num_cells/2 + 1
1381 :
1382 8172 : index_min_dist = MINLOC(abs_cell_vectors(1:num_cells/2 + 1), DIM=1)
1383 :
1384 1434 : xcell = index_to_cell_unsorted(1, index_min_dist)
1385 1434 : ycell = index_to_cell_unsorted(2, index_min_dist)
1386 1434 : zcell = index_to_cell_unsorted(3, index_min_dist)
1387 :
1388 1434 : index_to_cell(1, i_cell) = xcell
1389 1434 : index_to_cell(2, i_cell) = ycell
1390 1434 : index_to_cell(3, i_cell) = zcell
1391 :
1392 1434 : cell_to_index(xcell, ycell, zcell) = i_cell
1393 :
1394 1764 : abs_cell_vectors(index_min_dist) = 10000000000.0_dp
1395 :
1396 : END DO
1397 :
1398 : ! now all the remaining cells
1399 1434 : DO i_cell = num_cells/2 + 2, num_cells
1400 :
1401 10608 : index_min_dist = MINLOC(abs_cell_vectors(1:num_cells), DIM=1)
1402 :
1403 1104 : xcell = index_to_cell_unsorted(1, index_min_dist)
1404 1104 : ycell = index_to_cell_unsorted(2, index_min_dist)
1405 1104 : zcell = index_to_cell_unsorted(3, index_min_dist)
1406 :
1407 1104 : index_to_cell(1, i_cell) = xcell
1408 1104 : index_to_cell(2, i_cell) = ycell
1409 1104 : index_to_cell(3, i_cell) = zcell
1410 :
1411 1104 : cell_to_index(xcell, ycell, zcell) = i_cell
1412 :
1413 1434 : abs_cell_vectors(index_min_dist) = 10000000000.0_dp
1414 :
1415 : END DO
1416 :
1417 330 : DEALLOCATE (index_to_cell_unsorted, cell_to_index_unsorted, abs_cell_vectors)
1418 :
1419 330 : CALL timestop(handle)
1420 :
1421 330 : END SUBROUTINE init_cell_index_rpa
1422 :
1423 : ! **************************************************************************************************
1424 : !> \brief ...
1425 : !> \param i_cell_R ...
1426 : !> \param i_cell_S ...
1427 : !> \param i_cell_R_minus_S ...
1428 : !> \param index_to_cell_3c ...
1429 : !> \param cell_to_index_3c ...
1430 : !> \param index_to_cell_dm ...
1431 : !> \param R_minus_S_needed ...
1432 : !> \param do_kpoints_cubic_RPA ...
1433 : ! **************************************************************************************************
1434 13880 : SUBROUTINE get_diff_index_3c(i_cell_R, i_cell_S, i_cell_R_minus_S, index_to_cell_3c, &
1435 : cell_to_index_3c, index_to_cell_dm, R_minus_S_needed, &
1436 : do_kpoints_cubic_RPA)
1437 :
1438 : INTEGER, INTENT(IN) :: i_cell_R, i_cell_S
1439 : INTEGER, INTENT(OUT) :: i_cell_R_minus_S
1440 : INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(IN) :: index_to_cell_3c
1441 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
1442 : INTENT(IN) :: cell_to_index_3c
1443 : INTEGER, DIMENSION(:, :), INTENT(IN), POINTER :: index_to_cell_dm
1444 : LOGICAL, INTENT(OUT) :: R_minus_S_needed
1445 : LOGICAL, INTENT(IN) :: do_kpoints_cubic_RPA
1446 :
1447 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_diff_index_3c'
1448 :
1449 : INTEGER :: handle, x_cell_R, x_cell_R_minus_S, x_cell_S, y_cell_R, y_cell_R_minus_S, &
1450 : y_cell_S, z_cell_R, z_cell_R_minus_S, z_cell_S
1451 :
1452 13880 : CALL timeset(routineN, handle)
1453 :
1454 13880 : IF (do_kpoints_cubic_RPA) THEN
1455 :
1456 5040 : x_cell_R = index_to_cell_3c(1, i_cell_R)
1457 5040 : y_cell_R = index_to_cell_3c(2, i_cell_R)
1458 5040 : z_cell_R = index_to_cell_3c(3, i_cell_R)
1459 :
1460 5040 : x_cell_S = index_to_cell_dm(1, i_cell_S)
1461 5040 : y_cell_S = index_to_cell_dm(2, i_cell_S)
1462 5040 : z_cell_S = index_to_cell_dm(3, i_cell_S)
1463 :
1464 5040 : x_cell_R_minus_S = x_cell_R - x_cell_S
1465 5040 : y_cell_R_minus_S = y_cell_R - y_cell_S
1466 5040 : z_cell_R_minus_S = z_cell_R - z_cell_S
1467 :
1468 : IF (x_cell_R_minus_S >= LBOUND(cell_to_index_3c, 1) .AND. &
1469 : x_cell_R_minus_S <= UBOUND(cell_to_index_3c, 1) .AND. &
1470 : y_cell_R_minus_S >= LBOUND(cell_to_index_3c, 2) .AND. &
1471 : y_cell_R_minus_S <= UBOUND(cell_to_index_3c, 2) .AND. &
1472 35160 : z_cell_R_minus_S >= LBOUND(cell_to_index_3c, 3) .AND. &
1473 : z_cell_R_minus_S <= UBOUND(cell_to_index_3c, 3)) THEN
1474 :
1475 4740 : i_cell_R_minus_S = cell_to_index_3c(x_cell_R_minus_S, y_cell_R_minus_S, z_cell_R_minus_S)
1476 :
1477 : ! 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
1478 4740 : IF (i_cell_R_minus_S == 0) THEN
1479 :
1480 0 : R_minus_S_needed = .FALSE.
1481 0 : i_cell_R_minus_S = 0
1482 :
1483 : ELSE
1484 :
1485 4740 : R_minus_S_needed = .TRUE.
1486 :
1487 : END IF
1488 :
1489 : ELSE
1490 :
1491 300 : i_cell_R_minus_S = 0
1492 300 : R_minus_S_needed = .FALSE.
1493 :
1494 : END IF
1495 :
1496 : ELSE ! no k-points
1497 :
1498 8840 : R_minus_S_needed = .TRUE.
1499 8840 : i_cell_R_minus_S = 1
1500 :
1501 : END IF
1502 :
1503 13880 : CALL timestop(handle)
1504 :
1505 13880 : END SUBROUTINE get_diff_index_3c
1506 :
1507 : ! **************************************************************************************************
1508 : !> \brief ...
1509 : !> \param i_cell_R ...
1510 : !> \param i_cell_S ...
1511 : !> \param i_cell_T ...
1512 : !> \param i_cell_R_minus_S_minus_T ...
1513 : !> \param index_to_cell_3c ...
1514 : !> \param cell_to_index_3c ...
1515 : !> \param index_to_cell_dm ...
1516 : !> \param R_minus_S_minus_T_needed ...
1517 : !> \param do_kpoints_cubic_RPA ...
1518 : ! **************************************************************************************************
1519 16400 : SUBROUTINE get_diff_diff_index_3c(i_cell_R, i_cell_S, i_cell_T, i_cell_R_minus_S_minus_T, &
1520 8200 : index_to_cell_3c, cell_to_index_3c, index_to_cell_dm, &
1521 : R_minus_S_minus_T_needed, &
1522 : do_kpoints_cubic_RPA)
1523 :
1524 : INTEGER, INTENT(IN) :: i_cell_R, i_cell_S, i_cell_T
1525 : INTEGER, INTENT(OUT) :: i_cell_R_minus_S_minus_T
1526 : INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell_3c
1527 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
1528 : INTENT(IN) :: cell_to_index_3c
1529 : INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell_dm
1530 : LOGICAL, INTENT(OUT) :: R_minus_S_minus_T_needed
1531 : LOGICAL, INTENT(IN) :: do_kpoints_cubic_RPA
1532 :
1533 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_diff_diff_index_3c'
1534 :
1535 : INTEGER :: handle, x_cell_R, x_cell_R_minus_S_minus_T, x_cell_S, x_cell_T, y_cell_R, &
1536 : y_cell_R_minus_S_minus_T, y_cell_S, y_cell_T, z_cell_R, z_cell_R_minus_S_minus_T, &
1537 : z_cell_S, z_cell_T
1538 :
1539 8200 : CALL timeset(routineN, handle)
1540 :
1541 8200 : IF (do_kpoints_cubic_RPA) THEN
1542 :
1543 3780 : x_cell_R = index_to_cell_3c(1, i_cell_R)
1544 3780 : y_cell_R = index_to_cell_3c(2, i_cell_R)
1545 3780 : z_cell_R = index_to_cell_3c(3, i_cell_R)
1546 :
1547 3780 : x_cell_S = index_to_cell_dm(1, i_cell_S)
1548 3780 : y_cell_S = index_to_cell_dm(2, i_cell_S)
1549 3780 : z_cell_S = index_to_cell_dm(3, i_cell_S)
1550 :
1551 3780 : x_cell_T = index_to_cell_dm(1, i_cell_T)
1552 3780 : y_cell_T = index_to_cell_dm(2, i_cell_T)
1553 3780 : z_cell_T = index_to_cell_dm(3, i_cell_T)
1554 :
1555 3780 : x_cell_R_minus_S_minus_T = x_cell_R - x_cell_S - x_cell_T
1556 3780 : y_cell_R_minus_S_minus_T = y_cell_R - y_cell_S - y_cell_T
1557 3780 : z_cell_R_minus_S_minus_T = z_cell_R - z_cell_S - z_cell_T
1558 :
1559 : IF (x_cell_R_minus_S_minus_T >= LBOUND(cell_to_index_3c, 1) .AND. &
1560 : x_cell_R_minus_S_minus_T <= UBOUND(cell_to_index_3c, 1) .AND. &
1561 : y_cell_R_minus_S_minus_T >= LBOUND(cell_to_index_3c, 2) .AND. &
1562 : y_cell_R_minus_S_minus_T <= UBOUND(cell_to_index_3c, 2) .AND. &
1563 26400 : z_cell_R_minus_S_minus_T >= LBOUND(cell_to_index_3c, 3) .AND. &
1564 : z_cell_R_minus_S_minus_T <= UBOUND(cell_to_index_3c, 3)) THEN
1565 :
1566 : i_cell_R_minus_S_minus_T = cell_to_index_3c(x_cell_R_minus_S_minus_T, &
1567 : y_cell_R_minus_S_minus_T, &
1568 3480 : z_cell_R_minus_S_minus_T)
1569 :
1570 : ! index 0 means that there are only no 3c matrix elements because R-S-T is too big
1571 3480 : IF (i_cell_R_minus_S_minus_T == 0) THEN
1572 :
1573 0 : R_minus_S_minus_T_needed = .FALSE.
1574 :
1575 : ELSE
1576 :
1577 3480 : R_minus_S_minus_T_needed = .TRUE.
1578 :
1579 : END IF
1580 :
1581 : ELSE
1582 :
1583 300 : i_cell_R_minus_S_minus_T = 0
1584 300 : R_minus_S_minus_T_needed = .FALSE.
1585 :
1586 : END IF
1587 :
1588 : ! no k-kpoints
1589 : ELSE
1590 :
1591 4420 : R_minus_S_minus_T_needed = .TRUE.
1592 4420 : i_cell_R_minus_S_minus_T = 1
1593 :
1594 : END IF
1595 :
1596 8200 : CALL timestop(handle)
1597 :
1598 8200 : END SUBROUTINE get_diff_diff_index_3c
1599 :
1600 1482 : END MODULE rpa_im_time
|
