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