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