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