Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines to calculate EXX within GW
10 : !> \par History
11 : !> 07.2020 separated from mp2.F [F. Stein, code by Jan Wilhelm]
12 : !> 07.2024 determine number of corrected MOs from BSE cutoffs [Maximilian Graml]
13 : !> \author Jan Wilhelm, Frederick Stein
14 : ! **************************************************************************************************
15 : MODULE rpa_gw_sigma_x
16 : USE admm_methods, ONLY: admm_mo_merge_ks_matrix
17 : USE admm_types, ONLY: admm_type,&
18 : get_admm_env
19 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale_and_add_fm
20 : USE cp_cfm_types, ONLY: cp_cfm_create,&
21 : cp_cfm_get_info,&
22 : cp_cfm_release,&
23 : cp_cfm_type
24 : USE cp_control_types, ONLY: dft_control_type
25 : USE cp_dbcsr_api, ONLY: &
26 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_multiply, dbcsr_p_type, &
27 : dbcsr_release, dbcsr_release_p, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
28 : dbcsr_type_symmetric
29 : USE cp_dbcsr_contrib, ONLY: dbcsr_get_diag
30 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
31 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
32 : copy_fm_to_dbcsr,&
33 : dbcsr_allocate_matrix_set,&
34 : dbcsr_deallocate_matrix_set
35 : USE cp_files, ONLY: close_file,&
36 : open_file
37 : USE cp_fm_struct, ONLY: cp_fm_struct_type
38 : USE cp_fm_types, ONLY: cp_fm_create,&
39 : cp_fm_get_info,&
40 : cp_fm_release,&
41 : cp_fm_type
42 : USE hfx_energy_potential, ONLY: integrate_four_center
43 : USE hfx_exx, ONLY: calc_exx_admm_xc_contributions,&
44 : exx_post_hfx,&
45 : exx_pre_hfx
46 : USE hfx_ri, ONLY: hfx_ri_update_ks
47 : USE input_constants, ONLY: do_admm_basis_projection,&
48 : do_admm_purify_none,&
49 : gw_print_exx,&
50 : gw_read_exx,&
51 : xc_none
52 : USE input_section_types, ONLY: section_vals_get,&
53 : section_vals_get_subs_vals,&
54 : section_vals_type,&
55 : section_vals_val_get,&
56 : section_vals_val_set
57 : USE kinds, ONLY: dp
58 : USE kpoint_methods, ONLY: rskp_transform
59 : USE kpoint_types, ONLY: get_kpoint_info,&
60 : kpoint_env_type,&
61 : kpoint_type
62 : USE machine, ONLY: m_walltime
63 : USE mathconstants, ONLY: gaussi,&
64 : z_one,&
65 : z_zero
66 : USE message_passing, ONLY: mp_para_env_type
67 : USE mo_window, ONLY: determine_mo_window,&
68 : mo_window_type
69 : USE mp2_integrals, ONLY: compute_kpoints
70 : USE mp2_ri_2c, ONLY: trunc_coulomb_for_exchange
71 : USE mp2_types, ONLY: mp2_type
72 : USE parallel_gemm_api, ONLY: parallel_gemm
73 : USE physcon, ONLY: evolt
74 : USE qs_energy_types, ONLY: qs_energy_type
75 : USE qs_environment_types, ONLY: get_qs_env,&
76 : qs_environment_type
77 : USE qs_ks_methods, ONLY: qs_ks_build_kohn_sham_matrix
78 : USE qs_ks_types, ONLY: qs_ks_env_type
79 : USE qs_mo_types, ONLY: get_mo_set,&
80 : mo_set_type
81 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
82 : USE qs_rho_types, ONLY: qs_rho_get,&
83 : qs_rho_type
84 : USE rpa_gw, ONLY: compute_minus_vxc_kpoints,&
85 : trafo_to_mo_and_kpoints
86 : USE rpa_gw_kpoints_util, ONLY: get_bandstruc_and_k_dependent_MOs
87 :
88 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
89 :
90 : #include "./base/base_uses.f90"
91 :
92 : IMPLICIT NONE
93 :
94 : PRIVATE
95 :
96 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_gw_sigma_x'
97 :
98 : PUBLIC :: compute_vec_Sigma_x_minus_vxc_gw
99 :
100 : CONTAINS
101 :
102 : ! **************************************************************************************************
103 : !> \brief ...
104 : !> \param qs_env ...
105 : !> \param mp2_env ...
106 : !> \param mos_mp2 ...
107 : !> \param energy_ex ...
108 : !> \param energy_xc_admm ...
109 : !> \param t3 ...
110 : !> \param unit_nr ...
111 : !> \par History
112 : !> 04.2015 created
113 : !> \author Jan Wilhelm
114 : ! **************************************************************************************************
115 122 : SUBROUTINE compute_vec_Sigma_x_minus_vxc_gw(qs_env, mp2_env, mos_mp2, energy_ex, energy_xc_admm, t3, unit_nr)
116 : TYPE(qs_environment_type), POINTER :: qs_env
117 : TYPE(mp2_type) :: mp2_env
118 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos_mp2
119 : REAL(KIND=dp), INTENT(OUT) :: energy_ex, energy_xc_admm(2), t3
120 : INTEGER, INTENT(IN) :: unit_nr
121 :
122 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_vec_Sigma_x_minus_vxc_gw'
123 :
124 : CHARACTER(4) :: occ_virt
125 : CHARACTER(LEN=40) :: line
126 : INTEGER :: dimen, gw_corr_lev_occ, gw_corr_lev_tot, gw_corr_lev_virt, handle, homo, &
127 : homo_reduced_bse, i_img, ikp, irep, ispin, iunit, max_corr_lev_occ, max_corr_lev_virt, &
128 : myfun, myfun_aux, myfun_prim, n_level_gw, n_level_gw_ref, n_rep_hf, nkp, nkp_Sigma, nmo, &
129 : nspins, print_exx, virtual_reduced_bse
130 : LOGICAL :: calc_ints, charge_constrain_tmp, do_admm_rpa, do_hfx, do_kpoints_cubic_RPA, &
131 : do_kpoints_from_Gamma, do_ri_Sigma_x, really_read_line
132 : REAL(KIND=dp) :: E_GAP_GW, E_HOMO_GW, E_LUMO_GW, eh1, ehfx, eigval_dft, eigval_hf_at_dft, &
133 : energy_exc, energy_exc1, energy_exc1_aux_fit, energy_exc_aux_fit, energy_total, &
134 : exx_minus_vxc, hfx_fraction, min_direct_HF_at_DFT_gap, t1, t2, tmp
135 122 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: matrix_tmp_2_diag
136 122 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Eigenval_kp_HF_at_DFT, vec_Sigma_x
137 122 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Eigenval_kp, vec_Sigma_x_minus_vxc_gw, &
138 122 : vec_Sigma_x_minus_vxc_gw_im
139 122 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
140 : TYPE(admm_type), POINTER :: admm_env
141 : TYPE(cp_fm_type), POINTER :: mo_coeff
142 122 : TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:) :: mat_exchange_for_kp_from_gamma
143 122 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, &
144 122 : matrix_ks_aux_fit_hfx, rho_ao, &
145 122 : rho_ao_aux_fit
146 122 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_2d, matrix_ks_kp_im, &
147 122 : matrix_ks_kp_re, matrix_ks_transl, matrix_sigma_x_minus_vxc, matrix_sigma_x_minus_vxc_im, &
148 122 : rho_ao_2d
149 : TYPE(dbcsr_type) :: matrix_tmp, matrix_tmp_2, mo_coeff_b
150 : TYPE(dft_control_type), POINTER :: dft_control
151 : TYPE(kpoint_type), POINTER :: kpoints, kpoints_Sigma
152 : TYPE(mo_window_type) :: bse_window
153 : TYPE(mp_para_env_type), POINTER :: para_env
154 : TYPE(qs_energy_type), POINTER :: energy
155 : TYPE(qs_ks_env_type), POINTER :: ks_env
156 : TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
157 : TYPE(section_vals_type), POINTER :: hfx_sections, input, xc_section, &
158 : xc_section_admm_aux, &
159 : xc_section_admm_prim
160 :
161 122 : NULLIFY (admm_env, matrix_ks, matrix_ks_aux_fit, rho_ao, matrix_sigma_x_minus_vxc, input, &
162 122 : xc_section, xc_section_admm_aux, xc_section_admm_prim, hfx_sections, rho, &
163 122 : dft_control, para_env, ks_env, mo_coeff, matrix_sigma_x_minus_vxc_im, matrix_ks_aux_fit_hfx, &
164 122 : rho_aux_fit, rho_ao_aux_fit)
165 :
166 122 : CALL timeset(routineN, handle)
167 :
168 122 : t1 = m_walltime()
169 :
170 122 : do_admm_rpa = mp2_env%ri_rpa%do_admm
171 122 : do_ri_Sigma_x = mp2_env%ri_g0w0%do_ri_Sigma_x
172 122 : do_kpoints_cubic_RPA = qs_env%mp2_env%ri_rpa_im_time%do_im_time_kpoints
173 122 : do_kpoints_from_Gamma = qs_env%mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma
174 122 : print_exx = mp2_env%ri_g0w0%print_exx
175 :
176 122 : IF (do_kpoints_cubic_RPA) THEN
177 0 : CPASSERT(do_ri_Sigma_x)
178 : END IF
179 :
180 : IF (do_kpoints_cubic_RPA) THEN
181 :
182 : CALL get_qs_env(qs_env, &
183 : admm_env=admm_env, &
184 : matrix_ks_kp=matrix_ks_transl, &
185 : rho=rho, &
186 : input=input, &
187 : dft_control=dft_control, &
188 : para_env=para_env, &
189 : kpoints=kpoints, &
190 : ks_env=ks_env, &
191 0 : energy=energy)
192 0 : nkp = kpoints%nkp
193 :
194 : ELSE
195 :
196 : CALL get_qs_env(qs_env, &
197 : admm_env=admm_env, &
198 : matrix_ks=matrix_ks, &
199 : rho=rho, &
200 : input=input, &
201 : dft_control=dft_control, &
202 : para_env=para_env, &
203 : ks_env=ks_env, &
204 122 : energy=energy)
205 122 : nkp = 1
206 122 : CALL qs_rho_get(rho, rho_ao=rho_ao)
207 :
208 122 : IF (do_admm_rpa) THEN
209 : CALL get_admm_env(admm_env, matrix_ks_aux_fit=matrix_ks_aux_fit, rho_aux_fit=rho_aux_fit, &
210 8 : matrix_ks_aux_fit_hfx=matrix_ks_aux_fit_hfx)
211 8 : CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux_fit)
212 :
213 : ! RPA/GW with ADMM for EXX or the exchange self-energy only implemented for
214 : ! ADMM_PURIFICATION_METHOD NONE
215 : ! METHOD BASIS_PROJECTION
216 : ! in the admm section
217 8 : CPASSERT(admm_env%purification_method == do_admm_purify_none)
218 8 : CPASSERT(dft_control%admm_control%method == do_admm_basis_projection)
219 : END IF
220 : END IF
221 :
222 122 : nspins = dft_control%nspins
223 :
224 : ! safe ks matrix for later: we will transform matrix_ks
225 : ! to T-cell index and then to k-points for band structure calculation
226 122 : IF (do_kpoints_from_Gamma) THEN
227 : ! not yet there: open shell
228 66 : ALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_ks(nspins))
229 34 : DO ispin = 1, nspins
230 18 : NULLIFY (qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix)
231 18 : ALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix)
232 : CALL dbcsr_create(qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix, &
233 18 : template=matrix_ks(ispin)%matrix)
234 : CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, &
235 34 : qs_env%mp2_env%ri_g0w0%matrix_ks(ispin)%matrix)
236 :
237 : END DO
238 : END IF
239 :
240 122 : IF (do_kpoints_cubic_RPA) THEN
241 :
242 0 : CALL allocate_matrix_ks_kp(matrix_ks_transl, matrix_ks_kp_re, matrix_ks_kp_im, kpoints)
243 0 : CALL transform_matrix_ks_to_kp(matrix_ks_transl, matrix_ks_kp_re, matrix_ks_kp_im, kpoints)
244 :
245 0 : DO ispin = 1, nspins
246 0 : DO i_img = 1, SIZE(matrix_ks_transl, 2)
247 0 : CALL dbcsr_set(matrix_ks_transl(ispin, i_img)%matrix, 0.0_dp)
248 : END DO
249 : END DO
250 :
251 : END IF
252 :
253 : ! initialize matrix_sigma_x_minus_vxc
254 122 : NULLIFY (matrix_sigma_x_minus_vxc)
255 122 : CALL dbcsr_allocate_matrix_set(matrix_sigma_x_minus_vxc, nspins, nkp)
256 122 : IF (do_kpoints_cubic_RPA) THEN
257 0 : NULLIFY (matrix_sigma_x_minus_vxc_im)
258 0 : CALL dbcsr_allocate_matrix_set(matrix_sigma_x_minus_vxc_im, nspins, nkp)
259 : END IF
260 :
261 264 : DO ispin = 1, nspins
262 406 : DO ikp = 1, nkp
263 :
264 284 : IF (do_kpoints_cubic_RPA) THEN
265 :
266 0 : ALLOCATE (matrix_sigma_x_minus_vxc(ispin, ikp)%matrix)
267 : CALL dbcsr_create(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, &
268 : template=matrix_ks_kp_re(1, 1)%matrix, &
269 0 : matrix_type=dbcsr_type_symmetric)
270 :
271 0 : CALL dbcsr_copy(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, matrix_ks_kp_re(ispin, ikp)%matrix)
272 0 : CALL dbcsr_set(matrix_ks_kp_re(ispin, ikp)%matrix, 0.0_dp)
273 :
274 0 : ALLOCATE (matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix)
275 : CALL dbcsr_create(matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix, &
276 : template=matrix_ks_kp_im(1, 1)%matrix, &
277 0 : matrix_type=dbcsr_type_antisymmetric)
278 :
279 0 : CALL dbcsr_copy(matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix, matrix_ks_kp_im(ispin, ikp)%matrix)
280 0 : CALL dbcsr_set(matrix_ks_kp_im(ispin, ikp)%matrix, 0.0_dp)
281 :
282 : ELSE
283 :
284 142 : ALLOCATE (matrix_sigma_x_minus_vxc(ispin, ikp)%matrix)
285 : CALL dbcsr_create(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, &
286 142 : template=matrix_ks(1)%matrix)
287 :
288 142 : CALL dbcsr_copy(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, matrix_ks(ispin)%matrix)
289 142 : CALL dbcsr_set(matrix_ks(ispin)%matrix, 0.0_dp)
290 :
291 : END IF
292 :
293 : END DO
294 : END DO
295 :
296 : ! set DFT functional to none and hfx_fraction to zero
297 122 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
298 122 : CALL section_vals_get(hfx_sections, explicit=do_hfx)
299 :
300 122 : IF (do_hfx) THEN
301 18 : hfx_fraction = qs_env%x_data(1, 1)%general_parameter%fraction
302 54 : qs_env%x_data(:, :)%general_parameter%fraction = 0.0_dp
303 : END IF
304 122 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
305 : CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
306 122 : i_val=myfun)
307 : CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
308 122 : i_val=xc_none)
309 :
310 : ! in ADMM, also set the XC functional for ADMM correction to none
311 : ! do not do this if we do ADMM for Sigma_x
312 122 : IF (dft_control%do_admm) THEN
313 : xc_section_admm_aux => section_vals_get_subs_vals(admm_env%xc_section_aux, &
314 8 : "XC_FUNCTIONAL")
315 : CALL section_vals_val_get(xc_section_admm_aux, "_SECTION_PARAMETERS_", &
316 8 : i_val=myfun_aux)
317 : CALL section_vals_val_set(xc_section_admm_aux, "_SECTION_PARAMETERS_", &
318 8 : i_val=xc_none)
319 :
320 : ! the same for the primary basis
321 : xc_section_admm_prim => section_vals_get_subs_vals(admm_env%xc_section_primary, &
322 8 : "XC_FUNCTIONAL")
323 : CALL section_vals_val_get(xc_section_admm_prim, "_SECTION_PARAMETERS_", &
324 8 : i_val=myfun_prim)
325 : CALL section_vals_val_set(xc_section_admm_prim, "_SECTION_PARAMETERS_", &
326 8 : i_val=xc_none)
327 :
328 : ! for ADMMQ/S, set the charge_constrain to false (otherwise wrong results)
329 8 : charge_constrain_tmp = .FALSE.
330 8 : IF (admm_env%charge_constrain) THEN
331 0 : admm_env%charge_constrain = .FALSE.
332 0 : charge_constrain_tmp = .TRUE.
333 : END IF
334 :
335 : END IF
336 :
337 : ! if we do ADMM for Sigma_x, we write the ADMM correction into matrix_ks_aux_fit
338 : ! and therefore we should set it to zero
339 122 : IF (do_admm_rpa) THEN
340 18 : DO ispin = 1, nspins
341 18 : CALL dbcsr_set(matrix_ks_aux_fit(ispin)%matrix, 0.0_dp)
342 : END DO
343 : END IF
344 :
345 122 : IF (.NOT. mp2_env%ri_g0w0%update_xc_energy) THEN
346 96 : energy_total = energy%total
347 96 : energy_exc = energy%exc
348 96 : energy_exc1 = energy%exc1
349 96 : energy_exc_aux_fit = energy%ex
350 96 : energy_exc1_aux_fit = energy%exc_aux_fit
351 96 : energy_ex = energy%exc1_aux_fit
352 : END IF
353 :
354 : ! Remove the Exchange-correlation energy contributions from the total energy
355 : energy%total = energy%total - (energy%exc + energy%exc1 + energy%ex + &
356 122 : energy%exc_aux_fit + energy%exc1_aux_fit)
357 :
358 : ! calculate KS-matrix without XC and without HF
359 : CALL qs_ks_build_kohn_sham_matrix(qs_env=qs_env, calculate_forces=.FALSE., &
360 122 : just_energy=.FALSE.)
361 :
362 122 : IF (.NOT. mp2_env%ri_g0w0%update_xc_energy) THEN
363 96 : energy%exc = energy_exc
364 96 : energy%exc1 = energy_exc1
365 96 : energy%exc_aux_fit = energy_ex
366 96 : energy%exc1_aux_fit = energy_exc_aux_fit
367 96 : energy%ex = energy_exc1_aux_fit
368 96 : energy%total = energy_total
369 : END IF
370 :
371 : ! set the DFT functional and HF fraction back
372 : CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
373 122 : i_val=myfun)
374 122 : IF (do_hfx) THEN
375 54 : qs_env%x_data(:, :)%general_parameter%fraction = hfx_fraction
376 : END IF
377 :
378 122 : IF (dft_control%do_admm) THEN
379 : xc_section_admm_aux => section_vals_get_subs_vals(admm_env%xc_section_aux, &
380 8 : "XC_FUNCTIONAL")
381 : xc_section_admm_prim => section_vals_get_subs_vals(admm_env%xc_section_primary, &
382 8 : "XC_FUNCTIONAL")
383 :
384 : CALL section_vals_val_set(xc_section_admm_aux, "_SECTION_PARAMETERS_", &
385 8 : i_val=myfun_aux)
386 : CALL section_vals_val_set(xc_section_admm_prim, "_SECTION_PARAMETERS_", &
387 8 : i_val=myfun_prim)
388 8 : IF (charge_constrain_tmp) THEN
389 0 : admm_env%charge_constrain = .TRUE.
390 : END IF
391 : END IF
392 :
393 122 : IF (do_kpoints_cubic_RPA) THEN
394 0 : CALL transform_matrix_ks_to_kp(matrix_ks_transl, matrix_ks_kp_re, matrix_ks_kp_im, kpoints)
395 : END IF
396 :
397 : ! remove the single-particle part (kin. En + Hartree pot) and change the sign
398 264 : DO ispin = 1, nspins
399 264 : IF (do_kpoints_cubic_RPA) THEN
400 0 : DO ikp = 1, nkp
401 0 : CALL dbcsr_add(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, matrix_ks_kp_re(ispin, ikp)%matrix, -1.0_dp, 1.0_dp)
402 0 : 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)
403 : END DO
404 : ELSE
405 142 : CALL dbcsr_add(matrix_sigma_x_minus_vxc(ispin, 1)%matrix, matrix_ks(ispin)%matrix, -1.0_dp, 1.0_dp)
406 : END IF
407 : END DO
408 :
409 122 : IF (do_kpoints_cubic_RPA) THEN
410 :
411 : CALL transform_sigma_x_minus_vxc_to_MO_basis(kpoints, matrix_sigma_x_minus_vxc, &
412 : matrix_sigma_x_minus_vxc_im, &
413 : vec_Sigma_x_minus_vxc_gw, &
414 : vec_Sigma_x_minus_vxc_gw_im, &
415 0 : para_env, nmo, mp2_env)
416 :
417 : ELSE
418 :
419 264 : DO ispin = 1, nspins
420 142 : CALL dbcsr_set(matrix_ks(ispin)%matrix, 0.0_dp)
421 264 : IF (do_admm_rpa) THEN
422 10 : CALL dbcsr_set(matrix_ks_aux_fit(ispin)%matrix, 0.0_dp)
423 : END IF
424 : END DO
425 :
426 122 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
427 :
428 122 : CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
429 :
430 : ! in most cases, we calculate the exchange self-energy here. But if we do only RI for
431 : ! the exchange self-energy, we do not calculate exchange here
432 122 : ehfx = 0.0_dp
433 122 : IF (.NOT. do_ri_Sigma_x) THEN
434 :
435 48 : CALL exx_pre_hfx(hfx_sections, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx)
436 48 : calc_ints = .NOT. qs_env%mp2_env%ri_rpa%reuse_hfx
437 :
438 : ! add here HFX (=Sigma_exchange) to matrix_sigma_x_minus_vxc
439 96 : DO irep = 1, n_rep_hf
440 48 : IF (do_admm_rpa) THEN
441 8 : matrix_ks_2d(1:nspins, 1:1) => matrix_ks_aux_fit(1:nspins)
442 8 : rho_ao_2d(1:nspins, 1:1) => rho_ao_aux_fit(1:nspins)
443 : ELSE
444 40 : matrix_ks_2d(1:nspins, 1:1) => matrix_ks(1:nspins)
445 40 : rho_ao_2d(1:nspins, 1:1) => rho_ao(1:nspins)
446 : END IF
447 :
448 96 : IF (qs_env%mp2_env%ri_rpa%x_data(irep, 1)%do_hfx_ri) THEN
449 : CALL hfx_ri_update_ks(qs_env, qs_env%mp2_env%ri_rpa%x_data(irep, 1)%ri_data, matrix_ks_2d, ehfx, &
450 : rho_ao=rho_ao_2d, geometry_did_change=calc_ints, nspins=nspins, &
451 0 : hf_fraction=qs_env%mp2_env%ri_rpa%x_data(irep, 1)%general_parameter%fraction)
452 :
453 0 : IF (do_admm_rpa) THEN
454 : !for ADMMS, we need the exchange matrix k(d) for both spins
455 0 : DO ispin = 1, nspins
456 : CALL dbcsr_copy(matrix_ks_aux_fit_hfx(ispin)%matrix, matrix_ks_2d(ispin, 1)%matrix, &
457 0 : name="HF exch. part of matrix_ks_aux_fit for ADMMS")
458 : END DO
459 : END IF
460 : ELSE
461 : CALL integrate_four_center(qs_env, qs_env%mp2_env%ri_rpa%x_data, matrix_ks_2d, eh1, &
462 : rho_ao_2d, hfx_sections, &
463 : para_env, calc_ints, irep, .TRUE., &
464 48 : ispin=1)
465 48 : ehfx = ehfx + eh1
466 : END IF
467 : END DO
468 :
469 : !ADMM XC correction
470 48 : IF (do_admm_rpa) THEN
471 : CALL calc_exx_admm_xc_contributions(qs_env=qs_env, &
472 : matrix_prim=matrix_ks, &
473 : matrix_aux=matrix_ks_aux_fit, &
474 : x_data=qs_env%mp2_env%ri_rpa%x_data, &
475 : exc=energy_xc_admm(1), &
476 : exc_aux_fit=energy_xc_admm(2), &
477 : calc_forces=.FALSE., &
478 8 : use_virial=.FALSE.)
479 : END IF
480 :
481 48 : IF (do_kpoints_from_Gamma .AND. print_exx == gw_print_exx) THEN
482 0 : ALLOCATE (mat_exchange_for_kp_from_gamma(1))
483 :
484 0 : DO ispin = 1, 1
485 0 : NULLIFY (mat_exchange_for_kp_from_gamma(ispin)%matrix)
486 0 : ALLOCATE (mat_exchange_for_kp_from_gamma(ispin)%matrix)
487 0 : CALL dbcsr_create(mat_exchange_for_kp_from_gamma(ispin)%matrix, template=matrix_ks(ispin)%matrix)
488 0 : CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, mat_exchange_for_kp_from_gamma(ispin)%matrix)
489 : END DO
490 :
491 : END IF
492 :
493 48 : CALL exx_post_hfx(qs_env, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx)
494 : END IF
495 :
496 122 : energy_ex = ehfx
497 :
498 : ! transform Fock-Matrix (calculated in integrate_four_center, written in matrix_ks_aux_fit in case
499 : ! of ADMM) from ADMM basis to primary basis
500 122 : IF (do_admm_rpa) THEN
501 8 : CALL admm_mo_merge_ks_matrix(qs_env)
502 : END IF
503 :
504 264 : DO ispin = 1, nspins
505 264 : CALL dbcsr_add(matrix_sigma_x_minus_vxc(ispin, 1)%matrix, matrix_ks(ispin)%matrix, 1.0_dp, 1.0_dp)
506 : END DO
507 :
508 : ! safe matrix_sigma_x_minus_vxc for later: for example, we will transform matrix_sigma_x_minus_vxc
509 : ! to T-cell index and then to k-points for band structure calculation
510 122 : IF (do_kpoints_from_Gamma) THEN
511 : ! not yet there: open shell
512 66 : ALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(nspins))
513 34 : DO ispin = 1, nspins
514 18 : NULLIFY (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix)
515 18 : ALLOCATE (qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix)
516 : CALL dbcsr_create(qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix, &
517 18 : template=matrix_ks(ispin)%matrix)
518 :
519 : CALL dbcsr_desymmetrize(matrix_sigma_x_minus_vxc(ispin, 1)%matrix, &
520 34 : qs_env%mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix)
521 :
522 : END DO
523 : END IF
524 :
525 122 : CALL dbcsr_desymmetrize(matrix_ks(1)%matrix, mo_coeff_b)
526 122 : CALL dbcsr_set(mo_coeff_b, 0.0_dp)
527 :
528 : ! Transform matrix_sigma_x_minus_vxc to MO basis
529 264 : DO ispin = 1, nspins
530 :
531 : CALL get_mo_set(mo_set=mos_mp2(ispin), &
532 : mo_coeff=mo_coeff, &
533 : eigenvalues=mo_eigenvalues, &
534 : nmo=nmo, &
535 : homo=homo, &
536 142 : nao=dimen)
537 :
538 142 : IF (ispin == 1) THEN
539 :
540 610 : ALLOCATE (vec_Sigma_x_minus_vxc_gw(nmo, nspins, nkp))
541 122 : vec_Sigma_x_minus_vxc_gw = 0.0_dp
542 :
543 366 : ALLOCATE (matrix_tmp_2_diag(dimen))
544 : END IF
545 :
546 142 : CALL dbcsr_set(mo_coeff_b, 0.0_dp)
547 142 : CALL copy_fm_to_dbcsr(mo_coeff, mo_coeff_b, keep_sparsity=.FALSE.)
548 :
549 : ! initialize matrix_tmp and matrix_tmp2
550 142 : IF (ispin == 1) THEN
551 122 : CALL dbcsr_create(matrix_tmp, template=mo_coeff_b)
552 122 : CALL dbcsr_copy(matrix_tmp, mo_coeff_b)
553 122 : CALL dbcsr_set(matrix_tmp, 0.0_dp)
554 :
555 122 : CALL dbcsr_create(matrix_tmp_2, template=mo_coeff_b)
556 122 : CALL dbcsr_copy(matrix_tmp_2, mo_coeff_b)
557 122 : CALL dbcsr_set(matrix_tmp_2, 0.0_dp)
558 : END IF
559 :
560 142 : gw_corr_lev_occ = mp2_env%ri_g0w0%corr_mos_occ
561 142 : gw_corr_lev_virt = mp2_env%ri_g0w0%corr_mos_virt
562 :
563 : ! If SVD is used to invert overlap matrix (for CHOLESKY OFF), some MOs are removed
564 : ! Therefore, setting the number of gw_corr_lev_virt simply to dimen - homo leads to index problems
565 : ! Instead, we take into account the removed MOs
566 142 : max_corr_lev_occ = homo
567 142 : max_corr_lev_virt = nmo - homo
568 :
569 : ! If BSE is invoked, manipulate corrected MO number
570 142 : IF (mp2_env%bse%do_bse) THEN
571 : ! Logic: If cutoff is negative, all MOs are included in BSE, i.e. we need to correct them all
572 : ! If cutoff is positive, we can reduce the number of MOs to be corrected and force gw_corr_lev_...
573 : ! to a sufficiently large number by setting it to -2 and read indices afterwards
574 : ! Handling for occupied levels
575 56 : IF (mp2_env%bse%bse_cutoff_occ < 0) THEN
576 22 : gw_corr_lev_occ = -1
577 : ELSE
578 34 : IF (gw_corr_lev_occ > 0) THEN
579 34 : gw_corr_lev_occ = -2
580 : END IF
581 : END IF
582 : ! Handling for virtual levels
583 56 : IF (mp2_env%bse%bse_cutoff_empty < 0) THEN
584 12 : gw_corr_lev_virt = -1
585 : ELSE
586 44 : IF (gw_corr_lev_virt > 0) THEN
587 44 : gw_corr_lev_virt = -2
588 : END IF
589 : END IF
590 :
591 : ! Obtain indices from DFT if gw_corr... are set to -2
592 : CALL determine_mo_window(mo_eigenvalues, nmo, homo, &
593 56 : mp2_env%bse%bse_cutoff_occ, mp2_env%bse%bse_cutoff_empty, bse_window)
594 56 : homo_reduced_bse = homo - bse_window%first_mo + 1
595 56 : virtual_reduced_bse = bse_window%last_mo - homo
596 56 : IF (gw_corr_lev_occ == -2) THEN
597 34 : CPWARN("BSE cutoff overwrites user input for CORR_MOS_OCC")
598 34 : gw_corr_lev_occ = homo_reduced_bse
599 : END IF
600 56 : IF (gw_corr_lev_virt == -2) THEN
601 44 : CPWARN("BSE cutoff overwrites user input for CORR_MOS_VIRT")
602 44 : gw_corr_lev_virt = virtual_reduced_bse
603 : END IF
604 : END IF
605 :
606 : ! if requested number of occ/virt levels for correction either exceed the number of
607 : ! occ/virt levels or the requested number is negative, default to correct all
608 : ! occ/virt level energies
609 142 : IF (gw_corr_lev_occ > homo .OR. gw_corr_lev_occ < 0) gw_corr_lev_occ = max_corr_lev_occ
610 142 : IF (gw_corr_lev_virt > max_corr_lev_virt .OR. gw_corr_lev_virt < 0) gw_corr_lev_virt = max_corr_lev_virt
611 142 : IF (ispin == 1) THEN
612 122 : mp2_env%ri_g0w0%corr_mos_occ = gw_corr_lev_occ
613 122 : mp2_env%ri_g0w0%corr_mos_virt = gw_corr_lev_virt
614 20 : ELSE IF (ispin == 2) THEN
615 : ! ensure that the total number of corrected MOs is the same for alpha and beta, important
616 : ! for parallelization
617 20 : IF (mp2_env%ri_g0w0%corr_mos_occ + mp2_env%ri_g0w0%corr_mos_virt /= &
618 : gw_corr_lev_occ + gw_corr_lev_virt) THEN
619 10 : gw_corr_lev_virt = mp2_env%ri_g0w0%corr_mos_occ + mp2_env%ri_g0w0%corr_mos_virt - gw_corr_lev_occ
620 : END IF
621 20 : mp2_env%ri_g0w0%corr_mos_occ_beta = gw_corr_lev_occ
622 20 : mp2_env%ri_g0w0%corr_mos_virt_beta = gw_corr_lev_virt
623 :
624 : END IF
625 :
626 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_sigma_x_minus_vxc(ispin, 1)%matrix, &
627 : mo_coeff_b, 0.0_dp, matrix_tmp, first_column=homo + 1 - gw_corr_lev_occ, &
628 142 : last_column=homo + gw_corr_lev_virt)
629 :
630 : CALL dbcsr_multiply('T', 'N', 1.0_dp, mo_coeff_b, &
631 : matrix_tmp, 0.0_dp, matrix_tmp_2, first_row=homo + 1 - gw_corr_lev_occ, &
632 142 : last_row=homo + gw_corr_lev_virt)
633 :
634 142 : CALL dbcsr_get_diag(matrix_tmp_2, matrix_tmp_2_diag)
635 3742 : vec_Sigma_x_minus_vxc_gw(1:nmo, ispin, 1) = matrix_tmp_2_diag(1:nmo)
636 :
637 142 : CALL dbcsr_set(matrix_tmp, 0.0_dp)
638 406 : CALL dbcsr_set(matrix_tmp_2, 0.0_dp)
639 :
640 : END DO
641 :
642 122 : CALL para_env%sum(vec_Sigma_x_minus_vxc_gw)
643 :
644 : END IF
645 :
646 122 : CALL dbcsr_release(mo_coeff_b)
647 122 : CALL dbcsr_release(matrix_tmp)
648 122 : CALL dbcsr_release(matrix_tmp_2)
649 122 : IF (do_kpoints_cubic_RPA) THEN
650 0 : CALL dbcsr_deallocate_matrix_set(matrix_ks_kp_re)
651 0 : CALL dbcsr_deallocate_matrix_set(matrix_ks_kp_im)
652 : END IF
653 :
654 264 : DO ispin = 1, nspins
655 406 : DO ikp = 1, nkp
656 142 : CALL dbcsr_release_p(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix)
657 284 : IF (do_kpoints_cubic_RPA) THEN
658 0 : CALL dbcsr_release_p(matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix)
659 : END IF
660 : END DO
661 : END DO
662 :
663 610 : ALLOCATE (mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(nmo, nspins, nkp))
664 :
665 122 : IF (print_exx == gw_print_exx) THEN
666 :
667 0 : IF (do_kpoints_from_Gamma) THEN
668 :
669 0 : gw_corr_lev_tot = gw_corr_lev_occ + gw_corr_lev_virt
670 :
671 : CALL get_qs_env(qs_env=qs_env, &
672 0 : kpoints=kpoints)
673 :
674 0 : CALL trunc_coulomb_for_exchange(qs_env)
675 :
676 0 : CALL compute_kpoints(qs_env, kpoints, unit_nr)
677 :
678 0 : ALLOCATE (Eigenval_kp(nmo, 1, nspins))
679 :
680 0 : CALL get_bandstruc_and_k_dependent_MOs(qs_env, Eigenval_kp)
681 :
682 0 : CALL compute_minus_vxc_kpoints(qs_env)
683 :
684 0 : nkp_Sigma = SIZE(Eigenval_kp, 2)
685 :
686 0 : ALLOCATE (vec_Sigma_x(nmo, nkp_Sigma))
687 0 : vec_Sigma_x(:, :) = 0.0_dp
688 :
689 : CALL trafo_to_mo_and_kpoints(qs_env, &
690 : mat_exchange_for_kp_from_gamma(1)%matrix, &
691 : vec_Sigma_x(homo - gw_corr_lev_occ + 1:homo + gw_corr_lev_virt, :), &
692 0 : homo, gw_corr_lev_occ, gw_corr_lev_virt, 1)
693 :
694 0 : CALL dbcsr_release(mat_exchange_for_kp_from_gamma(1)%matrix)
695 0 : DEALLOCATE (mat_exchange_for_kp_from_gamma(1)%matrix)
696 0 : DEALLOCATE (mat_exchange_for_kp_from_gamma)
697 :
698 0 : DEALLOCATE (vec_Sigma_x_minus_vxc_gw)
699 :
700 0 : ALLOCATE (vec_Sigma_x_minus_vxc_gw(nmo, nspins, nkp_Sigma))
701 :
702 : vec_Sigma_x_minus_vxc_gw(:, 1, :) = vec_Sigma_x(:, :) + &
703 0 : qs_env%mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, :)
704 :
705 0 : kpoints_Sigma => qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma
706 :
707 : ELSE
708 :
709 0 : nkp_Sigma = 1
710 :
711 : END IF
712 :
713 0 : IF (unit_nr > 0) THEN
714 :
715 0 : ALLOCATE (Eigenval_kp_HF_at_DFT(nmo, nkp_Sigma))
716 0 : Eigenval_kp_HF_at_DFT(:, :) = Eigenval_kp(:, :, 1) + vec_Sigma_x_minus_vxc_gw(:, 1, :)
717 :
718 0 : min_direct_HF_at_DFT_gap = 100.0_dp
719 :
720 0 : WRITE (unit_nr, '(T3,A)') ''
721 0 : WRITE (unit_nr, '(T3,A)') 'Exchange energies'
722 0 : WRITE (unit_nr, '(T3,A)') '-----------------'
723 0 : WRITE (unit_nr, '(T3,A)') ''
724 0 : WRITE (unit_nr, '(T6,2A)') 'MO e_n^DFT Sigma_x-vxc e_n^HF@DFT'
725 0 : DO ikp = 1, nkp_Sigma
726 0 : IF (nkp_Sigma > 1) THEN
727 0 : WRITE (unit_nr, '(T3,A)') ''
728 0 : WRITE (unit_nr, '(T3,A7,I3,A3,I3,A8,3F7.3,A12,3F7.3)') 'Kpoint ', ikp, ' /', nkp_Sigma, &
729 0 : ' xkp =', kpoints_Sigma%xkp(1, ikp), kpoints_Sigma%xkp(2, ikp), &
730 0 : kpoints_Sigma%xkp(3, ikp), ' and xkp =', -kpoints_Sigma%xkp(1, ikp), &
731 0 : -kpoints_Sigma%xkp(2, ikp), -kpoints_Sigma%xkp(3, ikp)
732 0 : WRITE (unit_nr, '(T3,A)') ''
733 : END IF
734 0 : DO n_level_gw = 1, gw_corr_lev_occ + gw_corr_lev_virt
735 :
736 0 : n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
737 0 : IF (n_level_gw <= gw_corr_lev_occ) THEN
738 0 : occ_virt = 'occ'
739 : ELSE
740 0 : occ_virt = 'vir'
741 : END IF
742 :
743 0 : eigval_dft = Eigenval_kp(n_level_gw_ref, ikp, 1)*evolt
744 0 : exx_minus_vxc = REAL(vec_Sigma_x_minus_vxc_gw(n_level_gw_ref, 1, ikp)*evolt, kind=dp)
745 0 : eigval_hf_at_dft = Eigenval_kp_HF_at_DFT(n_level_gw_ref, ikp)*evolt
746 :
747 : WRITE (unit_nr, '(T4,I4,3A,3F21.3,3F21.3,3F21.3)') &
748 0 : n_level_gw_ref, ' ( ', occ_virt, ') ', eigval_dft, exx_minus_vxc, eigval_hf_at_dft
749 :
750 : END DO
751 0 : E_HOMO_GW = MAXVAL(Eigenval_kp_HF_at_DFT(homo - gw_corr_lev_occ + 1:homo, ikp))
752 0 : E_LUMO_GW = MINVAL(Eigenval_kp_HF_at_DFT(homo + 1:homo + gw_corr_lev_virt, ikp))
753 0 : E_GAP_GW = E_LUMO_GW - E_HOMO_GW
754 : IF (E_GAP_GW < min_direct_HF_at_DFT_gap) min_direct_HF_at_DFT_gap = E_GAP_GW
755 0 : WRITE (unit_nr, '(T3,A)') ''
756 0 : WRITE (unit_nr, '(T3,A,F53.2)') 'HF@DFT HOMO-LUMO gap (eV)', E_GAP_GW*evolt
757 0 : WRITE (unit_nr, '(T3,A)') ''
758 : END DO
759 :
760 0 : WRITE (unit_nr, '(T3,A)') ''
761 0 : WRITE (unit_nr, '(T3,A)') ''
762 0 : WRITE (unit_nr, '(T3,A,F63.3)') 'HF@DFT direct bandgap (eV)', min_direct_HF_at_DFT_gap*evolt
763 :
764 0 : WRITE (unit_nr, '(T3,A)') ''
765 0 : WRITE (unit_nr, '(T3,A)') 'End of exchange energies'
766 0 : WRITE (unit_nr, '(T3,A)') '------------------------'
767 0 : WRITE (unit_nr, '(T3,A)') ''
768 :
769 0 : CPABORT('Stop after printing exchange energies.')
770 :
771 : ELSE
772 0 : CALL para_env%sync()
773 : END IF
774 :
775 : END IF
776 :
777 122 : IF (print_exx == gw_read_exx) THEN
778 :
779 0 : CALL open_file(unit_number=iunit, file_name="exx.out")
780 :
781 0 : really_read_line = .FALSE.
782 :
783 : DO WHILE (.TRUE.)
784 :
785 0 : READ (iunit, '(A)') line
786 :
787 0 : IF (line == " End of exchange energies ") EXIT
788 :
789 0 : IF (really_read_line) THEN
790 :
791 0 : READ (line(1:7), *) n_level_gw_ref
792 0 : READ (line(17:40), *) tmp
793 :
794 0 : DO ikp = 1, SIZE(vec_Sigma_x_minus_vxc_gw, 3)
795 0 : vec_Sigma_x_minus_vxc_gw(n_level_gw_ref, 1, ikp) = tmp/evolt
796 : END DO
797 :
798 : END IF
799 :
800 0 : IF (line == " MO Sigma_x-vxc ") really_read_line = .TRUE.
801 :
802 : END DO
803 :
804 0 : CALL close_file(iunit)
805 :
806 : END IF
807 :
808 : ! store vec_Sigma_x_minus_vxc_gw in the mp2_environment
809 3986 : mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, :, :) = vec_Sigma_x_minus_vxc_gw(:, :, :)
810 :
811 : ! clean up
812 122 : DEALLOCATE (matrix_sigma_x_minus_vxc, vec_Sigma_x_minus_vxc_gw)
813 122 : IF (do_kpoints_cubic_RPA) THEN
814 0 : DEALLOCATE (matrix_sigma_x_minus_vxc_im)
815 : END IF
816 :
817 122 : t2 = m_walltime()
818 :
819 122 : t3 = t2 - t1
820 :
821 122 : CALL timestop(handle)
822 :
823 488 : END SUBROUTINE compute_vec_Sigma_x_minus_vxc_gw
824 :
825 : ! **************************************************************************************************
826 : !> \brief ...
827 : !> \param kpoints ...
828 : !> \param matrix_sigma_x_minus_vxc ...
829 : !> \param matrix_sigma_x_minus_vxc_im ...
830 : !> \param vec_Sigma_x_minus_vxc_gw ...
831 : !> \param vec_Sigma_x_minus_vxc_gw_im ...
832 : !> \param para_env ...
833 : !> \param nmo ...
834 : !> \param mp2_env ...
835 : ! **************************************************************************************************
836 0 : SUBROUTINE transform_sigma_x_minus_vxc_to_MO_basis(kpoints, matrix_sigma_x_minus_vxc, &
837 : matrix_sigma_x_minus_vxc_im, vec_Sigma_x_minus_vxc_gw, &
838 : vec_Sigma_x_minus_vxc_gw_im, para_env, nmo, mp2_env)
839 :
840 : TYPE(kpoint_type), POINTER :: kpoints
841 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_sigma_x_minus_vxc, &
842 : matrix_sigma_x_minus_vxc_im
843 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: vec_Sigma_x_minus_vxc_gw, &
844 : vec_Sigma_x_minus_vxc_gw_im
845 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
846 : INTEGER :: nmo
847 : TYPE(mp2_type) :: mp2_env
848 :
849 : CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_sigma_x_minus_vxc_to_MO_basis'
850 :
851 : INTEGER :: dimen, gw_corr_lev_occ, gw_corr_lev_virt, handle, homo, i_global, iiB, ikp, &
852 : ispin, j_global, jjB, max_corr_lev_occ, max_corr_lev_virt, ncol_local, nkp, nrow_local, &
853 : nspins
854 : INTEGER, DIMENSION(2) :: kp_range
855 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
856 : REAL(KIND=dp) :: imval, reval
857 : TYPE(cp_cfm_type) :: cfm_mos, cfm_sigma_x_minus_vxc, &
858 : cfm_sigma_x_minus_vxc_mo_basis, cfm_tmp
859 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
860 : TYPE(cp_fm_type) :: fwork_im, fwork_re
861 : TYPE(kpoint_env_type), POINTER :: kp
862 : TYPE(mo_set_type), POINTER :: mo_set, mo_set_im, mo_set_re
863 :
864 0 : CALL timeset(routineN, handle)
865 :
866 0 : mo_set => kpoints%kp_env(1)%kpoint_env%mos(1, 1)
867 0 : CALL get_mo_set(mo_set, nmo=nmo)
868 :
869 0 : nspins = SIZE(matrix_sigma_x_minus_vxc, 1)
870 0 : CALL get_kpoint_info(kpoints, nkp=nkp, kp_range=kp_range)
871 :
872 : ! if this CPASSERT is wrong, please make sure that the kpoint group size PARALLEL_GROUP_SIZE
873 : ! in the kpoint environment &DFT &KPOINTS is -1
874 0 : CPASSERT(kp_range(1) == 1 .AND. kp_range(2) == nkp)
875 :
876 0 : ALLOCATE (vec_Sigma_x_minus_vxc_gw(nmo, nspins, nkp))
877 0 : vec_Sigma_x_minus_vxc_gw = 0.0_dp
878 :
879 0 : ALLOCATE (vec_Sigma_x_minus_vxc_gw_im(nmo, nspins, nkp))
880 0 : vec_Sigma_x_minus_vxc_gw_im = 0.0_dp
881 :
882 0 : CALL cp_fm_get_info(mo_set%mo_coeff, matrix_struct=matrix_struct)
883 0 : CALL cp_fm_create(fwork_re, matrix_struct)
884 0 : CALL cp_fm_create(fwork_im, matrix_struct)
885 0 : CALL cp_cfm_create(cfm_mos, matrix_struct)
886 0 : CALL cp_cfm_create(cfm_sigma_x_minus_vxc, matrix_struct)
887 0 : CALL cp_cfm_create(cfm_sigma_x_minus_vxc_mo_basis, matrix_struct)
888 0 : CALL cp_cfm_create(cfm_tmp, matrix_struct)
889 :
890 : CALL cp_cfm_get_info(matrix=cfm_sigma_x_minus_vxc_mo_basis, &
891 : nrow_local=nrow_local, &
892 : ncol_local=ncol_local, &
893 : row_indices=row_indices, &
894 0 : col_indices=col_indices)
895 :
896 : ! Transform matrix_sigma_x_minus_vxc to MO basis
897 0 : DO ikp = 1, nkp
898 :
899 0 : kp => kpoints%kp_env(ikp)%kpoint_env
900 :
901 0 : DO ispin = 1, nspins
902 :
903 : ! v_xc_n to fm matrix
904 0 : CALL copy_dbcsr_to_fm(matrix_sigma_x_minus_vxc(ispin, ikp)%matrix, fwork_re)
905 0 : CALL copy_dbcsr_to_fm(matrix_sigma_x_minus_vxc_im(ispin, ikp)%matrix, fwork_im)
906 :
907 0 : CALL cp_cfm_scale_and_add_fm(z_zero, cfm_sigma_x_minus_vxc, z_one, fwork_re)
908 0 : CALL cp_cfm_scale_and_add_fm(z_one, cfm_sigma_x_minus_vxc, gaussi, fwork_im)
909 :
910 : ! get real part (1) and imag. part (2) of the mo coeffs
911 0 : mo_set_re => kp%mos(1, ispin)
912 0 : mo_set_im => kp%mos(2, ispin)
913 :
914 0 : CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mos, z_one, mo_set_re%mo_coeff)
915 0 : CALL cp_cfm_scale_and_add_fm(z_one, cfm_mos, gaussi, mo_set_im%mo_coeff)
916 :
917 : ! tmp = V(k)*C(k)
918 : CALL parallel_gemm('N', 'N', nmo, nmo, nmo, z_one, cfm_sigma_x_minus_vxc, &
919 0 : cfm_mos, z_zero, cfm_tmp)
920 :
921 : ! V_n(k) = C^H(k)*tmp
922 : CALL parallel_gemm('C', 'N', nmo, nmo, nmo, z_one, cfm_mos, cfm_tmp, &
923 0 : z_zero, cfm_sigma_x_minus_vxc_mo_basis)
924 :
925 0 : DO jjB = 1, ncol_local
926 :
927 0 : j_global = col_indices(jjB)
928 :
929 0 : DO iiB = 1, nrow_local
930 :
931 0 : i_global = row_indices(iiB)
932 :
933 0 : IF (j_global == i_global .AND. i_global <= nmo) THEN
934 :
935 0 : reval = REAL(cfm_sigma_x_minus_vxc_mo_basis%local_data(iiB, jjB), kind=dp)
936 0 : imval = AIMAG(cfm_sigma_x_minus_vxc_mo_basis%local_data(iiB, jjB))
937 :
938 0 : vec_Sigma_x_minus_vxc_gw(i_global, ispin, ikp) = reval
939 0 : vec_Sigma_x_minus_vxc_gw_im(i_global, ispin, ikp) = imval
940 :
941 : END IF
942 :
943 : END DO
944 :
945 : END DO
946 :
947 : END DO
948 :
949 : END DO
950 :
951 0 : CALL para_env%sum(vec_Sigma_x_minus_vxc_gw)
952 0 : CALL para_env%sum(vec_Sigma_x_minus_vxc_gw_im)
953 : ! also adjust in the case of kpoints too big gw_corr_lev_occ and gw_corr_lev_virt
954 0 : DO ispin = 1, nspins
955 : CALL get_mo_set(mo_set=kpoints%kp_env(1)%kpoint_env%mos(ispin, 1), &
956 0 : homo=homo, nao=dimen)
957 : ! If SVD is used to invert overlap matrix (for CHOLESKY OFF), some MOs are removed
958 : ! Therefore, setting the number of gw_corr_lev_virt simply to dimen - homo leads to index problems
959 : ! Instead, we take into account the removed MOs
960 0 : max_corr_lev_occ = homo
961 0 : max_corr_lev_virt = nmo - homo
962 :
963 0 : gw_corr_lev_occ = mp2_env%ri_g0w0%corr_mos_occ
964 0 : gw_corr_lev_virt = mp2_env%ri_g0w0%corr_mos_virt
965 : ! if corrected occ/virt levels exceed the number of occ/virt levels or are negative,
966 : ! correct all occ/virt level energies
967 0 : IF (gw_corr_lev_occ > homo .OR. gw_corr_lev_occ < 0) gw_corr_lev_occ = max_corr_lev_occ
968 0 : IF (gw_corr_lev_virt > max_corr_lev_virt .OR. gw_corr_lev_virt < 0) gw_corr_lev_virt = max_corr_lev_virt
969 0 : IF (ispin == 1) THEN
970 0 : mp2_env%ri_g0w0%corr_mos_occ = gw_corr_lev_occ
971 0 : mp2_env%ri_g0w0%corr_mos_virt = gw_corr_lev_virt
972 0 : ELSE IF (ispin == 2) THEN
973 : ! ensure that the total number of corrected MOs is the same for alpha and beta, important
974 : ! for parallelization
975 0 : IF (mp2_env%ri_g0w0%corr_mos_occ + mp2_env%ri_g0w0%corr_mos_virt /= &
976 : gw_corr_lev_occ + gw_corr_lev_virt) THEN
977 0 : gw_corr_lev_virt = mp2_env%ri_g0w0%corr_mos_occ + mp2_env%ri_g0w0%corr_mos_virt - gw_corr_lev_occ
978 : END IF
979 0 : mp2_env%ri_g0w0%corr_mos_occ_beta = gw_corr_lev_occ
980 0 : mp2_env%ri_g0w0%corr_mos_virt_beta = gw_corr_lev_virt
981 : END IF
982 : END DO
983 :
984 0 : CALL cp_fm_release(fwork_re)
985 0 : CALL cp_fm_release(fwork_im)
986 0 : CALL cp_cfm_release(cfm_mos)
987 0 : CALL cp_cfm_release(cfm_sigma_x_minus_vxc)
988 0 : CALL cp_cfm_release(cfm_sigma_x_minus_vxc_mo_basis)
989 0 : CALL cp_cfm_release(cfm_tmp)
990 :
991 0 : CALL timestop(handle)
992 :
993 0 : END SUBROUTINE transform_sigma_x_minus_vxc_to_MO_basis
994 :
995 : ! **************************************************************************************************
996 : !> \brief ...
997 : !> \param matrix_ks_transl ...
998 : !> \param matrix_ks_kp_re ...
999 : !> \param matrix_ks_kp_im ...
1000 : !> \param kpoints ...
1001 : ! **************************************************************************************************
1002 0 : SUBROUTINE transform_matrix_ks_to_kp(matrix_ks_transl, matrix_ks_kp_re, matrix_ks_kp_im, kpoints)
1003 :
1004 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_transl, matrix_ks_kp_re, &
1005 : matrix_ks_kp_im
1006 : TYPE(kpoint_type), POINTER :: kpoints
1007 :
1008 : CHARACTER(len=*), PARAMETER :: routineN = 'transform_matrix_ks_to_kp'
1009 :
1010 : INTEGER :: handle, ikp, ispin, nkp, nspin
1011 0 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1012 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1013 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1014 0 : POINTER :: sab_nl
1015 :
1016 0 : CALL timeset(routineN, handle)
1017 :
1018 0 : NULLIFY (sab_nl)
1019 0 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, sab_nl=sab_nl, cell_to_index=cell_to_index)
1020 :
1021 0 : CPASSERT(ASSOCIATED(sab_nl))
1022 :
1023 0 : nspin = SIZE(matrix_ks_transl, 1)
1024 :
1025 0 : DO ikp = 1, nkp
1026 0 : DO ispin = 1, nspin
1027 :
1028 0 : CALL dbcsr_set(matrix_ks_kp_re(ispin, ikp)%matrix, 0.0_dp)
1029 0 : CALL dbcsr_set(matrix_ks_kp_im(ispin, ikp)%matrix, 0.0_dp)
1030 : CALL rskp_transform(rmatrix=matrix_ks_kp_re(ispin, ikp)%matrix, &
1031 : cmatrix=matrix_ks_kp_im(ispin, ikp)%matrix, &
1032 : rsmat=matrix_ks_transl, ispin=ispin, &
1033 0 : xkp=xkp(1:3, ikp), cell_to_index=cell_to_index, sab_nl=sab_nl)
1034 :
1035 : END DO
1036 : END DO
1037 :
1038 0 : CALL timestop(handle)
1039 :
1040 0 : END SUBROUTINE transform_matrix_ks_to_kp
1041 :
1042 : ! **************************************************************************************************
1043 : !> \brief ...
1044 : !> \param matrix_ks_transl ...
1045 : !> \param matrix_ks_kp_re ...
1046 : !> \param matrix_ks_kp_im ...
1047 : !> \param kpoints ...
1048 : ! **************************************************************************************************
1049 0 : SUBROUTINE allocate_matrix_ks_kp(matrix_ks_transl, matrix_ks_kp_re, matrix_ks_kp_im, kpoints)
1050 :
1051 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_transl, matrix_ks_kp_re, &
1052 : matrix_ks_kp_im
1053 : TYPE(kpoint_type), POINTER :: kpoints
1054 :
1055 : CHARACTER(len=*), PARAMETER :: routineN = 'allocate_matrix_ks_kp'
1056 :
1057 : INTEGER :: handle, ikp, ispin, nkp, nspin
1058 0 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1059 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
1060 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1061 0 : POINTER :: sab_nl
1062 :
1063 0 : CALL timeset(routineN, handle)
1064 :
1065 0 : NULLIFY (sab_nl)
1066 0 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, sab_nl=sab_nl, cell_to_index=cell_to_index)
1067 :
1068 0 : CPASSERT(ASSOCIATED(sab_nl))
1069 :
1070 0 : nspin = SIZE(matrix_ks_transl, 1)
1071 :
1072 0 : NULLIFY (matrix_ks_kp_re, matrix_ks_kp_im)
1073 0 : CALL dbcsr_allocate_matrix_set(matrix_ks_kp_re, nspin, nkp)
1074 0 : CALL dbcsr_allocate_matrix_set(matrix_ks_kp_im, nspin, nkp)
1075 :
1076 0 : DO ikp = 1, nkp
1077 0 : DO ispin = 1, nspin
1078 :
1079 0 : ALLOCATE (matrix_ks_kp_re(ispin, ikp)%matrix)
1080 0 : ALLOCATE (matrix_ks_kp_im(ispin, ikp)%matrix)
1081 :
1082 : CALL dbcsr_create(matrix_ks_kp_re(ispin, ikp)%matrix, &
1083 : template=matrix_ks_transl(1, 1)%matrix, &
1084 0 : matrix_type=dbcsr_type_symmetric)
1085 : CALL dbcsr_create(matrix_ks_kp_im(ispin, ikp)%matrix, &
1086 : template=matrix_ks_transl(1, 1)%matrix, &
1087 0 : matrix_type=dbcsr_type_antisymmetric)
1088 :
1089 0 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_kp_re(ispin, ikp)%matrix, sab_nl)
1090 0 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_kp_im(ispin, ikp)%matrix, sab_nl)
1091 :
1092 0 : CALL dbcsr_set(matrix_ks_kp_re(ispin, ikp)%matrix, 0.0_dp)
1093 0 : CALL dbcsr_set(matrix_ks_kp_im(ispin, ikp)%matrix, 0.0_dp)
1094 :
1095 : END DO
1096 : END DO
1097 :
1098 0 : CALL timestop(handle)
1099 :
1100 0 : END SUBROUTINE allocate_matrix_ks_kp
1101 :
1102 : END MODULE rpa_gw_sigma_x
|
