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