(git:b279b6b)
rpa_gw_sigma_x.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 to calculate EXX within GW
10 !> \par History
11 !> 07.2020 separated from mp2.F [F. Stein, code by Jan Wilhelm]
12 !> \author Jan Wilhelm, Frederick Stein
13 ! **************************************************************************************************
14 MODULE rpa_gw_sigma_x
15  USE admm_methods, ONLY: admm_mo_merge_ks_matrix
16  USE admm_types, ONLY: admm_type,&
17  get_admm_env
18  USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale_and_add_fm
19  USE cp_cfm_types, ONLY: cp_cfm_create,&
20  cp_cfm_get_info,&
21  cp_cfm_release,&
22  cp_cfm_type
23  USE cp_control_types, ONLY: dft_control_type
24  USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
25  USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
26  copy_fm_to_dbcsr,&
27  dbcsr_allocate_matrix_set,&
28  dbcsr_deallocate_matrix_set
29  USE cp_files, ONLY: close_file,&
30  open_file
31  USE cp_fm_struct, ONLY: cp_fm_struct_type
32  USE cp_fm_types, ONLY: cp_fm_create,&
33  cp_fm_get_info,&
34  cp_fm_release,&
35  cp_fm_type
36  USE dbcsr_api, ONLY: &
37  dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_get_diag, dbcsr_multiply, &
38  dbcsr_p_type, dbcsr_release, dbcsr_release_p, dbcsr_set, dbcsr_type, &
39  dbcsr_type_antisymmetric, dbcsr_type_symmetric
40  USE hfx_energy_potential, ONLY: integrate_four_center
41  USE hfx_exx, ONLY: calc_exx_admm_xc_contributions,&
42  exx_post_hfx,&
43  exx_pre_hfx
44  USE hfx_ri, ONLY: hfx_ri_update_ks
45  USE input_constants, ONLY: do_admm_basis_projection,&
46  do_admm_purify_none,&
47  gw_print_exx,&
48  gw_read_exx,&
49  xc_none
50  USE input_section_types, ONLY: section_vals_get,&
51  section_vals_get_subs_vals,&
52  section_vals_type,&
53  section_vals_val_get,&
54  section_vals_val_set
55  USE kinds, ONLY: dp
56  USE kpoint_methods, ONLY: rskp_transform
57  USE kpoint_types, ONLY: get_kpoint_info,&
58  kpoint_env_type,&
59  kpoint_type
60  USE machine, ONLY: m_walltime
61  USE mathconstants, ONLY: gaussi,&
62  z_one,&
63  z_zero
64  USE message_passing, ONLY: mp_para_env_type
65  USE mp2_integrals, ONLY: compute_kpoints
66  USE mp2_ri_2c, ONLY: setup_trunc_coulomb_pot_for_exchange_self_energy
67  USE mp2_types, ONLY: mp2_type
68  USE parallel_gemm_api, ONLY: parallel_gemm
69  USE physcon, ONLY: evolt
70  USE qs_energy_types, ONLY: qs_energy_type
71  USE qs_environment_types, ONLY: get_qs_env,&
72  qs_environment_type
73  USE qs_ks_methods, ONLY: qs_ks_build_kohn_sham_matrix
74  USE qs_ks_types, ONLY: qs_ks_env_type
75  USE qs_mo_types, ONLY: get_mo_set,&
76  mo_set_type
77  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
78  USE qs_rho_types, ONLY: qs_rho_get,&
79  qs_rho_type
80  USE rpa_gw, ONLY: compute_minus_vxc_kpoints,&
81  trafo_to_mo_and_kpoints
82  USE rpa_gw_kpoints_util, ONLY: get_bandstruc_and_k_dependent_mos
83 
84 !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
85 
86 #include "./base/base_uses.f90"
87 
88  IMPLICIT NONE
89 
90  PRIVATE
91 
92  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_gw_sigma_x'
93 
94  PUBLIC :: compute_vec_sigma_x_minus_vxc_gw
95 
96 CONTAINS
97 
98 ! **************************************************************************************************
99 !> \brief ...
100 !> \param qs_env ...
101 !> \param mp2_env ...
102 !> \param mos_mp2 ...
103 !> \param energy_ex ...
104 !> \param energy_xc_admm ...
105 !> \param t3 ...
106 !> \param unit_nr ...
107 !> \par History
108 !> 04.2015 created
109 !> \author Jan Wilhelm
110 ! **************************************************************************************************
111  SUBROUTINE compute_vec_sigma_x_minus_vxc_gw(qs_env, mp2_env, mos_mp2, energy_ex, energy_xc_admm, t3, unit_nr)
112  TYPE(qs_environment_type), POINTER :: qs_env
113  TYPE(mp2_type) :: mp2_env
114  TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos_mp2
115  REAL(kind=dp), INTENT(OUT) :: energy_ex, energy_xc_admm(2), t3
116  INTEGER, INTENT(IN) :: unit_nr
117 
118  CHARACTER(len=*), PARAMETER :: routinen = 'compute_vec_Sigma_x_minus_vxc_gw'
119 
120  CHARACTER(4) :: occ_virt
121  CHARACTER(LEN=40) :: line
122  INTEGER :: dimen, gw_corr_lev_occ, gw_corr_lev_tot, gw_corr_lev_virt, handle, homo, i_img, &
123  ikp, irep, ispin, iunit, myfun, myfun_aux, myfun_prim, n_level_gw, n_level_gw_ref, &
124  n_rep_hf, nkp, nkp_sigma, nmo, nspins, print_exx
125  LOGICAL :: calc_ints, charge_constrain_tmp, do_admm_rpa, do_hfx, do_kpoints_cubic_rpa, &
126  do_kpoints_from_gamma, do_ri_sigma_x, really_read_line
127  REAL(kind=dp) :: e_gap_gw, e_homo_gw, e_lumo_gw, eh1, ehfx, eigval_dft, eigval_hf_at_dft, &
128  energy_exc, energy_exc1, energy_exc1_aux_fit, energy_exc_aux_fit, energy_total, &
129  exx_minus_vxc, hfx_fraction, min_direct_hf_at_dft_gap, t1, t2, tmp
130  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenval_kp_hf_at_dft, vec_sigma_x
131  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: eigenval_kp, vec_sigma_x_minus_vxc_gw, &
132  vec_sigma_x_minus_vxc_gw_im
133  TYPE(admm_type), POINTER :: admm_env
134  TYPE(cp_fm_type), POINTER :: mo_coeff
135  TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:) :: mat_exchange_for_kp_from_gamma
136  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, &
137  matrix_ks_aux_fit_hfx, rho_ao, &
138  rho_ao_aux_fit
139  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_2d, matrix_ks_kp_im, &
140  matrix_ks_kp_re, matrix_ks_transl, matrix_sigma_x_minus_vxc, matrix_sigma_x_minus_vxc_im, &
141  rho_ao_2d
142  TYPE(dbcsr_type) :: matrix_tmp, matrix_tmp_2, mo_coeff_b
143  TYPE(dft_control_type), POINTER :: dft_control
144  TYPE(kpoint_type), POINTER :: kpoints, kpoints_sigma
145  TYPE(mp_para_env_type), POINTER :: para_env
146  TYPE(qs_energy_type), POINTER :: energy
147  TYPE(qs_ks_env_type), POINTER :: ks_env
148  TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
149  TYPE(section_vals_type), POINTER :: hfx_sections, input, xc_section, &
150  xc_section_admm_aux, &
151  xc_section_admm_prim
152 
153  NULLIFY (admm_env, matrix_ks, matrix_ks_aux_fit, rho_ao, matrix_sigma_x_minus_vxc, input, &
154  xc_section, xc_section_admm_aux, xc_section_admm_prim, hfx_sections, rho, &
155  dft_control, para_env, ks_env, mo_coeff, matrix_sigma_x_minus_vxc_im, matrix_ks_aux_fit_hfx, &
156  rho_aux_fit, rho_ao_aux_fit)
157 
158  CALL timeset(routinen, handle)
159 
160  t1 = m_walltime()
161 
162  do_admm_rpa = mp2_env%ri_rpa%do_admm
163  do_ri_sigma_x = mp2_env%ri_g0w0%do_ri_Sigma_x
164  do_kpoints_cubic_rpa = qs_env%mp2_env%ri_rpa_im_time%do_im_time_kpoints
165  do_kpoints_from_gamma = qs_env%mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma
166  print_exx = mp2_env%ri_g0w0%print_exx
167 
168  IF (do_kpoints_cubic_rpa) THEN
169  cpassert(do_ri_sigma_x)
170  END IF
171 
172  IF (do_kpoints_cubic_rpa) THEN
173 
174  CALL get_qs_env(qs_env, &
175  admm_env=admm_env, &
176  matrix_ks_kp=matrix_ks_transl, &
177  rho=rho, &
178  input=input, &
179  dft_control=dft_control, &
180  para_env=para_env, &
181  kpoints=kpoints, &
182  ks_env=ks_env, &
183  energy=energy)
184  nkp = kpoints%nkp
185 
186  ELSE
187 
188  CALL get_qs_env(qs_env, &
189  admm_env=admm_env, &
190  matrix_ks=matrix_ks, &
191  rho=rho, &
192  input=input, &
193  dft_control=dft_control, &
194  para_env=para_env, &
195  ks_env=ks_env, &
196  energy=energy)
197  nkp = 1
198  CALL qs_rho_get(rho, rho_ao=rho_ao)
199 
200  IF (do_admm_rpa) THEN
201  CALL get_admm_env(admm_env, matrix_ks_aux_fit=matrix_ks_aux_fit, rho_aux_fit=rho_aux_fit, &
202  matrix_ks_aux_fit_hfx=matrix_ks_aux_fit_hfx)
203  CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux_fit)
204 
205  ! RPA/GW with ADMM for EXX or the exchange self-energy only implemented for
206  ! ADMM_PURIFICATION_METHOD NONE
207  ! METHOD BASIS_PROJECTION
208  ! in the admm section
209  cpassert(admm_env%purification_method == do_admm_purify_none)
210  cpassert(dft_control%admm_control%method == do_admm_basis_projection)
211  END IF
212  END IF
213 
214  nspins = dft_control%nspins
215 
216  ! safe ks matrix for later: we will transform matrix_ks
217  ! to T-cell index and then to k-points for band structure calculation
218  IF (do_kpoints_from_gamma) THEN
219  ! not yet there: open shell
220  ALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_ks(nspins))
221  DO ispin = 1, nspins
222  NULLIFY (qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix)
223  ALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix)
224  CALL dbcsr_create(qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix, &
225  template=matrix_ks(ispin)%matrix)
226  CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, &
227  qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix)
228 
229  END DO
230  END IF
231 
232  IF (do_kpoints_cubic_rpa) THEN
233 
234  CALL allocate_matrix_ks_kp(matrix_ks_transl, matrix_ks_kp_re, matrix_ks_kp_im, kpoints)
235  CALL transform_matrix_ks_to_kp(matrix_ks_transl, matrix_ks_kp_re, matrix_ks_kp_im, kpoints)
236 
237  DO ispin = 1, nspins
238  DO i_img = 1, SIZE(matrix_ks_transl, 2)
239  CALL dbcsr_set(matrix_ks_transl(ispin, i_img)%matrix, 0.0_dp)
240  END DO
241  END DO
242 
243  END IF
244 
245  ! initialize matrix_sigma_x_minus_vxc
246  NULLIFY (matrix_sigma_x_minus_vxc)
247  CALL dbcsr_allocate_matrix_set(matrix_sigma_x_minus_vxc, nspins, nkp)
248  IF (do_kpoints_cubic_rpa) THEN
249  NULLIFY (matrix_sigma_x_minus_vxc_im)
250  CALL dbcsr_allocate_matrix_set(matrix_sigma_x_minus_vxc_im, nspins, nkp)
251  END IF
252 
253  DO ispin = 1, nspins
254  DO ikp = 1, nkp
255 
256  IF (do_kpoints_cubic_rpa) THEN
257 
258  ALLOCATE (matrix_sigma_x_minus_vxc(ispin, ikp)%matrix)
259  CALL dbcsr_create(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, &
260  template=matrix_ks_kp_re(1, 1)%matrix, &
261  matrix_type=dbcsr_type_symmetric)
262 
263  CALL dbcsr_copy(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, matrix_ks_kp_re(ispin, ikp)%matrix)
264  CALL dbcsr_set(matrix_ks_kp_re(ispin, ikp)%matrix, 0.0_dp)
265 
266  ALLOCATE (matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix)
267  CALL dbcsr_create(matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix, &
268  template=matrix_ks_kp_im(1, 1)%matrix, &
269  matrix_type=dbcsr_type_antisymmetric)
270 
271  CALL dbcsr_copy(matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix, matrix_ks_kp_im(ispin, ikp)%matrix)
272  CALL dbcsr_set(matrix_ks_kp_im(ispin, ikp)%matrix, 0.0_dp)
273 
274  ELSE
275 
276  ALLOCATE (matrix_sigma_x_minus_vxc(ispin, ikp)%matrix)
277  CALL dbcsr_create(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, &
278  template=matrix_ks(1)%matrix)
279 
280  CALL dbcsr_copy(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, matrix_ks(ispin)%matrix)
281  CALL dbcsr_set(matrix_ks(ispin)%matrix, 0.0_dp)
282 
283  END IF
284 
285  END DO
286  END DO
287 
288  ! set DFT functional to none and hfx_fraction to zero
289  hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
290  CALL section_vals_get(hfx_sections, explicit=do_hfx)
291 
292  IF (do_hfx) THEN
293  hfx_fraction = qs_env%x_data(1, 1)%general_parameter%fraction
294  qs_env%x_data(:, :)%general_parameter%fraction = 0.0_dp
295  END IF
296  xc_section => section_vals_get_subs_vals(input, "DFT%XC")
297  CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
298  i_val=myfun)
299  CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
300  i_val=xc_none)
301 
302  ! in ADMM, also set the XC functional for ADMM correction to none
303  ! do not do this if we do ADMM for Sigma_x
304  IF (dft_control%do_admm) THEN
305  xc_section_admm_aux => section_vals_get_subs_vals(admm_env%xc_section_aux, &
306  "XC_FUNCTIONAL")
307  CALL section_vals_val_get(xc_section_admm_aux, "_SECTION_PARAMETERS_", &
308  i_val=myfun_aux)
309  CALL section_vals_val_set(xc_section_admm_aux, "_SECTION_PARAMETERS_", &
310  i_val=xc_none)
311 
312  ! the same for the primary basis
313  xc_section_admm_prim => section_vals_get_subs_vals(admm_env%xc_section_primary, &
314  "XC_FUNCTIONAL")
315  CALL section_vals_val_get(xc_section_admm_prim, "_SECTION_PARAMETERS_", &
316  i_val=myfun_prim)
317  CALL section_vals_val_set(xc_section_admm_prim, "_SECTION_PARAMETERS_", &
318  i_val=xc_none)
319 
320  ! for ADMMQ/S, set the charge_constrain to false (otherwise wrong results)
321  charge_constrain_tmp = .false.
322  IF (admm_env%charge_constrain) THEN
323  admm_env%charge_constrain = .false.
324  charge_constrain_tmp = .true.
325  END IF
326 
327  END IF
328 
329  ! if we do ADMM for Sigma_x, we write the ADMM correction into matrix_ks_aux_fit
330  ! and therefore we should set it to zero
331  IF (do_admm_rpa) THEN
332  DO ispin = 1, nspins
333  CALL dbcsr_set(matrix_ks_aux_fit(ispin)%matrix, 0.0_dp)
334  END DO
335  END IF
336 
337  IF (.NOT. mp2_env%ri_g0w0%update_xc_energy) THEN
338  energy_total = energy%total
339  energy_exc = energy%exc
340  energy_exc1 = energy%exc1
341  energy_exc_aux_fit = energy%ex
342  energy_exc1_aux_fit = energy%exc_aux_fit
343  energy_ex = energy%exc1_aux_fit
344  END IF
345 
346  ! Remove the Exchange-correlation energy contributions from the total energy
347  energy%total = energy%total - (energy%exc + energy%exc1 + energy%ex + &
348  energy%exc_aux_fit + energy%exc1_aux_fit)
349 
350  ! calculate KS-matrix without XC and without HF
351  CALL qs_ks_build_kohn_sham_matrix(qs_env=qs_env, calculate_forces=.false., &
352  just_energy=.false.)
353 
354  IF (.NOT. mp2_env%ri_g0w0%update_xc_energy) THEN
355  energy%exc = energy_exc
356  energy%exc1 = energy_exc1
357  energy%exc_aux_fit = energy_ex
358  energy%exc1_aux_fit = energy_exc_aux_fit
359  energy%ex = energy_exc1_aux_fit
360  energy%total = energy_total
361  END IF
362 
363  ! set the DFT functional and HF fraction back
364  CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
365  i_val=myfun)
366  IF (do_hfx) THEN
367  qs_env%x_data(:, :)%general_parameter%fraction = hfx_fraction
368  END IF
369 
370  IF (dft_control%do_admm) THEN
371  xc_section_admm_aux => section_vals_get_subs_vals(admm_env%xc_section_aux, &
372  "XC_FUNCTIONAL")
373  xc_section_admm_prim => section_vals_get_subs_vals(admm_env%xc_section_primary, &
374  "XC_FUNCTIONAL")
375 
376  CALL section_vals_val_set(xc_section_admm_aux, "_SECTION_PARAMETERS_", &
377  i_val=myfun_aux)
378  CALL section_vals_val_set(xc_section_admm_prim, "_SECTION_PARAMETERS_", &
379  i_val=myfun_prim)
380  IF (charge_constrain_tmp) THEN
381  admm_env%charge_constrain = .true.
382  END IF
383  END IF
384 
385  IF (do_kpoints_cubic_rpa) THEN
386  CALL transform_matrix_ks_to_kp(matrix_ks_transl, matrix_ks_kp_re, matrix_ks_kp_im, kpoints)
387  END IF
388 
389  ! remove the single-particle part (kin. En + Hartree pot) and change the sign
390  DO ispin = 1, nspins
391  IF (do_kpoints_cubic_rpa) THEN
392  DO ikp = 1, nkp
393  CALL dbcsr_add(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, matrix_ks_kp_re(ispin, ikp)%matrix, -1.0_dp, 1.0_dp)
394  CALL dbcsr_add(matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix, matrix_ks_kp_im(ispin, ikp)%matrix, -1.0_dp, 1.0_dp)
395  END DO
396  ELSE
397  CALL dbcsr_add(matrix_sigma_x_minus_vxc(ispin, 1)%matrix, matrix_ks(ispin)%matrix, -1.0_dp, 1.0_dp)
398  END IF
399  END DO
400 
401  IF (do_kpoints_cubic_rpa) THEN
402 
403  CALL transform_sigma_x_minus_vxc_to_mo_basis(kpoints, matrix_sigma_x_minus_vxc, &
404  matrix_sigma_x_minus_vxc_im, &
405  vec_sigma_x_minus_vxc_gw, &
406  vec_sigma_x_minus_vxc_gw_im, &
407  para_env, nmo, mp2_env)
408 
409  ELSE
410 
411  DO ispin = 1, nspins
412  CALL dbcsr_set(matrix_ks(ispin)%matrix, 0.0_dp)
413  IF (do_admm_rpa) THEN
414  CALL dbcsr_set(matrix_ks_aux_fit(ispin)%matrix, 0.0_dp)
415  END IF
416  END DO
417 
418  hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
419 
420  CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
421 
422  ! in most cases, we calculate the exchange self-energy here. But if we do only RI for
423  ! the exchange self-energy, we do not calculate exchange here
424  ehfx = 0.0_dp
425  IF (.NOT. do_ri_sigma_x) THEN
426 
427  CALL exx_pre_hfx(hfx_sections, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx)
428  calc_ints = .NOT. qs_env%mp2_env%ri_rpa%reuse_hfx
429 
430  ! add here HFX (=Sigma_exchange) to matrix_sigma_x_minus_vxc
431  DO irep = 1, n_rep_hf
432  IF (do_admm_rpa) THEN
433  matrix_ks_2d(1:nspins, 1:1) => matrix_ks_aux_fit(1:nspins)
434  rho_ao_2d(1:nspins, 1:1) => rho_ao_aux_fit(1:nspins)
435  ELSE
436  matrix_ks_2d(1:nspins, 1:1) => matrix_ks(1:nspins)
437  rho_ao_2d(1:nspins, 1:1) => rho_ao(1:nspins)
438  END IF
439 
440  IF (qs_env%mp2_env%ri_rpa%x_data(irep, 1)%do_hfx_ri) THEN
441  CALL hfx_ri_update_ks(qs_env, qs_env%mp2_env%ri_rpa%x_data(irep, 1)%ri_data, matrix_ks_2d, ehfx, &
442  rho_ao=rho_ao_2d, geometry_did_change=calc_ints, nspins=nspins, &
443  hf_fraction=qs_env%mp2_env%ri_rpa%x_data(irep, 1)%general_parameter%fraction)
444 
445  IF (do_admm_rpa) THEN
446  !for ADMMS, we need the exchange matrix k(d) for both spins
447  DO ispin = 1, nspins
448  CALL dbcsr_copy(matrix_ks_aux_fit_hfx(ispin)%matrix, matrix_ks_2d(ispin, 1)%matrix, &
449  name="HF exch. part of matrix_ks_aux_fit for ADMMS")
450  END DO
451  END IF
452  ELSE
453  CALL integrate_four_center(qs_env, qs_env%mp2_env%ri_rpa%x_data, matrix_ks_2d, eh1, &
454  rho_ao_2d, hfx_sections, &
455  para_env, calc_ints, irep, .true., &
456  ispin=1)
457  ehfx = ehfx + eh1
458  END IF
459  END DO
460 
461  !ADMM XC correction
462  IF (do_admm_rpa) THEN
463  CALL calc_exx_admm_xc_contributions(qs_env=qs_env, &
464  matrix_prim=matrix_ks, &
465  matrix_aux=matrix_ks_aux_fit, &
466  x_data=qs_env%mp2_env%ri_rpa%x_data, &
467  exc=energy_xc_admm(1), &
468  exc_aux_fit=energy_xc_admm(2), &
469  calc_forces=.false., &
470  use_virial=.false.)
471  END IF
472 
473  IF (do_kpoints_from_gamma .AND. print_exx == gw_print_exx) THEN
474  ! JW not yet there: open shell
475  ALLOCATE (mat_exchange_for_kp_from_gamma(1))
476 
477  DO ispin = 1, 1
478  NULLIFY (mat_exchange_for_kp_from_gamma(ispin)%matrix)
479  ALLOCATE (mat_exchange_for_kp_from_gamma(ispin)%matrix)
480  CALL dbcsr_create(mat_exchange_for_kp_from_gamma(ispin)%matrix, template=matrix_ks(ispin)%matrix)
481  CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, mat_exchange_for_kp_from_gamma(ispin)%matrix)
482  END DO
483 
484  END IF
485 
486  CALL exx_post_hfx(qs_env, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx)
487  END IF
488 
489  energy_ex = ehfx
490 
491  ! transform Fock-Matrix (calculated in integrate_four_center, written in matrix_ks_aux_fit in case
492  ! of ADMM) from ADMM basis to primary basis
493  IF (do_admm_rpa) THEN
494  CALL admm_mo_merge_ks_matrix(qs_env)
495  END IF
496 
497  DO ispin = 1, nspins
498  CALL dbcsr_add(matrix_sigma_x_minus_vxc(ispin, 1)%matrix, matrix_ks(ispin)%matrix, 1.0_dp, 1.0_dp)
499  END DO
500 
501  ! safe matrix_sigma_x_minus_vxc for later: for example, we will transform matrix_sigma_x_minus_vxc
502  ! to T-cell index and then to k-points for band structure calculation
503  IF (do_kpoints_from_gamma) THEN
504  ! not yet there: open shell
505  ALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(nspins))
506  DO ispin = 1, nspins
507  NULLIFY (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix)
508  ALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix)
509  CALL dbcsr_create(qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix, &
510  template=matrix_ks(ispin)%matrix)
511 
512  CALL dbcsr_desymmetrize(matrix_sigma_x_minus_vxc(ispin, 1)%matrix, &
513  qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix)
514 
515  END DO
516  END IF
517 
518  CALL dbcsr_desymmetrize(matrix_ks(1)%matrix, mo_coeff_b)
519  CALL dbcsr_set(mo_coeff_b, 0.0_dp)
520 
521  ! Transform matrix_sigma_x_minus_vxc to MO basis
522  DO ispin = 1, nspins
523 
524  CALL get_mo_set(mo_set=mos_mp2(ispin), &
525  mo_coeff=mo_coeff, &
526  nmo=nmo, &
527  homo=homo, &
528  nao=dimen)
529 
530  IF (ispin == 1) THEN
531 
532  ALLOCATE (vec_sigma_x_minus_vxc_gw(nmo, nspins, nkp))
533  vec_sigma_x_minus_vxc_gw = 0.0_dp
534  END IF
535 
536  CALL dbcsr_set(mo_coeff_b, 0.0_dp)
537  CALL copy_fm_to_dbcsr(mo_coeff, mo_coeff_b, keep_sparsity=.false.)
538 
539  ! initialize matrix_tmp and matrix_tmp2
540  IF (ispin == 1) THEN
541  CALL dbcsr_create(matrix_tmp, template=mo_coeff_b)
542  CALL dbcsr_copy(matrix_tmp, mo_coeff_b)
543  CALL dbcsr_set(matrix_tmp, 0.0_dp)
544 
545  CALL dbcsr_create(matrix_tmp_2, template=mo_coeff_b)
546  CALL dbcsr_copy(matrix_tmp_2, mo_coeff_b)
547  CALL dbcsr_set(matrix_tmp_2, 0.0_dp)
548  END IF
549 
550  gw_corr_lev_occ = mp2_env%ri_g0w0%corr_mos_occ
551  gw_corr_lev_virt = mp2_env%ri_g0w0%corr_mos_virt
552  ! if requested number of occ/virt levels for correction either exceed the number of
553  ! occ/virt levels or the requested number is negative, default to correct all
554  ! occ/virt level energies
555  IF (gw_corr_lev_occ > homo .OR. gw_corr_lev_occ < 0) gw_corr_lev_occ = homo
556  IF (gw_corr_lev_virt > dimen - homo .OR. gw_corr_lev_virt < 0) gw_corr_lev_virt = dimen - homo
557  IF (ispin == 1) THEN
558  mp2_env%ri_g0w0%corr_mos_occ = gw_corr_lev_occ
559  mp2_env%ri_g0w0%corr_mos_virt = gw_corr_lev_virt
560  ELSE IF (ispin == 2) THEN
561  ! ensure that the total number of corrected MOs is the same for alpha and beta, important
562  ! for parallelization
563  IF (mp2_env%ri_g0w0%corr_mos_occ + mp2_env%ri_g0w0%corr_mos_virt /= &
564  gw_corr_lev_occ + gw_corr_lev_virt) THEN
565  gw_corr_lev_virt = mp2_env%ri_g0w0%corr_mos_occ + mp2_env%ri_g0w0%corr_mos_virt - gw_corr_lev_occ
566  END IF
567  mp2_env%ri_g0w0%corr_mos_occ_beta = gw_corr_lev_occ
568  mp2_env%ri_g0w0%corr_mos_virt_beta = gw_corr_lev_virt
569 
570  END IF
571 
572  CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_sigma_x_minus_vxc(ispin, 1)%matrix, &
573  mo_coeff_b, 0.0_dp, matrix_tmp, first_column=homo + 1 - gw_corr_lev_occ, &
574  last_column=homo + gw_corr_lev_virt)
575 
576  CALL dbcsr_multiply('T', 'N', 1.0_dp, mo_coeff_b, &
577  matrix_tmp, 0.0_dp, matrix_tmp_2, first_row=homo + 1 - gw_corr_lev_occ, &
578  last_row=homo + gw_corr_lev_virt)
579 
580  CALL dbcsr_get_diag(matrix_tmp_2, vec_sigma_x_minus_vxc_gw(:, ispin, 1))
581 
582  CALL dbcsr_set(matrix_tmp, 0.0_dp)
583  CALL dbcsr_set(matrix_tmp_2, 0.0_dp)
584 
585  END DO
586 
587  CALL para_env%sum(vec_sigma_x_minus_vxc_gw)
588 
589  END IF
590 
591  CALL dbcsr_release(mo_coeff_b)
592  CALL dbcsr_release(matrix_tmp)
593  CALL dbcsr_release(matrix_tmp_2)
594  IF (do_kpoints_cubic_rpa) THEN
595  CALL dbcsr_deallocate_matrix_set(matrix_ks_kp_re)
596  CALL dbcsr_deallocate_matrix_set(matrix_ks_kp_im)
597  END IF
598 
599  DO ispin = 1, nspins
600  DO ikp = 1, nkp
601  CALL dbcsr_release_p(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix)
602  IF (do_kpoints_cubic_rpa) THEN
603  CALL dbcsr_release_p(matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix)
604  END IF
605  END DO
606  END DO
607 
608  ALLOCATE (mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(nmo, nspins, nkp))
609 
610  IF (print_exx == gw_print_exx) THEN
611 
612  IF (do_kpoints_from_gamma) THEN
613 
614  gw_corr_lev_tot = gw_corr_lev_occ + gw_corr_lev_virt
615 
616  CALL get_qs_env(qs_env=qs_env, &
617  kpoints=kpoints)
618 
619  CALL setup_trunc_coulomb_pot_for_exchange_self_energy(qs_env)
620 
621  CALL compute_kpoints(qs_env, kpoints, unit_nr)
622 
623  ALLOCATE (eigenval_kp(nmo, 1, nspins))
624 
625  CALL get_bandstruc_and_k_dependent_mos(qs_env, eigenval_kp)
626 
627  CALL compute_minus_vxc_kpoints(qs_env)
628 
629  nkp_sigma = SIZE(eigenval_kp, 2)
630 
631  ALLOCATE (vec_sigma_x(nmo, nkp_sigma))
632  vec_sigma_x(:, :) = 0.0_dp
633 
634  CALL trafo_to_mo_and_kpoints(qs_env, &
635  mat_exchange_for_kp_from_gamma(1)%matrix, &
636  vec_sigma_x(homo - gw_corr_lev_occ + 1:homo + gw_corr_lev_virt, :), &
637  homo, gw_corr_lev_occ, gw_corr_lev_virt, 1)
638 
639  CALL dbcsr_release(mat_exchange_for_kp_from_gamma(1)%matrix)
640  DEALLOCATE (mat_exchange_for_kp_from_gamma(1)%matrix)
641  DEALLOCATE (mat_exchange_for_kp_from_gamma)
642 
643  DEALLOCATE (vec_sigma_x_minus_vxc_gw)
644 
645  ALLOCATE (vec_sigma_x_minus_vxc_gw(nmo, nspins, nkp_sigma))
646 
647  vec_sigma_x_minus_vxc_gw(:, 1, :) = vec_sigma_x(:, :) + &
648  qs_env%mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, :)
649 
650  kpoints_sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma
651 
652  ELSE
653 
654  nkp_sigma = 1
655 
656  END IF
657 
658  IF (unit_nr > 0) THEN
659 
660  ALLOCATE (eigenval_kp_hf_at_dft(nmo, nkp_sigma))
661  eigenval_kp_hf_at_dft(:, :) = eigenval_kp(:, :, 1) + vec_sigma_x_minus_vxc_gw(:, 1, :)
662 
663  min_direct_hf_at_dft_gap = 100.0_dp
664 
665  WRITE (unit_nr, '(T3,A)') ''
666  WRITE (unit_nr, '(T3,A)') 'Exchange energies'
667  WRITE (unit_nr, '(T3,A)') '-----------------'
668  WRITE (unit_nr, '(T3,A)') ''
669  WRITE (unit_nr, '(T6,2A)') 'MO e_n^DFT Sigma_x-vxc e_n^HF@DFT'
670  DO ikp = 1, nkp_sigma
671  IF (nkp_sigma > 1) THEN
672  WRITE (unit_nr, '(T3,A)') ''
673  WRITE (unit_nr, '(T3,A7,I3,A3,I3,A8,3F7.3,A12,3F7.3)') 'Kpoint ', ikp, ' /', nkp_sigma, &
674  ' xkp =', kpoints_sigma%xkp(1, ikp), kpoints_sigma%xkp(2, ikp), &
675  kpoints_sigma%xkp(3, ikp), ' and xkp =', -kpoints_sigma%xkp(1, ikp), &
676  -kpoints_sigma%xkp(2, ikp), -kpoints_sigma%xkp(3, ikp)
677  WRITE (unit_nr, '(T3,A)') ''
678  END IF
679  DO n_level_gw = 1, gw_corr_lev_occ + gw_corr_lev_virt
680 
681  n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
682  IF (n_level_gw <= gw_corr_lev_occ) THEN
683  occ_virt = 'occ'
684  ELSE
685  occ_virt = 'vir'
686  END IF
687 
688  eigval_dft = eigenval_kp(n_level_gw_ref, ikp, 1)*evolt
689  exx_minus_vxc = real(vec_sigma_x_minus_vxc_gw(n_level_gw_ref, 1, ikp)*evolt, kind=dp)
690  eigval_hf_at_dft = eigenval_kp_hf_at_dft(n_level_gw_ref, ikp)*evolt
691 
692  WRITE (unit_nr, '(T4,I4,3A,3F21.3,3F21.3,3F21.3)') &
693  n_level_gw_ref, ' ( ', occ_virt, ') ', eigval_dft, exx_minus_vxc, eigval_hf_at_dft
694 
695  END DO
696  e_homo_gw = maxval(eigenval_kp_hf_at_dft(homo - gw_corr_lev_occ + 1:homo, ikp))
697  e_lumo_gw = minval(eigenval_kp_hf_at_dft(homo + 1:homo + gw_corr_lev_virt, ikp))
698  e_gap_gw = e_lumo_gw - e_homo_gw
699  IF (e_gap_gw < min_direct_hf_at_dft_gap) min_direct_hf_at_dft_gap = e_gap_gw
700  WRITE (unit_nr, '(T3,A)') ''
701  WRITE (unit_nr, '(T3,A,F53.2)') 'HF@DFT HOMO-LUMO gap (eV)', e_gap_gw*evolt
702  WRITE (unit_nr, '(T3,A)') ''
703  END DO
704 
705  WRITE (unit_nr, '(T3,A)') ''
706  WRITE (unit_nr, '(T3,A)') ''
707  WRITE (unit_nr, '(T3,A,F63.3)') 'HF@DFT direct bandgap (eV)', min_direct_hf_at_dft_gap*evolt
708 
709  WRITE (unit_nr, '(T3,A)') ''
710  WRITE (unit_nr, '(T3,A)') 'End of exchange energies'
711  WRITE (unit_nr, '(T3,A)') '------------------------'
712  WRITE (unit_nr, '(T3,A)') ''
713 
714  cpabort('Stop after printing exchange energies.')
715 
716  ELSE
717  CALL para_env%sync()
718  END IF
719 
720  END IF
721 
722  IF (print_exx == gw_read_exx) THEN
723 
724  CALL open_file(unit_number=iunit, file_name="exx.out")
725 
726  really_read_line = .false.
727 
728  DO WHILE (.true.)
729 
730  READ (iunit, '(A)') line
731 
732  IF (line == " End of exchange energies ") EXIT
733 
734  IF (really_read_line) THEN
735 
736  READ (line(1:7), *) n_level_gw_ref
737  READ (line(17:40), *) tmp
738 
739  DO ikp = 1, SIZE(vec_sigma_x_minus_vxc_gw, 3)
740  vec_sigma_x_minus_vxc_gw(n_level_gw_ref, 1, ikp) = tmp/evolt
741  END DO
742 
743  END IF
744 
745  IF (line == " MO Sigma_x-vxc ") really_read_line = .true.
746 
747  END DO
748 
749  CALL close_file(iunit)
750 
751  END IF
752 
753  ! store vec_Sigma_x_minus_vxc_gw in the mp2_environment
754  mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, :, :) = vec_sigma_x_minus_vxc_gw(:, :, :)
755 
756  ! clean up
757  DEALLOCATE (matrix_sigma_x_minus_vxc, vec_sigma_x_minus_vxc_gw)
758  IF (do_kpoints_cubic_rpa) THEN
759  DEALLOCATE (matrix_sigma_x_minus_vxc_im)
760  END IF
761 
762  t2 = m_walltime()
763 
764  t3 = t2 - t1
765 
766  CALL timestop(handle)
767 
768  END SUBROUTINE compute_vec_sigma_x_minus_vxc_gw
769 
770 ! **************************************************************************************************
771 !> \brief ...
772 !> \param kpoints ...
773 !> \param matrix_sigma_x_minus_vxc ...
774 !> \param matrix_sigma_x_minus_vxc_im ...
775 !> \param vec_Sigma_x_minus_vxc_gw ...
776 !> \param vec_Sigma_x_minus_vxc_gw_im ...
777 !> \param para_env ...
778 !> \param nmo ...
779 !> \param mp2_env ...
780 ! **************************************************************************************************
781  SUBROUTINE transform_sigma_x_minus_vxc_to_mo_basis(kpoints, matrix_sigma_x_minus_vxc, &
782  matrix_sigma_x_minus_vxc_im, vec_Sigma_x_minus_vxc_gw, &
783  vec_Sigma_x_minus_vxc_gw_im, para_env, nmo, mp2_env)
784 
785  TYPE(kpoint_type), POINTER :: kpoints
786  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_sigma_x_minus_vxc, &
787  matrix_sigma_x_minus_vxc_im
788  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: vec_sigma_x_minus_vxc_gw, &
789  vec_sigma_x_minus_vxc_gw_im
790  TYPE(mp_para_env_type), INTENT(IN) :: para_env
791  INTEGER :: nmo
792  TYPE(mp2_type) :: mp2_env
793 
794  CHARACTER(LEN=*), PARAMETER :: routinen = 'transform_sigma_x_minus_vxc_to_MO_basis'
795 
796  INTEGER :: dimen, gw_corr_lev_occ, gw_corr_lev_virt, handle, homo, i_global, iib, ikp, &
797  ispin, j_global, jjb, ncol_local, nkp, nrow_local, nspins
798  INTEGER, DIMENSION(2) :: kp_range
799  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
800  REAL(kind=dp) :: imval, reval
801  TYPE(cp_cfm_type) :: cfm_mos, cfm_sigma_x_minus_vxc, &
802  cfm_sigma_x_minus_vxc_mo_basis, cfm_tmp
803  TYPE(cp_fm_struct_type), POINTER :: matrix_struct
804  TYPE(cp_fm_type) :: fwork_im, fwork_re
805  TYPE(kpoint_env_type), POINTER :: kp
806  TYPE(mo_set_type), POINTER :: mo_set, mo_set_im, mo_set_re
807 
808  CALL timeset(routinen, handle)
809 
810  mo_set => kpoints%kp_env(1)%kpoint_env%mos(1, 1)
811  CALL get_mo_set(mo_set, nmo=nmo)
812 
813  nspins = SIZE(matrix_sigma_x_minus_vxc, 1)
814  CALL get_kpoint_info(kpoints, nkp=nkp, kp_range=kp_range)
815 
816  ! if this CPASSERT is wrong, please make sure that the kpoint group size PARALLEL_GROUP_SIZE
817  ! in the kpoint environment &DFT &KPOINTS is -1
818  cpassert(kp_range(1) == 1 .AND. kp_range(2) == nkp)
819 
820  ALLOCATE (vec_sigma_x_minus_vxc_gw(nmo, nspins, nkp))
821  vec_sigma_x_minus_vxc_gw = 0.0_dp
822 
823  ALLOCATE (vec_sigma_x_minus_vxc_gw_im(nmo, nspins, nkp))
824  vec_sigma_x_minus_vxc_gw_im = 0.0_dp
825 
826  CALL cp_fm_get_info(mo_set%mo_coeff, matrix_struct=matrix_struct)
827  CALL cp_fm_create(fwork_re, matrix_struct)
828  CALL cp_fm_create(fwork_im, matrix_struct)
829  CALL cp_cfm_create(cfm_mos, matrix_struct)
830  CALL cp_cfm_create(cfm_sigma_x_minus_vxc, matrix_struct)
831  CALL cp_cfm_create(cfm_sigma_x_minus_vxc_mo_basis, matrix_struct)
832  CALL cp_cfm_create(cfm_tmp, matrix_struct)
833 
834  CALL cp_cfm_get_info(matrix=cfm_sigma_x_minus_vxc_mo_basis, &
835  nrow_local=nrow_local, &
836  ncol_local=ncol_local, &
837  row_indices=row_indices, &
838  col_indices=col_indices)
839 
840  ! Transform matrix_sigma_x_minus_vxc to MO basis
841  DO ikp = 1, nkp
842 
843  kp => kpoints%kp_env(ikp)%kpoint_env
844 
845  DO ispin = 1, nspins
846 
847  ! v_xc_n to fm matrix
848  CALL copy_dbcsr_to_fm(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, fwork_re)
849  CALL copy_dbcsr_to_fm(matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix, fwork_im)
850 
851  CALL cp_cfm_scale_and_add_fm(z_zero, cfm_sigma_x_minus_vxc, z_one, fwork_re)
852  CALL cp_cfm_scale_and_add_fm(z_one, cfm_sigma_x_minus_vxc, gaussi, fwork_im)
853 
854  ! get real part (1) and imag. part (2) of the mo coeffs
855  mo_set_re => kp%mos(1, ispin)
856  mo_set_im => kp%mos(2, ispin)
857 
858  CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mos, z_one, mo_set_re%mo_coeff)
859  CALL cp_cfm_scale_and_add_fm(z_one, cfm_mos, gaussi, mo_set_im%mo_coeff)
860 
861  ! tmp = V(k)*C(k)
862  CALL parallel_gemm('N', 'N', nmo, nmo, nmo, z_one, cfm_sigma_x_minus_vxc, &
863  cfm_mos, z_zero, cfm_tmp)
864 
865  ! V_n(k) = C^H(k)*tmp
866  CALL parallel_gemm('C', 'N', nmo, nmo, nmo, z_one, cfm_mos, cfm_tmp, &
867  z_zero, cfm_sigma_x_minus_vxc_mo_basis)
868 
869  DO jjb = 1, ncol_local
870 
871  j_global = col_indices(jjb)
872 
873  DO iib = 1, nrow_local
874 
875  i_global = row_indices(iib)
876 
877  IF (j_global == i_global .AND. i_global <= nmo) THEN
878 
879  reval = real(cfm_sigma_x_minus_vxc_mo_basis%local_data(iib, jjb), kind=dp)
880  imval = aimag(cfm_sigma_x_minus_vxc_mo_basis%local_data(iib, jjb))
881 
882  vec_sigma_x_minus_vxc_gw(i_global, ispin, ikp) = reval
883  vec_sigma_x_minus_vxc_gw_im(i_global, ispin, ikp) = imval
884 
885  END IF
886 
887  END DO
888 
889  END DO
890 
891  END DO
892 
893  END DO
894 
895  CALL para_env%sum(vec_sigma_x_minus_vxc_gw)
896  CALL para_env%sum(vec_sigma_x_minus_vxc_gw_im)
897 
898  ! also adjust in the case of kpoints too big gw_corr_lev_occ and gw_corr_lev_virt
899  DO ispin = 1, nspins
900  CALL get_mo_set(mo_set=kpoints%kp_env(1)%kpoint_env%mos(ispin, 1), &
901  homo=homo, nao=dimen)
902  gw_corr_lev_occ = mp2_env%ri_g0w0%corr_mos_occ
903  gw_corr_lev_virt = mp2_env%ri_g0w0%corr_mos_virt
904  ! if corrected occ/virt levels exceed the number of occ/virt levels,
905  ! correct all occ/virt level energies
906  IF (gw_corr_lev_occ > homo) gw_corr_lev_occ = homo
907  IF (gw_corr_lev_virt > dimen - homo) gw_corr_lev_virt = dimen - homo
908  IF (ispin == 1) THEN
909  mp2_env%ri_g0w0%corr_mos_occ = gw_corr_lev_occ
910  mp2_env%ri_g0w0%corr_mos_virt = gw_corr_lev_virt
911  ELSE IF (ispin == 2) THEN
912  ! ensure that the total number of corrected MOs is the same for alpha and beta, important
913  ! for parallelization
914  IF (mp2_env%ri_g0w0%corr_mos_occ + mp2_env%ri_g0w0%corr_mos_virt /= &
915  gw_corr_lev_occ + gw_corr_lev_virt) THEN
916  gw_corr_lev_virt = mp2_env%ri_g0w0%corr_mos_occ + mp2_env%ri_g0w0%corr_mos_virt - gw_corr_lev_occ
917  END IF
918  mp2_env%ri_g0w0%corr_mos_occ_beta = gw_corr_lev_occ
919  mp2_env%ri_g0w0%corr_mos_virt_beta = gw_corr_lev_virt
920  END IF
921  END DO
922 
923  CALL cp_fm_release(fwork_re)
924  CALL cp_fm_release(fwork_im)
925  CALL cp_cfm_release(cfm_mos)
926  CALL cp_cfm_release(cfm_sigma_x_minus_vxc)
927  CALL cp_cfm_release(cfm_sigma_x_minus_vxc_mo_basis)
928  CALL cp_cfm_release(cfm_tmp)
929 
930  CALL timestop(handle)
931 
932  END SUBROUTINE
933 
934 ! **************************************************************************************************
935 !> \brief ...
936 !> \param matrix_ks_transl ...
937 !> \param matrix_ks_kp_re ...
938 !> \param matrix_ks_kp_im ...
939 !> \param kpoints ...
940 ! **************************************************************************************************
941  SUBROUTINE transform_matrix_ks_to_kp(matrix_ks_transl, matrix_ks_kp_re, matrix_ks_kp_im, kpoints)
942 
943  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_transl, matrix_ks_kp_re, &
944  matrix_ks_kp_im
945  TYPE(kpoint_type), POINTER :: kpoints
946 
947  CHARACTER(len=*), PARAMETER :: routinen = 'transform_matrix_ks_to_kp'
948 
949  INTEGER :: handle, ikp, ispin, nkp, nspin
950  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
951  REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
952  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
953  POINTER :: sab_nl
954 
955  CALL timeset(routinen, handle)
956 
957  NULLIFY (sab_nl)
958  CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, sab_nl=sab_nl, cell_to_index=cell_to_index)
959 
960  cpassert(ASSOCIATED(sab_nl))
961 
962  nspin = SIZE(matrix_ks_transl, 1)
963 
964  DO ikp = 1, nkp
965  DO ispin = 1, nspin
966 
967  CALL dbcsr_set(matrix_ks_kp_re(ispin, ikp)%matrix, 0.0_dp)
968  CALL dbcsr_set(matrix_ks_kp_im(ispin, ikp)%matrix, 0.0_dp)
969  CALL rskp_transform(rmatrix=matrix_ks_kp_re(ispin, ikp)%matrix, &
970  cmatrix=matrix_ks_kp_im(ispin, ikp)%matrix, &
971  rsmat=matrix_ks_transl, ispin=ispin, &
972  xkp=xkp(1:3, ikp), cell_to_index=cell_to_index, sab_nl=sab_nl)
973 
974  END DO
975  END DO
976 
977  CALL timestop(handle)
978 
979  END SUBROUTINE
980 
981 ! **************************************************************************************************
982 !> \brief ...
983 !> \param matrix_ks_transl ...
984 !> \param matrix_ks_kp_re ...
985 !> \param matrix_ks_kp_im ...
986 !> \param kpoints ...
987 ! **************************************************************************************************
988  SUBROUTINE allocate_matrix_ks_kp(matrix_ks_transl, matrix_ks_kp_re, matrix_ks_kp_im, kpoints)
989 
990  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_transl, matrix_ks_kp_re, &
991  matrix_ks_kp_im
992  TYPE(kpoint_type), POINTER :: kpoints
993 
994  CHARACTER(len=*), PARAMETER :: routinen = 'allocate_matrix_ks_kp'
995 
996  INTEGER :: handle, ikp, ispin, nkp, nspin
997  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
998  REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
999  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1000  POINTER :: sab_nl
1001 
1002  CALL timeset(routinen, handle)
1003 
1004  NULLIFY (sab_nl)
1005  CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, sab_nl=sab_nl, cell_to_index=cell_to_index)
1006 
1007  cpassert(ASSOCIATED(sab_nl))
1008 
1009  nspin = SIZE(matrix_ks_transl, 1)
1010 
1011  NULLIFY (matrix_ks_kp_re, matrix_ks_kp_im)
1012  CALL dbcsr_allocate_matrix_set(matrix_ks_kp_re, nspin, nkp)
1013  CALL dbcsr_allocate_matrix_set(matrix_ks_kp_im, nspin, nkp)
1014 
1015  DO ikp = 1, nkp
1016  DO ispin = 1, nspin
1017 
1018  ALLOCATE (matrix_ks_kp_re(ispin, ikp)%matrix)
1019  ALLOCATE (matrix_ks_kp_im(ispin, ikp)%matrix)
1020 
1021  CALL dbcsr_create(matrix_ks_kp_re(ispin, ikp)%matrix, &
1022  template=matrix_ks_transl(1, 1)%matrix, &
1023  matrix_type=dbcsr_type_symmetric)
1024  CALL dbcsr_create(matrix_ks_kp_im(ispin, ikp)%matrix, &
1025  template=matrix_ks_transl(1, 1)%matrix, &
1026  matrix_type=dbcsr_type_antisymmetric)
1027 
1028  CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_kp_re(ispin, ikp)%matrix, sab_nl)
1029  CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_kp_im(ispin, ikp)%matrix, sab_nl)
1030 
1031  CALL dbcsr_set(matrix_ks_kp_re(ispin, ikp)%matrix, 0.0_dp)
1032  CALL dbcsr_set(matrix_ks_kp_im(ispin, ikp)%matrix, 0.0_dp)
1033 
1034  END DO
1035  END DO
1036 
1037  CALL timestop(handle)
1038 
1039  END SUBROUTINE
1040 
1041 END MODULE rpa_gw_sigma_x
1042 
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
admm_methods
Contains ADMM methods which require molecular orbitals.
Definition: admm_methods.F:15
admm_methods::admm_mo_merge_ks_matrix
subroutine, public admm_mo_merge_ks_matrix(qs_env)
...
Definition: admm_methods.F:778
admm_types
Types and set/get functions for auxiliary density matrix methods.
Definition: admm_types.F:15
admm_types::get_admm_env
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition: admm_types.F:593
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:355
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_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_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition: cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition: cp_dbcsr_cp2k_link.F:445
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_files
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition: cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition: cp_files.F:119
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_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
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_energy_potential
Routines to calculate HFX energy and potential.
Definition: hfx_energy_potential.F:14
hfx_energy_potential::integrate_four_center
subroutine, public integrate_four_center(qs_env, x_data, ks_matrix, ehfx, rho_ao, hfx_section, para_env, geometry_did_change, irep, distribute_fock_matrix, ispin)
computes four center integrals for a full basis set and updates the Kohn-Sham-Matrix and energy....
Definition: hfx_energy_potential.F:137
hfx_exx
Routines to calculate EXX in RPA and energy correction methods.
Definition: hfx_exx.F:16
hfx_exx::calc_exx_admm_xc_contributions
subroutine, public calc_exx_admm_xc_contributions(qs_env, matrix_prim, matrix_aux, x_data, exc, exc_aux_fit, calc_forces, use_virial)
Calculate the RI_RPAHF / EC_ENVHF ADMM XC contributions to the KS matrices and the respective energie...
Definition: hfx_exx.F:629
hfx_exx::exx_pre_hfx
subroutine, public exx_pre_hfx(ext_hfx_section, x_data, reuse_hfx)
Prepare the external x_data for integration. Simply change the HFX fraction in case the qs_envx_data ...
Definition: hfx_exx.F:735
hfx_exx::exx_post_hfx
subroutine, public exx_post_hfx(qs_env, x_data, reuse_hfx)
Revert back to the proper HFX fraction in case qs_envx_data is reused.
Definition: hfx_exx.F:760
hfx_ri
RI-methods for HFX.
Definition: hfx_ri.F:12
hfx_ri::hfx_ri_update_ks
subroutine, public hfx_ri_update_ks(qs_env, ri_data, ks_matrix, ehfx, mos, rho_ao, geometry_did_change, nspins, hf_fraction)
...
Definition: hfx_ri.F:1033
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_admm_purify_none
integer, parameter, public do_admm_purify_none
Definition: input_constants.F:848
input_constants::gw_print_exx
integer, parameter, public gw_print_exx
Definition: input_constants.F:1063
input_constants::do_admm_basis_projection
integer, parameter, public do_admm_basis_projection
Definition: input_constants.F:857
input_constants::gw_read_exx
integer, parameter, public gw_read_exx
Definition: input_constants.F:1063
input_constants::xc_none
integer, parameter, public xc_none
Definition: input_constants.F:551
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_val_set
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
Definition: input_section_types.F:1216
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition: input_section_types.F:697
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
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::rskp_transform
subroutine, public rskp_transform(rmatrix, cmatrix, rsmat, ispin, xkp, cell_to_index, sab_nl, is_complex, rs_sign)
Transformation of real space matrices to a kpoint.
Definition: kpoint_methods.F:780
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::z_zero
complex(kind=dp), parameter, public z_zero
Definition: mathconstants.F:143
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
mp2_integrals
Routines to calculate and distribute 2c- and 3c- integrals for RI.
Definition: mp2_integrals.F:14
mp2_integrals::compute_kpoints
subroutine, public compute_kpoints(qs_env, kpoints, unit_nr)
...
Definition: mp2_integrals.F:1213
mp2_ri_2c
Framework for 2c-integrals for RI.
Definition: mp2_ri_2c.F:14
mp2_ri_2c::setup_trunc_coulomb_pot_for_exchange_self_energy
subroutine, public setup_trunc_coulomb_pot_for_exchange_self_energy(qs_env, trunc_coulomb, rel_cutoff_trunc_coulomb_ri_x)
...
Definition: mp2_ri_2c.F:1606
mp2_types
Types needed for MP2 calculations.
Definition: mp2_types.F:14
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition: physcon.F:183
qs_energy_types
Definition: qs_energy_types.F:14
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_ks_methods
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition: qs_ks_methods.F:22
qs_ks_methods::qs_ks_build_kohn_sham_matrix
subroutine, public qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, print_active, ext_ks_matrix)
routine where the real calculations are made: the KS matrix is calculated
Definition: qs_ks_methods.F:173
qs_ks_types
Definition: qs_ks_types.F:15
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_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
rpa_gw_kpoints_util
Routines treating GW and RPA calculations with kpoints.
Definition: rpa_gw_kpoints_util.F:13
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_sigma_x
Routines to calculate EXX within GW.
Definition: rpa_gw_sigma_x.F:14
rpa_gw_sigma_x::compute_vec_sigma_x_minus_vxc_gw
subroutine, public compute_vec_sigma_x_minus_vxc_gw(qs_env, mp2_env, mos_mp2, energy_ex, energy_xc_admm, t3, unit_nr)
...
Definition: rpa_gw_sigma_x.F:112
rpa_gw
Routines for GW, continuous development [Jan Wilhelm].
Definition: rpa_gw.F:14
rpa_gw::trafo_to_mo_and_kpoints
subroutine, public trafo_to_mo_and_kpoints(qs_env, mat_self_energy_ao_ao, vec_Sigma, homo, gw_corr_lev_occ, gw_corr_lev_virt, ispin)
...
Definition: rpa_gw.F:6313
rpa_gw::compute_minus_vxc_kpoints
subroutine, public compute_minus_vxc_kpoints(qs_env)
...
Definition: rpa_gw.F:6188