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