Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 : minimum_image_weights
17 : USE cp_blacs_env, ONLY: cp_blacs_env_type
18 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_add_on_diag,&
19 : cp_cfm_get_diag,&
20 : cp_cfm_scale_and_add_fm,&
21 : cp_cfm_uplo_to_full
22 : USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose,&
23 : cp_cfm_cholesky_invert
24 : USE cp_cfm_diag, ONLY: cp_cfm_geeig,&
25 : cp_cfm_geeig_canon,&
26 : cp_cfm_power
27 : USE cp_cfm_types, ONLY: cp_cfm_copy_general,&
28 : cp_cfm_create,&
29 : cp_cfm_get_info,&
30 : cp_cfm_release,&
31 : cp_cfm_set_all,&
32 : cp_cfm_to_cfm,&
33 : cp_cfm_type,&
34 : cp_fm_to_cfm
35 : USE cp_control_types, ONLY: dft_control_type
36 : USE cp_dbcsr_api, ONLY: &
37 : dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_filter, &
38 : dbcsr_get_block_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
39 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
40 : dbcsr_release, dbcsr_set, dbcsr_transposed, dbcsr_type, dbcsr_type_no_symmetry
41 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_all_blocks
42 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
43 : copy_fm_to_dbcsr,&
44 : dbcsr_allocate_matrix_set
45 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
46 : USE cp_fm_struct, ONLY: cp_fm_struct_type
47 : USE cp_fm_types, ONLY: cp_fm_copy_general,&
48 : cp_fm_create,&
49 : cp_fm_release,&
50 : cp_fm_set_all,&
51 : cp_fm_type
52 : USE gw_utils_fm, ONLY: cfm_contract_ABA
53 : USE hfx_types, ONLY: hfx_release
54 : USE input_constants, ONLY: cholesky_off,&
55 : kp_weights_W_auto,&
56 : kp_weights_W_tailored,&
57 : kp_weights_W_uniform
58 : USE kinds, ONLY: dp
59 : USE kpoint_methods, ONLY: kpoint_env_initialize,&
60 : kpoint_initialize_mo_set,&
61 : kpoint_initialize_mos
62 : USE kpoint_types, ONLY: get_kpoint_info,&
63 : kpoint_env_type,&
64 : kpoint_type
65 : USE machine, ONLY: m_walltime
66 : USE mathconstants, ONLY: gaussi,&
67 : twopi,&
68 : z_one,&
69 : z_zero
70 : USE mathlib, ONLY: invmat
71 : USE message_passing, ONLY: mp_para_env_type
72 : USE parallel_gemm_api, ONLY: parallel_gemm
73 : USE particle_types, ONLY: particle_type
74 : USE qs_band_structure, ONLY: calculate_kpoints_for_bs
75 : USE qs_environment_types, ONLY: get_qs_env,&
76 : qs_environment_type
77 : USE qs_mo_types, ONLY: get_mo_set
78 : USE qs_scf_types, ONLY: qs_scf_env_type
79 : USE rpa_gw_im_time_util, ONLY: compute_weight_re_im,&
80 : get_atom_index_from_basis_function_index
81 : USE rpa_im_time, ONLY: init_cell_index_rpa
82 : USE scf_control_types, ONLY: scf_control_type
83 : USE time_frequency_grids, ONLY: time_frequency_grid_type
84 : #include "./base/base_uses.f90"
85 :
86 : IMPLICIT NONE
87 :
88 : PRIVATE
89 :
90 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_gw_kpoints_util'
91 :
92 : PUBLIC :: invert_eps_compute_W_and_Erpa_kp, real_space_to_kpoint_transform_rpa, &
93 : get_mat_cell_T_from_mat_gamma, get_bandstruc_and_k_dependent_MOs, &
94 : compute_wkp_W, mat_kp_from_mat_gamma, cfm_screened_interaction_from_eps_cholesky
95 :
96 : CONTAINS
97 :
98 : ! **************************************************************************************************
99 : !> \brief Form the screened interaction from a Cholesky-factorized dielectric matrix.
100 : !> \param eps_chol input Cholesky factor; overwritten by epsilon^-1-Id or Wc
101 : !> \param B Coulomb factor satisfying V=B^H B
102 : !> \param Wc optional separate output; if absent, eps_chol is overwritten by Wc
103 : !> \param n_active optional number of physical rows/columns
104 : ! **************************************************************************************************
105 5842 : SUBROUTINE cfm_screened_interaction_from_eps_cholesky(eps_chol, B, Wc, n_active)
106 : TYPE(cp_cfm_type), INTENT(INOUT) :: eps_chol
107 : TYPE(cp_cfm_type), INTENT(IN) :: B
108 : TYPE(cp_cfm_type), INTENT(INOUT), OPTIONAL :: Wc
109 : INTEGER, INTENT(IN), OPTIONAL :: n_active
110 :
111 5842 : CALL cp_cfm_cholesky_invert(eps_chol)
112 5842 : CALL cp_cfm_uplo_to_full(eps_chol)
113 5842 : CALL cp_cfm_add_on_diag(eps_chol, -z_one, n_active)
114 5842 : CALL cfm_contract_ABA(B, eps_chol, Wc)
115 :
116 5842 : END SUBROUTINE cfm_screened_interaction_from_eps_cholesky
117 :
118 : ! **************************************************************************************************
119 : !> \brief ...
120 : !> \param dimen_RI ...
121 : !> \param jquad ...
122 : !> \param nkp ...
123 : !> \param count_ev_sc_GW ...
124 : !> \param para_env ...
125 : !> \param Erpa ...
126 : !> \param grid ...
127 : !> \param wkp_W ...
128 : !> \param do_gw_im_time ...
129 : !> \param do_ri_Sigma_x ...
130 : !> \param do_kpoints_from_Gamma ...
131 : !> \param cfm_mat_Q ...
132 : !> \param ikp_local ...
133 : !> \param mat_P_omega ...
134 : !> \param mat_P_omega_kp ...
135 : !> \param qs_env ...
136 : !> \param eps_filter_im_time ...
137 : !> \param unit_nr ...
138 : !> \param kpoints ...
139 : !> \param cfm_mat_Minv_L_kpoints ...
140 : !> \param cfm_matrix_L_kpoints ...
141 : !> \param fm_mat_W ...
142 : !> \param cfm_mat_RI_global_work ...
143 : !> \param mat_MinvVMinv ...
144 : !> \param fm_matrix_Minv ...
145 : !> \param fm_matrix_Minv_Vtrunc_Minv ...
146 : ! **************************************************************************************************
147 132 : SUBROUTINE invert_eps_compute_W_and_Erpa_kp(dimen_RI, jquad, nkp, count_ev_sc_GW, para_env, &
148 264 : Erpa, grid, wkp_W, do_gw_im_time, &
149 : do_ri_Sigma_x, do_kpoints_from_Gamma, &
150 132 : cfm_mat_Q, ikp_local, mat_P_omega, mat_P_omega_kp, &
151 : qs_env, eps_filter_im_time, unit_nr, kpoints, cfm_mat_Minv_L_kpoints, &
152 168 : cfm_matrix_L_kpoints, fm_mat_W, &
153 : cfm_mat_RI_global_work, mat_MinvVMinv, fm_matrix_Minv, &
154 : fm_matrix_Minv_Vtrunc_Minv)
155 :
156 : INTEGER, INTENT(IN) :: dimen_RI, jquad, nkp, count_ev_sc_GW
157 : TYPE(mp_para_env_type), POINTER :: para_env
158 : REAL(KIND=dp), INTENT(INOUT) :: Erpa
159 : TYPE(time_frequency_grid_type), INTENT(IN) :: grid
160 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wkp_W
161 : LOGICAL, INTENT(IN) :: do_gw_im_time, do_ri_Sigma_x, &
162 : do_kpoints_from_Gamma
163 : TYPE(cp_cfm_type), INTENT(INOUT) :: cfm_mat_Q
164 : INTEGER, INTENT(IN) :: ikp_local
165 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: mat_P_omega, mat_P_omega_kp
166 : TYPE(qs_environment_type), POINTER :: qs_env
167 : REAL(KIND=dp), INTENT(IN) :: eps_filter_im_time
168 : INTEGER, INTENT(IN) :: unit_nr
169 : TYPE(kpoint_type), POINTER :: kpoints
170 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_mat_Minv_L_kpoints, &
171 : cfm_matrix_L_kpoints
172 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_W
173 : TYPE(cp_cfm_type) :: cfm_mat_RI_global_work
174 : TYPE(dbcsr_p_type), INTENT(IN) :: mat_MinvVMinv
175 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_Minv, &
176 : fm_matrix_Minv_Vtrunc_Minv
177 :
178 : CHARACTER(LEN=*), PARAMETER :: routineN = 'invert_eps_compute_W_and_Erpa_kp'
179 :
180 : INTEGER :: handle, ikp, num_integ_points
181 : LOGICAL :: do_this_ikp
182 : REAL(KIND=dp) :: t1, t2
183 132 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag_Qomega
184 :
185 132 : CALL timeset(routineN, handle)
186 :
187 132 : num_integ_points = SIZE(grid%frequency)
188 :
189 132 : t1 = m_walltime()
190 :
191 132 : IF (do_kpoints_from_Gamma) THEN
192 96 : CALL get_mat_cell_T_from_mat_gamma(mat_P_omega(jquad, :), qs_env, kpoints, jquad, unit_nr)
193 : END IF
194 :
195 : CALL transform_P_from_real_space_to_kpoints(mat_P_omega, mat_P_omega_kp, &
196 132 : kpoints, eps_filter_im_time, jquad)
197 :
198 396 : ALLOCATE (diag_Qomega(dimen_RI))
199 :
200 132 : IF (unit_nr > 0) WRITE (unit_nr, '(/T3,A,1X,I3)') &
201 66 : 'GW_INFO| Computing chi and W frequency point', jquad
202 :
203 2700 : DO ikp = 1, nkp
204 :
205 : ! parallization, we either have all kpoints on all processors or a single kpoint per group
206 2568 : do_this_ikp = (ikp_local == -1) .OR. (ikp_local == 0 .AND. ikp == 1) .OR. (ikp_local == ikp)
207 : IF (.NOT. do_this_ikp) CYCLE
208 :
209 : ! 1. remove all spurious negative eigenvalues from P(iw,k), multiplication Q(iw,k) = K^H(k)P(iw,k)K(k)
210 : CALL compute_Q_kp_RPA(cfm_mat_Q, &
211 : mat_P_omega_kp, &
212 : cfm_mat_Minv_L_kpoints(ikp), &
213 : cfm_mat_RI_global_work, &
214 : ikp, nkp, ikp_local, para_env, &
215 2568 : qs_env%mp2_env%ri_rpa_im_time%make_chi_pos_definite)
216 :
217 : ! 2. Cholesky decomposition of Id + Q(iw,k)
218 2568 : CALL cholesky_decomp_Q(cfm_mat_Q, diag_Qomega, dimen_RI)
219 :
220 : ! 3. Computing E_c^RPA = E_c^RPA + a_w/N_k*sum_k ln[det(1+Q(iw,k))-Tr(Q(iw,k))]
221 : CALL frequency_and_kpoint_integration(Erpa, cfm_mat_Q, para_env, diag_Qomega, &
222 2568 : dimen_RI, grid%frequency_weights(jquad), kpoints%wkp(ikp))
223 :
224 2700 : IF (do_gw_im_time) THEN
225 :
226 : ! compute S^-1*V*S^-1 for exchange part of the self-energy in real space as W in real space
227 2496 : IF (do_ri_Sigma_x .AND. jquad == 1 .AND. count_ev_sc_GW == 1 .AND. do_kpoints_from_Gamma) THEN
228 :
229 312 : CALL dbcsr_set(mat_MinvVMinv%matrix, 0.0_dp)
230 312 : CALL copy_fm_to_dbcsr(fm_matrix_Minv_Vtrunc_Minv(1, 1), mat_MinvVMinv%matrix, keep_sparsity=.FALSE.)
231 :
232 : END IF
233 2496 : IF (do_kpoints_from_Gamma) THEN
234 : CALL compute_Wc_real_space_tau_GW(fm_mat_W, cfm_mat_Q, cfm_matrix_L_kpoints(ikp), &
235 : dimen_RI, jquad, &
236 : ikp, grid, &
237 2496 : ikp_local, para_env, kpoints, qs_env, wkp_W)
238 : END IF
239 :
240 : END IF
241 : END DO
242 :
243 : ! after the transform of (eps(iw)-1)^-1 from iw to it is done, multiply with V^1/2 to obtain W(it)
244 132 : IF (do_gw_im_time .AND. do_kpoints_from_Gamma .AND. jquad == num_integ_points) THEN
245 16 : CALL Wc_to_Minv_Wc_Minv(fm_mat_W, fm_matrix_Minv, para_env, dimen_RI, num_integ_points)
246 16 : CALL deallocate_kp_matrices(cfm_matrix_L_kpoints, cfm_mat_Minv_L_kpoints)
247 : END IF
248 :
249 132 : DEALLOCATE (diag_Qomega)
250 :
251 132 : t2 = m_walltime()
252 :
253 132 : IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,T56,F25.1)') 'Execution time (s):', t2 - t1
254 :
255 132 : CALL timestop(handle)
256 :
257 132 : END SUBROUTINE invert_eps_compute_W_and_Erpa_kp
258 :
259 : ! **************************************************************************************************
260 : !> \brief ...
261 : !> \param cfm_matrix_L_kpoints ...
262 : !> \param cfm_mat_Minv_L_kpoints ...
263 : ! **************************************************************************************************
264 16 : SUBROUTINE deallocate_kp_matrices(cfm_matrix_L_kpoints, cfm_mat_Minv_L_kpoints)
265 :
266 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_matrix_L_kpoints, &
267 : cfm_mat_Minv_L_kpoints
268 :
269 : CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_kp_matrices'
270 :
271 : INTEGER :: handle, ikp
272 :
273 16 : CALL timeset(routineN, handle)
274 :
275 16 : IF (ALLOCATED(cfm_mat_Minv_L_kpoints)) THEN
276 432 : DO ikp = 1, SIZE(cfm_mat_Minv_L_kpoints)
277 432 : IF (ASSOCIATED(cfm_mat_Minv_L_kpoints(ikp)%matrix_struct)) CALL cp_cfm_release(cfm_mat_Minv_L_kpoints(ikp))
278 : END DO
279 16 : DEALLOCATE (cfm_mat_Minv_L_kpoints)
280 : END IF
281 16 : IF (ALLOCATED(cfm_matrix_L_kpoints)) THEN
282 432 : DO ikp = 1, SIZE(cfm_matrix_L_kpoints)
283 432 : IF (ASSOCIATED(cfm_matrix_L_kpoints(ikp)%matrix_struct)) CALL cp_cfm_release(cfm_matrix_L_kpoints(ikp))
284 : END DO
285 16 : DEALLOCATE (cfm_matrix_L_kpoints)
286 : END IF
287 :
288 16 : CALL timestop(handle)
289 :
290 16 : END SUBROUTINE deallocate_kp_matrices
291 :
292 : ! **************************************************************************************************
293 : !> \brief ...
294 : !> \param cfm_mat_Q ...
295 : !> \param mat_P_omega_kp ...
296 : !> \param cfm_mat_L ...
297 : !> \param cfm_mat_RI_global_work ...
298 : !> \param ikp ...
299 : !> \param nkp ...
300 : !> \param ikp_local ...
301 : !> \param para_env ...
302 : !> \param make_chi_pos_definite ...
303 : ! **************************************************************************************************
304 2568 : SUBROUTINE compute_Q_kp_RPA(cfm_mat_Q, mat_P_omega_kp, cfm_mat_L, &
305 : cfm_mat_RI_global_work, ikp, nkp, ikp_local, para_env, &
306 : make_chi_pos_definite)
307 :
308 : TYPE(cp_cfm_type) :: cfm_mat_Q
309 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: mat_P_omega_kp
310 : TYPE(cp_cfm_type) :: cfm_mat_L, cfm_mat_RI_global_work
311 : INTEGER, INTENT(IN) :: ikp, nkp, ikp_local
312 : TYPE(mp_para_env_type), POINTER :: para_env
313 : LOGICAL, INTENT(IN) :: make_chi_pos_definite
314 :
315 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Q_kp_RPA'
316 :
317 : INTEGER :: handle
318 :
319 2568 : CALL timeset(routineN, handle)
320 :
321 : ! 1. Convert the dbcsr matrix mat_P_omega_kp (that is chi(k,iw)) to a full matrix and
322 : ! distribute it to subgroups
323 2568 : CALL mat_P_to_subgroup(mat_P_omega_kp, cfm_mat_RI_global_work, cfm_mat_Q, ikp, nkp, ikp_local, para_env)
324 :
325 : ! 2. Remove all negative eigenvalues from chi(k,iw)
326 2568 : IF (make_chi_pos_definite) THEN
327 2568 : CALL cp_cfm_power(cfm_mat_Q, threshold=0.0_dp, exponent=1.0_dp)
328 : END IF
329 :
330 : ! Q(iw,k) = L^H(k)*P(iw,k)*L(k)
331 2568 : CALL cfm_contract_ABA(cfm_mat_L, cfm_mat_Q)
332 :
333 2568 : CALL timestop(handle)
334 :
335 2568 : END SUBROUTINE compute_Q_kp_RPA
336 :
337 : ! **************************************************************************************************
338 : !> \brief ...
339 : !> \param mat_P_omega_kp ...
340 : !> \param cfm_mat_RI_global_work ...
341 : !> \param cfm_mat_Q ...
342 : !> \param ikp ...
343 : !> \param nkp ...
344 : !> \param ikp_local ...
345 : !> \param para_env ...
346 : ! **************************************************************************************************
347 2568 : SUBROUTINE mat_P_to_subgroup(mat_P_omega_kp, cfm_mat_RI_global_work, cfm_mat_Q, ikp, nkp, ikp_local, para_env)
348 :
349 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: mat_P_omega_kp
350 : TYPE(cp_cfm_type), INTENT(INOUT) :: cfm_mat_RI_global_work, cfm_mat_Q
351 : INTEGER, INTENT(IN) :: ikp, nkp, ikp_local
352 : TYPE(mp_para_env_type), POINTER :: para_env
353 :
354 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mat_P_to_subgroup'
355 :
356 : INTEGER :: handle, jkp
357 : TYPE(cp_cfm_type) :: cfm_dummy
358 : TYPE(cp_fm_type) :: fm_mat_work
359 : TYPE(dbcsr_type), POINTER :: mat_P_omega_im, mat_P_omega_re
360 :
361 2568 : CALL timeset(routineN, handle)
362 :
363 2568 : CALL cp_fm_create(fm_mat_work, cfm_mat_RI_global_work%matrix_struct)
364 :
365 67608 : DO jkp = 1, nkp
366 65040 : IF (ikp_local == -1 .AND. jkp /= ikp) CYCLE
367 :
368 2568 : mat_P_omega_re => mat_P_omega_kp(1, jkp)%matrix
369 2568 : CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
370 2568 : CALL copy_dbcsr_to_fm(mat_P_omega_re, fm_mat_work)
371 2568 : CALL cp_fm_to_cfm(msourcer=fm_mat_work, mtarget=cfm_mat_RI_global_work)
372 :
373 2568 : mat_P_omega_im => mat_P_omega_kp(2, jkp)%matrix
374 2568 : CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
375 2568 : CALL copy_dbcsr_to_fm(mat_P_omega_im, fm_mat_work)
376 2568 : CALL cp_cfm_scale_and_add_fm(z_one, cfm_mat_RI_global_work, gaussi, fm_mat_work)
377 :
378 2568 : IF (ikp_local /= -1) CALL para_env%sync()
379 :
380 2568 : IF (ikp_local == -1 .OR. ikp_local == jkp) THEN
381 2568 : CALL cp_cfm_copy_general(cfm_mat_RI_global_work, cfm_mat_Q, para_env)
382 : ELSE
383 0 : CALL cp_cfm_copy_general(cfm_mat_RI_global_work, cfm_dummy, para_env)
384 : END IF
385 :
386 5136 : IF (ikp_local /= -1) CALL para_env%sync()
387 : END DO
388 :
389 2568 : CALL cp_fm_release(fm_mat_work)
390 2568 : CALL para_env%sync()
391 :
392 2568 : CALL timestop(handle)
393 :
394 2568 : END SUBROUTINE mat_P_to_subgroup
395 :
396 : ! **************************************************************************************************
397 : !> \brief ...
398 : !> \param cfm_mat_Q ...
399 : !> \param diag_Qomega ...
400 : !> \param dimen_RI ...
401 : ! **************************************************************************************************
402 2568 : SUBROUTINE cholesky_decomp_Q(cfm_mat_Q, diag_Qomega, dimen_RI)
403 :
404 : TYPE(cp_cfm_type), INTENT(INOUT) :: cfm_mat_Q
405 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: diag_Qomega
406 : INTEGER, INTENT(IN) :: dimen_RI
407 :
408 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cholesky_decomp_Q'
409 :
410 : INTEGER :: handle, info_chol
411 5136 : COMPLEX(KIND=dp), DIMENSION(dimen_RI) :: diag_Q
412 :
413 2568 : CALL timeset(routineN, handle)
414 :
415 : ! epsilon = I + Q; the Hermitian Q has a real diagonal for the RPA energy.
416 2568 : CALL cp_cfm_get_diag(cfm_mat_Q, diag_Q, n_active=dimen_RI)
417 182640 : diag_Qomega = REAL(diag_Q, KIND=dp)
418 2568 : CALL cp_cfm_add_on_diag(cfm_mat_Q, z_one, n_active=dimen_RI)
419 :
420 2568 : CALL cp_cfm_cholesky_decompose(matrix=cfm_mat_Q, n=dimen_RI, info_out=info_chol)
421 :
422 2568 : CPASSERT(info_chol == 0)
423 :
424 2568 : CALL timestop(handle)
425 :
426 2568 : END SUBROUTINE cholesky_decomp_Q
427 :
428 : ! **************************************************************************************************
429 : !> \brief ...
430 : !> \param Erpa ...
431 : !> \param cfm_mat_Q ...
432 : !> \param para_env ...
433 : !> \param diag_Qomega ...
434 : !> \param dimen_RI ...
435 : !> \param freq_weight ...
436 : !> \param kp_weight ...
437 : ! **************************************************************************************************
438 2568 : SUBROUTINE frequency_and_kpoint_integration(Erpa, cfm_mat_Q, para_env, diag_Qomega, &
439 : dimen_RI, freq_weight, kp_weight)
440 :
441 : REAL(KIND=dp), INTENT(INOUT) :: Erpa
442 : TYPE(cp_cfm_type), INTENT(IN) :: cfm_mat_Q
443 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
444 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: diag_Qomega
445 : INTEGER, INTENT(IN) :: dimen_RI
446 : REAL(KIND=dp), INTENT(IN) :: freq_weight, kp_weight
447 :
448 : CHARACTER(LEN=*), PARAMETER :: routineN = 'frequency_and_kpoint_integration'
449 :
450 : INTEGER :: handle, iiB
451 5136 : COMPLEX(KIND=dp), DIMENSION(dimen_RI) :: diag_chol
452 : REAL(KIND=dp) :: FComega
453 2568 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Q_log
454 :
455 2568 : CALL timeset(routineN, handle)
456 :
457 7704 : ALLOCATE (Q_log(dimen_RI))
458 : ! epsilon = U^H U; log(det(epsilon)) = 2*sum(log(real(U_ii))).
459 2568 : CALL cp_cfm_get_diag(cfm_mat_Q, diag_chol, n_active=dimen_RI)
460 182640 : Q_log(:) = 2.0_dp*LOG(REAL(diag_chol, KIND=dp))
461 :
462 : FComega = 0.0_dp
463 182640 : DO iiB = 1, dimen_RI
464 180072 : IF (MODULO(iiB, para_env%num_pe) /= para_env%mepos) CYCLE
465 182640 : FComega = FComega + (Q_log(iiB) - diag_Qomega(iiB))/2.0_dp
466 : END DO
467 :
468 2568 : Erpa = Erpa + FComega*freq_weight*kp_weight
469 :
470 2568 : DEALLOCATE (Q_log)
471 :
472 2568 : CALL timestop(handle)
473 :
474 2568 : END SUBROUTINE frequency_and_kpoint_integration
475 :
476 : ! **************************************************************************************************
477 : !> \brief Split a real Gamma matrix among nearest images using normalized equal weights.
478 : !> \param mat_P_omega ...
479 : !> \param qs_env ...
480 : !> \param kpoints ...
481 : !> \param jquad ...
482 : !> \param unit_nr ...
483 : ! **************************************************************************************************
484 210 : SUBROUTINE get_mat_cell_T_from_mat_gamma(mat_P_omega, qs_env, kpoints, jquad, unit_nr)
485 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: mat_P_omega
486 : TYPE(qs_environment_type), POINTER :: qs_env
487 : TYPE(kpoint_type), POINTER :: kpoints
488 : INTEGER, INTENT(IN) :: jquad, unit_nr
489 :
490 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mat_cell_T_from_mat_gamma'
491 :
492 : INTEGER :: col, handle, i_cell, num_cells_P, row
493 : INTEGER, DIMENSION(3) :: periodic
494 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_P
495 : LOGICAL :: found
496 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: weights
497 210 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block, image_block
498 : TYPE(cell_type), POINTER :: cell
499 : TYPE(dbcsr_iterator_type) :: iter
500 210 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
501 :
502 210 : CALL timeset(routineN, handle)
503 :
504 210 : NULLIFY (cell, particle_set)
505 210 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
506 210 : CALL get_cell(cell=cell, periodic=periodic)
507 :
508 : ! Nearest-neighbour candidate cells depend on periodicity, not on the target k-point mesh.
509 840 : CALL init_cell_index_rpa(2*periodic + 1, kpoints%cell_to_index, kpoints%index_to_cell, cell)
510 210 : index_to_cell_P => kpoints%index_to_cell
511 210 : num_cells_P = SIZE(index_to_cell_P, 2)
512 210 : CPASSERT(SIZE(mat_P_omega) >= num_cells_P)
513 630 : ALLOCATE (weights(num_cells_P))
514 :
515 1890 : DO i_cell = 2, num_cells_P
516 1890 : CALL dbcsr_copy(mat_P_omega(i_cell)%matrix, mat_P_omega(1)%matrix)
517 : END DO
518 :
519 210 : IF (jquad == 1 .AND. unit_nr > 0) THEN
520 8 : WRITE (unit_nr, '(T3,A,T66,ES15.2)') 'GW_INFO| RI regularization parameter: ', &
521 16 : qs_env%mp2_env%ri_rpa_im_time%regularization_RI
522 8 : WRITE (unit_nr, '(T3,A,T66,ES15.2)') 'GW_INFO| eps_eigval_S: ', &
523 16 : qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S
524 8 : IF (qs_env%mp2_env%ri_rpa_im_time%make_chi_pos_definite) THEN
525 : WRITE (unit_nr, '(T3,A,T81)') &
526 8 : 'GW_INFO| Make chi(iw,k) positive definite? TRUE'
527 : ELSE
528 : WRITE (unit_nr, '(T3,A,T81)') &
529 0 : 'GW_INFO| Make chi(iw,k) positive definite? FALSE'
530 : END IF
531 : END IF
532 :
533 210 : CALL dbcsr_iterator_start(iter, mat_P_omega(1)%matrix)
534 1065 : DO WHILE (dbcsr_iterator_blocks_left(iter))
535 855 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
536 : CALL minimum_image_weights(particle_set(row)%r, particle_set(col)%r, &
537 855 : cell, index_to_cell_P, weights)
538 : ! Keep the Gamma block unchanged until its contributions to all other images are copied.
539 8760 : DO i_cell = num_cells_P, 1, -1
540 : CALL dbcsr_get_block_p(matrix=mat_P_omega(i_cell)%matrix, row=row, col=col, &
541 7695 : block=image_block, found=found)
542 7695 : CPASSERT(found)
543 5344992 : image_block = weights(i_cell)*data_block
544 : END DO
545 : END DO
546 210 : CALL dbcsr_iterator_stop(iter)
547 :
548 210 : CALL timestop(handle)
549 :
550 630 : END SUBROUTINE get_mat_cell_T_from_mat_gamma
551 :
552 : ! **************************************************************************************************
553 : !> \brief ...
554 : !> \param mat_P_omega ...
555 : !> \param mat_P_omega_kp ...
556 : !> \param kpoints ...
557 : !> \param eps_filter_im_time ...
558 : !> \param jquad ...
559 : ! **************************************************************************************************
560 132 : SUBROUTINE transform_P_from_real_space_to_kpoints(mat_P_omega, mat_P_omega_kp, &
561 : kpoints, eps_filter_im_time, jquad)
562 :
563 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: mat_P_omega, mat_P_omega_kp
564 : TYPE(kpoint_type), POINTER :: kpoints
565 : REAL(kind=dp), INTENT(IN) :: eps_filter_im_time
566 : INTEGER, INTENT(IN) :: jquad
567 :
568 : CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_P_from_real_space_to_kpoints'
569 :
570 : INTEGER :: handle, icell, nkp, num_integ_points
571 :
572 132 : CALL timeset(routineN, handle)
573 :
574 132 : num_integ_points = SIZE(mat_P_omega, 1)
575 132 : nkp = SIZE(mat_P_omega, 2)
576 :
577 : CALL real_space_to_kpoint_transform_rpa(mat_P_omega_kp(1, :), mat_P_omega_kp(2, :), mat_P_omega(jquad, :), &
578 132 : kpoints, eps_filter_im_time)
579 :
580 2700 : DO icell = 1, SIZE(mat_P_omega, 2)
581 2568 : CALL dbcsr_set(mat_P_omega(jquad, icell)%matrix, 0.0_dp)
582 2700 : CALL dbcsr_filter(mat_P_omega(jquad, icell)%matrix, 1.0_dp)
583 : END DO
584 :
585 132 : CALL timestop(handle)
586 :
587 132 : END SUBROUTINE transform_P_from_real_space_to_kpoints
588 :
589 : ! **************************************************************************************************
590 : !> \brief ...
591 : !> \param real_mat_kp ...
592 : !> \param imag_mat_kp ...
593 : !> \param mat_real_space ...
594 : !> \param kpoints ...
595 : !> \param eps_filter_im_time ...
596 : !> \param real_mat_real_space ...
597 : ! **************************************************************************************************
598 246 : SUBROUTINE real_space_to_kpoint_transform_rpa(real_mat_kp, imag_mat_kp, mat_real_space, &
599 : kpoints, eps_filter_im_time, real_mat_real_space)
600 :
601 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: real_mat_kp, imag_mat_kp, mat_real_space
602 : TYPE(kpoint_type), POINTER :: kpoints
603 : REAL(KIND=dp), INTENT(IN) :: eps_filter_im_time
604 : LOGICAL, INTENT(IN), OPTIONAL :: real_mat_real_space
605 :
606 : CHARACTER(LEN=*), PARAMETER :: routineN = 'real_space_to_kpoint_transform_rpa'
607 :
608 : INTEGER :: handle, i_cell, ik, nkp, num_cells
609 : INTEGER, DIMENSION(3) :: cell
610 246 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
611 : LOGICAL :: my_real_mat_real_space
612 : REAL(KIND=dp) :: arg, coskl, sinkl
613 246 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
614 : TYPE(dbcsr_type) :: mat_work
615 :
616 246 : CALL timeset(routineN, handle)
617 :
618 246 : my_real_mat_real_space = .TRUE.
619 246 : IF (PRESENT(real_mat_real_space)) my_real_mat_real_space = real_mat_real_space
620 :
621 : CALL dbcsr_create(matrix=mat_work, &
622 : template=real_mat_kp(1)%matrix, &
623 246 : matrix_type=dbcsr_type_no_symmetry)
624 246 : CALL dbcsr_reserve_all_blocks(mat_work)
625 246 : CALL dbcsr_set(mat_work, 0.0_dp)
626 :
627 : ! this kpoint environme t should be the kpoints for D(it) and X(it) created in init_cell_index_rpa
628 246 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp)
629 :
630 246 : NULLIFY (index_to_cell)
631 246 : index_to_cell => kpoints%index_to_cell
632 :
633 246 : num_cells = SIZE(index_to_cell, 2)
634 :
635 246 : CPASSERT(SIZE(mat_real_space) >= num_cells/2 + 1)
636 :
637 3408 : DO ik = 1, nkp
638 :
639 3162 : CALL dbcsr_reserve_all_blocks(real_mat_kp(ik)%matrix)
640 3162 : CALL dbcsr_reserve_all_blocks(imag_mat_kp(ik)%matrix)
641 :
642 3162 : CALL dbcsr_set(real_mat_kp(ik)%matrix, 0.0_dp)
643 3162 : CALL dbcsr_set(imag_mat_kp(ik)%matrix, 0.0_dp)
644 :
645 18756 : DO i_cell = 1, num_cells/2 + 1
646 :
647 62376 : cell(:) = index_to_cell(:, i_cell)
648 :
649 15594 : arg = REAL(cell(1), dp)*xkp(1, ik) + REAL(cell(2), dp)*xkp(2, ik) + REAL(cell(3), dp)*xkp(3, ik)
650 15594 : coskl = COS(twopi*arg)
651 15594 : sinkl = SIN(twopi*arg)
652 :
653 15594 : IF (my_real_mat_real_space) THEN
654 15354 : CALL dbcsr_add_local(real_mat_kp(ik)%matrix, mat_real_space(i_cell)%matrix, 1.0_dp, coskl)
655 15354 : CALL dbcsr_add_local(imag_mat_kp(ik)%matrix, mat_real_space(i_cell)%matrix, 1.0_dp, sinkl)
656 : ELSE
657 240 : CALL dbcsr_add_local(real_mat_kp(ik)%matrix, mat_real_space(i_cell)%matrix, 1.0_dp, -sinkl)
658 240 : CALL dbcsr_add_local(imag_mat_kp(ik)%matrix, mat_real_space(i_cell)%matrix, 1.0_dp, coskl)
659 : END IF
660 :
661 18756 : IF (.NOT. (cell(1) == 0 .AND. cell(2) == 0 .AND. cell(3) == 0)) THEN
662 :
663 12432 : CALL dbcsr_transposed(mat_work, mat_real_space(i_cell)%matrix)
664 :
665 12432 : IF (my_real_mat_real_space) THEN
666 12240 : CALL dbcsr_add_local(real_mat_kp(ik)%matrix, mat_work, 1.0_dp, coskl)
667 12240 : CALL dbcsr_add_local(imag_mat_kp(ik)%matrix, mat_work, 1.0_dp, -sinkl)
668 : ELSE
669 : ! for an imaginary real-space matrix, we need to consider the imaginary unit
670 : ! and we need to take into account that the transposed gives an extra "-" sign
671 : ! because the transposed is actually Hermitian conjugate
672 192 : CALL dbcsr_add_local(real_mat_kp(ik)%matrix, mat_work, 1.0_dp, -sinkl)
673 192 : CALL dbcsr_add_local(imag_mat_kp(ik)%matrix, mat_work, 1.0_dp, -coskl)
674 : END IF
675 :
676 12432 : CALL dbcsr_set(mat_work, 0.0_dp)
677 :
678 : END IF
679 :
680 : END DO
681 :
682 3162 : CALL dbcsr_filter(real_mat_kp(ik)%matrix, eps_filter_im_time)
683 3408 : CALL dbcsr_filter(imag_mat_kp(ik)%matrix, eps_filter_im_time)
684 :
685 : END DO
686 :
687 246 : CALL dbcsr_release(mat_work)
688 :
689 246 : CALL timestop(handle)
690 :
691 246 : END SUBROUTINE real_space_to_kpoint_transform_rpa
692 :
693 : ! **************************************************************************************************
694 : !> \brief ...
695 : !> \param mat_a ...
696 : !> \param mat_b ...
697 : !> \param alpha ...
698 : !> \param beta ...
699 : ! **************************************************************************************************
700 56052 : SUBROUTINE dbcsr_add_local(mat_a, mat_b, alpha, beta)
701 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_a, mat_b
702 : REAL(kind=dp), INTENT(IN) :: alpha, beta
703 :
704 : INTEGER :: col, row
705 : LOGICAL :: found
706 56052 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_to_compute, data_block
707 : TYPE(dbcsr_iterator_type) :: iter
708 :
709 56052 : CALL dbcsr_iterator_start(iter, mat_b)
710 289170 : DO WHILE (dbcsr_iterator_blocks_left(iter))
711 233118 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
712 :
713 233118 : NULLIFY (block_to_compute)
714 : CALL dbcsr_get_block_p(matrix=mat_a, &
715 233118 : row=row, col=col, block=block_to_compute, found=found)
716 :
717 233118 : CPASSERT(found)
718 :
719 259745328 : block_to_compute(:, :) = alpha*block_to_compute(:, :) + beta*data_block(:, :)
720 :
721 : END DO
722 56052 : CALL dbcsr_iterator_stop(iter)
723 :
724 56052 : END SUBROUTINE dbcsr_add_local
725 :
726 : ! **************************************************************************************************
727 : !> \brief ...
728 : !> \param fm_mat_W_tau ...
729 : !> \param cfm_mat_Q ...
730 : !> \param cfm_mat_L ...
731 : !> \param dimen_RI ...
732 : !> \param jquad ...
733 : !> \param ikp ...
734 : !> \param grid ...
735 : !> \param ikp_local ...
736 : !> \param para_env ...
737 : !> \param kpoints ...
738 : !> \param qs_env ...
739 : !> \param wkp_W ...
740 : ! **************************************************************************************************
741 2496 : SUBROUTINE compute_Wc_real_space_tau_GW(fm_mat_W_tau, cfm_mat_Q, cfm_mat_L, &
742 : dimen_RI, jquad, &
743 : ikp, grid, ikp_local, &
744 2496 : para_env, kpoints, qs_env, wkp_W)
745 :
746 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_W_tau
747 : TYPE(cp_cfm_type), INTENT(INOUT) :: cfm_mat_Q
748 : TYPE(cp_cfm_type), INTENT(IN) :: cfm_mat_L
749 : INTEGER, INTENT(IN) :: dimen_RI, jquad, ikp
750 : TYPE(time_frequency_grid_type), INTENT(IN) :: grid
751 : INTEGER, INTENT(IN) :: ikp_local
752 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
753 : TYPE(kpoint_type), INTENT(IN), POINTER :: kpoints
754 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
755 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wkp_W
756 :
757 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Wc_real_space_tau_GW'
758 :
759 : INTEGER :: handle, handle2, iatom, iatom_old, iquad, irow, jatom, jatom_old, jcol, jkp, &
760 : ncol_local, nkp, nrow_local, num_cells, num_integ_points
761 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_RI_index
762 2496 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
763 2496 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
764 : REAL(KIND=dp) :: contribution, weight, weight_im, &
765 : weight_re
766 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
767 2496 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
768 2496 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
769 : TYPE(cell_type), POINTER :: cell
770 : TYPE(cp_cfm_type) :: cfm_mat_work_2
771 : TYPE(cp_fm_type) :: fm_dummy, fm_mat_work_global, &
772 : fm_mat_work_local
773 2496 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
774 :
775 2496 : CALL timeset(routineN, handle)
776 :
777 2496 : num_integ_points = SIZE(grid%imaginary_time)
778 :
779 2496 : CALL timeset(routineN//"_1", handle2)
780 :
781 2496 : CALL cp_cfm_create(cfm_mat_work_2, cfm_mat_Q%matrix_struct)
782 2496 : CALL cp_cfm_set_all(cfm_mat_work_2, z_zero)
783 :
784 2496 : CALL cp_fm_create(fm_mat_work_global, fm_mat_W_tau(1)%matrix_struct)
785 2496 : CALL cp_fm_set_all(fm_mat_work_global, 0.0_dp)
786 :
787 2496 : CALL cp_fm_create(fm_mat_work_local, cfm_mat_Q%matrix_struct)
788 2496 : CALL cp_fm_set_all(fm_mat_work_local, 0.0_dp)
789 :
790 2496 : CALL timestop(handle2)
791 :
792 2496 : CALL timeset(routineN//"_2", handle2)
793 :
794 2496 : CALL cfm_screened_interaction_from_eps_cholesky(cfm_mat_Q, cfm_mat_L, cfm_mat_work_2, dimen_RI)
795 :
796 2496 : CALL timestop(handle2)
797 :
798 2496 : CALL timeset(routineN//"_4", handle2)
799 :
800 2496 : CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp, nkp=nkp)
801 2496 : index_to_cell => kpoints%index_to_cell
802 2496 : num_cells = SIZE(index_to_cell, 2)
803 :
804 7488 : ALLOCATE (atom_from_RI_index(dimen_RI))
805 :
806 2496 : CALL get_atom_index_from_basis_function_index(qs_env, atom_from_RI_index, dimen_RI, "RI_AUX")
807 :
808 2496 : NULLIFY (cell, particle_set)
809 2496 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
810 2496 : CALL get_cell(cell=cell, h=hmat)
811 2496 : iatom_old = 0
812 2496 : jatom_old = 0
813 :
814 : CALL cp_cfm_get_info(matrix=cfm_mat_Q, &
815 : nrow_local=nrow_local, &
816 : ncol_local=ncol_local, &
817 : row_indices=row_indices, &
818 2496 : col_indices=col_indices)
819 :
820 89544 : DO irow = 1, nrow_local
821 6863376 : DO jcol = 1, ncol_local
822 :
823 6773832 : iatom = atom_from_RI_index(row_indices(irow))
824 6773832 : jatom = atom_from_RI_index(col_indices(jcol))
825 :
826 6773832 : IF (iatom /= iatom_old .OR. jatom /= jatom_old) THEN
827 :
828 : ! symmetrize=.FALSE. necessary since we already have a symmetrized index_to_cell
829 : CALL compute_weight_re_im(weight_re, weight_im, &
830 : num_cells, iatom, jatom, xkp(1:3, ikp), wkp_W(ikp), &
831 258336 : cell, index_to_cell, hmat, particle_set)
832 :
833 258336 : iatom_old = iatom
834 258336 : jatom_old = jatom
835 :
836 : END IF
837 :
838 : contribution = weight_re*REAL(cfm_mat_work_2%local_data(irow, jcol)) + &
839 6773832 : weight_im*AIMAG(cfm_mat_work_2%local_data(irow, jcol))
840 :
841 6860880 : fm_mat_work_local%local_data(irow, jcol) = fm_mat_work_local%local_data(irow, jcol) + contribution
842 :
843 : END DO
844 : END DO
845 :
846 2496 : CALL timestop(handle2)
847 :
848 2496 : CALL timeset(routineN//"_5", handle2)
849 :
850 2496 : IF (ikp_local == -1) THEN
851 :
852 2496 : CALL cp_fm_copy_general(fm_mat_work_local, fm_mat_work_global, para_env)
853 :
854 17472 : DO iquad = 1, num_integ_points
855 :
856 14976 : weight = grid%cosine_frequency_to_time(iquad, jquad)
857 :
858 14976 : IF (jquad == 1 .AND. ikp == 1) THEN
859 96 : CALL cp_fm_set_all(matrix=fm_mat_W_tau(iquad), alpha=0.0_dp)
860 : END IF
861 :
862 17472 : CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_mat_W_tau(iquad), beta=weight, matrix_b=fm_mat_work_global)
863 :
864 : END DO
865 :
866 : ELSE
867 :
868 0 : DO jkp = 1, nkp
869 :
870 0 : CALL para_env%sync()
871 :
872 0 : IF (ikp_local == jkp) THEN
873 0 : CALL cp_fm_copy_general(fm_mat_work_local, fm_mat_work_global, para_env)
874 : ELSE
875 0 : CALL cp_fm_copy_general(fm_dummy, fm_mat_work_global, para_env)
876 : END IF
877 :
878 0 : CALL para_env%sync()
879 :
880 0 : DO iquad = 1, num_integ_points
881 :
882 0 : weight = grid%cosine_frequency_to_time(iquad, jquad)
883 :
884 0 : IF (jquad == 1 .AND. jkp == 1) THEN
885 0 : CALL cp_fm_set_all(matrix=fm_mat_W_tau(iquad), alpha=0.0_dp)
886 : END IF
887 :
888 : CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_mat_W_tau(iquad), beta=weight, &
889 0 : matrix_b=fm_mat_work_global)
890 :
891 : END DO
892 :
893 : END DO
894 :
895 : END IF
896 :
897 2496 : CALL cp_cfm_release(cfm_mat_work_2)
898 2496 : CALL cp_fm_release(fm_mat_work_global)
899 2496 : CALL cp_fm_release(fm_mat_work_local)
900 :
901 2496 : DEALLOCATE (atom_from_RI_index)
902 :
903 2496 : CALL timestop(handle2)
904 :
905 2496 : CALL timestop(handle)
906 :
907 17472 : END SUBROUTINE compute_Wc_real_space_tau_GW
908 :
909 : ! **************************************************************************************************
910 : !> \brief ...
911 : !> \param fm_mat_W ...
912 : !> \param fm_matrix_Minv ...
913 : !> \param para_env ...
914 : !> \param dimen_RI ...
915 : !> \param num_integ_points ...
916 : ! **************************************************************************************************
917 16 : SUBROUTINE Wc_to_Minv_Wc_Minv(fm_mat_W, fm_matrix_Minv, para_env, dimen_RI, num_integ_points)
918 : TYPE(cp_fm_type), DIMENSION(:) :: fm_mat_W
919 : TYPE(cp_fm_type), DIMENSION(:, :) :: fm_matrix_Minv
920 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
921 : INTEGER :: dimen_RI, num_integ_points
922 :
923 : CHARACTER(LEN=*), PARAMETER :: routineN = 'Wc_to_Minv_Wc_Minv'
924 :
925 : INTEGER :: handle, jquad
926 : TYPE(cp_fm_type) :: fm_work_Minv, fm_work_Minv_W
927 :
928 16 : CALL timeset(routineN, handle)
929 :
930 16 : CALL cp_fm_create(fm_work_Minv, fm_mat_W(1)%matrix_struct)
931 16 : CALL cp_fm_copy_general(fm_matrix_Minv(1, 1), fm_work_Minv, para_env)
932 :
933 16 : CALL cp_fm_create(fm_work_Minv_W, fm_mat_W(1)%matrix_struct)
934 :
935 112 : DO jquad = 1, num_integ_points
936 :
937 : CALL parallel_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_work_Minv, fm_mat_W(jquad), &
938 96 : 0.0_dp, fm_work_Minv_W)
939 : CALL parallel_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_work_Minv_W, fm_work_Minv, &
940 112 : 0.0_dp, fm_mat_W(jquad))
941 :
942 : END DO
943 :
944 16 : CALL cp_fm_release(fm_work_Minv)
945 :
946 16 : CALL cp_fm_release(fm_work_Minv_W)
947 :
948 16 : CALL timestop(handle)
949 :
950 16 : END SUBROUTINE Wc_to_Minv_Wc_Minv
951 :
952 : ! **************************************************************************************************
953 : !> \brief ...
954 : !> \param qs_env ...
955 : !> \param wkp_W ...
956 : !> \param wkp_V ...
957 : !> \param kpoints ...
958 : !> \param h_inv ...
959 : !> \param periodic ...
960 : ! **************************************************************************************************
961 22 : SUBROUTINE compute_wkp_W(qs_env, wkp_W, wkp_V, kpoints, h_inv, periodic)
962 :
963 : TYPE(qs_environment_type), POINTER :: qs_env
964 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
965 : INTENT(OUT) :: wkp_W, wkp_V
966 : TYPE(kpoint_type), POINTER :: kpoints
967 : REAL(KIND=dp), DIMENSION(3, 3) :: h_inv
968 : INTEGER, DIMENSION(3) :: periodic
969 :
970 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_wkp_W'
971 :
972 : INTEGER :: handle, i_x, ikp, info, j_y, k_z, &
973 : kpoint_weights_W_method, n_x, n_y, &
974 : n_z, nkp, nsuperfine, num_lin_eqs
975 : REAL(KIND=dp) :: exp_kpoints, integral, k_sq, weight
976 : REAL(KIND=dp), DIMENSION(3) :: k_vec, x_vec
977 22 : REAL(KIND=dp), DIMENSION(:), POINTER :: right_side, wkp, wkp_tmp
978 22 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: matrix_lin_eqs, xkp
979 :
980 22 : CALL timeset(routineN, handle)
981 :
982 22 : kpoint_weights_W_method = qs_env%mp2_env%ri_rpa_im_time%kpoint_weights_W_method
983 :
984 22 : CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp, nkp=nkp)
985 :
986 : ! we determine the kpoint weights of the Monkhors Pack mesh new
987 : ! such that the functions 1/k^2, 1/k and const are integrated exactly
988 : ! in the Brillouin zone
989 : ! this is done by minimizing sum_i |w_i|^2 where w_i are the weights of
990 : ! the i-th kpoint under the following constraints:
991 : ! 1) 1/k^2, 1/k and const are integrated exactly
992 : ! 2) the kpoint weights of kpoints with identical absolute value are
993 : ! the same, of e.g. (1/8,3/8,3/8) same weight as for (3/8,1/8,3/8)
994 : ! for 1d and 2d materials: we use ordinary Monkhorst-Pack weights, checked
995 : ! by SUM(periodic) == 3
996 88 : ALLOCATE (wkp_V(nkp), wkp_W(nkp))
997 :
998 : ! for exchange part of self-energy, we use truncated Coulomb operator that should be fine
999 : ! with uniform weights (without k-point extrapolation)
1000 22 : IF (ALLOCATED(qs_env%mp2_env%ri_rpa_im_time%wkp_V)) THEN
1001 432 : wkp_V(:) = qs_env%mp2_env%ri_rpa_im_time%wkp_V(:)
1002 : ELSE
1003 18 : wkp_V(:) = wkp(:)
1004 : END IF
1005 :
1006 22 : IF (kpoint_weights_W_method == kp_weights_W_uniform) THEN
1007 :
1008 : ! in the k-point weights wkp, there might be k-point extrapolation included
1009 450 : wkp_W(:) = wkp(:)
1010 :
1011 0 : ELSE IF (kpoint_weights_W_method == kp_weights_W_tailored .OR. &
1012 : kpoint_weights_W_method == kp_weights_W_auto) THEN
1013 :
1014 0 : IF (kpoint_weights_W_method == kp_weights_W_tailored) THEN
1015 0 : exp_kpoints = qs_env%mp2_env%ri_rpa_im_time%exp_tailored_weights
1016 : END IF
1017 :
1018 0 : IF (kpoint_weights_W_method == kp_weights_W_auto) THEN
1019 0 : IF (SUM(periodic) == 2) exp_kpoints = -1.0_dp
1020 : END IF
1021 :
1022 : ! first, compute the integral of f(k)=1/k^2 and 1/k on super fine grid
1023 0 : nsuperfine = 500
1024 0 : integral = 0.0_dp
1025 :
1026 0 : IF (periodic(1) == 1) THEN
1027 : n_x = nsuperfine
1028 : ELSE
1029 0 : n_x = 1
1030 : END IF
1031 0 : IF (periodic(2) == 1) THEN
1032 : n_y = nsuperfine
1033 : ELSE
1034 0 : n_y = 1
1035 : END IF
1036 0 : IF (periodic(3) == 1) THEN
1037 : n_z = nsuperfine
1038 : ELSE
1039 0 : n_z = 1
1040 : END IF
1041 :
1042 : ! actually, there is the factor *det_3x3(h_inv) missing to account for the
1043 : ! integration volume but for wkp det_3x3(h_inv) is needed
1044 0 : weight = 1.0_dp/(REAL(n_x, dp)*REAL(n_y, dp)*REAL(n_z, dp))
1045 0 : DO i_x = 1, n_x
1046 0 : DO j_y = 1, n_y
1047 0 : DO k_z = 1, n_z
1048 :
1049 0 : IF (periodic(1) == 1) THEN
1050 0 : x_vec(1) = (REAL(i_x - nsuperfine/2, dp) - 0.5_dp)/REAL(nsuperfine, dp)
1051 : ELSE
1052 0 : x_vec(1) = 0.0_dp
1053 : END IF
1054 0 : IF (periodic(2) == 1) THEN
1055 0 : x_vec(2) = (REAL(j_y - nsuperfine/2, dp) - 0.5_dp)/REAL(nsuperfine, dp)
1056 : ELSE
1057 0 : x_vec(2) = 0.0_dp
1058 : END IF
1059 0 : IF (periodic(3) == 1) THEN
1060 0 : x_vec(3) = (REAL(k_z - nsuperfine/2, dp) - 0.5_dp)/REAL(nsuperfine, dp)
1061 : ELSE
1062 0 : x_vec(3) = 0.0_dp
1063 : END IF
1064 :
1065 0 : k_vec = MATMUL(h_inv(1:3, 1:3), x_vec)
1066 0 : k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2
1067 0 : integral = integral + weight*k_sq**(exp_kpoints*0.5_dp)
1068 :
1069 : END DO
1070 : END DO
1071 : END DO
1072 :
1073 0 : num_lin_eqs = nkp + 2
1074 :
1075 0 : ALLOCATE (matrix_lin_eqs(num_lin_eqs, num_lin_eqs))
1076 0 : matrix_lin_eqs(:, :) = 0.0_dp
1077 :
1078 0 : DO ikp = 1, nkp
1079 :
1080 0 : k_vec = MATMUL(h_inv(1:3, 1:3), xkp(1:3, ikp))
1081 0 : k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2
1082 :
1083 0 : matrix_lin_eqs(ikp, ikp) = 2.0_dp
1084 0 : matrix_lin_eqs(ikp, nkp + 1) = 1.0_dp
1085 0 : matrix_lin_eqs(nkp + 1, ikp) = 1.0_dp
1086 :
1087 0 : matrix_lin_eqs(ikp, nkp + 2) = k_sq**(exp_kpoints*0.5_dp)
1088 0 : matrix_lin_eqs(nkp + 2, ikp) = k_sq**(exp_kpoints*0.5_dp)
1089 :
1090 : END DO
1091 :
1092 0 : CALL invmat(matrix_lin_eqs, info)
1093 : ! check whether inversion was successful
1094 0 : CPASSERT(info == 0)
1095 :
1096 0 : ALLOCATE (right_side(num_lin_eqs))
1097 0 : right_side = 0.0_dp
1098 0 : right_side(nkp + 1) = 1.0_dp
1099 : ! divide integral by two because CP2K k-mesh already considers symmetry k <-> -k
1100 0 : right_side(nkp + 2) = integral
1101 :
1102 0 : ALLOCATE (wkp_tmp(num_lin_eqs))
1103 :
1104 0 : wkp_tmp(1:num_lin_eqs) = MATMUL(matrix_lin_eqs, right_side)
1105 :
1106 0 : wkp_W(1:nkp) = wkp_tmp(1:nkp)
1107 :
1108 0 : DEALLOCATE (matrix_lin_eqs, right_side, wkp_tmp)
1109 :
1110 : END IF
1111 :
1112 22 : CALL timestop(handle)
1113 :
1114 22 : END SUBROUTINE compute_wkp_W
1115 :
1116 : ! **************************************************************************************************
1117 : !> \brief ...
1118 : !> \param qs_env ...
1119 : !> \param Eigenval_kp ...
1120 : ! **************************************************************************************************
1121 16 : SUBROUTINE get_bandstruc_and_k_dependent_MOs(qs_env, Eigenval_kp)
1122 : TYPE(qs_environment_type), POINTER :: qs_env
1123 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Eigenval_kp
1124 :
1125 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_bandstruc_and_k_dependent_MOs'
1126 :
1127 : INTEGER :: handle, ikp, ispin, nmo, nspins
1128 : INTEGER, DIMENSION(3) :: nkp_grid_G
1129 16 : REAL(KIND=dp), DIMENSION(:), POINTER :: ev
1130 16 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: kpgeneral
1131 : TYPE(kpoint_type), POINTER :: kpoints_Sigma
1132 : TYPE(mp_para_env_type), POINTER :: para_env
1133 :
1134 16 : CALL timeset(routineN, handle)
1135 :
1136 : NULLIFY (qs_env%mp2_env%ri_rpa_im_time%kpoints_G, &
1137 16 : qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma, &
1138 16 : qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma_no_xc, &
1139 16 : para_env)
1140 :
1141 16 : nkp_grid_G(1:3) = [1, 1, 1]
1142 :
1143 16 : CALL get_qs_env(qs_env=qs_env, para_env=para_env)
1144 :
1145 : CALL create_kp_and_calc_kp_orbitals(qs_env, qs_env%mp2_env%ri_rpa_im_time%kpoints_G, &
1146 : "MONKHORST-PACK", para_env%num_pe, &
1147 16 : mp_grid=nkp_grid_G(1:3))
1148 :
1149 16 : IF (qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
1150 :
1151 : ! set up k-points for GW band structure calculation, will be completed later
1152 16 : CALL get_kpgeneral_for_Sigma_kpoints(qs_env, kpgeneral)
1153 :
1154 : CALL create_kp_and_calc_kp_orbitals(qs_env, qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma, &
1155 : "GENERAL", para_env%num_pe, &
1156 16 : kpgeneral=kpgeneral)
1157 :
1158 : CALL create_kp_and_calc_kp_orbitals(qs_env, qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma_no_xc, &
1159 : "GENERAL", para_env%num_pe, &
1160 16 : kpgeneral=kpgeneral, with_xc_terms=.FALSE.)
1161 :
1162 16 : kpoints_Sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma
1163 16 : nmo = SIZE(Eigenval_kp, 1)
1164 16 : nspins = SIZE(Eigenval_kp, 3)
1165 :
1166 48 : ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%Eigenval_Gamma(nmo))
1167 340 : qs_env%mp2_env%ri_rpa_im_time%Eigenval_Gamma(:) = Eigenval_kp(:, 1, 1)
1168 :
1169 16 : DEALLOCATE (Eigenval_kp)
1170 :
1171 80 : ALLOCATE (Eigenval_kp(nmo, kpoints_Sigma%nkp, nspins))
1172 :
1173 136 : DO ikp = 1, kpoints_Sigma%nkp
1174 :
1175 272 : DO ispin = 1, nspins
1176 :
1177 136 : ev => kpoints_Sigma%kp_env(ikp)%kpoint_env%mos(ispin)%eigenvalues
1178 :
1179 3016 : Eigenval_kp(:, ikp, ispin) = ev(:)
1180 :
1181 : END DO
1182 :
1183 : END DO
1184 :
1185 16 : DEALLOCATE (kpgeneral)
1186 :
1187 : END IF
1188 :
1189 16 : CALL release_hfx_stuff(qs_env)
1190 :
1191 16 : CALL timestop(handle)
1192 :
1193 16 : END SUBROUTINE get_bandstruc_and_k_dependent_MOs
1194 :
1195 : ! **************************************************************************************************
1196 : !> \brief releases part of the given qs_env in order to save memory
1197 : !> \param qs_env the object to release
1198 : ! **************************************************************************************************
1199 16 : SUBROUTINE release_hfx_stuff(qs_env)
1200 : TYPE(qs_environment_type), POINTER :: qs_env
1201 :
1202 16 : IF (ASSOCIATED(qs_env%x_data) .AND. .NOT. qs_env%mp2_env%ri_g0w0%do_ri_Sigma_x) THEN
1203 2 : CALL hfx_release(qs_env%x_data)
1204 : END IF
1205 :
1206 16 : END SUBROUTINE release_hfx_stuff
1207 :
1208 : ! **************************************************************************************************
1209 : !> \brief ...
1210 : !> \param qs_env ...
1211 : !> \param kpoints ...
1212 : !> \param scheme ...
1213 : !> \param group_size_ext ...
1214 : !> \param mp_grid ...
1215 : !> \param kpgeneral ...
1216 : !> \param with_xc_terms ...
1217 : !> \param kp_shift ...
1218 : !> \param gamma_centered ...
1219 : ! **************************************************************************************************
1220 336 : SUBROUTINE create_kp_and_calc_kp_orbitals(qs_env, kpoints, scheme, &
1221 48 : group_size_ext, mp_grid, kpgeneral, with_xc_terms, &
1222 : kp_shift, gamma_centered)
1223 :
1224 : TYPE(qs_environment_type), POINTER :: qs_env
1225 : TYPE(kpoint_type), POINTER :: kpoints
1226 : CHARACTER(LEN=*), INTENT(IN) :: scheme
1227 : INTEGER :: group_size_ext
1228 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: mp_grid
1229 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
1230 : OPTIONAL :: kpgeneral
1231 : LOGICAL, OPTIONAL :: with_xc_terms
1232 : REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: kp_shift
1233 : LOGICAL, INTENT(IN), OPTIONAL :: gamma_centered
1234 :
1235 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_kp_and_calc_kp_orbitals'
1236 :
1237 : INTEGER :: handle, i_re_im, ikp, ispin, nkp, nspins
1238 : INTEGER, DIMENSION(3) :: periodic
1239 : LOGICAL :: my_with_xc_terms
1240 48 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1241 : TYPE(cell_type), POINTER :: cell
1242 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1243 : TYPE(cp_cfm_type) :: cksmat, cmos, csmat, cwork
1244 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1245 : TYPE(cp_fm_type) :: fm_work
1246 48 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_s_desymm
1247 48 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_ks_kp, mat_s_kp
1248 : TYPE(dft_control_type), POINTER :: dft_control
1249 : TYPE(kpoint_env_type), POINTER :: kp
1250 : TYPE(mp_para_env_type), POINTER :: para_env
1251 : TYPE(qs_scf_env_type), POINTER :: scf_env
1252 : TYPE(scf_control_type), POINTER :: scf_control
1253 :
1254 48 : CALL timeset(routineN, handle)
1255 :
1256 48 : my_with_xc_terms = .TRUE.
1257 48 : IF (PRESENT(with_xc_terms)) my_with_xc_terms = with_xc_terms
1258 :
1259 : CALL get_qs_env(qs_env, &
1260 : para_env=para_env, &
1261 : blacs_env=blacs_env, &
1262 : matrix_s=matrix_s, &
1263 : scf_env=scf_env, &
1264 : scf_control=scf_control, &
1265 48 : cell=cell)
1266 :
1267 : ! get kpoints
1268 : CALL calculate_kpoints_for_bs(kpoints, scheme, kpgeneral=kpgeneral, mp_grid=mp_grid, &
1269 : group_size_ext=group_size_ext, kp_shift=kp_shift, &
1270 64 : gamma_centered=gamma_centered)
1271 :
1272 48 : CALL kpoint_env_initialize(kpoints, para_env, blacs_env)
1273 :
1274 : ! calculate all MOs that are accessible in the given
1275 : ! Gaussian AO basis, therefore nadd=1E10
1276 48 : CALL kpoint_initialize_mos(kpoints, qs_env%mos, 2000000000)
1277 48 : CALL kpoint_initialize_mo_set(kpoints)
1278 :
1279 48 : CALL get_cell(cell=cell, periodic=periodic)
1280 :
1281 192 : CALL init_cell_index_rpa(2*periodic + 1, kpoints%cell_to_index, kpoints%index_to_cell, cell)
1282 :
1283 : ! get S(k)
1284 48 : CALL get_qs_env(qs_env, matrix_s=matrix_s, scf_env=scf_env, scf_control=scf_control, dft_control=dft_control)
1285 :
1286 48 : NULLIFY (matrix_s_desymm)
1287 48 : CALL dbcsr_allocate_matrix_set(matrix_s_desymm, 1)
1288 48 : ALLOCATE (matrix_s_desymm(1)%matrix)
1289 : CALL dbcsr_create(matrix=matrix_s_desymm(1)%matrix, template=matrix_s(1)%matrix, &
1290 48 : matrix_type=dbcsr_type_no_symmetry)
1291 48 : CALL dbcsr_desymmetrize(matrix_s(1)%matrix, matrix_s_desymm(1)%matrix)
1292 :
1293 48 : CALL mat_kp_from_mat_gamma(qs_env, mat_s_kp, matrix_s_desymm(1)%matrix, kpoints, 1)
1294 :
1295 48 : CALL get_kpoint_info(kpoints, nkp=nkp)
1296 :
1297 48 : matrix_struct => kpoints%kp_env(1)%kpoint_env%wmat(1, 1)%matrix_struct
1298 :
1299 48 : CALL cp_cfm_create(cksmat, matrix_struct)
1300 48 : CALL cp_cfm_create(csmat, matrix_struct)
1301 48 : CALL cp_cfm_create(cmos, matrix_struct)
1302 48 : CALL cp_cfm_create(cwork, matrix_struct)
1303 48 : CALL cp_fm_create(fm_work, matrix_struct)
1304 :
1305 48 : nspins = dft_control%nspins
1306 :
1307 102 : DO ispin = 1, nspins
1308 :
1309 : ! get H(k)
1310 54 : IF (my_with_xc_terms) THEN
1311 36 : CALL mat_kp_from_mat_gamma(qs_env, mat_ks_kp, qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix, kpoints, ispin)
1312 : ELSE
1313 : CALL mat_kp_from_mat_gamma(qs_env, mat_ks_kp, qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix, &
1314 18 : kpoints, ispin)
1315 : END IF
1316 :
1317 392 : DO ikp = 1, nkp
1318 :
1319 290 : CALL copy_dbcsr_to_fm(mat_ks_kp(ikp, 1)%matrix, kpoints%kp_env(ikp)%kpoint_env%wmat(1, ispin))
1320 290 : CALL cp_cfm_scale_and_add_fm(z_zero, cksmat, z_one, kpoints%kp_env(ikp)%kpoint_env%wmat(1, ispin))
1321 :
1322 290 : CALL copy_dbcsr_to_fm(mat_ks_kp(ikp, 2)%matrix, kpoints%kp_env(ikp)%kpoint_env%wmat(2, ispin))
1323 290 : CALL cp_cfm_scale_and_add_fm(z_one, cksmat, gaussi, kpoints%kp_env(ikp)%kpoint_env%wmat(2, ispin))
1324 :
1325 290 : CALL copy_dbcsr_to_fm(mat_s_kp(ikp, 1)%matrix, fm_work)
1326 290 : CALL cp_cfm_scale_and_add_fm(z_zero, csmat, z_one, fm_work)
1327 :
1328 290 : CALL copy_dbcsr_to_fm(mat_s_kp(ikp, 2)%matrix, fm_work)
1329 290 : CALL cp_cfm_scale_and_add_fm(z_one, csmat, gaussi, fm_work)
1330 :
1331 290 : kp => kpoints%kp_env(ikp)%kpoint_env
1332 :
1333 290 : CALL get_mo_set(kp%mos(ispin), eigenvalues=eigenvalues)
1334 :
1335 290 : IF (scf_env%cholesky_method == cholesky_off .OR. &
1336 : qs_env%mp2_env%ri_rpa_im_time%make_overlap_mat_ao_pos_definite) THEN
1337 0 : CALL cp_cfm_geeig_canon(cksmat, csmat, cmos, eigenvalues, cwork, scf_control%eps_eigval)
1338 : ELSE
1339 290 : CALL cp_cfm_geeig(cksmat, csmat, cmos, eigenvalues, cwork)
1340 : END IF
1341 :
1342 290 : CALL cp_cfm_to_cfm(cmos, kp%mos(ispin)%cmo_coeff)
1343 :
1344 12382 : kp%mos(ispin)%eigenvalues = eigenvalues
1345 :
1346 : END DO
1347 :
1348 : END DO
1349 :
1350 304 : DO ikp = 1, nkp
1351 816 : DO i_re_im = 1, 2
1352 768 : CALL dbcsr_deallocate_matrix(mat_ks_kp(ikp, i_re_im)%matrix)
1353 : END DO
1354 : END DO
1355 48 : DEALLOCATE (mat_ks_kp)
1356 :
1357 304 : DO ikp = 1, nkp
1358 816 : DO i_re_im = 1, 2
1359 768 : CALL dbcsr_deallocate_matrix(mat_s_kp(ikp, i_re_im)%matrix)
1360 : END DO
1361 : END DO
1362 48 : DEALLOCATE (mat_s_kp)
1363 :
1364 48 : CALL dbcsr_deallocate_matrix(matrix_s_desymm(1)%matrix)
1365 48 : DEALLOCATE (matrix_s_desymm)
1366 :
1367 48 : CALL cp_cfm_release(cksmat)
1368 48 : CALL cp_cfm_release(csmat)
1369 48 : CALL cp_cfm_release(cwork)
1370 48 : CALL cp_cfm_release(cmos)
1371 48 : CALL cp_fm_release(fm_work)
1372 :
1373 48 : CALL timestop(handle)
1374 :
1375 48 : END SUBROUTINE create_kp_and_calc_kp_orbitals
1376 :
1377 : ! **************************************************************************************************
1378 : !> \brief ...
1379 : !> \param qs_env ...
1380 : !> \param mat_kp ...
1381 : !> \param mat_gamma ...
1382 : !> \param kpoints ...
1383 : !> \param ispin ...
1384 : !> \param real_mat_real_space ...
1385 : ! **************************************************************************************************
1386 114 : SUBROUTINE mat_kp_from_mat_gamma(qs_env, mat_kp, mat_gamma, kpoints, ispin, real_mat_real_space)
1387 :
1388 : TYPE(qs_environment_type), POINTER :: qs_env
1389 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_kp
1390 : TYPE(dbcsr_type) :: mat_gamma
1391 : TYPE(kpoint_type), POINTER :: kpoints
1392 : INTEGER :: ispin
1393 : LOGICAL, INTENT(IN), OPTIONAL :: real_mat_real_space
1394 :
1395 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mat_kp_from_mat_gamma'
1396 :
1397 : INTEGER :: handle, i_cell, i_re_im, ikp, nkp, &
1398 : num_cells
1399 : INTEGER, DIMENSION(3) :: periodic
1400 114 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1401 114 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1402 : TYPE(cell_type), POINTER :: cell
1403 114 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_real_space
1404 :
1405 114 : CALL timeset(routineN, handle)
1406 :
1407 114 : CALL get_qs_env(qs_env, cell=cell)
1408 114 : CALL get_cell(cell=cell, periodic=periodic)
1409 114 : num_cells = 3**(periodic(1) + periodic(2) + periodic(3))
1410 :
1411 114 : NULLIFY (mat_real_space)
1412 114 : CALL dbcsr_allocate_matrix_set(mat_real_space, num_cells)
1413 1140 : DO i_cell = 1, num_cells
1414 1026 : ALLOCATE (mat_real_space(i_cell)%matrix)
1415 : CALL dbcsr_create(matrix=mat_real_space(i_cell)%matrix, &
1416 1026 : template=mat_gamma)
1417 1026 : CALL dbcsr_reserve_all_blocks(mat_real_space(i_cell)%matrix)
1418 1140 : CALL dbcsr_set(mat_real_space(i_cell)%matrix, 0.0_dp)
1419 : END DO
1420 :
1421 114 : CALL dbcsr_copy(mat_real_space(1)%matrix, mat_gamma)
1422 :
1423 114 : CALL get_mat_cell_T_from_mat_gamma(mat_real_space, qs_env, kpoints, 2, 0)
1424 :
1425 114 : NULLIFY (xkp, cell_to_index)
1426 114 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, cell_to_index=cell_to_index)
1427 :
1428 114 : IF (ispin == 1) THEN
1429 108 : NULLIFY (mat_kp)
1430 108 : CALL dbcsr_allocate_matrix_set(mat_kp, nkp, 2)
1431 668 : DO ikp = 1, nkp
1432 1788 : DO i_re_im = 1, 2
1433 1120 : ALLOCATE (mat_kp(ikp, i_re_im)%matrix)
1434 1120 : CALL dbcsr_create(matrix=mat_kp(ikp, i_re_im)%matrix, template=mat_gamma)
1435 1120 : CALL dbcsr_reserve_all_blocks(mat_kp(ikp, i_re_im)%matrix)
1436 1680 : CALL dbcsr_set(mat_kp(ikp, i_re_im)%matrix, 0.0_dp)
1437 : END DO
1438 : END DO
1439 : END IF
1440 :
1441 114 : IF (PRESENT(real_mat_real_space)) THEN
1442 : CALL real_space_to_kpoint_transform_rpa(mat_kp(:, 1), mat_kp(:, 2), mat_real_space, kpoints, 0.0_dp, &
1443 12 : real_mat_real_space)
1444 : ELSE
1445 102 : CALL real_space_to_kpoint_transform_rpa(mat_kp(:, 1), mat_kp(:, 2), mat_real_space, kpoints, 0.0_dp)
1446 : END IF
1447 :
1448 1140 : DO i_cell = 1, num_cells
1449 1140 : CALL dbcsr_deallocate_matrix(mat_real_space(i_cell)%matrix)
1450 : END DO
1451 114 : DEALLOCATE (mat_real_space)
1452 :
1453 114 : CALL timestop(handle)
1454 :
1455 114 : END SUBROUTINE mat_kp_from_mat_gamma
1456 :
1457 : ! **************************************************************************************************
1458 : !> \brief ...
1459 : !> \param qs_env ...
1460 : !> \param kpgeneral ...
1461 : ! **************************************************************************************************
1462 16 : SUBROUTINE get_kpgeneral_for_Sigma_kpoints(qs_env, kpgeneral)
1463 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1464 : REAL(kind=dp), DIMENSION(:, :), POINTER :: kpgeneral
1465 :
1466 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_kpgeneral_for_Sigma_kpoints'
1467 :
1468 : INTEGER :: handle, i_kp_in_kp_line, i_special_kp, &
1469 : i_x, ikk, j_y, k_z, n_kp_in_kp_line, &
1470 : n_special_kp
1471 16 : INTEGER, DIMENSION(:), POINTER :: nkp_grid
1472 :
1473 16 : CALL timeset(routineN, handle)
1474 :
1475 16 : n_special_kp = qs_env%mp2_env%ri_g0w0%n_special_kp
1476 16 : n_kp_in_kp_line = qs_env%mp2_env%ri_g0w0%n_kp_in_kp_line
1477 16 : IF (n_special_kp > 0) THEN
1478 14 : qs_env%mp2_env%ri_g0w0%nkp_self_energy_special_kp = n_kp_in_kp_line*(n_special_kp - 1) + 1
1479 : ELSE
1480 2 : qs_env%mp2_env%ri_g0w0%nkp_self_energy_special_kp = 0
1481 : END IF
1482 :
1483 : qs_env%mp2_env%ri_g0w0%nkp_self_energy_monkh_pack = qs_env%mp2_env%ri_g0w0%kp_grid_Sigma(1)* &
1484 : qs_env%mp2_env%ri_g0w0%kp_grid_Sigma(2)* &
1485 16 : qs_env%mp2_env%ri_g0w0%kp_grid_Sigma(3)
1486 :
1487 : qs_env%mp2_env%ri_g0w0%nkp_self_energy = qs_env%mp2_env%ri_g0w0%nkp_self_energy_special_kp + &
1488 16 : qs_env%mp2_env%ri_g0w0%nkp_self_energy_monkh_pack
1489 :
1490 48 : ALLOCATE (kpgeneral(3, qs_env%mp2_env%ri_g0w0%nkp_self_energy))
1491 :
1492 16 : IF (n_special_kp > 0) THEN
1493 :
1494 112 : kpgeneral(1:3, 1) = qs_env%mp2_env%ri_g0w0%xkp_special_kp(1:3, 1)
1495 :
1496 14 : ikk = 1
1497 :
1498 28 : DO i_special_kp = 2, n_special_kp
1499 70 : DO i_kp_in_kp_line = 1, n_kp_in_kp_line
1500 :
1501 42 : ikk = ikk + 1
1502 : kpgeneral(1:3, ikk) = qs_env%mp2_env%ri_g0w0%xkp_special_kp(1:3, i_special_kp - 1) + &
1503 : REAL(i_kp_in_kp_line, KIND=dp)/REAL(n_kp_in_kp_line, KIND=dp)* &
1504 : (qs_env%mp2_env%ri_g0w0%xkp_special_kp(1:3, i_special_kp) - &
1505 350 : qs_env%mp2_env%ri_g0w0%xkp_special_kp(1:3, i_special_kp - 1))
1506 :
1507 : END DO
1508 : END DO
1509 :
1510 : ELSE
1511 :
1512 : ikk = 0
1513 :
1514 : END IF
1515 :
1516 16 : nkp_grid => qs_env%mp2_env%ri_g0w0%kp_grid_Sigma
1517 :
1518 48 : DO i_x = 1, nkp_grid(1)
1519 112 : DO j_y = 1, nkp_grid(2)
1520 160 : DO k_z = 1, nkp_grid(3)
1521 64 : ikk = ikk + 1
1522 64 : kpgeneral(1, ikk) = REAL(2*i_x - nkp_grid(1) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(1), KIND=dp))
1523 64 : kpgeneral(2, ikk) = REAL(2*j_y - nkp_grid(2) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(2), KIND=dp))
1524 128 : kpgeneral(3, ikk) = REAL(2*k_z - nkp_grid(3) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(3), KIND=dp))
1525 : END DO
1526 : END DO
1527 : END DO
1528 :
1529 16 : CALL timestop(handle)
1530 :
1531 16 : END SUBROUTINE get_kpgeneral_for_Sigma_kpoints
1532 :
1533 0 : END MODULE rpa_gw_kpoints_util
|
