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