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