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