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