(git:b279b6b)
rpa_gw_kpoints_util.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 treating GW and RPA calculations with kpoints
10 !> \par History
11 !> since 2018 continuous development [J. Wilhelm]
12 ! **************************************************************************************************
13 MODULE rpa_gw_kpoints_util
14  USE cell_types, ONLY: cell_type,&
15  get_cell,&
16  pbc
17  USE cp_blacs_env, ONLY: cp_blacs_env_type
18  USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_decompose,&
19  cp_cfm_cholesky_invert,&
20  cp_cfm_column_scale,&
21  cp_cfm_scale_and_add,&
22  cp_cfm_scale_and_add_fm,&
23  cp_cfm_transpose
24  USE cp_cfm_diag, ONLY: cp_cfm_geeig,&
25  cp_cfm_geeig_canon,&
26  cp_cfm_heevd
27  USE cp_cfm_types, ONLY: cp_cfm_create,&
28  cp_cfm_get_info,&
29  cp_cfm_release,&
30  cp_cfm_set_all,&
31  cp_cfm_to_cfm,&
32  cp_cfm_to_fm,&
33  cp_cfm_type
34  USE cp_control_types, ONLY: dft_control_type
35  USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
36  copy_fm_to_dbcsr,&
37  dbcsr_allocate_matrix_set
38  USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
39  USE cp_fm_struct, ONLY: cp_fm_struct_type
40  USE cp_fm_types, ONLY: cp_fm_copy_general,&
41  cp_fm_create,&
42  cp_fm_release,&
43  cp_fm_set_all,&
44  cp_fm_type
45  USE dbcsr_api, ONLY: &
46  dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_filter, &
47  dbcsr_get_block_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
48  dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
49  dbcsr_release, dbcsr_reserve_all_blocks, dbcsr_set, dbcsr_transposed, dbcsr_type, &
50  dbcsr_type_no_symmetry
51  USE hfx_types, ONLY: hfx_release
52  USE input_constants, ONLY: cholesky_off,&
53  kp_weights_w_auto,&
54  kp_weights_w_tailored,&
55  kp_weights_w_uniform
56  USE kinds, ONLY: dp
57  USE kpoint_methods, ONLY: kpoint_env_initialize,&
58  kpoint_initialize_mo_set,&
59  kpoint_initialize_mos
60  USE kpoint_types, ONLY: get_kpoint_info,&
61  kpoint_env_type,&
62  kpoint_type
63  USE machine, ONLY: m_walltime
64  USE mathconstants, ONLY: gaussi,&
65  twopi,&
66  z_one,&
67  z_zero
68  USE mathlib, ONLY: invmat
69  USE message_passing, ONLY: mp_para_env_type
70  USE parallel_gemm_api, ONLY: parallel_gemm
71  USE particle_types, ONLY: particle_type
72  USE qs_band_structure, ONLY: calculate_kpoints_for_bs
73  USE qs_environment_types, ONLY: get_qs_env,&
74  qs_environment_type
75  USE qs_mo_types, ONLY: get_mo_set
76  USE qs_scf_types, ONLY: qs_scf_env_type
77  USE rpa_gw_im_time_util, ONLY: compute_weight_re_im,&
78  get_atom_index_from_basis_function_index
79  USE rpa_im_time, ONLY: init_cell_index_rpa
80  USE scf_control_types, ONLY: scf_control_type
81 #include "./base/base_uses.f90"
82 
83  IMPLICIT NONE
84 
85  PRIVATE
86 
87  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_gw_kpoints_util'
88 
89  PUBLIC :: invert_eps_compute_w_and_erpa_kp, cp_cfm_power, real_space_to_kpoint_transform_rpa, &
90  get_mat_cell_t_from_mat_gamma, get_bandstruc_and_k_dependent_mos, &
91  compute_wkp_w, mat_kp_from_mat_gamma, cp_cfm_upper_to_full
92 
93 CONTAINS
94 
95 ! **************************************************************************************************
96 !> \brief ...
97 !> \param dimen_RI ...
98 !> \param num_integ_points ...
99 !> \param jquad ...
100 !> \param nkp ...
101 !> \param count_ev_sc_GW ...
102 !> \param para_env ...
103 !> \param Erpa ...
104 !> \param tau_tj ...
105 !> \param tj ...
106 !> \param wj ...
107 !> \param weights_cos_tf_w_to_t ...
108 !> \param wkp_W ...
109 !> \param do_gw_im_time ...
110 !> \param do_ri_Sigma_x ...
111 !> \param do_kpoints_from_Gamma ...
112 !> \param cfm_mat_Q ...
113 !> \param ikp_local ...
114 !> \param mat_P_omega ...
115 !> \param mat_P_omega_kp ...
116 !> \param qs_env ...
117 !> \param eps_filter_im_time ...
118 !> \param unit_nr ...
119 !> \param kpoints ...
120 !> \param fm_mat_Minv_L_kpoints ...
121 !> \param fm_matrix_L_kpoints ...
122 !> \param fm_mat_W ...
123 !> \param fm_mat_RI_global_work ...
124 !> \param mat_MinvVMinv ...
125 !> \param fm_matrix_Minv ...
126 !> \param fm_matrix_Minv_Vtrunc_Minv ...
127 ! **************************************************************************************************
128  SUBROUTINE invert_eps_compute_w_and_erpa_kp(dimen_RI, num_integ_points, jquad, nkp, count_ev_sc_GW, para_env, &
129  Erpa, tau_tj, tj, wj, weights_cos_tf_w_to_t, wkp_W, do_gw_im_time, &
130  do_ri_Sigma_x, do_kpoints_from_Gamma, &
131  cfm_mat_Q, ikp_local, mat_P_omega, mat_P_omega_kp, &
132  qs_env, eps_filter_im_time, unit_nr, kpoints, fm_mat_Minv_L_kpoints, &
133  fm_matrix_L_kpoints, fm_mat_W, &
134  fm_mat_RI_global_work, mat_MinvVMinv, fm_matrix_Minv, &
135  fm_matrix_Minv_Vtrunc_Minv)
136 
137  INTEGER, INTENT(IN) :: dimen_ri, num_integ_points, jquad, nkp, &
138  count_ev_sc_gw
139  TYPE(mp_para_env_type), POINTER :: para_env
140  REAL(kind=dp), INTENT(INOUT) :: erpa
141  REAL(kind=dp), DIMENSION(0:num_integ_points), &
142  INTENT(IN) :: tau_tj
143  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: tj, wj
144  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
145  INTENT(IN) :: weights_cos_tf_w_to_t
146  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wkp_w
147  LOGICAL, INTENT(IN) :: do_gw_im_time, do_ri_sigma_x, &
148  do_kpoints_from_gamma
149  TYPE(cp_cfm_type), INTENT(IN) :: cfm_mat_q
150  INTEGER, INTENT(IN) :: ikp_local
151  TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: mat_p_omega, mat_p_omega_kp
152  TYPE(qs_environment_type), POINTER :: qs_env
153  REAL(kind=dp), INTENT(IN) :: eps_filter_im_time
154  INTEGER, INTENT(IN) :: unit_nr
155  TYPE(kpoint_type), POINTER :: kpoints
156  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_mat_minv_l_kpoints, &
157  fm_matrix_l_kpoints
158  TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_w
159  TYPE(cp_fm_type) :: fm_mat_ri_global_work
160  TYPE(dbcsr_p_type), INTENT(IN) :: mat_minvvminv
161  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_minv, &
162  fm_matrix_minv_vtrunc_minv
163 
164  CHARACTER(LEN=*), PARAMETER :: routinen = 'invert_eps_compute_W_and_Erpa_kp'
165 
166  INTEGER :: handle, ikp
167  LOGICAL :: do_this_ikp
168  REAL(kind=dp) :: t1, t2
169  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: trace_qomega
170 
171  CALL timeset(routinen, handle)
172 
173  t1 = m_walltime()
174 
175  IF (do_kpoints_from_gamma) THEN
176  CALL get_mat_cell_t_from_mat_gamma(mat_p_omega(jquad, :), qs_env, kpoints, jquad, unit_nr)
177  END IF
178 
179  CALL transform_p_from_real_space_to_kpoints(mat_p_omega, mat_p_omega_kp, &
180  kpoints, eps_filter_im_time, jquad)
181 
182  ALLOCATE (trace_qomega(dimen_ri))
183 
184  IF (unit_nr > 0) WRITE (unit_nr, '(/T3,A,1X,I3)') &
185  'GW_INFO| Computing chi and W frequency point', jquad
186 
187  DO ikp = 1, nkp
188 
189  ! parallization, we either have all kpoints on all processors or a single kpoint per group
190  do_this_ikp = (ikp_local == -1) .OR. (ikp_local == 0 .AND. ikp == 1) .OR. (ikp_local == ikp)
191  IF (.NOT. do_this_ikp) cycle
192 
193  ! 1. remove all spurious negative eigenvalues from P(iw,k), multiplication Q(iw,k) = K^H(k)P(iw,k)K(k)
194  CALL compute_q_kp_rpa(cfm_mat_q, &
195  mat_p_omega_kp, &
196  fm_mat_minv_l_kpoints(ikp, 1), &
197  fm_mat_minv_l_kpoints(ikp, 2), &
198  fm_mat_ri_global_work, &
199  dimen_ri, ikp, nkp, ikp_local, para_env, &
200  qs_env%mp2_env%ri_rpa_im_time%make_chi_pos_definite)
201 
202  ! 2. Cholesky decomposition of Id + Q(iw,k)
203  CALL cholesky_decomp_q(cfm_mat_q, para_env, trace_qomega, dimen_ri)
204 
205  ! 3. Computing E_c^RPA = E_c^RPA + a_w/N_k*sum_k ln[det(1+Q(iw,k))-Tr(Q(iw,k))]
206  CALL frequency_and_kpoint_integration(erpa, cfm_mat_q, para_env, trace_qomega, &
207  dimen_ri, wj(jquad), kpoints%wkp(ikp))
208 
209  IF (do_gw_im_time) THEN
210 
211  ! compute S^-1*V*S^-1 for exchange part of the self-energy in real space as W in real space
212  IF (do_ri_sigma_x .AND. jquad == 1 .AND. count_ev_sc_gw == 1 .AND. do_kpoints_from_gamma) THEN
213 
214  CALL dbcsr_set(mat_minvvminv%matrix, 0.0_dp)
215  CALL copy_fm_to_dbcsr(fm_matrix_minv_vtrunc_minv(1, 1), mat_minvvminv%matrix, keep_sparsity=.false.)
216 
217  END IF
218  IF (do_kpoints_from_gamma) THEN
219  CALL compute_wc_real_space_tau_gw(fm_mat_w, cfm_mat_q, &
220  fm_matrix_l_kpoints(ikp, 1), &
221  fm_matrix_l_kpoints(ikp, 2), &
222  dimen_ri, num_integ_points, jquad, &
223  ikp, tj, tau_tj, weights_cos_tf_w_to_t, &
224  ikp_local, para_env, kpoints, qs_env, wkp_w)
225  END IF
226 
227  END IF
228  END DO
229 
230  ! after the transform of (eps(iw)-1)^-1 from iw to it is done, multiply with V^1/2 to obtain W(it)
231  IF (do_gw_im_time .AND. do_kpoints_from_gamma .AND. jquad == num_integ_points) THEN
232  CALL wc_to_minv_wc_minv(fm_mat_w, fm_matrix_minv, para_env, dimen_ri, num_integ_points)
233  CALL deallocate_kp_matrices(fm_matrix_l_kpoints, fm_mat_minv_l_kpoints)
234  END IF
235 
236  DEALLOCATE (trace_qomega)
237 
238  t2 = m_walltime()
239 
240  IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,T56,F25.1)') 'Execution time (s):', t2 - t1
241 
242  CALL timestop(handle)
243 
244  END SUBROUTINE invert_eps_compute_w_and_erpa_kp
245 
246 ! **************************************************************************************************
247 !> \brief ...
248 !> \param fm_matrix_L_kpoints ...
249 !> \param fm_mat_Minv_L_kpoints ...
250 ! **************************************************************************************************
251  SUBROUTINE deallocate_kp_matrices(fm_matrix_L_kpoints, fm_mat_Minv_L_kpoints)
252 
253  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_l_kpoints, &
254  fm_mat_minv_l_kpoints
255 
256  CHARACTER(LEN=*), PARAMETER :: routinen = 'deallocate_kp_matrices'
257 
258  INTEGER :: handle
259 
260  CALL timeset(routinen, handle)
261 
262  CALL cp_fm_release(fm_mat_minv_l_kpoints)
263  CALL cp_fm_release(fm_matrix_l_kpoints)
264 
265  CALL timestop(handle)
266 
267  END SUBROUTINE deallocate_kp_matrices
268 
269 ! **************************************************************************************************
270 !> \brief ...
271 !> \param matrix ...
272 !> \param threshold ...
273 !> \param exponent ...
274 !> \param min_eigval ...
275 ! **************************************************************************************************
276  SUBROUTINE cp_cfm_power(matrix, threshold, exponent, min_eigval)
277  TYPE(cp_cfm_type), INTENT(INOUT) :: matrix
278  REAL(kind=dp) :: threshold, exponent
279  REAL(kind=dp), OPTIONAL :: min_eigval
280 
281  CHARACTER(LEN=*), PARAMETER :: routinen = 'cp_cfm_power'
282  COMPLEX(KIND=dp), PARAMETER :: czero = cmplx(0.0_dp, 0.0_dp, kind=dp)
283 
284  COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_exponent
285  INTEGER :: handle, i, ncol_global, nrow_global
286  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
287  TYPE(cp_cfm_type) :: cfm_work
288 
289  CALL timeset(routinen, handle)
290 
291  CALL cp_cfm_create(cfm_work, matrix%matrix_struct)
292  CALL cp_cfm_set_all(cfm_work, z_zero)
293 
294  ! Test that matrix is square
295  CALL cp_cfm_get_info(matrix, nrow_global=nrow_global, ncol_global=ncol_global)
296  cpassert(nrow_global == ncol_global)
297  ALLOCATE (eigenvalues(nrow_global))
298  eigenvalues(:) = 0.0_dp
299  ALLOCATE (eigenvalues_exponent(nrow_global))
300  eigenvalues_exponent(:) = czero
301 
302  ! Diagonalize matrix: get eigenvectors and eigenvalues
303  CALL cp_cfm_heevd(matrix, cfm_work, eigenvalues)
304 
305  DO i = 1, nrow_global
306  IF (eigenvalues(i) > threshold) THEN
307  eigenvalues_exponent(i) = cmplx((eigenvalues(i))**(0.5_dp*exponent), threshold, kind=dp)
308  ELSE
309  IF (PRESENT(min_eigval)) THEN
310  eigenvalues_exponent(i) = cmplx(min_eigval, 0.0_dp, kind=dp)
311  ELSE
312  eigenvalues_exponent(i) = czero
313  END IF
314  END IF
315  END DO
316 
317  CALL cp_cfm_column_scale(cfm_work, eigenvalues_exponent)
318 
319  CALL parallel_gemm("N", "C", nrow_global, nrow_global, nrow_global, z_one, &
320  cfm_work, cfm_work, z_zero, matrix)
321 
322  DEALLOCATE (eigenvalues, eigenvalues_exponent)
323 
324  CALL cp_cfm_release(cfm_work)
325 
326  CALL timestop(handle)
327 
328  END SUBROUTINE cp_cfm_power
329 
330 ! **************************************************************************************************
331 !> \brief ...
332 !> \param cfm_mat_Q ...
333 !> \param mat_P_omega_kp ...
334 !> \param fm_mat_L_re ...
335 !> \param fm_mat_L_im ...
336 !> \param fm_mat_RI_global_work ...
337 !> \param dimen_RI ...
338 !> \param ikp ...
339 !> \param nkp ...
340 !> \param ikp_local ...
341 !> \param para_env ...
342 !> \param make_chi_pos_definite ...
343 ! **************************************************************************************************
344  SUBROUTINE compute_q_kp_rpa(cfm_mat_Q, mat_P_omega_kp, fm_mat_L_re, fm_mat_L_im, &
345  fm_mat_RI_global_work, dimen_RI, ikp, nkp, ikp_local, para_env, &
346  make_chi_pos_definite)
347 
348  TYPE(cp_cfm_type) :: cfm_mat_q
349  TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: mat_p_omega_kp
350  TYPE(cp_fm_type) :: fm_mat_l_re, fm_mat_l_im, &
351  fm_mat_ri_global_work
352  INTEGER, INTENT(IN) :: dimen_ri, ikp, nkp, ikp_local
353  TYPE(mp_para_env_type), POINTER :: para_env
354  LOGICAL, INTENT(IN) :: make_chi_pos_definite
355 
356  CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_Q_kp_RPA'
357 
358  INTEGER :: handle
359  TYPE(cp_cfm_type) :: cfm_mat_l, cfm_mat_work
360  TYPE(cp_fm_type) :: fm_mat_work
361 
362  CALL timeset(routinen, handle)
363 
364  CALL cp_cfm_create(cfm_mat_work, fm_mat_l_re%matrix_struct)
365  CALL cp_cfm_set_all(cfm_mat_work, z_zero)
366 
367  CALL cp_cfm_create(cfm_mat_l, fm_mat_l_re%matrix_struct)
368  CALL cp_cfm_set_all(cfm_mat_l, z_zero)
369 
370  CALL cp_fm_create(fm_mat_work, fm_mat_l_re%matrix_struct)
371  CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
372 
373  ! 1. Convert the dbcsr matrix mat_P_omega_kp (that is chi(k,iw)) to a full matrix and
374  ! distribute it to subgroups
375  CALL mat_p_to_subgroup(mat_p_omega_kp, fm_mat_ri_global_work, &
376  fm_mat_work, cfm_mat_q, ikp, nkp, ikp_local, para_env)
377 
378  ! 2. Remove all negative eigenvalues from chi(k,iw)
379  IF (make_chi_pos_definite) THEN
380  CALL cp_cfm_power(cfm_mat_q, threshold=0.0_dp, exponent=1.0_dp)
381  END IF
382 
383  ! 3. Copy fm_mat_L_re and fm_mat_L_re to cfm_mat_L
384  CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mat_l, z_one, fm_mat_l_re)
385  CALL cp_cfm_scale_and_add_fm(z_one, cfm_mat_l, gaussi, fm_mat_l_im)
386 
387  ! 4. work = P(iw,k)*L(k)
388  CALL parallel_gemm('N', 'N', dimen_ri, dimen_ri, dimen_ri, z_one, cfm_mat_q, cfm_mat_l, &
389  z_zero, cfm_mat_work)
390 
391  ! 5. Q(iw,k) = L^H(k)*work
392  CALL parallel_gemm('C', 'N', dimen_ri, dimen_ri, dimen_ri, z_one, cfm_mat_l, cfm_mat_work, &
393  z_zero, cfm_mat_q)
394 
395  CALL cp_cfm_release(cfm_mat_work)
396  CALL cp_cfm_release(cfm_mat_l)
397  CALL cp_fm_release(fm_mat_work)
398 
399  CALL timestop(handle)
400 
401  END SUBROUTINE compute_q_kp_rpa
402 
403 ! **************************************************************************************************
404 !> \brief ...
405 !> \param mat_P_omega_kp ...
406 !> \param fm_mat_RI_global_work ...
407 !> \param fm_mat_work ...
408 !> \param cfm_mat_Q ...
409 !> \param ikp ...
410 !> \param nkp ...
411 !> \param ikp_local ...
412 !> \param para_env ...
413 ! **************************************************************************************************
414  SUBROUTINE mat_p_to_subgroup(mat_P_omega_kp, fm_mat_RI_global_work, &
415  fm_mat_work, cfm_mat_Q, ikp, nkp, ikp_local, para_env)
416 
417  TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: mat_p_omega_kp
418  TYPE(cp_fm_type), INTENT(IN) :: fm_mat_ri_global_work, fm_mat_work
419  TYPE(cp_cfm_type), INTENT(IN) :: cfm_mat_q
420  INTEGER, INTENT(IN) :: ikp, nkp, ikp_local
421  TYPE(mp_para_env_type), POINTER :: para_env
422 
423  CHARACTER(LEN=*), PARAMETER :: routinen = 'mat_P_to_subgroup'
424 
425  INTEGER :: handle, jkp
426  TYPE(cp_fm_type) :: fm_dummy
427  TYPE(dbcsr_type), POINTER :: mat_p_omega_im, mat_p_omega_re
428 
429  CALL timeset(routinen, handle)
430 
431  IF (ikp_local == -1) THEN
432 
433  mat_p_omega_re => mat_p_omega_kp(1, ikp)%matrix
434  CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
435  CALL copy_dbcsr_to_fm(mat_p_omega_re, fm_mat_work)
436  CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mat_q, z_one, fm_mat_work)
437 
438  mat_p_omega_im => mat_p_omega_kp(2, ikp)%matrix
439  CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
440  CALL copy_dbcsr_to_fm(mat_p_omega_im, fm_mat_work)
441  CALL cp_cfm_scale_and_add_fm(z_one, cfm_mat_q, gaussi, fm_mat_work)
442 
443  ELSE
444 
445  CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
446 
447  DO jkp = 1, nkp
448 
449  mat_p_omega_re => mat_p_omega_kp(1, jkp)%matrix
450 
451  CALL cp_fm_set_all(fm_mat_ri_global_work, 0.0_dp)
452  CALL copy_dbcsr_to_fm(mat_p_omega_re, fm_mat_ri_global_work)
453 
454  CALL para_env%sync()
455 
456  IF (ikp_local == jkp) THEN
457  CALL cp_fm_copy_general(fm_mat_ri_global_work, fm_mat_work, para_env)
458  ELSE
459  CALL cp_fm_copy_general(fm_mat_ri_global_work, fm_dummy, para_env)
460  END IF
461 
462  CALL para_env%sync()
463 
464  END DO
465 
466  CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mat_q, z_one, fm_mat_work)
467 
468  CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
469 
470  DO jkp = 1, nkp
471 
472  mat_p_omega_im => mat_p_omega_kp(2, jkp)%matrix
473 
474  CALL cp_fm_set_all(fm_mat_ri_global_work, 0.0_dp)
475  CALL copy_dbcsr_to_fm(mat_p_omega_im, fm_mat_ri_global_work)
476 
477  CALL para_env%sync()
478 
479  IF (ikp_local == jkp) THEN
480  CALL cp_fm_copy_general(fm_mat_ri_global_work, fm_mat_work, para_env)
481  ELSE
482  CALL cp_fm_copy_general(fm_mat_ri_global_work, fm_dummy, para_env)
483  END IF
484 
485  CALL para_env%sync()
486 
487  END DO
488 
489  CALL cp_cfm_scale_and_add_fm(z_one, cfm_mat_q, gaussi, fm_mat_work)
490 
491  CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
492 
493  END IF
494 
495  CALL para_env%sync()
496 
497  CALL timestop(handle)
498 
499  END SUBROUTINE mat_p_to_subgroup
500 
501 ! **************************************************************************************************
502 !> \brief ...
503 !> \param cfm_mat_Q ...
504 !> \param para_env ...
505 !> \param trace_Qomega ...
506 !> \param dimen_RI ...
507 ! **************************************************************************************************
508  SUBROUTINE cholesky_decomp_q(cfm_mat_Q, para_env, trace_Qomega, dimen_RI)
509 
510  TYPE(cp_cfm_type), INTENT(IN) :: cfm_mat_q
511  TYPE(mp_para_env_type), INTENT(IN) :: para_env
512  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: trace_qomega
513  INTEGER, INTENT(IN) :: dimen_ri
514 
515  CHARACTER(LEN=*), PARAMETER :: routinen = 'cholesky_decomp_Q'
516 
517  INTEGER :: handle, i_global, iib, info_chol, &
518  j_global, jjb, ncol_local, nrow_local
519  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
520  TYPE(cp_cfm_type) :: cfm_mat_q_tmp, cfm_mat_work
521 
522  CALL timeset(routinen, handle)
523 
524  CALL cp_cfm_create(cfm_mat_work, cfm_mat_q%matrix_struct)
525  CALL cp_cfm_set_all(cfm_mat_work, z_zero)
526 
527  CALL cp_cfm_create(cfm_mat_q_tmp, cfm_mat_q%matrix_struct)
528  CALL cp_cfm_set_all(cfm_mat_q_tmp, z_zero)
529 
530  ! get info of fm_mat_Q
531  CALL cp_cfm_get_info(matrix=cfm_mat_q, &
532  nrow_local=nrow_local, &
533  ncol_local=ncol_local, &
534  row_indices=row_indices, &
535  col_indices=col_indices)
536 
537  ! calculate the trace of Q and add 1 on the diagonal
538  trace_qomega = 0.0_dp
539 !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
540 !$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,trace_Qomega,cfm_mat_Q,dimen_RI)
541  DO jjb = 1, ncol_local
542  j_global = col_indices(jjb)
543  DO iib = 1, nrow_local
544  i_global = row_indices(iib)
545  IF (j_global == i_global .AND. i_global <= dimen_ri) THEN
546  trace_qomega(i_global) = real(cfm_mat_q%local_data(iib, jjb))
547  cfm_mat_q%local_data(iib, jjb) = cfm_mat_q%local_data(iib, jjb) + z_one
548  END IF
549  END DO
550  END DO
551  CALL para_env%sum(trace_qomega)
552 
553  CALL cp_cfm_to_cfm(cfm_mat_q, cfm_mat_q_tmp)
554 
555  CALL cp_cfm_cholesky_decompose(matrix=cfm_mat_q, n=dimen_ri, info_out=info_chol)
556 
557  cpassert(info_chol == 0)
558 
559  CALL cp_cfm_release(cfm_mat_work)
560  CALL cp_cfm_release(cfm_mat_q_tmp)
561 
562  CALL timestop(handle)
563 
564  END SUBROUTINE cholesky_decomp_q
565 
566 ! **************************************************************************************************
567 !> \brief ...
568 !> \param Erpa ...
569 !> \param cfm_mat_Q ...
570 !> \param para_env ...
571 !> \param trace_Qomega ...
572 !> \param dimen_RI ...
573 !> \param freq_weight ...
574 !> \param kp_weight ...
575 ! **************************************************************************************************
576  SUBROUTINE frequency_and_kpoint_integration(Erpa, cfm_mat_Q, para_env, trace_Qomega, &
577  dimen_RI, freq_weight, kp_weight)
578 
579  REAL(kind=dp), INTENT(INOUT) :: erpa
580  TYPE(cp_cfm_type), INTENT(IN) :: cfm_mat_q
581  TYPE(mp_para_env_type), INTENT(IN) :: para_env
582  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: trace_qomega
583  INTEGER, INTENT(IN) :: dimen_ri
584  REAL(kind=dp), INTENT(IN) :: freq_weight, kp_weight
585 
586  CHARACTER(LEN=*), PARAMETER :: routinen = 'frequency_and_kpoint_integration'
587 
588  INTEGER :: handle, i_global, iib, j_global, jjb, &
589  ncol_local, nrow_local
590  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
591  REAL(kind=dp) :: fcomega
592  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: q_log
593 
594  CALL timeset(routinen, handle)
595 
596  ! get info of cholesky_decomposed(fm_mat_Q)
597  CALL cp_cfm_get_info(matrix=cfm_mat_q, &
598  nrow_local=nrow_local, &
599  ncol_local=ncol_local, &
600  row_indices=row_indices, &
601  col_indices=col_indices)
602 
603  ALLOCATE (q_log(dimen_ri))
604  q_log = 0.0_dp
605 !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
606 !$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,Q_log,cfm_mat_Q,dimen_RI)
607  DO jjb = 1, ncol_local
608  j_global = col_indices(jjb)
609  DO iib = 1, nrow_local
610  i_global = row_indices(iib)
611  IF (j_global == i_global .AND. i_global <= dimen_ri) THEN
612  q_log(i_global) = 2.0_dp*log(real(cfm_mat_q%local_data(iib, jjb)))
613  END IF
614  END DO
615  END DO
616  CALL para_env%sum(q_log)
617 
618  fcomega = 0.0_dp
619  DO iib = 1, dimen_ri
620  IF (modulo(iib, para_env%num_pe) /= para_env%mepos) cycle
621  ! FComega=FComega+(LOG(Q_log(iiB))-trace_Qomega(iiB))/2.0_dp
622  fcomega = fcomega + (q_log(iib) - trace_qomega(iib))/2.0_dp
623  END DO
624 
625  erpa = erpa + fcomega*freq_weight*kp_weight
626 
627  DEALLOCATE (q_log)
628 
629  CALL timestop(handle)
630 
631  END SUBROUTINE frequency_and_kpoint_integration
632 
633 ! **************************************************************************************************
634 !> \brief ...
635 !> \param tj_dummy ...
636 !> \param tau_tj_dummy ...
637 !> \param weights_cos_tf_w_to_t_dummy ...
638 ! **************************************************************************************************
639  SUBROUTINE get_dummys(tj_dummy, tau_tj_dummy, weights_cos_tf_w_to_t_dummy)
640 
641  REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
642  INTENT(INOUT) :: tj_dummy, tau_tj_dummy
643  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
644  INTENT(INOUT) :: weights_cos_tf_w_to_t_dummy
645 
646  CHARACTER(LEN=*), PARAMETER :: routinen = 'get_dummys'
647 
648  INTEGER :: handle
649 
650  CALL timeset(routinen, handle)
651 
652  ALLOCATE (weights_cos_tf_w_to_t_dummy(1, 1))
653  ALLOCATE (tj_dummy(1))
654  ALLOCATE (tau_tj_dummy(1))
655 
656  tj_dummy(1) = 0.0_dp
657  tau_tj_dummy(1) = 0.0_dp
658  weights_cos_tf_w_to_t_dummy(1, 1) = 1.0_dp
659 
660  CALL timestop(handle)
661 
662  END SUBROUTINE
663 
664 ! **************************************************************************************************
665 !> \brief ...
666 !> \param tj_dummy ...
667 !> \param tau_tj_dummy ...
668 !> \param weights_cos_tf_w_to_t_dummy ...
669 ! **************************************************************************************************
670  SUBROUTINE release_dummys(tj_dummy, tau_tj_dummy, weights_cos_tf_w_to_t_dummy)
671 
672  REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
673  INTENT(INOUT) :: tj_dummy, tau_tj_dummy
674  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
675  INTENT(INOUT) :: weights_cos_tf_w_to_t_dummy
676 
677  CHARACTER(LEN=*), PARAMETER :: routinen = 'release_dummys'
678 
679  INTEGER :: handle
680 
681  CALL timeset(routinen, handle)
682 
683  DEALLOCATE (weights_cos_tf_w_to_t_dummy)
684  DEALLOCATE (tj_dummy)
685  DEALLOCATE (tau_tj_dummy)
686 
687  CALL timestop(handle)
688 
689  END SUBROUTINE
690 
691 ! **************************************************************************************************
692 !> \brief ...
693 !> \param mat_P_omega ...
694 !> \param qs_env ...
695 !> \param kpoints ...
696 !> \param jquad ...
697 !> \param unit_nr ...
698 ! **************************************************************************************************
699  SUBROUTINE get_mat_cell_t_from_mat_gamma(mat_P_omega, qs_env, kpoints, jquad, unit_nr)
700  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: mat_p_omega
701  TYPE(qs_environment_type), POINTER :: qs_env
702  TYPE(kpoint_type), POINTER :: kpoints
703  INTEGER, INTENT(IN) :: jquad, unit_nr
704 
705  CHARACTER(LEN=*), PARAMETER :: routinen = 'get_mat_cell_T_from_mat_gamma'
706 
707  INTEGER :: col, handle, i_cell, i_dim, j_cell, &
708  num_cells_p, num_integ_points, row
709  INTEGER, DIMENSION(3) :: cell_grid_p, periodic
710  INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_p
711  LOGICAL :: i_cell_is_the_minimum_image_cell
712  REAL(kind=dp) :: abs_rab_cell_i, abs_rab_cell_j
713  REAL(kind=dp), DIMENSION(3) :: cell_vector, cell_vector_j, rab_cell_i, &
714  rab_cell_j
715  REAL(kind=dp), DIMENSION(3, 3) :: hmat
716  REAL(kind=dp), DIMENSION(:, :), POINTER :: data_block
717  TYPE(cell_type), POINTER :: cell
718  TYPE(dbcsr_iterator_type) :: iter
719  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
720 
721  CALL timeset(routinen, handle)
722 
723  NULLIFY (cell, particle_set)
724  CALL get_qs_env(qs_env, cell=cell, &
725  particle_set=particle_set)
726  CALL get_cell(cell=cell, h=hmat, periodic=periodic)
727 
728  DO i_dim = 1, 3
729  ! we have at most 3 neigboring cells per dimension and at least one because
730  ! the density response at Gamma is only divided to neighboring
731  IF (periodic(i_dim) == 1) THEN
732  cell_grid_p(i_dim) = max(min((kpoints%nkp_grid(i_dim)/2)*2 - 1, 1), 3)
733  ELSE
734  cell_grid_p(i_dim) = 1
735  END IF
736  END DO
737 
738  ! overwrite the cell indices in kpoints
739  CALL init_cell_index_rpa(cell_grid_p, kpoints%cell_to_index, kpoints%index_to_cell, cell)
740 
741  index_to_cell_p => kpoints%index_to_cell
742 
743  num_cells_p = SIZE(index_to_cell_p, 2)
744 
745  num_integ_points = SIZE(mat_p_omega, 1)
746 
747  ! first, copy the Gamma-only result from mat_P_omega(1) into all other matrices and
748  ! remove the blocks later which do not belong to the cell index
749  DO i_cell = 2, num_cells_p
750  CALL dbcsr_copy(mat_p_omega(i_cell)%matrix, &
751  mat_p_omega(1)%matrix)
752  END DO
753 
754  IF (jquad == 1 .AND. unit_nr > 0) THEN
755  WRITE (unit_nr, '(T3,A,T66,ES15.2)') 'GW_INFO| RI regularization parameter: ', &
756  qs_env%mp2_env%ri_rpa_im_time%regularization_RI
757  WRITE (unit_nr, '(T3,A,T66,ES15.2)') 'GW_INFO| eps_eigval_S: ', &
758  qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S
759  IF (qs_env%mp2_env%ri_rpa_im_time%make_chi_pos_definite) THEN
760  WRITE (unit_nr, '(T3,A,T81)') &
761  'GW_INFO| Make chi(iw,k) positive definite? TRUE'
762  ELSE
763  WRITE (unit_nr, '(T3,A,T81)') &
764  'GW_INFO| Make chi(iw,k) positive definite? FALSE'
765  END IF
766 
767  END IF
768 
769  DO i_cell = 1, num_cells_p
770 
771  CALL dbcsr_iterator_start(iter, mat_p_omega(i_cell)%matrix)
772  DO WHILE (dbcsr_iterator_blocks_left(iter))
773  CALL dbcsr_iterator_next_block(iter, row, col, data_block)
774 
775  cell_vector(1:3) = matmul(hmat, real(index_to_cell_p(1:3, i_cell), dp))
776  rab_cell_i(1:3) = pbc(particle_set(row)%r(1:3), cell) - &
777  (pbc(particle_set(col)%r(1:3), cell) + cell_vector(1:3))
778  abs_rab_cell_i = sqrt(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2)
779 
780  ! minimum image convention
781  i_cell_is_the_minimum_image_cell = .true.
782  DO j_cell = 1, num_cells_p
783  cell_vector_j(1:3) = matmul(hmat, real(index_to_cell_p(1:3, j_cell), dp))
784  rab_cell_j(1:3) = pbc(particle_set(row)%r(1:3), cell) - &
785  (pbc(particle_set(col)%r(1:3), cell) + cell_vector_j(1:3))
786  abs_rab_cell_j = sqrt(rab_cell_j(1)**2 + rab_cell_j(2)**2 + rab_cell_j(3)**2)
787 
788  IF (abs_rab_cell_i > abs_rab_cell_j + 1.0e-6_dp) THEN
789  i_cell_is_the_minimum_image_cell = .false.
790  END IF
791  END DO
792 
793  IF (.NOT. i_cell_is_the_minimum_image_cell) THEN
794  data_block(:, :) = data_block(:, :)*0.0_dp
795  END IF
796 
797  END DO
798  CALL dbcsr_iterator_stop(iter)
799 
800  END DO
801 
802  CALL timestop(handle)
803 
804  END SUBROUTINE get_mat_cell_t_from_mat_gamma
805 
806 ! **************************************************************************************************
807 !> \brief ...
808 !> \param mat_P_omega ...
809 !> \param mat_P_omega_kp ...
810 !> \param kpoints ...
811 !> \param eps_filter_im_time ...
812 !> \param jquad ...
813 ! **************************************************************************************************
814  SUBROUTINE transform_p_from_real_space_to_kpoints(mat_P_omega, mat_P_omega_kp, &
815  kpoints, eps_filter_im_time, jquad)
816 
817  TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: mat_p_omega, mat_p_omega_kp
818  TYPE(kpoint_type), POINTER :: kpoints
819  REAL(kind=dp), INTENT(IN) :: eps_filter_im_time
820  INTEGER, INTENT(IN) :: jquad
821 
822  CHARACTER(LEN=*), PARAMETER :: routinen = 'transform_P_from_real_space_to_kpoints'
823 
824  INTEGER :: handle, icell, nkp, num_integ_points
825 
826  CALL timeset(routinen, handle)
827 
828  num_integ_points = SIZE(mat_p_omega, 1)
829  nkp = SIZE(mat_p_omega, 2)
830 
831  CALL real_space_to_kpoint_transform_rpa(mat_p_omega_kp(1, :), mat_p_omega_kp(2, :), mat_p_omega(jquad, :), &
832  kpoints, eps_filter_im_time)
833 
834  DO icell = 1, SIZE(mat_p_omega, 2)
835  CALL dbcsr_set(mat_p_omega(jquad, icell)%matrix, 0.0_dp)
836  CALL dbcsr_filter(mat_p_omega(jquad, icell)%matrix, 1.0_dp)
837  END DO
838 
839  CALL timestop(handle)
840 
841  END SUBROUTINE transform_p_from_real_space_to_kpoints
842 
843 ! **************************************************************************************************
844 !> \brief ...
845 !> \param real_mat_kp ...
846 !> \param imag_mat_kp ...
847 !> \param mat_real_space ...
848 !> \param kpoints ...
849 !> \param eps_filter_im_time ...
850 !> \param real_mat_real_space ...
851 ! **************************************************************************************************
852  SUBROUTINE real_space_to_kpoint_transform_rpa(real_mat_kp, imag_mat_kp, mat_real_space, &
853  kpoints, eps_filter_im_time, real_mat_real_space)
854 
855  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: real_mat_kp, imag_mat_kp, mat_real_space
856  TYPE(kpoint_type), POINTER :: kpoints
857  REAL(kind=dp), INTENT(IN) :: eps_filter_im_time
858  LOGICAL, INTENT(IN), OPTIONAL :: real_mat_real_space
859 
860  CHARACTER(LEN=*), PARAMETER :: routinen = 'real_space_to_kpoint_transform_rpa'
861 
862  INTEGER :: handle, i_cell, ik, nkp, num_cells
863  INTEGER, DIMENSION(3) :: cell
864  INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
865  LOGICAL :: my_real_mat_real_space
866  REAL(kind=dp) :: arg, coskl, sinkl
867  REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
868  TYPE(dbcsr_type) :: mat_work
869 
870  CALL timeset(routinen, handle)
871 
872  my_real_mat_real_space = .true.
873  IF (PRESENT(real_mat_real_space)) my_real_mat_real_space = real_mat_real_space
874 
875  CALL dbcsr_create(matrix=mat_work, &
876  template=real_mat_kp(1)%matrix, &
877  matrix_type=dbcsr_type_no_symmetry)
878  CALL dbcsr_reserve_all_blocks(mat_work)
879  CALL dbcsr_set(mat_work, 0.0_dp)
880 
881  ! this kpoint environme t should be the kpoints for D(it) and X(it) created in init_cell_index_rpa
882  CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp)
883 
884  NULLIFY (index_to_cell)
885  index_to_cell => kpoints%index_to_cell
886 
887  num_cells = SIZE(index_to_cell, 2)
888 
889  cpassert(SIZE(mat_real_space) >= num_cells/2 + 1)
890 
891  DO ik = 1, nkp
892 
893  CALL dbcsr_reserve_all_blocks(real_mat_kp(ik)%matrix)
894  CALL dbcsr_reserve_all_blocks(imag_mat_kp(ik)%matrix)
895 
896  CALL dbcsr_set(real_mat_kp(ik)%matrix, 0.0_dp)
897  CALL dbcsr_set(imag_mat_kp(ik)%matrix, 0.0_dp)
898 
899  DO i_cell = 1, num_cells/2 + 1
900 
901  cell(:) = index_to_cell(:, i_cell)
902 
903  arg = real(cell(1), dp)*xkp(1, ik) + real(cell(2), dp)*xkp(2, ik) + real(cell(3), dp)*xkp(3, ik)
904  coskl = cos(twopi*arg)
905  sinkl = sin(twopi*arg)
906 
907  IF (my_real_mat_real_space) THEN
908  CALL dbcsr_add_local(real_mat_kp(ik)%matrix, mat_real_space(i_cell)%matrix, 1.0_dp, coskl)
909  CALL dbcsr_add_local(imag_mat_kp(ik)%matrix, mat_real_space(i_cell)%matrix, 1.0_dp, sinkl)
910  ELSE
911  CALL dbcsr_add_local(real_mat_kp(ik)%matrix, mat_real_space(i_cell)%matrix, 1.0_dp, -sinkl)
912  CALL dbcsr_add_local(imag_mat_kp(ik)%matrix, mat_real_space(i_cell)%matrix, 1.0_dp, coskl)
913  END IF
914 
915  IF (.NOT. (cell(1) == 0 .AND. cell(2) == 0 .AND. cell(3) == 0)) THEN
916 
917  CALL dbcsr_transposed(mat_work, mat_real_space(i_cell)%matrix)
918 
919  IF (my_real_mat_real_space) THEN
920  CALL dbcsr_add_local(real_mat_kp(ik)%matrix, mat_work, 1.0_dp, coskl)
921  CALL dbcsr_add_local(imag_mat_kp(ik)%matrix, mat_work, 1.0_dp, -sinkl)
922  ELSE
923  ! for an imaginary real-space matrix, we need to consider the imaginary unit
924  ! and we need to take into account that the transposed gives an extra "-" sign
925  ! because the transposed is actually Hermitian conjugate
926  CALL dbcsr_add_local(real_mat_kp(ik)%matrix, mat_work, 1.0_dp, -sinkl)
927  CALL dbcsr_add_local(imag_mat_kp(ik)%matrix, mat_work, 1.0_dp, -coskl)
928  END IF
929 
930  CALL dbcsr_set(mat_work, 0.0_dp)
931 
932  END IF
933 
934  END DO
935 
936  CALL dbcsr_filter(real_mat_kp(ik)%matrix, eps_filter_im_time)
937  CALL dbcsr_filter(imag_mat_kp(ik)%matrix, eps_filter_im_time)
938 
939  END DO
940 
941  CALL dbcsr_release(mat_work)
942 
943  CALL timestop(handle)
944 
945  END SUBROUTINE real_space_to_kpoint_transform_rpa
946 
947 ! **************************************************************************************************
948 !> \brief ...
949 !> \param mat_a ...
950 !> \param mat_b ...
951 !> \param alpha ...
952 !> \param beta ...
953 ! **************************************************************************************************
954  SUBROUTINE dbcsr_add_local(mat_a, mat_b, alpha, beta)
955  TYPE(dbcsr_type), INTENT(INOUT) :: mat_a, mat_b
956  REAL(kind=dp), INTENT(IN) :: alpha, beta
957 
958  INTEGER :: col, row
959  LOGICAL :: found
960  REAL(kind=dp), DIMENSION(:, :), POINTER :: block_to_compute, data_block
961  TYPE(dbcsr_iterator_type) :: iter
962 
963  CALL dbcsr_iterator_start(iter, mat_b)
964  DO WHILE (dbcsr_iterator_blocks_left(iter))
965  CALL dbcsr_iterator_next_block(iter, row, col, data_block)
966 
967  NULLIFY (block_to_compute)
968  CALL dbcsr_get_block_p(matrix=mat_a, &
969  row=row, col=col, block=block_to_compute, found=found)
970 
971  cpassert(found)
972 
973  block_to_compute(:, :) = alpha*block_to_compute(:, :) + beta*data_block(:, :)
974 
975  END DO
976  CALL dbcsr_iterator_stop(iter)
977 
978  END SUBROUTINE dbcsr_add_local
979 
980 ! **************************************************************************************************
981 !> \brief ...
982 !> \param fm_mat_W_tau ...
983 !> \param cfm_mat_Q ...
984 !> \param fm_mat_L_re ...
985 !> \param fm_mat_L_im ...
986 !> \param dimen_RI ...
987 !> \param num_integ_points ...
988 !> \param jquad ...
989 !> \param ikp ...
990 !> \param tj ...
991 !> \param tau_tj ...
992 !> \param weights_cos_tf_w_to_t ...
993 !> \param ikp_local ...
994 !> \param para_env ...
995 !> \param kpoints ...
996 !> \param qs_env ...
997 !> \param wkp_W ...
998 ! **************************************************************************************************
999  SUBROUTINE compute_wc_real_space_tau_gw(fm_mat_W_tau, cfm_mat_Q, fm_mat_L_re, fm_mat_L_im, &
1000  dimen_RI, num_integ_points, jquad, &
1001  ikp, tj, tau_tj, weights_cos_tf_w_to_t, ikp_local, &
1002  para_env, kpoints, qs_env, wkp_W)
1003 
1004  TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_w_tau
1005  TYPE(cp_cfm_type), INTENT(IN) :: cfm_mat_q
1006  TYPE(cp_fm_type), INTENT(IN) :: fm_mat_l_re, fm_mat_l_im
1007  INTEGER, INTENT(IN) :: dimen_ri, num_integ_points, jquad, ikp
1008  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: tj
1009  REAL(kind=dp), DIMENSION(0:num_integ_points), &
1010  INTENT(IN) :: tau_tj
1011  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: weights_cos_tf_w_to_t
1012  INTEGER, INTENT(IN) :: ikp_local
1013  TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
1014  TYPE(kpoint_type), INTENT(IN), POINTER :: kpoints
1015  TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1016  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wkp_w
1017 
1018  CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_Wc_real_space_tau_GW'
1019 
1020  INTEGER :: handle, handle2, i_global, iatom, iatom_old, iib, iquad, irow, j_global, jatom, &
1021  jatom_old, jcol, jjb, jkp, ncol_local, nkp, nrow_local, num_cells
1022  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_ri_index
1023  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1024  INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
1025  REAL(kind=dp) :: contribution, omega, tau, weight, &
1026  weight_im, weight_re
1027  REAL(kind=dp), DIMENSION(3, 3) :: hmat
1028  REAL(kind=dp), DIMENSION(:), POINTER :: wkp
1029  REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
1030  TYPE(cell_type), POINTER :: cell
1031  TYPE(cp_cfm_type) :: cfm_mat_l, cfm_mat_work, cfm_mat_work_2
1032  TYPE(cp_fm_type) :: fm_dummy, fm_mat_work_global, &
1033  fm_mat_work_local
1034  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1035 
1036  CALL timeset(routinen, handle)
1037 
1038  CALL timeset(routinen//"_1", handle2)
1039 
1040  CALL cp_cfm_create(cfm_mat_work, cfm_mat_q%matrix_struct)
1041  CALL cp_cfm_set_all(cfm_mat_work, z_zero)
1042 
1043  CALL cp_cfm_create(cfm_mat_work_2, cfm_mat_q%matrix_struct)
1044  CALL cp_cfm_set_all(cfm_mat_work_2, z_zero)
1045 
1046  CALL cp_cfm_create(cfm_mat_l, cfm_mat_q%matrix_struct)
1047  CALL cp_cfm_set_all(cfm_mat_l, z_zero)
1048 
1049  ! Copy fm_mat_L_re and fm_mat_L_re to cfm_mat_L
1050  CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mat_l, z_one, fm_mat_l_re)
1051  CALL cp_cfm_scale_and_add_fm(z_one, cfm_mat_l, gaussi, fm_mat_l_im)
1052 
1053  CALL cp_fm_create(fm_mat_work_global, fm_mat_w_tau(1)%matrix_struct)
1054  CALL cp_fm_set_all(fm_mat_work_global, 0.0_dp)
1055 
1056  CALL cp_fm_create(fm_mat_work_local, cfm_mat_q%matrix_struct)
1057  CALL cp_fm_set_all(fm_mat_work_local, 0.0_dp)
1058 
1059  CALL timestop(handle2)
1060 
1061  CALL timeset(routinen//"_2", handle2)
1062 
1063  ! calculate [1+Q(iw')]^-1
1064  CALL cp_cfm_cholesky_invert(cfm_mat_q)
1065 
1066  ! symmetrize the result
1067  CALL cp_cfm_upper_to_full(cfm_mat_q)
1068 
1069  ! subtract exchange part by subtracing identity matrix from epsilon
1070  CALL cp_cfm_get_info(matrix=cfm_mat_q, &
1071  nrow_local=nrow_local, &
1072  ncol_local=ncol_local, &
1073  row_indices=row_indices, &
1074  col_indices=col_indices)
1075 
1076  DO jjb = 1, ncol_local
1077  j_global = col_indices(jjb)
1078  DO iib = 1, nrow_local
1079  i_global = row_indices(iib)
1080  IF (j_global == i_global .AND. i_global <= dimen_ri) THEN
1081  cfm_mat_q%local_data(iib, jjb) = cfm_mat_q%local_data(iib, jjb) - z_one
1082  END IF
1083  END DO
1084  END DO
1085 
1086  CALL timestop(handle2)
1087 
1088  CALL timeset(routinen//"_3", handle2)
1089 
1090  ! work = epsilon(iw,k)*V^1/2(k)
1091  CALL parallel_gemm('N', 'N', dimen_ri, dimen_ri, dimen_ri, z_one, cfm_mat_q, cfm_mat_l, &
1092  z_zero, cfm_mat_work)
1093 
1094  ! W(iw,k) = V^1/2(k)*work
1095  CALL parallel_gemm('N', 'N', dimen_ri, dimen_ri, dimen_ri, z_one, cfm_mat_l, cfm_mat_work, &
1096  z_zero, cfm_mat_work_2)
1097 
1098  CALL timestop(handle2)
1099 
1100  CALL timeset(routinen//"_4", handle2)
1101 
1102  CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp, nkp=nkp)
1103  index_to_cell => kpoints%index_to_cell
1104  num_cells = SIZE(index_to_cell, 2)
1105 
1106  CALL cp_cfm_set_all(cfm_mat_work, z_zero)
1107 
1108  ALLOCATE (atom_from_ri_index(dimen_ri))
1109 
1110  CALL get_atom_index_from_basis_function_index(qs_env, atom_from_ri_index, dimen_ri, "RI_AUX")
1111 
1112  NULLIFY (cell, particle_set)
1113  CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1114  CALL get_cell(cell=cell, h=hmat)
1115  iatom_old = 0
1116  jatom_old = 0
1117 
1118  CALL cp_cfm_get_info(matrix=cfm_mat_q, &
1119  nrow_local=nrow_local, &
1120  ncol_local=ncol_local, &
1121  row_indices=row_indices, &
1122  col_indices=col_indices)
1123 
1124  DO irow = 1, nrow_local
1125  DO jcol = 1, ncol_local
1126 
1127  iatom = atom_from_ri_index(row_indices(irow))
1128  jatom = atom_from_ri_index(col_indices(jcol))
1129 
1130  IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old) THEN
1131 
1132  ! symmetrize=.FALSE. necessary since we already have a symmetrized index_to_cell
1133  CALL compute_weight_re_im(weight_re, weight_im, &
1134  num_cells, iatom, jatom, xkp(1:3, ikp), wkp_w(ikp), &
1135  cell, index_to_cell, hmat, particle_set)
1136 
1137  iatom_old = iatom
1138  jatom_old = jatom
1139 
1140  END IF
1141 
1142  contribution = weight_re*real(cfm_mat_work_2%local_data(irow, jcol)) + &
1143  weight_im*aimag(cfm_mat_work_2%local_data(irow, jcol))
1144 
1145  fm_mat_work_local%local_data(irow, jcol) = fm_mat_work_local%local_data(irow, jcol) + contribution
1146 
1147  END DO
1148  END DO
1149 
1150  CALL timestop(handle2)
1151 
1152  CALL timeset(routinen//"_5", handle2)
1153 
1154  IF (ikp_local == -1) THEN
1155 
1156  CALL cp_fm_copy_general(fm_mat_work_local, fm_mat_work_global, para_env)
1157 
1158  DO iquad = 1, num_integ_points
1159 
1160  omega = tj(jquad)
1161  tau = tau_tj(iquad)
1162  weight = weights_cos_tf_w_to_t(iquad, jquad)*cos(tau*omega)
1163 
1164  IF (jquad == 1 .AND. ikp == 1) THEN
1165  CALL cp_fm_set_all(matrix=fm_mat_w_tau(iquad), alpha=0.0_dp)
1166  END IF
1167 
1168  CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_mat_w_tau(iquad), beta=weight, matrix_b=fm_mat_work_global)
1169 
1170  END DO
1171 
1172  ELSE
1173 
1174  DO jkp = 1, nkp
1175 
1176  CALL para_env%sync()
1177 
1178  IF (ikp_local == jkp) THEN
1179  CALL cp_fm_copy_general(fm_mat_work_local, fm_mat_work_global, para_env)
1180  ELSE
1181  CALL cp_fm_copy_general(fm_dummy, fm_mat_work_global, para_env)
1182  END IF
1183 
1184  CALL para_env%sync()
1185 
1186  DO iquad = 1, num_integ_points
1187 
1188  omega = tj(jquad)
1189  tau = tau_tj(iquad)
1190  weight = weights_cos_tf_w_to_t(iquad, jquad)*cos(tau*omega)
1191 
1192  IF (jquad == 1 .AND. jkp == 1) THEN
1193  CALL cp_fm_set_all(matrix=fm_mat_w_tau(iquad), alpha=0.0_dp)
1194  END IF
1195 
1196  CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_mat_w_tau(iquad), beta=weight, &
1197  matrix_b=fm_mat_work_global)
1198 
1199  END DO
1200 
1201  END DO
1202 
1203  END IF
1204 
1205  CALL cp_cfm_release(cfm_mat_work)
1206  CALL cp_cfm_release(cfm_mat_work_2)
1207  CALL cp_cfm_release(cfm_mat_l)
1208  CALL cp_fm_release(fm_mat_work_global)
1209  CALL cp_fm_release(fm_mat_work_local)
1210 
1211  DEALLOCATE (atom_from_ri_index)
1212 
1213  CALL timestop(handle2)
1214 
1215  CALL timestop(handle)
1216 
1217  END SUBROUTINE compute_wc_real_space_tau_gw
1218 
1219 ! **************************************************************************************************
1220 !> \brief ...
1221 !> \param fm_mat_W ...
1222 !> \param fm_matrix_Minv ...
1223 !> \param para_env ...
1224 !> \param dimen_RI ...
1225 !> \param num_integ_points ...
1226 ! **************************************************************************************************
1227  SUBROUTINE wc_to_minv_wc_minv(fm_mat_W, fm_matrix_Minv, para_env, dimen_RI, num_integ_points)
1228  TYPE(cp_fm_type), DIMENSION(:) :: fm_mat_w
1229  TYPE(cp_fm_type), DIMENSION(:, :) :: fm_matrix_minv
1230  TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
1231  INTEGER :: dimen_ri, num_integ_points
1232 
1233  CHARACTER(LEN=*), PARAMETER :: routinen = 'Wc_to_Minv_Wc_Minv'
1234 
1235  INTEGER :: handle, jquad
1236  TYPE(cp_fm_type) :: fm_work_minv, fm_work_minv_w
1237 
1238  CALL timeset(routinen, handle)
1239 
1240  CALL cp_fm_create(fm_work_minv, fm_mat_w(1)%matrix_struct)
1241  CALL cp_fm_copy_general(fm_matrix_minv(1, 1), fm_work_minv, para_env)
1242 
1243  CALL cp_fm_create(fm_work_minv_w, fm_mat_w(1)%matrix_struct)
1244 
1245  DO jquad = 1, num_integ_points
1246 
1247  CALL parallel_gemm('N', 'N', dimen_ri, dimen_ri, dimen_ri, 1.0_dp, fm_work_minv, fm_mat_w(jquad), &
1248  0.0_dp, fm_work_minv_w)
1249  CALL parallel_gemm('N', 'N', dimen_ri, dimen_ri, dimen_ri, 1.0_dp, fm_work_minv_w, fm_work_minv, &
1250  0.0_dp, fm_mat_w(jquad))
1251 
1252  END DO
1253 
1254  CALL cp_fm_release(fm_work_minv)
1255 
1256  CALL cp_fm_release(fm_work_minv_w)
1257 
1258  CALL timestop(handle)
1259 
1260  END SUBROUTINE wc_to_minv_wc_minv
1261 
1262 ! **************************************************************************************************
1263 !> \brief ...
1264 !> \param qs_env ...
1265 !> \param wkp_W ...
1266 !> \param wkp_V ...
1267 !> \param kpoints ...
1268 !> \param h_inv ...
1269 !> \param periodic ...
1270 ! **************************************************************************************************
1271  SUBROUTINE compute_wkp_w(qs_env, wkp_W, wkp_V, kpoints, h_inv, periodic)
1272 
1273  TYPE(qs_environment_type), POINTER :: qs_env
1274  REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
1275  INTENT(OUT) :: wkp_w, wkp_v
1276  TYPE(kpoint_type), POINTER :: kpoints
1277  REAL(kind=dp), DIMENSION(3, 3) :: h_inv
1278  INTEGER, DIMENSION(3) :: periodic
1279 
1280  CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_wkp_W'
1281 
1282  INTEGER :: handle, i_x, ikp, info, j_y, k_z, &
1283  kpoint_weights_w_method, n_x, n_y, &
1284  n_z, nkp, nsuperfine, num_lin_eqs
1285  REAL(kind=dp) :: exp_kpoints, integral, k_sq, weight
1286  REAL(kind=dp), DIMENSION(3) :: k_vec, x_vec
1287  REAL(kind=dp), DIMENSION(:), POINTER :: right_side, wkp, wkp_tmp
1288  REAL(kind=dp), DIMENSION(:, :), POINTER :: matrix_lin_eqs, xkp
1289 
1290  CALL timeset(routinen, handle)
1291 
1292  kpoint_weights_w_method = qs_env%mp2_env%ri_rpa_im_time%kpoint_weights_W_method
1293 
1294  CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp, nkp=nkp)
1295 
1296  ! we determine the kpoint weights of the Monkhors Pack mesh new
1297  ! such that the functions 1/k^2, 1/k and const are integrated exactly
1298  ! in the Brillouin zone
1299  ! this is done by minimizing sum_i |w_i|^2 where w_i are the weights of
1300  ! the i-th kpoint under the following constraints:
1301  ! 1) 1/k^2, 1/k and const are integrated exactly
1302  ! 2) the kpoint weights of kpoints with identical absolute value are
1303  ! the same, of e.g. (1/8,3/8,3/8) same weight as for (3/8,1/8,3/8)
1304  ! for 1d and 2d materials: we use ordinary Monkhorst-Pack weights, checked
1305  ! by SUM(periodic) == 3
1306  ALLOCATE (wkp_v(nkp), wkp_w(nkp))
1307 
1308  ! for exchange part of self-energy, we use truncated Coulomb operator that should be fine
1309  ! with uniform weights (without k-point extrapolation)
1310  IF (ALLOCATED(qs_env%mp2_env%ri_rpa_im_time%wkp_V)) THEN
1311  wkp_v(:) = qs_env%mp2_env%ri_rpa_im_time%wkp_V(:)
1312  ELSE
1313  wkp_v(:) = wkp(:)
1314  END IF
1315 
1316  IF (kpoint_weights_w_method == kp_weights_w_uniform) THEN
1317 
1318  ! in the k-point weights wkp, there might be k-point extrapolation included
1319  wkp_w(:) = wkp(:)
1320 
1321  ELSE IF (kpoint_weights_w_method == kp_weights_w_tailored .OR. &
1322  kpoint_weights_w_method == kp_weights_w_auto) THEN
1323 
1324  IF (kpoint_weights_w_method == kp_weights_w_tailored) &
1325  exp_kpoints = qs_env%mp2_env%ri_rpa_im_time%exp_tailored_weights
1326 
1327  IF (kpoint_weights_w_method == kp_weights_w_auto) THEN
1328  IF (sum(periodic) == 2) exp_kpoints = -1.0_dp
1329  END IF
1330 
1331  ! first, compute the integral of f(k)=1/k^2 and 1/k on super fine grid
1332  nsuperfine = 500
1333  integral = 0.0_dp
1334 
1335  IF (periodic(1) == 1) THEN
1336  n_x = nsuperfine
1337  ELSE
1338  n_x = 1
1339  END IF
1340  IF (periodic(2) == 1) THEN
1341  n_y = nsuperfine
1342  ELSE
1343  n_y = 1
1344  END IF
1345  IF (periodic(3) == 1) THEN
1346  n_z = nsuperfine
1347  ELSE
1348  n_z = 1
1349  END IF
1350 
1351  ! actually, there is the factor *det_3x3(h_inv) missing to account for the
1352  ! integration volume but for wkp det_3x3(h_inv) is needed
1353  weight = 1.0_dp/(real(n_x, dp)*real(n_y, dp)*real(n_z, dp))
1354  DO i_x = 1, n_x
1355  DO j_y = 1, n_y
1356  DO k_z = 1, n_z
1357 
1358  IF (periodic(1) == 1) THEN
1359  x_vec(1) = (real(i_x - nsuperfine/2, dp) - 0.5_dp)/real(nsuperfine, dp)
1360  ELSE
1361  x_vec(1) = 0.0_dp
1362  END IF
1363  IF (periodic(2) == 1) THEN
1364  x_vec(2) = (real(j_y - nsuperfine/2, dp) - 0.5_dp)/real(nsuperfine, dp)
1365  ELSE
1366  x_vec(2) = 0.0_dp
1367  END IF
1368  IF (periodic(3) == 1) THEN
1369  x_vec(3) = (real(k_z - nsuperfine/2, dp) - 0.5_dp)/real(nsuperfine, dp)
1370  ELSE
1371  x_vec(3) = 0.0_dp
1372  END IF
1373 
1374  k_vec = matmul(h_inv(1:3, 1:3), x_vec)
1375  k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2
1376  integral = integral + weight*k_sq**(exp_kpoints*0.5_dp)
1377 
1378  END DO
1379  END DO
1380  END DO
1381 
1382  num_lin_eqs = nkp + 2
1383 
1384  ALLOCATE (matrix_lin_eqs(num_lin_eqs, num_lin_eqs))
1385  matrix_lin_eqs(:, :) = 0.0_dp
1386 
1387  DO ikp = 1, nkp
1388 
1389  k_vec = matmul(h_inv(1:3, 1:3), xkp(1:3, ikp))
1390  k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2
1391 
1392  matrix_lin_eqs(ikp, ikp) = 2.0_dp
1393  matrix_lin_eqs(ikp, nkp + 1) = 1.0_dp
1394  matrix_lin_eqs(nkp + 1, ikp) = 1.0_dp
1395 
1396  matrix_lin_eqs(ikp, nkp + 2) = k_sq**(exp_kpoints*0.5_dp)
1397  matrix_lin_eqs(nkp + 2, ikp) = k_sq**(exp_kpoints*0.5_dp)
1398 
1399  END DO
1400 
1401  CALL invmat(matrix_lin_eqs, info)
1402  ! check whether inversion was successful
1403  cpassert(info == 0)
1404 
1405  ALLOCATE (right_side(num_lin_eqs))
1406  right_side = 0.0_dp
1407  right_side(nkp + 1) = 1.0_dp
1408  ! divide integral by two because CP2K k-mesh already considers symmetry k <-> -k
1409  right_side(nkp + 2) = integral
1410 
1411  ALLOCATE (wkp_tmp(num_lin_eqs))
1412 
1413  wkp_tmp(1:num_lin_eqs) = matmul(matrix_lin_eqs, right_side)
1414 
1415  wkp_w(1:nkp) = wkp_tmp(1:nkp)
1416 
1417  DEALLOCATE (matrix_lin_eqs, right_side, wkp_tmp)
1418 
1419  END IF
1420 
1421  CALL timestop(handle)
1422 
1423  END SUBROUTINE
1424 
1425 ! **************************************************************************************************
1426 !> \brief ...
1427 !> \param cfm_mat_Q ...
1428 ! **************************************************************************************************
1429  SUBROUTINE cp_cfm_upper_to_full(cfm_mat_Q)
1430 
1431  TYPE(cp_cfm_type), INTENT(IN) :: cfm_mat_q
1432 
1433  CHARACTER(LEN=*), PARAMETER :: routinen = 'cp_cfm_upper_to_full'
1434 
1435  INTEGER :: handle, i_global, iib, j_global, jjb, &
1436  ncol_local, nrow_local
1437  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1438  TYPE(cp_cfm_type) :: cfm_mat_work
1439 
1440  CALL timeset(routinen, handle)
1441 
1442  CALL cp_cfm_create(cfm_mat_work, cfm_mat_q%matrix_struct)
1443 
1444  ! get info of fm_mat_Q
1445  CALL cp_cfm_get_info(matrix=cfm_mat_q, &
1446  nrow_local=nrow_local, &
1447  ncol_local=ncol_local, &
1448  row_indices=row_indices, &
1449  col_indices=col_indices)
1450 
1451  DO jjb = 1, ncol_local
1452  j_global = col_indices(jjb)
1453  DO iib = 1, nrow_local
1454  i_global = row_indices(iib)
1455  IF (j_global < i_global) THEN
1456  cfm_mat_q%local_data(iib, jjb) = z_zero
1457  END IF
1458  IF (j_global == i_global) THEN
1459  cfm_mat_q%local_data(iib, jjb) = cfm_mat_q%local_data(iib, jjb)/(2.0_dp, 0.0_dp)
1460  END IF
1461  END DO
1462  END DO
1463 
1464  CALL cp_cfm_transpose(cfm_mat_q, 'C', cfm_mat_work)
1465 
1466  CALL cp_cfm_scale_and_add(z_one, cfm_mat_q, z_one, cfm_mat_work)
1467 
1468  CALL cp_cfm_release(cfm_mat_work)
1469 
1470  CALL timestop(handle)
1471 
1472  END SUBROUTINE cp_cfm_upper_to_full
1473 
1474 ! **************************************************************************************************
1475 !> \brief ...
1476 !> \param qs_env ...
1477 !> \param Eigenval_kp ...
1478 ! **************************************************************************************************
1479  SUBROUTINE get_bandstruc_and_k_dependent_mos(qs_env, Eigenval_kp)
1480  TYPE(qs_environment_type), POINTER :: qs_env
1481  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: eigenval_kp
1482 
1483  CHARACTER(LEN=*), PARAMETER :: routinen = 'get_bandstruc_and_k_dependent_MOs'
1484 
1485  INTEGER :: handle, ikp, ispin, nmo, nspins
1486  INTEGER, DIMENSION(3) :: nkp_grid_g
1487  REAL(kind=dp), DIMENSION(:), POINTER :: ev
1488  REAL(kind=dp), DIMENSION(:, :), POINTER :: kpgeneral
1489  TYPE(kpoint_type), POINTER :: kpoints_sigma
1490  TYPE(mp_para_env_type), POINTER :: para_env
1491 
1492  CALL timeset(routinen, handle)
1493 
1494  NULLIFY (qs_env%mp2_env%ri_rpa_im_time%kpoints_G, &
1495  qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma, &
1496  qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma_no_xc, &
1497  para_env)
1498 
1499  nkp_grid_g(1:3) = (/1, 1, 1/)
1500 
1501  CALL get_qs_env(qs_env=qs_env, para_env=para_env)
1502 
1503  CALL create_kp_and_calc_kp_orbitals(qs_env, qs_env%mp2_env%ri_rpa_im_time%kpoints_G, &
1504  "MONKHORST-PACK", para_env%num_pe, &
1505  mp_grid=nkp_grid_g(1:3))
1506 
1507  IF (qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
1508 
1509  ! set up k-points for GW band structure calculation, will be completed later
1510  CALL get_kpgeneral_for_sigma_kpoints(qs_env, kpgeneral)
1511 
1512  CALL create_kp_and_calc_kp_orbitals(qs_env, qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma, &
1513  "GENERAL", para_env%num_pe, &
1514  kpgeneral=kpgeneral)
1515 
1516  CALL create_kp_and_calc_kp_orbitals(qs_env, qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma_no_xc, &
1517  "GENERAL", para_env%num_pe, &
1518  kpgeneral=kpgeneral, with_xc_terms=.false.)
1519 
1520  kpoints_sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma
1521  nmo = SIZE(eigenval_kp, 1)
1522  nspins = SIZE(eigenval_kp, 3)
1523 
1524  ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%Eigenval_Gamma(nmo))
1525  qs_env%mp2_env%ri_rpa_im_time%Eigenval_Gamma(:) = eigenval_kp(:, 1, 1)
1526 
1527  DEALLOCATE (eigenval_kp)
1528 
1529  ALLOCATE (eigenval_kp(nmo, kpoints_sigma%nkp, nspins))
1530 
1531  DO ikp = 1, kpoints_sigma%nkp
1532 
1533  DO ispin = 1, nspins
1534 
1535  ev => kpoints_sigma%kp_env(ikp)%kpoint_env%mos(1, ispin)%eigenvalues
1536 
1537  eigenval_kp(:, ikp, ispin) = ev(:)
1538 
1539  END DO
1540 
1541  END DO
1542 
1543  DEALLOCATE (kpgeneral)
1544 
1545  END IF
1546 
1547  CALL release_hfx_stuff(qs_env)
1548 
1549  CALL timestop(handle)
1550 
1551  END SUBROUTINE get_bandstruc_and_k_dependent_mos
1552 
1553 ! **************************************************************************************************
1554 !> \brief releases part of the given qs_env in order to save memory
1555 !> \param qs_env the object to release
1556 ! **************************************************************************************************
1557  SUBROUTINE release_hfx_stuff(qs_env)
1558  TYPE(qs_environment_type), POINTER :: qs_env
1559 
1560  IF (ASSOCIATED(qs_env%x_data) .AND. .NOT. qs_env%mp2_env%ri_g0w0%do_ri_Sigma_x) THEN
1561  CALL hfx_release(qs_env%x_data)
1562  END IF
1563 
1564  END SUBROUTINE release_hfx_stuff
1565 
1566 ! **************************************************************************************************
1567 !> \brief ...
1568 !> \param qs_env ...
1569 !> \param kpoints ...
1570 !> \param scheme ...
1571 !> \param group_size_ext ...
1572 !> \param mp_grid ...
1573 !> \param kpgeneral ...
1574 !> \param with_xc_terms ...
1575 ! **************************************************************************************************
1576  SUBROUTINE create_kp_and_calc_kp_orbitals(qs_env, kpoints, scheme, &
1577  group_size_ext, mp_grid, kpgeneral, with_xc_terms)
1578 
1579  TYPE(qs_environment_type), POINTER :: qs_env
1580  TYPE(kpoint_type), POINTER :: kpoints
1581  CHARACTER(LEN=*), INTENT(IN) :: scheme
1582  INTEGER :: group_size_ext
1583  INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: mp_grid
1584  REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
1585  OPTIONAL :: kpgeneral
1586  LOGICAL, OPTIONAL :: with_xc_terms
1587 
1588  CHARACTER(LEN=*), PARAMETER :: routinen = 'create_kp_and_calc_kp_orbitals'
1589  COMPLEX(KIND=dp), PARAMETER :: cone = cmplx(1.0_dp, 0.0_dp, kind=dp), &
1590  czero = cmplx(0.0_dp, 0.0_dp, kind=dp), ione = cmplx(0.0_dp, 1.0_dp, kind=dp)
1591 
1592  INTEGER :: handle, i_dim, i_re_im, ikp, ispin, nkp, &
1593  nspins
1594  INTEGER, DIMENSION(3) :: cell_grid, periodic
1595  LOGICAL :: my_with_xc_terms
1596  REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
1597  TYPE(cell_type), POINTER :: cell
1598  TYPE(cp_blacs_env_type), POINTER :: blacs_env
1599  TYPE(cp_cfm_type) :: cksmat, cmos, csmat, cwork
1600  TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1601  TYPE(cp_fm_type) :: fm_work
1602  TYPE(cp_fm_type), POINTER :: imos, rmos
1603  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_s_desymm
1604  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_ks_kp, mat_s_kp
1605  TYPE(dft_control_type), POINTER :: dft_control
1606  TYPE(kpoint_env_type), POINTER :: kp
1607  TYPE(mp_para_env_type), POINTER :: para_env
1608  TYPE(qs_scf_env_type), POINTER :: scf_env
1609  TYPE(scf_control_type), POINTER :: scf_control
1610 
1611  CALL timeset(routinen, handle)
1612 
1613  my_with_xc_terms = .true.
1614  IF (PRESENT(with_xc_terms)) my_with_xc_terms = with_xc_terms
1615 
1616  CALL get_qs_env(qs_env, &
1617  para_env=para_env, &
1618  blacs_env=blacs_env, &
1619  matrix_s=matrix_s, &
1620  scf_env=scf_env, &
1621  scf_control=scf_control, &
1622  cell=cell)
1623 
1624  ! get kpoints
1625  CALL calculate_kpoints_for_bs(kpoints, scheme, kpgeneral=kpgeneral, mp_grid=mp_grid, &
1626  group_size_ext=group_size_ext)
1627 
1628  CALL kpoint_env_initialize(kpoints, para_env, blacs_env)
1629 
1630  ! calculate all MOs that are accessible in the given
1631  ! Gaussian AO basis, therefore nadd=1E10
1632  CALL kpoint_initialize_mos(kpoints, qs_env%mos, 2000000000)
1633  CALL kpoint_initialize_mo_set(kpoints)
1634 
1635  CALL get_cell(cell=cell, periodic=periodic)
1636 
1637  DO i_dim = 1, 3
1638  ! we have at most 3 neigboring cells per dimension and at least one because
1639  ! the density response at Gamma is only divided to neighboring
1640  IF (periodic(i_dim) == 1) THEN
1641  cell_grid(i_dim) = max(min((kpoints%nkp_grid(i_dim)/2)*2 - 1, 1), 3)
1642  ELSE
1643  cell_grid(i_dim) = 1
1644  END IF
1645  END DO
1646  CALL init_cell_index_rpa(cell_grid, kpoints%cell_to_index, kpoints%index_to_cell, cell)
1647 
1648  ! get S(k)
1649  CALL get_qs_env(qs_env, matrix_s=matrix_s, scf_env=scf_env, scf_control=scf_control, dft_control=dft_control)
1650 
1651  NULLIFY (matrix_s_desymm)
1652  CALL dbcsr_allocate_matrix_set(matrix_s_desymm, 1)
1653  ALLOCATE (matrix_s_desymm(1)%matrix)
1654  CALL dbcsr_create(matrix=matrix_s_desymm(1)%matrix, template=matrix_s(1)%matrix, &
1655  matrix_type=dbcsr_type_no_symmetry)
1656  CALL dbcsr_desymmetrize(matrix_s(1)%matrix, matrix_s_desymm(1)%matrix)
1657 
1658  CALL mat_kp_from_mat_gamma(qs_env, mat_s_kp, matrix_s_desymm(1)%matrix, kpoints, 1)
1659 
1660  CALL get_kpoint_info(kpoints, nkp=nkp)
1661 
1662  matrix_struct => kpoints%kp_env(1)%kpoint_env%wmat(1, 1)%matrix_struct
1663 
1664  CALL cp_cfm_create(cksmat, matrix_struct)
1665  CALL cp_cfm_create(csmat, matrix_struct)
1666  CALL cp_cfm_create(cmos, matrix_struct)
1667  CALL cp_cfm_create(cwork, matrix_struct)
1668  CALL cp_fm_create(fm_work, matrix_struct)
1669 
1670  nspins = dft_control%nspins
1671 
1672  DO ispin = 1, nspins
1673 
1674  ! get H(k)
1675  IF (my_with_xc_terms) THEN
1676  CALL mat_kp_from_mat_gamma(qs_env, mat_ks_kp, qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix, kpoints, ispin)
1677  ELSE
1678  CALL mat_kp_from_mat_gamma(qs_env, mat_ks_kp, qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix, &
1679  kpoints, ispin)
1680  END IF
1681 
1682  DO ikp = 1, nkp
1683 
1684  CALL copy_dbcsr_to_fm(mat_ks_kp(ikp, 1)%matrix, kpoints%kp_env(ikp)%kpoint_env%wmat(1, ispin))
1685  CALL cp_cfm_scale_and_add_fm(czero, cksmat, cone, kpoints%kp_env(ikp)%kpoint_env%wmat(1, ispin))
1686 
1687  CALL copy_dbcsr_to_fm(mat_ks_kp(ikp, 2)%matrix, kpoints%kp_env(ikp)%kpoint_env%wmat(2, ispin))
1688  CALL cp_cfm_scale_and_add_fm(cone, cksmat, ione, kpoints%kp_env(ikp)%kpoint_env%wmat(2, ispin))
1689 
1690  CALL copy_dbcsr_to_fm(mat_s_kp(ikp, 1)%matrix, fm_work)
1691  CALL cp_cfm_scale_and_add_fm(czero, csmat, cone, fm_work)
1692 
1693  CALL copy_dbcsr_to_fm(mat_s_kp(ikp, 2)%matrix, fm_work)
1694  CALL cp_cfm_scale_and_add_fm(cone, csmat, ione, fm_work)
1695 
1696  kp => kpoints%kp_env(ikp)%kpoint_env
1697 
1698  CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
1699  CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
1700 
1701  IF (scf_env%cholesky_method == cholesky_off .OR. &
1702  qs_env%mp2_env%ri_rpa_im_time%make_overlap_mat_ao_pos_definite) THEN
1703  CALL cp_cfm_geeig_canon(cksmat, csmat, cmos, eigenvalues, cwork, scf_control%eps_eigval)
1704  ELSE
1705  CALL cp_cfm_geeig(cksmat, csmat, cmos, eigenvalues, cwork)
1706  END IF
1707 
1708  CALL cp_cfm_to_fm(cmos, rmos, imos)
1709 
1710  kp%mos(2, ispin)%eigenvalues = eigenvalues
1711 
1712  END DO
1713 
1714  END DO
1715 
1716  DO ikp = 1, nkp
1717  DO i_re_im = 1, 2
1718  CALL dbcsr_deallocate_matrix(mat_ks_kp(ikp, i_re_im)%matrix)
1719  END DO
1720  END DO
1721  DEALLOCATE (mat_ks_kp)
1722 
1723  DO ikp = 1, nkp
1724  DO i_re_im = 1, 2
1725  CALL dbcsr_deallocate_matrix(mat_s_kp(ikp, i_re_im)%matrix)
1726  END DO
1727  END DO
1728  DEALLOCATE (mat_s_kp)
1729 
1730  CALL dbcsr_deallocate_matrix(matrix_s_desymm(1)%matrix)
1731  DEALLOCATE (matrix_s_desymm)
1732 
1733  CALL cp_cfm_release(cksmat)
1734  CALL cp_cfm_release(csmat)
1735  CALL cp_cfm_release(cwork)
1736  CALL cp_cfm_release(cmos)
1737  CALL cp_fm_release(fm_work)
1738 
1739  CALL timestop(handle)
1740 
1741  END SUBROUTINE create_kp_and_calc_kp_orbitals
1742 
1743 ! **************************************************************************************************
1744 !> \brief ...
1745 !> \param qs_env ...
1746 !> \param mat_kp ...
1747 !> \param mat_gamma ...
1748 !> \param kpoints ...
1749 !> \param ispin ...
1750 !> \param real_mat_real_space ...
1751 ! **************************************************************************************************
1752  SUBROUTINE mat_kp_from_mat_gamma(qs_env, mat_kp, mat_gamma, kpoints, ispin, real_mat_real_space)
1753 
1754  TYPE(qs_environment_type), POINTER :: qs_env
1755  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_kp
1756  TYPE(dbcsr_type) :: mat_gamma
1757  TYPE(kpoint_type), POINTER :: kpoints
1758  INTEGER :: ispin
1759  LOGICAL, INTENT(IN), OPTIONAL :: real_mat_real_space
1760 
1761  CHARACTER(LEN=*), PARAMETER :: routinen = 'mat_kp_from_mat_gamma'
1762 
1763  INTEGER :: handle, i_cell, i_re_im, ikp, nkp, &
1764  num_cells
1765  INTEGER, DIMENSION(3) :: periodic
1766  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1767  REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
1768  TYPE(cell_type), POINTER :: cell
1769  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_real_space
1770 
1771  CALL timeset(routinen, handle)
1772 
1773  CALL get_qs_env(qs_env, cell=cell)
1774  CALL get_cell(cell=cell, periodic=periodic)
1775  num_cells = 3**(periodic(1) + periodic(2) + periodic(3))
1776 
1777  NULLIFY (mat_real_space)
1778  CALL dbcsr_allocate_matrix_set(mat_real_space, num_cells)
1779  DO i_cell = 1, num_cells
1780  ALLOCATE (mat_real_space(i_cell)%matrix)
1781  CALL dbcsr_create(matrix=mat_real_space(i_cell)%matrix, &
1782  template=mat_gamma)
1783  CALL dbcsr_reserve_all_blocks(mat_real_space(i_cell)%matrix)
1784  CALL dbcsr_set(mat_real_space(i_cell)%matrix, 0.0_dp)
1785  END DO
1786 
1787  CALL dbcsr_copy(mat_real_space(1)%matrix, mat_gamma)
1788 
1789  CALL get_mat_cell_t_from_mat_gamma(mat_real_space, qs_env, kpoints, 2, 0)
1790 
1791  NULLIFY (xkp, cell_to_index)
1792  CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, cell_to_index=cell_to_index)
1793 
1794  IF (ispin == 1) THEN
1795  NULLIFY (mat_kp)
1796  CALL dbcsr_allocate_matrix_set(mat_kp, nkp, 2)
1797  DO ikp = 1, nkp
1798  DO i_re_im = 1, 2
1799  ALLOCATE (mat_kp(ikp, i_re_im)%matrix)
1800  CALL dbcsr_create(matrix=mat_kp(ikp, i_re_im)%matrix, template=mat_gamma)
1801  CALL dbcsr_reserve_all_blocks(mat_kp(ikp, i_re_im)%matrix)
1802  CALL dbcsr_set(mat_kp(ikp, i_re_im)%matrix, 0.0_dp)
1803  END DO
1804  END DO
1805  END IF
1806 
1807  IF (PRESENT(real_mat_real_space)) THEN
1808  CALL real_space_to_kpoint_transform_rpa(mat_kp(:, 1), mat_kp(:, 2), mat_real_space, kpoints, 0.0_dp, &
1809  real_mat_real_space)
1810  ELSE
1811  CALL real_space_to_kpoint_transform_rpa(mat_kp(:, 1), mat_kp(:, 2), mat_real_space, kpoints, 0.0_dp)
1812  END IF
1813 
1814  DO i_cell = 1, num_cells
1815  CALL dbcsr_deallocate_matrix(mat_real_space(i_cell)%matrix)
1816  END DO
1817  DEALLOCATE (mat_real_space)
1818 
1819  CALL timestop(handle)
1820 
1821  END SUBROUTINE mat_kp_from_mat_gamma
1822 
1823 ! **************************************************************************************************
1824 !> \brief ...
1825 !> \param qs_env ...
1826 !> \param kpgeneral ...
1827 ! **************************************************************************************************
1828  SUBROUTINE get_kpgeneral_for_sigma_kpoints(qs_env, kpgeneral)
1829  TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1830  REAL(kind=dp), DIMENSION(:, :), POINTER :: kpgeneral
1831 
1832  CHARACTER(LEN=*), PARAMETER :: routinen = 'get_kpgeneral_for_Sigma_kpoints'
1833 
1834  INTEGER :: handle, i_kp_in_kp_line, i_special_kp, &
1835  i_x, ikk, j_y, k_z, n_kp_in_kp_line, &
1836  n_special_kp
1837  INTEGER, DIMENSION(:), POINTER :: nkp_grid
1838 
1839  CALL timeset(routinen, handle)
1840 
1841  n_special_kp = qs_env%mp2_env%ri_g0w0%n_special_kp
1842  n_kp_in_kp_line = qs_env%mp2_env%ri_g0w0%n_kp_in_kp_line
1843  IF (n_special_kp > 0) THEN
1844  qs_env%mp2_env%ri_g0w0%nkp_self_energy_special_kp = n_kp_in_kp_line*(n_special_kp - 1) + 1
1845  ELSE
1846  qs_env%mp2_env%ri_g0w0%nkp_self_energy_special_kp = 0
1847  END IF
1848 
1849  qs_env%mp2_env%ri_g0w0%nkp_self_energy_monkh_pack = qs_env%mp2_env%ri_g0w0%kp_grid_Sigma(1)* &
1850  qs_env%mp2_env%ri_g0w0%kp_grid_Sigma(2)* &
1851  qs_env%mp2_env%ri_g0w0%kp_grid_Sigma(3)
1852 
1853  qs_env%mp2_env%ri_g0w0%nkp_self_energy = qs_env%mp2_env%ri_g0w0%nkp_self_energy_special_kp + &
1854  qs_env%mp2_env%ri_g0w0%nkp_self_energy_monkh_pack
1855 
1856  ALLOCATE (kpgeneral(3, qs_env%mp2_env%ri_g0w0%nkp_self_energy))
1857 
1858  IF (n_special_kp > 0) THEN
1859 
1860  kpgeneral(1:3, 1) = qs_env%mp2_env%ri_g0w0%xkp_special_kp(1:3, 1)
1861 
1862  ikk = 1
1863 
1864  DO i_special_kp = 2, n_special_kp
1865  DO i_kp_in_kp_line = 1, n_kp_in_kp_line
1866 
1867  ikk = ikk + 1
1868  kpgeneral(1:3, ikk) = qs_env%mp2_env%ri_g0w0%xkp_special_kp(1:3, i_special_kp - 1) + &
1869  REAL(i_kp_in_kp_line, kind=dp)/real(n_kp_in_kp_line, kind=dp)* &
1870  (qs_env%mp2_env%ri_g0w0%xkp_special_kp(1:3, i_special_kp) - &
1871  qs_env%mp2_env%ri_g0w0%xkp_special_kp(1:3, i_special_kp - 1))
1872 
1873  END DO
1874  END DO
1875 
1876  ELSE
1877 
1878  ikk = 0
1879 
1880  END IF
1881 
1882  nkp_grid => qs_env%mp2_env%ri_g0w0%kp_grid_Sigma
1883 
1884  DO i_x = 1, nkp_grid(1)
1885  DO j_y = 1, nkp_grid(2)
1886  DO k_z = 1, nkp_grid(3)
1887  ikk = ikk + 1
1888  kpgeneral(1, ikk) = real(2*i_x - nkp_grid(1) - 1, kind=dp)/(2._dp*real(nkp_grid(1), kind=dp))
1889  kpgeneral(2, ikk) = real(2*j_y - nkp_grid(2) - 1, kind=dp)/(2._dp*real(nkp_grid(2), kind=dp))
1890  kpgeneral(3, ikk) = real(2*k_z - nkp_grid(3) - 1, kind=dp)/(2._dp*real(nkp_grid(3), kind=dp))
1891  END DO
1892  END DO
1893  END DO
1894 
1895  CALL timestop(handle)
1896 
1897  END SUBROUTINE get_kpgeneral_for_sigma_kpoints
1898 
1899 END MODULE rpa_gw_kpoints_util
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
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_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition: cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_scale_and_add
subroutine, public cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b).
Definition: cp_cfm_basic_linalg.F:244
cp_cfm_basic_linalg::cp_cfm_transpose
subroutine, public cp_cfm_transpose(matrix, trans, matrixt)
Transposes a BLACS distributed complex matrix.
Definition: cp_cfm_basic_linalg.F:1239
cp_cfm_basic_linalg::cp_cfm_scale_and_add_fm
subroutine, public cp_cfm_scale_and_add_fm(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b). where b is a real matrix (adapted from cp_cf...
Definition: cp_cfm_basic_linalg.F:355
cp_cfm_basic_linalg::cp_cfm_cholesky_invert
subroutine, public cp_cfm_cholesky_invert(matrix, n, info_out)
Used to replace Cholesky decomposition by the inverse.
Definition: cp_cfm_basic_linalg.F:982
cp_cfm_basic_linalg::cp_cfm_cholesky_decompose
subroutine, public cp_cfm_cholesky_decompose(matrix, n, info_out)
Used to replace a symmetric positive definite matrix M with its Cholesky decomposition U: M = U^T * U...
Definition: cp_cfm_basic_linalg.F:926
cp_cfm_basic_linalg::cp_cfm_column_scale
subroutine, public cp_cfm_column_scale(matrix_a, scaling)
Scales columns of the full matrix by corresponding factors.
Definition: cp_cfm_basic_linalg.F:649
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition: cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_geeig
subroutine, public cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition: cp_cfm_diag.F:145
cp_cfm_diag::cp_cfm_heevd
subroutine, public cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
Perform a diagonalisation of a complex matrix.
Definition: cp_cfm_diag.F:52
cp_cfm_diag::cp_cfm_geeig_canon
subroutine, public cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
General Eigenvalue Problem AX = BXE Use canonical orthogonalization.
Definition: cp_cfm_diag.F:191
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition: cp_cfm_types.F:12
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
Definition: cp_cfm_types.F:121
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition: cp_cfm_types.F:159
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition: cp_cfm_types.F:607
cp_cfm_types::cp_cfm_set_all
subroutine, public cp_cfm_set_all(matrix, alpha, beta)
Set all elements of the full matrix to alpha. Besides, set all diagonal matrix elements to beta (if g...
Definition: cp_cfm_types.F:179
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition: cp_cfm_types.F:765
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition: cp_dbcsr_operations.F:208
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_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition: cp_fm_basic_linalg.F:167
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_copy_general
subroutine, public cp_fm_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
Definition: cp_fm_types.F:1538
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
hfx_types
Types and set/get functions for HFX.
Definition: hfx_types.F:15
hfx_types::hfx_release
subroutine, public hfx_release(x_data)
This routine deallocates all data structures
Definition: hfx_types.F:1905
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::kp_weights_w_auto
integer, parameter, public kp_weights_w_auto
Definition: input_constants.F:1063
input_constants::kp_weights_w_uniform
integer, parameter, public kp_weights_w_uniform
Definition: input_constants.F:1063
input_constants::cholesky_off
integer, parameter, public cholesky_off
Definition: input_constants.F:195
input_constants::kp_weights_w_tailored
integer, parameter, public kp_weights_w_tailored
Definition: input_constants.F:1063
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kpoint_methods
Routines needed for kpoint calculation.
Definition: kpoint_methods.F:15
kpoint_methods::kpoint_initialize_mo_set
subroutine, public kpoint_initialize_mo_set(kpoint)
...
Definition: kpoint_methods.F:607
kpoint_methods::kpoint_initialize_mos
subroutine, public kpoint_initialize_mos(kpoint, mos, added_mos, for_aux_fit)
Initialize a set of MOs and density matrix for each kpoint (kpoint group)
Definition: kpoint_methods.F:454
kpoint_methods::kpoint_env_initialize
subroutine, public kpoint_env_initialize(kpoint, para_env, blacs_env, with_aux_fit)
Initialize the kpoint environment.
Definition: kpoint_methods.F:266
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_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::z_one
complex(kind=dp), parameter, public z_one
Definition: mathconstants.F:143
mathconstants::gaussi
complex(kind=dp), parameter, public gaussi
Definition: mathconstants.F:142
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition: mathconstants.F:112
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition: mathconstants.F:143
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
mathlib::invmat
subroutine, public invmat(a, info)
returns inverse of matrix using the lapack routines DGETRF and DGETRI
Definition: mathlib.F:543
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
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_band_structure
Calculation of band structures.
Definition: qs_band_structure.F:14
qs_band_structure::calculate_kpoints_for_bs
subroutine, public calculate_kpoints_for_bs(kpoint, scheme, group_size_ext, mp_grid, kpgeneral)
...
Definition: qs_band_structure.F:348
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_scf_types
module that contains the definitions of the scf types
Definition: qs_scf_types.F:14
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_gw_kpoints_util
Routines treating GW and RPA calculations with kpoints.
Definition: rpa_gw_kpoints_util.F:13
rpa_gw_kpoints_util::get_mat_cell_t_from_mat_gamma
subroutine, public get_mat_cell_t_from_mat_gamma(mat_P_omega, qs_env, kpoints, jquad, unit_nr)
...
Definition: rpa_gw_kpoints_util.F:700
rpa_gw_kpoints_util::invert_eps_compute_w_and_erpa_kp
subroutine, public invert_eps_compute_w_and_erpa_kp(dimen_RI, num_integ_points, jquad, nkp, count_ev_sc_GW, para_env, Erpa, tau_tj, tj, wj, weights_cos_tf_w_to_t, wkp_W, do_gw_im_time, do_ri_Sigma_x, do_kpoints_from_Gamma, cfm_mat_Q, ikp_local, mat_P_omega, mat_P_omega_kp, qs_env, eps_filter_im_time, unit_nr, kpoints, fm_mat_Minv_L_kpoints, fm_matrix_L_kpoints, fm_mat_W, fm_mat_RI_global_work, mat_MinvVMinv, fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv)
...
Definition: rpa_gw_kpoints_util.F:136
rpa_gw_kpoints_util::real_space_to_kpoint_transform_rpa
subroutine, public real_space_to_kpoint_transform_rpa(real_mat_kp, imag_mat_kp, mat_real_space, kpoints, eps_filter_im_time, real_mat_real_space)
...
Definition: rpa_gw_kpoints_util.F:854
rpa_gw_kpoints_util::compute_wkp_w
subroutine, public compute_wkp_w(qs_env, wkp_W, wkp_V, kpoints, h_inv, periodic)
...
Definition: rpa_gw_kpoints_util.F:1272
rpa_gw_kpoints_util::get_bandstruc_and_k_dependent_mos
subroutine, public get_bandstruc_and_k_dependent_mos(qs_env, Eigenval_kp)
...
Definition: rpa_gw_kpoints_util.F:1480
rpa_gw_kpoints_util::cp_cfm_upper_to_full
subroutine, public cp_cfm_upper_to_full(cfm_mat_Q)
...
Definition: rpa_gw_kpoints_util.F:1430
rpa_gw_kpoints_util::cp_cfm_power
subroutine, public cp_cfm_power(matrix, threshold, exponent, min_eigval)
...
Definition: rpa_gw_kpoints_util.F:277
rpa_gw_kpoints_util::mat_kp_from_mat_gamma
subroutine, public mat_kp_from_mat_gamma(qs_env, mat_kp, mat_gamma, kpoints, ispin, real_mat_real_space)
...
Definition: rpa_gw_kpoints_util.F:1753
rpa_im_time
Routines for low-scaling RPA/GW with imaginary time.
Definition: rpa_im_time.F:13
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
scf_control_types
parameters that control an scf iteration
Definition: scf_control_types.F:17