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