(git:6344619)
Loading...
Searching...
No Matches
post_scf_bandstructure_utils.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief
10!> \author Jan Wilhelm
11!> \date 07.2023
12! **************************************************************************************************
13MODULE post_scf_bandstructure_utils
14 USE atomic_kind_types, ONLY: atomic_kind_type,&
15 get_atomic_kind,&
16 get_atomic_kind_set
17 USE bse_input, ONLY: read_bse_section
18 USE bse_main, ONLY: bse_fulldiag_memory_check
19 USE bse_types, ONLY: bse_env_type
20 USE bse_util, ONLY: bse_window_from_dft
21 USE cell_types, ONLY: cell_type,&
22 get_cell
23 USE cp_blacs_env, ONLY: cp_blacs_env_type
24 USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale
25 USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose
26 USE cp_cfm_diag, ONLY: cp_cfm_geeig,&
27 cp_cfm_geeig_canon,&
28 cp_cfm_heevd
29 USE cp_cfm_types, ONLY: cp_cfm_create,&
30 cp_cfm_get_info,&
31 cp_cfm_release,&
32 cp_cfm_set_all,&
33 cp_cfm_to_cfm,&
34 cp_cfm_to_fm,&
35 cp_cfm_type,&
36 cp_fm_to_cfm
37 USE cp_control_types, ONLY: dft_control_type
38 USE cp_dbcsr_api, ONLY: &
39 dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_p_type, dbcsr_set, &
40 dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
41 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
42 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
43 copy_fm_to_dbcsr,&
44 dbcsr_allocate_matrix_set,&
45 dbcsr_deallocate_matrix_set
46 USE cp_files, ONLY: close_file,&
47 open_file
48 USE cp_fm_diag, ONLY: cp_fm_geeig_canon
49 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
50 cp_fm_struct_release,&
51 cp_fm_struct_type
52 USE cp_fm_types, ONLY: cp_fm_create,&
53 cp_fm_get_diag,&
54 cp_fm_get_info,&
55 cp_fm_release,&
56 cp_fm_set_all,&
57 cp_fm_to_fm,&
58 cp_fm_to_fm_submat,&
59 cp_fm_type
60 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit,&
61 cp_to_string
62 USE cp_parser_methods, ONLY: read_float_object
63 USE dft_plus_u, ONLY: lowdin_kp_apply,&
64 lowdin_kp_prepare
65 USE input_constants, ONLY: &
66 g0w0, bse_fulldiag, cohsex_off, evgw0, int_ldos_z, ri_rs_large_cell_gamma, &
67 ri_rs_non_periodic, smear_fermi_dirac, smear_gaussian, smear_mp, smear_mv, &
68 tensor_large_cell_gamma, tensor_small_cell_full_kp, tensor_static_cohsex
69 USE input_section_types, ONLY: section_vals_get,&
70 section_vals_get_subs_vals,&
71 section_vals_type,&
72 section_vals_val_get
73 USE kinds, ONLY: default_string_length,&
74 dp,&
75 max_line_length
76 USE kpoint_k_r_trafo_simple, ONLY: fm_gamma_to_kp
77 USE kpoint_methods, ONLY: kpoint_init_cell_index,&
78 rskp_transform
79 USE kpoint_types, ONLY: get_kpoint_info,&
80 kpoint_create,&
81 kpoint_env_type,&
82 kpoint_type
83 USE machine, ONLY: m_walltime
84 USE mathconstants, ONLY: gaussi,&
85 twopi,&
86 z_one,&
87 z_zero
88 USE message_passing, ONLY: mp_para_env_type
89 USE mo_window, ONLY: mo_window_type
90 USE parallel_gemm_api, ONLY: parallel_gemm
91 USE particle_types, ONLY: particle_type
92 USE physcon, ONLY: angstrom,&
93 evolt
94 USE post_scf_bandstructure_types, ONLY: band_edges_type,&
95 eps_qp_gap,&
96 max_qp_gap,&
97 post_scf_bandstructure_type
98 USE pw_env_types, ONLY: pw_env_get,&
99 pw_env_type
100 USE pw_pool_types, ONLY: pw_pool_type
101 USE pw_types, ONLY: pw_c1d_gs_type,&
102 pw_r3d_rs_type
103 USE qs_collocate_density, ONLY: calculate_rho_elec
104 USE qs_environment_types, ONLY: get_qs_env,&
105 qs_environment_type
106 USE qs_ks_types, ONLY: qs_ks_env_type
107 USE qs_mo_types, ONLY: get_mo_set,&
108 mo_set_type
109 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
110 USE rpa_gw_im_time_util, ONLY: compute_weight_re_im,&
111 get_atom_index_from_basis_function_index
112 USE scf_control_types, ONLY: scf_control_type
113 USE smearing_utils, ONLY: smearkp,&
114 smearkp2
115 USE soc_pseudopotential_methods, ONLY: v_soc_xyz_from_pseudopotential,&
116 remove_soc_outside_energy_window_mo
117 USE soc_pseudopotential_utils, ONLY: add_cfm_submat,&
118 add_dbcsr_submat,&
119 cfm_add_on_diag,&
120 create_cfm_double,&
121 get_cfm_submat
122 USE string_utilities, ONLY: uppercase
123#include "base/base_uses.f90"
124
125 IMPLICIT NONE
126
127 PRIVATE
128
129 PUBLIC :: create_and_init_bs_env, &
130 eval_bandstructure_properties, cfm_ikp_from_fm_gamma, &
131 mic_contribution_from_ikp, compute_xkp, kpoint_init_cell_index_simple, &
132 rsmat_to_kp, soc, get_vbm_cbm_bandgaps, get_all_vbm_cbm_bandgaps, &
133 allocate_gw_eigenvalues, gw_flavour_label
134
135 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'post_scf_bandstructure_utils'
136
137CONTAINS
138
139! **************************************************************************************************
140!> \brief Allocate the arrays holding the GW quasiparticle energies
141!> \param bs_env ...
142! **************************************************************************************************
143 SUBROUTINE allocate_gw_eigenvalues(bs_env)
144 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
145
146 CHARACTER(LEN=*), PARAMETER :: routinen = 'allocate_GW_eigenvalues'
147
148 INTEGER :: handle, n_ao, n_spin, nkp
149
150 CALL timeset(routinen, handle)
151
152 n_ao = bs_env%n_ao
153 nkp = bs_env%nkp_bs_and_DOS
154 n_spin = bs_env%n_spin
155
156 IF (.NOT. ALLOCATED(bs_env%eigenval_GW)) THEN
157 ALLOCATE (bs_env%eigenval_GW(n_ao, nkp, n_spin))
158 END IF
159 IF (.NOT. ALLOCATED(bs_env%eigenval_HF)) THEN
160 ALLOCATE (bs_env%eigenval_HF(n_ao, nkp, n_spin))
161 END IF
162
163 SELECT CASE (bs_env%gw_flavour)
164 CASE (g0w0, evgw0)
165 IF (.NOT. ALLOCATED(bs_env%eigenval_G0W0)) THEN
166 ALLOCATE (bs_env%eigenval_G0W0(n_ao, nkp, n_spin))
167 END IF
168 END SELECT
169
170 IF (bs_env%gw_flavour == evgw0 .AND. .NOT. ALLOCATED(bs_env%eigenval_evGW0)) THEN
171 ALLOCATE (bs_env%eigenval_evGW0(n_ao, nkp, n_spin))
172 END IF
173
174 CALL timestop(handle)
175
176 END SUBROUTINE allocate_gw_eigenvalues
177
178! **************************************************************************************************
179!> \brief Name of the GW flavour that was requested, for printing
180!> \param bs_env ...
181!> \return ...
182! **************************************************************************************************
183 FUNCTION gw_flavour_label(bs_env) RESULT(label)
184 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
185 CHARACTER(LEN=default_string_length) :: label
186
187 SELECT CASE (bs_env%gw_flavour)
188 CASE (evgw0)
189 label = "evGW0"
190 CASE DEFAULT
191 label = "G0W0"
192 END SELECT
193 IF (bs_env%do_static_cohsex) label = "static COHSEX"
194
195 END FUNCTION gw_flavour_label
196
197! **************************************************************************************************
198!> \brief ...
199!> \param qs_env ...
200!> \param bs_env ...
201!> \param post_scf_bandstructure_section ...
202! **************************************************************************************************
203 SUBROUTINE create_and_init_bs_env(qs_env, bs_env, post_scf_bandstructure_section)
204 TYPE(qs_environment_type), POINTER :: qs_env
205 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
206 TYPE(section_vals_type), POINTER :: post_scf_bandstructure_section
207
208 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_and_init_bs_env'
209
210 INTEGER :: handle
211
212 CALL timeset(routinen, handle)
213
214 ALLOCATE (bs_env)
215
216 CALL print_header(bs_env)
217
218 CALL read_bandstructure_input_parameters(bs_env, post_scf_bandstructure_section, qs_env)
219
220 CALL get_parameters_from_qs_env(qs_env, bs_env)
221
222 CALL set_heuristic_parameters(bs_env)
223
224 SELECT CASE (bs_env%gw_implementation)
225 CASE (tensor_large_cell_gamma, ri_rs_large_cell_gamma, ri_rs_non_periodic, tensor_static_cohsex)
226
227 CALL setup_kpoints_dos_large_cell_gamma(qs_env, bs_env, bs_env%kpoints_DOS)
228
229 CALL allocate_and_fill_fm_ks_fm_s(qs_env, bs_env)
230
231 CALL diagonalize_ks_matrix(bs_env)
232
233 CALL check_positive_definite_overlap_mat(bs_env, qs_env)
234
235 CASE (tensor_small_cell_full_kp)
236
237 CALL setup_kpoints_scf_desymm(qs_env, bs_env, bs_env%kpoints_scf_desymm, .true.)
238
239 CALL setup_kpoints_dos_small_cell_full_kp(bs_env, bs_env%kpoints_DOS)
240
241 CALL allocate_and_fill_fm_ks_fm_s(qs_env, bs_env)
242
243 CALL compute_cfm_mo_coeff_kp_and_eigenval_scf_kp(qs_env, bs_env)
244
245 END SELECT
246
247 CALL read_bse_input(qs_env, bs_env)
248
249 CALL timestop(handle)
250
251 END SUBROUTINE create_and_init_bs_env
252
253! **************************************************************************************************
254!> \brief What PROPERTIES%BSE needs the GW to keep and, with FULLDIAG, the early memory check of the
255!> full diagonalization; check_cp2k_input refused the input combinations the BSE cannot run
256!> \param qs_env source of the input
257!> \param bs_env SCF eigenvalues and occupations of the memory check; keep_W_freq_zero and keep_3c_nl are set
258! **************************************************************************************************
259 SUBROUTINE read_bse_input(qs_env, bs_env)
260 TYPE(qs_environment_type), POINTER :: qs_env
261 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
262
263 CHARACTER(LEN=*), PARAMETER :: routinen = 'read_bse_input'
264
265 INTEGER :: handle, n_spin
266 INTEGER, DIMENSION(2) :: n_mo
267 LOGICAL :: do_bse_props
268 TYPE(bse_env_type) :: bse_env
269 TYPE(mo_window_type) :: window
270 TYPE(section_vals_type), POINTER :: bse_sec, input
271
272 CALL timeset(routinen, handle)
273
274 NULLIFY (bse_sec, input)
275 CALL get_qs_env(qs_env, input=input)
276 bse_sec => section_vals_get_subs_vals(input, "PROPERTIES%BSE")
277 CALL section_vals_get(bse_sec, explicit=do_bse_props)
278
279 IF (do_bse_props) THEN
280 ! the full diagonalization's footprint is known here, before any GW work
281 CALL read_bse_section(bse_sec, bse_env)
282 IF (bse_env%bse_diag_method == bse_fulldiag) THEN
283 n_spin = bs_env%n_spin
284 ! the canonical orthogonalisation parks the linearly dependent directions above n_mo_retained
285 n_mo(:) = min(bs_env%n_mo_retained, bs_env%n_ao)
286 ! the BSE warns about a combined window when it starts; the estimate stays silent
287 CALL bse_window_from_dft(bs_env%eigenval_scf_Gamma(1:n_mo(1), 1:n_spin), n_mo(1:n_spin), &
288 bs_env%n_occ(1:n_spin), bse_env%bse_cutoff_occ, &
289 bse_env%bse_cutoff_empty, window, warn=.false.)
290 CALL bse_fulldiag_memory_check(bse_env, window, bs_env%n_occ(1:n_spin), bs_env%para_env, &
291 bs_env%unit_nr)
292 END IF
293 END IF
294 bs_env%keep_W_freq_zero = do_bse_props
295 bs_env%keep_3c_nl = do_bse_props
296
297 CALL timestop(handle)
298
299 END SUBROUTINE read_bse_input
300
301! **************************************************************************************************
302!> \brief ...
303!> \param bs_env ...
304!> \param bs_sec ...
305!> \param qs_env ...
306! **************************************************************************************************
307 SUBROUTINE read_bandstructure_input_parameters(bs_env, bs_sec, qs_env)
308 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
309 TYPE(section_vals_type), POINTER :: bs_sec
310 TYPE(qs_environment_type), POINTER :: qs_env
311
312 CHARACTER(LEN=*), PARAMETER :: routinen = 'read_bandstructure_input_parameters'
313
314 CHARACTER(LEN=default_string_length) :: ustr
315 CHARACTER(LEN=default_string_length), &
316 DIMENSION(:), POINTER :: string_ptr
317 CHARACTER(LEN=max_line_length) :: error_msg
318 INTEGER :: handle, i, ikp
319 REAL(kind=dp), DIMENSION(3) :: kpptr
320 REAL(kind=dp), DIMENSION(3, 3) :: cart_hmat
321 TYPE(cell_type), POINTER :: cell
322 TYPE(section_vals_type), POINTER :: dos_pdos_sec, floquet_sec, gw_sec, &
323 kp_bs_sec, ldos_sec, soc_sec
324
325 CALL timeset(routinen, handle)
326 NULLIFY (cell)
327 CALL get_qs_env(qs_env=qs_env, cell=cell)
328 cart_hmat(:, :) = cell%hmat(:, :)
329 IF (cell%input_cell_canonicalized) cart_hmat(:, :) = cell%input_hmat(:, :)
330
331 NULLIFY (gw_sec)
332 gw_sec => section_vals_get_subs_vals(bs_sec, "GW")
333 CALL section_vals_get(gw_sec, explicit=bs_env%do_gw)
334 CALL section_vals_val_get(gw_sec, "COHSEX", i_val=bs_env%cohsex_flavor)
335 bs_env%do_static_cohsex = (bs_env%cohsex_flavor /= cohsex_off)
336 CALL section_vals_val_get(gw_sec, "RI_RS%_SECTION_PARAMETERS_", l_val=bs_env%do_gw_ri_rs)
337
338 NULLIFY (soc_sec)
339 soc_sec => section_vals_get_subs_vals(bs_sec, "SOC")
340 CALL section_vals_get(soc_sec, explicit=bs_env%do_soc)
341
342 CALL section_vals_val_get(soc_sec, "SOC_WINDOW_OCC", r_val=bs_env%soc_window_occ)
343 CALL section_vals_val_get(soc_sec, "SOC_WINDOW_VIRT", r_val=bs_env%soc_window_virt)
344 CALL section_vals_val_get(soc_sec, "SOC_WINDOW_SMEARING", r_val=bs_env%soc_window_smearing)
345
346 CALL section_vals_val_get(soc_sec, "SOC_NBANDS_OCC", i_val=bs_env%soc_nbands_occ)
347 CALL section_vals_val_get(soc_sec, "SOC_NBANDS_VIRT", i_val=bs_env%soc_nbands_virt)
348
349 NULLIFY (dos_pdos_sec)
350 dos_pdos_sec => section_vals_get_subs_vals(bs_sec, "DOS")
351 CALL section_vals_get(dos_pdos_sec, explicit=bs_env%do_dos_pdos)
352
353 CALL section_vals_val_get(bs_sec, "DOS%KPOINTS", i_vals=bs_env%nkp_grid_DOS_input)
354 CALL section_vals_val_get(bs_sec, "DOS%ENERGY_WINDOW", r_val=bs_env%energy_window_DOS)
355 CALL section_vals_val_get(bs_sec, "DOS%ENERGY_STEP", r_val=bs_env%energy_step_DOS)
356 CALL section_vals_val_get(bs_sec, "DOS%BROADENING", r_val=bs_env%broadening_DOS)
357
358 NULLIFY (ldos_sec)
359 ldos_sec => section_vals_get_subs_vals(bs_sec, "DOS%LDOS")
360 CALL section_vals_get(ldos_sec, explicit=bs_env%do_ldos)
361
362 CALL section_vals_val_get(ldos_sec, "INTEGRATION", i_val=bs_env%int_ldos_xyz)
363 CALL section_vals_val_get(ldos_sec, "BIN_MESH", i_vals=bs_env%bin_mesh)
364
365 NULLIFY (kp_bs_sec)
366 kp_bs_sec => section_vals_get_subs_vals(bs_sec, "BANDSTRUCTURE_PATH")
367 CALL section_vals_val_get(kp_bs_sec, "NPOINTS", i_val=bs_env%input_kp_bs_npoints)
368 CALL section_vals_val_get(kp_bs_sec, "UNITS", c_val=ustr)
369 CALL uppercase(ustr)
370 CALL section_vals_val_get(kp_bs_sec, "SPECIAL_POINT", n_rep_val=bs_env%input_kp_bs_n_sp_pts)
371
372 NULLIFY (floquet_sec)
373 floquet_sec => section_vals_get_subs_vals(bs_sec, "FLOQUET")
374 CALL section_vals_get(floquet_sec, explicit=bs_env%do_floquet)
375 CALL section_vals_val_get(floquet_sec, "AMPLITUDE", r_val=bs_env%floquet_amplitude)
376 CALL section_vals_val_get(floquet_sec, "FREQUENCY", r_val=bs_env%floquet_omega)
377 CALL section_vals_val_get(floquet_sec, "POLARISATION", r_vals=bs_env%floquet_polarisation)
378 CALL section_vals_val_get(floquet_sec, "PHASE_OFFSETS", r_vals=bs_env%floquet_phi)
379 CALL section_vals_val_get(floquet_sec, "MAX_FLOQUET_INDEX", i_val=bs_env%max_floquet_index)
380 CALL section_vals_val_get(floquet_sec, "MAX_NMO", i_val=bs_env%max_nmo_floquet)
381 CALL section_vals_val_get(floquet_sec, "EPS_FLOQUET", r_val=bs_env%eps_floquet)
382 CALL section_vals_val_get(floquet_sec, "ENERGY_WINDOW", r_val=bs_env%energy_window_floquet)
383 CALL section_vals_val_get(floquet_sec, "ENERGY_STEP", r_val=bs_env%energy_step_floquet)
384 CALL section_vals_val_get(floquet_sec, "BROADENING", r_val=bs_env%broadening_floquet)
385 CALL section_vals_val_get(floquet_sec, "MEM_FILL_FRACTION", r_val=bs_env%floquet_mem_fill_fraction)
386 CALL section_vals_val_get(floquet_sec, "TEMPERATURE", r_val=bs_env%floquet_temperature)
387 CALL section_vals_val_get(floquet_sec, "FLOQUET_DOS_FILE_NAME", c_val=bs_env%floquet_dos_file)
388 CALL section_vals_val_get(floquet_sec, "QUASI_ENERGIES_FILE_NAME", c_val=bs_env%floquet_qe_file)
389 CALL section_vals_val_get(floquet_sec, "FLOQUET_BS_FILE_NAME", c_val=bs_env%floquet_bs_file)
390
391 ! read special points for band structure
392 ALLOCATE (bs_env%xkp_special(3, bs_env%input_kp_bs_n_sp_pts))
393 DO ikp = 1, bs_env%input_kp_bs_n_sp_pts
394 CALL section_vals_val_get(kp_bs_sec, "SPECIAL_POINT", i_rep_val=ikp, c_vals=string_ptr)
395 cpassert(SIZE(string_ptr(:), 1) == 4)
396 DO i = 1, 3
397 CALL read_float_object(string_ptr(i + 1), kpptr(i), error_msg)
398 IF (len_trim(error_msg) > 0) cpabort(trim(error_msg))
399 END DO
400 SELECT CASE (ustr)
401 CASE ("B_VECTOR")
402 bs_env%xkp_special(1:3, ikp) = kpptr(1:3)
403 CASE ("CART_ANGSTROM")
404 bs_env%xkp_special(1:3, ikp) = (kpptr(1)*cart_hmat(1, 1:3) + &
405 kpptr(2)*cart_hmat(2, 1:3) + &
406 kpptr(3)*cart_hmat(3, 1:3))/twopi*angstrom
407 CASE ("CART_BOHR")
408 bs_env%xkp_special(1:3, ikp) = (kpptr(1)*cart_hmat(1, 1:3) + &
409 kpptr(2)*cart_hmat(2, 1:3) + &
410 kpptr(3)*cart_hmat(3, 1:3))/twopi
411 CASE DEFAULT
412 cpabort("Unknown unit <"//trim(ustr)//"> specified for k-point definition")
413 END SELECT
414 END DO
415
416 CALL timestop(handle)
417
418 END SUBROUTINE read_bandstructure_input_parameters
419
420! **************************************************************************************************
421!> \brief ...
422!> \param bs_env ...
423! **************************************************************************************************
424 SUBROUTINE print_header(bs_env)
425
426 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
427
428 CHARACTER(LEN=*), PARAMETER :: routinen = 'print_header'
429
430 INTEGER :: handle, u
431
432 CALL timeset(routinen, handle)
433
434 bs_env%unit_nr = cp_logger_get_default_io_unit()
435
436 u = bs_env%unit_nr
437
438 IF (u > 0) THEN
439 WRITE (u, '(T2,A)') ' '
440 WRITE (u, '(T2,A)') repeat('-', 79)
441 WRITE (u, '(T2,A,A78)') '-', '-'
442 WRITE (u, '(T2,A,A51,A27)') '-', 'BANDSTRUCTURE CALCULATION', '-'
443 WRITE (u, '(T2,A,A78)') '-', '-'
444 WRITE (u, '(T2,A)') repeat('-', 79)
445 WRITE (u, '(T2,A)') ' '
446 END IF
447
448 CALL timestop(handle)
449
450 END SUBROUTINE print_header
451
452! **************************************************************************************************
453!> \brief ...
454!> \param qs_env ...
455!> \param bs_env ...
456!> \param kpoints ...
457! **************************************************************************************************
458 SUBROUTINE setup_kpoints_dos_large_cell_gamma(qs_env, bs_env, kpoints)
459
460 TYPE(qs_environment_type), POINTER :: qs_env
461 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
462 TYPE(kpoint_type), POINTER :: kpoints
463
464 CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_kpoints_DOS_large_cell_Gamma'
465
466 INTEGER :: handle, i_dim, i_kp_in_line, &
467 i_special_kp, ikk, n_kp_in_line, &
468 n_special_kp, nkp, nkp_only_bs, &
469 nkp_only_dos, u
470 INTEGER, DIMENSION(3) :: nkp_grid, periodic
471
472 CALL timeset(routinen, handle)
473
474 ! routine adapted from mp2_integrals.F
475 NULLIFY (kpoints)
476 CALL kpoint_create(kpoints)
477
478 kpoints%kp_scheme = "GENERAL"
479
480 n_special_kp = bs_env%input_kp_bs_n_sp_pts
481 n_kp_in_line = bs_env%input_kp_bs_npoints
482
483 periodic(1:3) = bs_env%periodic(1:3)
484
485 DO i_dim = 1, 3
486
487 cpassert(periodic(i_dim) == 0 .OR. periodic(i_dim) == 1)
488
489 IF (bs_env%nkp_grid_DOS_input(i_dim) < 0) THEN
490 IF (periodic(i_dim) == 1) nkp_grid(i_dim) = 2
491 IF (periodic(i_dim) == 0) nkp_grid(i_dim) = 1
492 ELSE
493 nkp_grid(i_dim) = bs_env%nkp_grid_DOS_input(i_dim)
494 END IF
495
496 END DO
497
498 ! use the k <-> -k symmetry to reduce the number of kpoints
499 IF (nkp_grid(1) > 1) THEN
500 nkp_only_dos = (nkp_grid(1) + 1)/2*nkp_grid(2)*nkp_grid(3)
501 ELSE IF (nkp_grid(2) > 1) THEN
502 nkp_only_dos = nkp_grid(1)*(nkp_grid(2) + 1)/2*nkp_grid(3)
503 ELSE IF (nkp_grid(3) > 1) THEN
504 nkp_only_dos = nkp_grid(1)*nkp_grid(2)*(nkp_grid(3) + 1)/2
505 ELSE
506 nkp_only_dos = 1
507 END IF
508
509 ! we will compute the GW QP levels for all k's in the bandstructure path but also
510 ! for all k-points from the SCF (e.g. for DOS or for self-consistent GW)
511 IF (n_special_kp > 0) THEN
512 nkp_only_bs = n_kp_in_line*(n_special_kp - 1) + 1
513 ELSE
514 nkp_only_bs = 0
515 END IF
516
517 nkp = nkp_only_dos + nkp_only_bs
518
519 kpoints%nkp_grid(1:3) = nkp_grid(1:3)
520 kpoints%nkp = nkp
521
522 bs_env%nkp_bs_and_DOS = nkp
523 bs_env%nkp_only_bs = nkp_only_bs
524 bs_env%nkp_only_DOS = nkp_only_dos
525
526 ALLOCATE (kpoints%xkp(3, nkp), kpoints%wkp(nkp))
527 kpoints%wkp(:) = 0.0_dp
528 kpoints%wkp(1:nkp_only_dos) = 1.0_dp/real(nkp_only_dos, kind=dp)
529
530 CALL compute_xkp(kpoints%xkp, 1, nkp_only_dos, nkp_grid)
531
532 IF (n_special_kp > 0) THEN
533 kpoints%xkp(1:3, nkp_only_dos + 1) = bs_env%xkp_special(1:3, 1)
534 ikk = nkp_only_dos + 1
535 DO i_special_kp = 2, n_special_kp
536 DO i_kp_in_line = 1, n_kp_in_line
537 ikk = ikk + 1
538 kpoints%xkp(1:3, ikk) = bs_env%xkp_special(1:3, i_special_kp - 1) + &
539 REAL(i_kp_in_line, kind=dp)/real(n_kp_in_line, kind=dp)* &
540 (bs_env%xkp_special(1:3, i_special_kp) - &
541 bs_env%xkp_special(1:3, i_special_kp - 1))
542 kpoints%wkp(ikk) = 0.0_dp
543 END DO
544 END DO
545 END IF
546
547 CALL kpoint_init_cell_index_simple(kpoints, qs_env)
548
549 u = bs_env%unit_nr
550
551 IF (u > 0 .AND. bs_env%do_periodic) THEN
552 IF (nkp_only_bs > 0) THEN
553 WRITE (u, fmt="(T2,1A,T77,I4)") &
554 "Number of special k-points for the bandstructure", n_special_kp
555 WRITE (u, fmt="(T2,1A,T77,I4)") "Number of k-points for the bandstructure", nkp
556 WRITE (u, fmt="(T2,1A,T69,3I4)") &
557 "K-point mesh for the density of states (DOS)", nkp_grid(1:3)
558 ELSE
559 WRITE (u, fmt="(T2,1A,T69,3I4)") &
560 "K-point mesh for the density of states (DOS) and the self-energy", nkp_grid(1:3)
561 END IF
562 END IF
563
564 CALL timestop(handle)
565
566 END SUBROUTINE setup_kpoints_dos_large_cell_gamma
567
568! **************************************************************************************************
569!> \brief ...
570!> \param qs_env ...
571!> \param bs_env ...
572!> \param kpoints ...
573!> \param do_print ...
574! **************************************************************************************************
575 SUBROUTINE setup_kpoints_scf_desymm(qs_env, bs_env, kpoints, do_print)
576 TYPE(qs_environment_type), POINTER :: qs_env
577 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
578 TYPE(kpoint_type), POINTER :: kpoints
579
580 CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_kpoints_scf_desymm'
581
582 INTEGER :: handle, i_cell_x, i_dim, img, j_cell_y, &
583 k_cell_z, nimages, nkp, u
584 INTEGER, DIMENSION(3) :: cell_grid, cixd, nkp_grid
585 TYPE(kpoint_type), POINTER :: kpoints_scf
586
587 LOGICAL:: do_print
588
589 CALL timeset(routinen, handle)
590
591 NULLIFY (kpoints)
592 CALL kpoint_create(kpoints)
593
594 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints_scf)
595
596 nkp_grid(1:3) = kpoints_scf%nkp_grid(1:3)
597 nkp = nkp_grid(1)*nkp_grid(2)*nkp_grid(3)
598
599 ! we need in periodic directions at least 4 k-points in the SCF
600 DO i_dim = 1, 3
601 IF (bs_env%periodic(i_dim) == 1) THEN
602 cpassert(nkp_grid(i_dim) >= 4)
603 END IF
604 END DO
605
606 kpoints%kp_scheme = "GENERAL"
607 kpoints%nkp_grid(1:3) = nkp_grid(1:3)
608 kpoints%nkp = nkp
609 bs_env%nkp_scf_desymm = nkp
610
611 ALLOCATE (kpoints%xkp(1:3, nkp))
612 CALL compute_xkp(kpoints%xkp, 1, nkp, nkp_grid)
613
614 ALLOCATE (kpoints%wkp(nkp))
615 kpoints%wkp(:) = 1.0_dp/real(nkp, kind=dp)
616
617 ! for example 4x3x6 kpoint grid -> 3x3x5 cell grid because we need the same number of
618 ! neighbor cells on both sides of the unit cell
619 cell_grid(1:3) = nkp_grid(1:3) - modulo(nkp_grid(1:3) + 1, 2)
620
621 ! cell index: for example for x: from -n_x/2 to +n_x/2, n_x: number of cells in x direction
622 cixd(1:3) = cell_grid(1:3)/2
623
624 nimages = cell_grid(1)*cell_grid(2)*cell_grid(3)
625
626 bs_env%nimages_scf_desymm = nimages
627 bs_env%cell_grid_scf_desymm(1:3) = cell_grid(1:3)
628
629 IF (ASSOCIATED(kpoints%index_to_cell)) DEALLOCATE (kpoints%index_to_cell)
630 IF (ASSOCIATED(kpoints%cell_to_index)) DEALLOCATE (kpoints%cell_to_index)
631
632 ALLOCATE (kpoints%cell_to_index(-cixd(1):cixd(1), -cixd(2):cixd(2), -cixd(3):cixd(3)))
633 ALLOCATE (kpoints%index_to_cell(3, nimages))
634
635 img = 0
636 DO i_cell_x = -cixd(1), cixd(1)
637 DO j_cell_y = -cixd(2), cixd(2)
638 DO k_cell_z = -cixd(3), cixd(3)
639 img = img + 1
640 kpoints%cell_to_index(i_cell_x, j_cell_y, k_cell_z) = img
641 kpoints%index_to_cell(1:3, img) = [i_cell_x, j_cell_y, k_cell_z]
642 END DO
643 END DO
644 END DO
645
646 u = bs_env%unit_nr
647 IF (u > 0 .AND. do_print) THEN
648 WRITE (u, fmt="(T2,A,I49)") "Number of cells for G, χ, W, Σ", nimages
649 END IF
650
651 CALL timestop(handle)
652
653 END SUBROUTINE setup_kpoints_scf_desymm
654
655! **************************************************************************************************
656!> \brief ...
657!> \param bs_env ...
658!> \param kpoints ...
659! **************************************************************************************************
660 SUBROUTINE setup_kpoints_dos_small_cell_full_kp(bs_env, kpoints)
661
662 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
663 TYPE(kpoint_type), POINTER :: kpoints
664
665 CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_kpoints_DOS_small_cell_full_kp'
666
667 INTEGER :: handle, i_kp_in_line, i_special_kp, ikk, &
668 n_kp_in_line, n_special_kp, nkp, &
669 nkp_only_bs, nkp_scf_desymm, u
670
671 CALL timeset(routinen, handle)
672
673 ! routine adapted from mp2_integrals.F
674 NULLIFY (kpoints)
675 CALL kpoint_create(kpoints)
676
677 n_special_kp = bs_env%input_kp_bs_n_sp_pts
678 n_kp_in_line = bs_env%input_kp_bs_npoints
679 nkp_scf_desymm = bs_env%nkp_scf_desymm
680
681 ! we will compute the GW QP levels for all k's in the bandstructure path but also
682 ! for all k-points from the SCF (e.g. for DOS or for self-consistent GW)
683 IF (n_special_kp > 0) THEN
684 nkp_only_bs = n_kp_in_line*(n_special_kp - 1) + 1
685 ELSE
686 nkp_only_bs = 0
687 END IF
688 nkp = nkp_only_bs + nkp_scf_desymm
689
690 ALLOCATE (kpoints%xkp(3, nkp))
691 ALLOCATE (kpoints%wkp(nkp))
692
693 kpoints%nkp = nkp
694
695 bs_env%nkp_bs_and_DOS = nkp
696 bs_env%nkp_only_bs = nkp_only_bs
697 bs_env%nkp_only_DOS = nkp_scf_desymm
698
699 kpoints%xkp(1:3, 1:nkp_scf_desymm) = bs_env%kpoints_scf_desymm%xkp(1:3, 1:nkp_scf_desymm)
700 kpoints%wkp(1:nkp_scf_desymm) = 1.0_dp/real(nkp_scf_desymm, kind=dp)
701
702 IF (n_special_kp > 0) THEN
703 kpoints%xkp(1:3, nkp_scf_desymm + 1) = bs_env%xkp_special(1:3, 1)
704 ikk = nkp_scf_desymm + 1
705 DO i_special_kp = 2, n_special_kp
706 DO i_kp_in_line = 1, n_kp_in_line
707 ikk = ikk + 1
708 kpoints%xkp(1:3, ikk) = bs_env%xkp_special(1:3, i_special_kp - 1) + &
709 REAL(i_kp_in_line, kind=dp)/real(n_kp_in_line, kind=dp)* &
710 (bs_env%xkp_special(1:3, i_special_kp) - &
711 bs_env%xkp_special(1:3, i_special_kp - 1))
712 kpoints%wkp(ikk) = 0.0_dp
713 END DO
714 END DO
715 END IF
716
717 IF (ASSOCIATED(kpoints%index_to_cell)) DEALLOCATE (kpoints%index_to_cell)
718
719 ALLOCATE (kpoints%index_to_cell(3, bs_env%nimages_scf_desymm))
720 kpoints%index_to_cell(:, :) = bs_env%kpoints_scf_desymm%index_to_cell(:, :)
721
722 u = bs_env%unit_nr
723
724 IF (u > 0) THEN
725 WRITE (u, fmt="(T2,1A,T77,I4)") "Number of special k-points for the bandstructure", &
726 n_special_kp
727 WRITE (u, fmt="(T2,1A,T77,I4)") "Number of k-points for the bandstructure", nkp
728 END IF
729
730 CALL timestop(handle)
731
732 END SUBROUTINE setup_kpoints_dos_small_cell_full_kp
733
734! **************************************************************************************************
735!> \brief ...
736!> \param qs_env ...
737!> \param bs_env ...
738! **************************************************************************************************
739 SUBROUTINE compute_cfm_mo_coeff_kp_and_eigenval_scf_kp(qs_env, bs_env)
740 TYPE(qs_environment_type), POINTER :: qs_env
741 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
742
743 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_cfm_mo_coeff_kp_and_eigenval_scf_kp'
744
745 INTEGER :: handle, ikp, ispin, nkp_bs_and_dos, &
746 nmo_retained
747 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_scf
748 REAL(kind=dp) :: cbm, vbm
749 REAL(kind=dp), DIMENSION(3) :: xkp
750 TYPE(cp_cfm_type) :: cfm_ks, cfm_mos, cfm_s
751 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
752 TYPE(kpoint_env_type), POINTER :: lowdin_kp
753 TYPE(kpoint_type), POINTER :: kpoints_scf
754 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
755 POINTER :: sab_nl
756
757 CALL timeset(routinen, handle)
758
759 CALL get_qs_env(qs_env, &
760 matrix_ks_kp=matrix_ks, &
761 matrix_s_kp=matrix_s, &
762 kpoints=kpoints_scf)
763
764 NULLIFY (sab_nl, lowdin_kp)
765 IF (ALLOCATED(kpoints_scf%lowdin_v)) THEN
766 IF (any(kpoints_scf%lowdin_v /= 0.0_dp)) ALLOCATE (lowdin_kp)
767 END IF
768 CALL get_kpoint_info(kpoints_scf, sab_nl=sab_nl, cell_to_index=cell_to_index_scf)
769
770 CALL cp_cfm_create(cfm_ks, bs_env%cfm_work_mo%matrix_struct)
771 CALL cp_cfm_create(cfm_s, bs_env%cfm_work_mo%matrix_struct)
772 CALL cp_cfm_create(cfm_mos, bs_env%cfm_work_mo%matrix_struct)
773
774 ! nkp_bs_and_DOS contains desymmetrized k-point mesh from SCF and k-points from GW bandstructure
775 nkp_bs_and_dos = bs_env%nkp_bs_and_DOS
776
777 CALL allocate_gw_eigenvalues(bs_env)
778
779 ALLOCATE (bs_env%cfm_mo_coeff_kp(nkp_bs_and_dos, bs_env%n_spin))
780 ALLOCATE (bs_env%cfm_ks_kp(nkp_bs_and_dos, bs_env%n_spin))
781 ALLOCATE (bs_env%cfm_s_kp(nkp_bs_and_dos))
782 DO ikp = 1, nkp_bs_and_dos
783 DO ispin = 1, bs_env%n_spin
784 CALL cp_cfm_create(bs_env%cfm_mo_coeff_kp(ikp, ispin), bs_env%cfm_work_mo%matrix_struct)
785 CALL cp_cfm_create(bs_env%cfm_ks_kp(ikp, ispin), bs_env%cfm_work_mo%matrix_struct)
786 END DO
787 CALL cp_cfm_create(bs_env%cfm_s_kp(ikp), bs_env%cfm_work_mo%matrix_struct)
788 END DO
789
790 DO ikp = 1, nkp_bs_and_dos
791
792 xkp(1:3) = bs_env%kpoints_DOS%xkp(1:3, ikp)
793
794 ! S(k) and its Lowdin projectors are shared by both spin channels.
795 CALL rsmat_to_kp(matrix_s, 1, xkp, cell_to_index_scf, sab_nl, bs_env, cfm_s)
796 CALL cp_cfm_to_cfm(cfm_s, bs_env%cfm_s_kp(ikp))
797 IF (ASSOCIATED(lowdin_kp)) THEN
798 CALL lowdin_kp_prepare(lowdin_kp, kpoints_scf%lowdin_orbitals, .false., cfm_s)
799 END IF
800 DO ispin = 1, bs_env%n_spin
801 ! h^KS^R -> h^KS(k)
802 CALL rsmat_to_kp(matrix_ks, ispin, xkp, cell_to_index_scf, sab_nl, bs_env, cfm_ks)
803 IF (ASSOCIATED(lowdin_kp)) THEN
804 CALL lowdin_kp_apply(lowdin_kp, kpoints_scf%lowdin_v(:, ispin), cmat=cfm_ks)
805 END IF
806
807 ! we store the complex KS matrix as fm matrix because the infrastructure for fm is
808 ! much nicer compared to cfm
809 CALL cp_cfm_to_cfm(cfm_ks, bs_env%cfm_ks_kp(ikp, ispin))
810 CALL cp_cfm_to_cfm(bs_env%cfm_s_kp(ikp), cfm_s)
811
812 ! Diagonalize KS-matrix via Rothaan-Hall equation:
813 ! H^KS(k) C(k) = S(k) C(k) ε(k)
814 CALL cp_cfm_geeig_canon(cfm_ks, cfm_s, cfm_mos, &
815 bs_env%eigenval_scf(:, ikp, ispin), &
816 bs_env%cfm_work_mo, bs_env%eps_eigval_mat_s, &
817 nmo_retained=nmo_retained)
818 bs_env%n_mo_retained = min(bs_env%n_mo_retained, nmo_retained)
819
820 ! we store the complex MO coeff as fm matrix because the infrastructure for fm is
821 ! much nicer compared to cfm
822 CALL cp_cfm_to_cfm(cfm_mos, bs_env%cfm_mo_coeff_kp(ikp, ispin))
823
824 END DO
825 END DO
826
827 DO ispin = 1, bs_env%n_spin
828 vbm = maxval(bs_env%eigenval_scf(bs_env%n_occ(ispin), :, ispin))
829 cbm = minval(bs_env%eigenval_scf(bs_env%n_occ(ispin) + 1, :, ispin))
830
831 bs_env%e_fermi(ispin) = 0.5_dp*(vbm + cbm)
832
833 END DO
834
835 CALL get_vbm_cbm_bandgaps(bs_env%band_edges_scf, bs_env%eigenval_scf, bs_env)
836
837 CALL cp_cfm_release(cfm_ks)
838 CALL cp_cfm_release(cfm_s)
839 CALL cp_cfm_release(cfm_mos)
840 IF (ASSOCIATED(lowdin_kp)) THEN
841 CALL cp_cfm_release(lowdin_kp%lowdin_cb)
842 CALL cp_cfm_release(lowdin_kp%lowdin_cwork)
843 CALL lowdin_kp%lowdin_gemm%destroy()
844 DEALLOCATE (lowdin_kp)
845 END IF
846
847 CALL timestop(handle)
848
849 END SUBROUTINE compute_cfm_mo_coeff_kp_and_eigenval_scf_kp
850
851! **************************************************************************************************
852!> \brief ...
853!> \param mat_rs ...
854!> \param ispin ...
855!> \param xkp ...
856!> \param cell_to_index_scf ...
857!> \param sab_nl ...
858!> \param bs_env ...
859!> \param cfm_kp ...
860!> \param imag_rs_mat ...
861! **************************************************************************************************
862 SUBROUTINE rsmat_to_kp(mat_rs, ispin, xkp, cell_to_index_scf, sab_nl, bs_env, cfm_kp, imag_rs_mat)
863 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_rs
864 INTEGER :: ispin
865 REAL(kind=dp), DIMENSION(3) :: xkp
866 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_scf
867 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
868 POINTER :: sab_nl
869 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
870 TYPE(cp_cfm_type) :: cfm_kp
871 LOGICAL, OPTIONAL :: imag_rs_mat
872
873 CHARACTER(LEN=*), PARAMETER :: routinen = 'rsmat_to_kp'
874
875 INTEGER :: handle
876 LOGICAL :: imag_rs_mat_private
877 TYPE(dbcsr_type), POINTER :: cmat, nsmat, rmat
878
879 CALL timeset(routinen, handle)
880
881 ALLOCATE (rmat, cmat, nsmat)
882
883 imag_rs_mat_private = .false.
884 IF (PRESENT(imag_rs_mat)) imag_rs_mat_private = imag_rs_mat
885
886 IF (imag_rs_mat_private) THEN
887 CALL dbcsr_create(rmat, template=mat_rs(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
888 CALL dbcsr_create(cmat, template=mat_rs(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
889 ELSE
890 CALL dbcsr_create(rmat, template=mat_rs(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
891 CALL dbcsr_create(cmat, template=mat_rs(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
892 END IF
893 CALL dbcsr_create(nsmat, template=mat_rs(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
894 CALL cp_dbcsr_alloc_block_from_nbl(rmat, sab_nl)
895 CALL cp_dbcsr_alloc_block_from_nbl(cmat, sab_nl)
896
897 CALL dbcsr_set(rmat, 0.0_dp)
898 CALL dbcsr_set(cmat, 0.0_dp)
899 CALL rskp_transform(rmatrix=rmat, cmatrix=cmat, rsmat=mat_rs, ispin=ispin, &
900 xkp=xkp, cell_to_index=cell_to_index_scf, sab_nl=sab_nl)
901
902 CALL dbcsr_desymmetrize(rmat, nsmat)
903 CALL copy_dbcsr_to_fm(nsmat, bs_env%fm_work_mo(1))
904 CALL dbcsr_desymmetrize(cmat, nsmat)
905 CALL copy_dbcsr_to_fm(nsmat, bs_env%fm_work_mo(2))
906 CALL cp_fm_to_cfm(bs_env%fm_work_mo(1), bs_env%fm_work_mo(2), cfm_kp)
907
908 CALL dbcsr_deallocate_matrix(rmat)
909 CALL dbcsr_deallocate_matrix(cmat)
910 CALL dbcsr_deallocate_matrix(nsmat)
911
912 CALL timestop(handle)
913
914 END SUBROUTINE rsmat_to_kp
915
916! **************************************************************************************************
917!> \brief ...
918!> \param bs_env ...
919! **************************************************************************************************
920 SUBROUTINE diagonalize_ks_matrix(bs_env)
921 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
922
923 CHARACTER(LEN=*), PARAMETER :: routinen = 'diagonalize_ks_matrix'
924
925 INTEGER :: handle, ikp, ispin, nmo_retained
926 REAL(kind=dp) :: cbm, vbm
927
928 CALL timeset(routinen, handle)
929
930 ALLOCATE (bs_env%eigenval_scf_Gamma(bs_env%n_ao, bs_env%n_spin))
931
932 DO ispin = 1, bs_env%n_spin
933
934 ! use work matrices because the matrices are overwritten in cp_fm_geeig_canon
935 CALL cp_fm_to_fm(bs_env%fm_ks_Gamma(ispin), bs_env%fm_work_mo(1))
936 CALL cp_fm_to_fm(bs_env%fm_s_Gamma, bs_env%fm_work_mo(2))
937
938 ! diagonalize the Kohn-Sham matrix to obtain MO coefficients and SCF eigenvalues
939 ! (at the Gamma-point)
940 CALL cp_fm_geeig_canon(bs_env%fm_work_mo(1), &
941 bs_env%fm_work_mo(2), &
942 bs_env%fm_mo_coeff_Gamma(ispin), &
943 bs_env%eigenval_scf_Gamma(:, ispin), &
944 bs_env%fm_work_mo(3), &
945 bs_env%eps_eigval_mat_s, &
946 nmo_retained=nmo_retained)
947 bs_env%n_mo_retained = min(bs_env%n_mo_retained, nmo_retained)
948
949 vbm = bs_env%eigenval_scf_Gamma(bs_env%n_occ(ispin), ispin)
950 cbm = bs_env%eigenval_scf_Gamma(bs_env%n_occ(ispin) + 1, ispin)
951
952 bs_env%e_fermi(ispin) = 0.5_dp*(vbm + cbm)
953
954 END DO
955
956 CALL timestop(handle)
957
958 ! Gamma-only path for molecules: eigenval_scf is filled here from the Gamma eigenvalues
959 DO ispin = 1, bs_env%n_spin
960 DO ikp = 1, bs_env%nkp_bs_and_DOS
961 bs_env%eigenval_scf(:, ikp, ispin) = bs_env%eigenval_scf_Gamma(:, ispin)
962 END DO
963 END DO
964
965 END SUBROUTINE diagonalize_ks_matrix
966
967! **************************************************************************************************
968!> \brief ...
969!> \param bs_env ...
970!> \param qs_env ...
971! **************************************************************************************************
972 SUBROUTINE check_positive_definite_overlap_mat(bs_env, qs_env)
973 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
974 TYPE(qs_environment_type), POINTER :: qs_env
975
976 CHARACTER(LEN=*), PARAMETER :: routinen = 'check_positive_definite_overlap_mat'
977
978 INTEGER :: handle, ikp, info, u
979 TYPE(cp_cfm_type) :: cfm_s_ikp
980
981 CALL timeset(routinen, handle)
982
983 DO ikp = 1, bs_env%kpoints_DOS%nkp
984
985 ! get S_µν(k_i) from S_µν(k=0)
986 CALL cfm_ikp_from_fm_gamma(cfm_s_ikp, bs_env%fm_s_Gamma, &
987 ikp, qs_env, bs_env%kpoints_DOS, "ORB")
988
989 ! check whether S_µν(k_i) is positive definite
990 CALL cp_cfm_cholesky_decompose(matrix=cfm_s_ikp, n=bs_env%n_ao, info_out=info)
991
992 ! check if Cholesky decomposition failed (Cholesky decomposition only works for
993 ! positive definite matrices
994 IF (info /= 0) THEN
995 u = bs_env%unit_nr
996
997 IF (u > 0) THEN
998 WRITE (u, fmt="(T2,A)") ""
999 WRITE (u, fmt="(T2,A)") "ERROR: The Cholesky decomposition "// &
1000 "of the k-point overlap matrix failed. This is"
1001 WRITE (u, fmt="(T2,A)") "because the algorithm is "// &
1002 "only correct in the limit of large cells. The cell of "
1003 WRITE (u, fmt="(T2,A)") "the calculation is too small. "// &
1004 "Use MULTIPLE_UNIT_CELL to create a larger cell "
1005 WRITE (u, fmt="(T2,A)") "and to prevent this error."
1006 END IF
1007
1008 CALL bs_env%para_env%sync()
1009 cpabort("Please see information on the error above.")
1010
1011 END IF ! Cholesky decomposition failed
1012
1013 END DO ! ikp
1014
1015 CALL cp_cfm_release(cfm_s_ikp)
1016
1017 CALL timestop(handle)
1018
1019 END SUBROUTINE check_positive_definite_overlap_mat
1020
1021! **************************************************************************************************
1022!> \brief ...
1023!> \param qs_env ...
1024!> \param bs_env ...
1025! **************************************************************************************************
1026 SUBROUTINE get_parameters_from_qs_env(qs_env, bs_env)
1027 TYPE(qs_environment_type), POINTER :: qs_env
1028 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1029
1030 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_parameters_from_qs_env'
1031
1032 INTEGER :: color_sub, handle, homo, n_ao, n_atom, &
1033 nelectron, u
1034 INTEGER, DIMENSION(3) :: periodic
1035 REAL(kind=dp) :: maxocc
1036 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1037 TYPE(cell_type), POINTER :: cell
1038 TYPE(dft_control_type), POINTER :: dft_control
1039 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1040 TYPE(mp_para_env_type), POINTER :: para_env
1041 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1042 TYPE(scf_control_type), POINTER :: scf_control
1043 TYPE(section_vals_type), POINTER :: input
1044
1045 CALL timeset(routinen, handle)
1046
1047 CALL get_qs_env(qs_env, &
1048 dft_control=dft_control, &
1049 scf_control=scf_control, &
1050 mos=mos)
1051
1052 bs_env%n_spin = dft_control%nspins
1053 IF (bs_env%n_spin == 1) bs_env%spin_degeneracy = 2.0_dp
1054 IF (bs_env%n_spin == 2) bs_env%spin_degeneracy = 1.0_dp
1055
1056 CALL get_mo_set(mo_set=mos(1), nao=n_ao, homo=homo, nelectron=nelectron, maxocc=maxocc)
1057 bs_env%n_ao = n_ao
1058 bs_env%n_occ_nosmear(1:2) = nint(nelectron/maxocc) !fix for moire large-cell Gamma
1059 bs_env%n_electron(1:2) = nelectron
1060 ! lowered by the diagonalization below if the canonical orthogonalization removes modes
1061 bs_env%n_mo_retained = n_ao
1062 bs_env%n_occ(1:2) = homo
1063 bs_env%n_vir(1:2) = n_ao - homo
1064 IF (bs_env%n_spin == 2) THEN
1065 CALL get_mo_set(mo_set=mos(2), homo=homo, nelectron=nelectron, maxocc=maxocc)
1066 bs_env%n_occ_nosmear(2) = nint(nelectron/maxocc)
1067 bs_env%n_occ(2) = homo
1068 bs_env%n_vir(2) = n_ao - homo
1069 bs_env%n_electron(2) = nelectron
1070 END IF
1071 bs_env%eps_eigval_mat_s = scf_control%eps_eigval
1072 ! reuse the ground-state SCF's &SMEAR settings (if any) for the metal-capable Fermi
1073 ! level search (get_fermi_level_1c/2c); SMEAR_ENERGY_WINDOW and SMEAR_LIST are not
1074 ! supported by Smearkp/Smearkp2, so fall back to the T=0 bisection for those
1075 bs_env%do_smear = .false.
1076 bs_env%smear_method = smear_fermi_dirac
1077 bs_env%smear_width = 0.0_dp
1078 IF (scf_control%smear%do_smear) THEN
1079 SELECT CASE (scf_control%smear%method)
1080 CASE (smear_fermi_dirac)
1081 bs_env%do_smear = .true.
1082 bs_env%smear_method = smear_fermi_dirac
1083 bs_env%smear_width = scf_control%smear%electronic_temperature
1084 CASE (smear_gaussian, smear_mp, smear_mv)
1085 bs_env%do_smear = .true.
1086 bs_env%smear_method = scf_control%smear%method
1087 bs_env%smear_width = scf_control%smear%smearing_width
1088 CASE DEFAULT
1089 IF (bs_env%unit_nr > 0) THEN
1090 WRITE (bs_env%unit_nr, '(T2,A)') &
1091 'SMEAR method not supported for the post-SCF Fermi level search, '// &
1092 'falling back to a T=0 bisection.'
1093 END IF
1094 END SELECT
1095 END IF
1096 cpassert(scf_control%smear%fixed_mag_mom <= 0.0_dp)
1097 ! get para_env from qs_env (bs_env%para_env is identical to para_env in qs_env)
1098 CALL get_qs_env(qs_env, para_env=para_env)
1099 color_sub = 0
1100 ALLOCATE (bs_env%para_env)
1101 CALL bs_env%para_env%from_split(para_env, color_sub)
1102
1103 CALL get_qs_env(qs_env, particle_set=particle_set)
1104
1105 n_atom = SIZE(particle_set)
1106 bs_env%n_atom = n_atom
1107
1108 CALL get_qs_env(qs_env=qs_env, cell=cell)
1109 CALL get_cell(cell=cell, periodic=periodic, h=hmat)
1110 bs_env%periodic(1:3) = periodic(1:3)
1111 bs_env%do_periodic = any(bs_env%periodic /= 0)
1112 bs_env%hmat(1:3, 1:3) = hmat
1113 bs_env%nimages_scf = dft_control%nimages
1114 IF (dft_control%nimages == 1) THEN
1115 IF (bs_env%do_gw_ri_rs) THEN
1116 IF (bs_env%do_periodic) THEN
1117 bs_env%gw_implementation = ri_rs_large_cell_gamma
1118 ELSE
1119 bs_env%gw_implementation = ri_rs_non_periodic
1120 END IF
1121 ELSE IF (bs_env%do_static_cohsex) THEN
1122 bs_env%gw_implementation = tensor_static_cohsex
1123 ELSE
1124 bs_env%gw_implementation = tensor_large_cell_gamma
1125 END IF
1126 ELSE IF (dft_control%nimages > 1) THEN
1127 IF (bs_env%do_gw_ri_rs) THEN
1128 cpabort("RI-RS Not Implemented for K-point Calculations")
1129 ELSE IF (bs_env%do_static_cohsex) THEN
1130 cpabort("COHSEX Not Implemented for K-point Calculations")
1131 ELSE
1132 bs_env%gw_implementation = tensor_small_cell_full_kp
1133 END IF
1134 ELSE
1135 cpabort("Wrong number of cells from DFT calculation.")
1136 END IF
1137
1138 u = bs_env%unit_nr
1139
1140 ! Marek : Get and save the rtp method
1141 CALL get_qs_env(qs_env=qs_env, input=input)
1142 CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%_SECTION_PARAMETERS_", i_val=bs_env%rtp_method)
1143
1144 IF (u > 0) THEN
1145 WRITE (u, fmt="(T2,2A,T73,I8)") "Number of occupied molecular orbitals (MOs) ", &
1146 "= Number of occupied bands", homo
1147 WRITE (u, fmt="(T2,2A,T73,I8)") "Number of unoccupied (= virtual) MOs ", &
1148 "= Number of unoccupied bands", n_ao - homo
1149 WRITE (u, fmt="(T2,A,T73,I8)") "Number of Gaussian basis functions for MOs", n_ao
1150 IF (bs_env%gw_implementation == tensor_small_cell_full_kp) THEN
1151 WRITE (u, fmt="(T2,2A,T73,I8)") "Number of cells considered in the DFT ", &
1152 "calculation", bs_env%nimages_scf
1153 END IF
1154 END IF
1155
1156 CALL timestop(handle)
1157
1158 END SUBROUTINE get_parameters_from_qs_env
1159
1160! **************************************************************************************************
1161!> \brief ...
1162!> \param bs_env ...
1163! **************************************************************************************************
1164 SUBROUTINE set_heuristic_parameters(bs_env)
1165 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1166
1167 CHARACTER(LEN=*), PARAMETER :: routinen = 'set_heuristic_parameters'
1168
1169 INTEGER :: handle
1170
1171 CALL timeset(routinen, handle)
1172
1173 bs_env%n_bins_max_for_printing = 5000
1174
1175 CALL timestop(handle)
1176
1177 END SUBROUTINE set_heuristic_parameters
1178
1179! **************************************************************************************************
1180!> \brief ...
1181!> \param qs_env ...
1182!> \param bs_env ...
1183! **************************************************************************************************
1184 SUBROUTINE allocate_and_fill_fm_ks_fm_s(qs_env, bs_env)
1185 TYPE(qs_environment_type), POINTER :: qs_env
1186 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1187
1188 CHARACTER(LEN=*), PARAMETER :: routinen = 'allocate_and_fill_fm_ks_fm_s'
1189
1190 INTEGER :: handle, i_work, ispin
1191 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1192 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1193 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
1194 TYPE(mp_para_env_type), POINTER :: para_env
1195
1196 CALL timeset(routinen, handle)
1197
1198 CALL get_qs_env(qs_env, &
1199 para_env=para_env, &
1200 blacs_env=blacs_env, &
1201 matrix_ks_kp=matrix_ks, &
1202 matrix_s_kp=matrix_s)
1203
1204 NULLIFY (fm_struct)
1205 CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=bs_env%n_ao, &
1206 ncol_global=bs_env%n_ao, para_env=para_env)
1207
1208 DO i_work = 1, SIZE(bs_env%fm_work_mo)
1209 CALL cp_fm_create(bs_env%fm_work_mo(i_work), fm_struct)
1210 END DO
1211
1212 CALL cp_cfm_create(bs_env%cfm_work_mo, fm_struct)
1213 CALL cp_cfm_create(bs_env%cfm_work_mo_2, fm_struct)
1214
1215 CALL cp_fm_create(bs_env%fm_s_Gamma, fm_struct)
1216 CALL copy_dbcsr_to_fm(matrix_s(1, 1)%matrix, bs_env%fm_s_Gamma)
1217
1218 DO ispin = 1, bs_env%n_spin
1219 CALL cp_fm_create(bs_env%fm_ks_Gamma(ispin), fm_struct)
1220 CALL copy_dbcsr_to_fm(matrix_ks(ispin, 1)%matrix, bs_env%fm_ks_Gamma(ispin))
1221 CALL cp_fm_create(bs_env%fm_mo_coeff_Gamma(ispin), fm_struct)
1222 END DO
1223
1224 CALL cp_fm_struct_release(fm_struct)
1225
1226 NULLIFY (bs_env%mat_ao_ao%matrix)
1227 ALLOCATE (bs_env%mat_ao_ao%matrix)
1228 CALL dbcsr_create(bs_env%mat_ao_ao%matrix, template=matrix_s(1, 1)%matrix, &
1229 matrix_type=dbcsr_type_no_symmetry)
1230
1231 ALLOCATE (bs_env%eigenval_scf(bs_env%n_ao, bs_env%nkp_bs_and_DOS, bs_env%n_spin))
1232
1233 CALL timestop(handle)
1234
1235 END SUBROUTINE allocate_and_fill_fm_ks_fm_s
1236
1237! **************************************************************************************************
1238!> \brief ...
1239!> \param qs_env ...
1240!> \param bs_env ...
1241! **************************************************************************************************
1242 SUBROUTINE eval_bandstructure_properties(qs_env, bs_env)
1243 TYPE(qs_environment_type), POINTER :: qs_env
1244 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1245
1246 CHARACTER(LEN=*), PARAMETER :: routinen = 'eval_bandstructure_properties'
1247
1248 CHARACTER(LEN=default_string_length) :: gw_label, gw_soc_label
1249 INTEGER :: handle, homo, homo_1, homo_2, &
1250 homo_spinor, ikp, ikp_for_file, ispin, &
1251 n_ao, n_e, nkind, nkp
1252 LOGICAL :: is_bandstruc_kpoint, print_dos_kpoints, &
1253 print_ikp
1254 REAL(kind=dp) :: broadening, e_max, e_max_g0w0, e_min, e_min_g0w0, e_total_window, &
1255 energy_step_dos, energy_window_dos, n_el_spinor, t1
1256 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dos_g0w0, dos_g0w0_soc, dos_scf, dos_scf_soc, &
1257 eigenval, eigenval_spinor, eigenval_spinor_g0w0, eigenval_spinor_no_soc
1258 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenval_spinor_all_gw, &
1259 eigenval_spinor_all_scf, pdos_g0w0, pdos_g0w0_soc, pdos_scf, pdos_scf_soc, proj_mo_on_kind
1260 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ldos_g0w0_2d, ldos_scf_2d, &
1261 ldos_scf_2d_soc
1262 TYPE(band_edges_type) :: band_edges_g0w0, band_edges_g0w0_soc, &
1263 band_edges_scf, band_edges_scf_guess, &
1264 band_edges_scf_soc
1265 TYPE(cp_cfm_type) :: cfm_ks_ikp, cfm_ks_ikp_spinor, cfm_mos_ikp_spinor, cfm_s_ikp, &
1266 cfm_s_ikp_copy, cfm_s_ikp_spinor, cfm_s_ikp_spinor_copy, cfm_soc_ikp_spinor, &
1267 cfm_spinor_wf_ikp, cfm_work_ikp, cfm_work_ikp_spinor
1268 TYPE(cp_cfm_type), DIMENSION(2) :: cfm_mos_ikp
1269
1270 CALL timeset(routinen, handle)
1271
1272 n_ao = bs_env%n_ao
1273
1274 energy_window_dos = bs_env%energy_window_DOS
1275 energy_step_dos = bs_env%energy_step_DOS
1276 broadening = bs_env%broadening_DOS
1277
1278 ! if we have done GW or a full kpoint SCF, we already have the band edges
1279 IF (bs_env%do_gw .OR. &
1280 bs_env%gw_implementation == tensor_small_cell_full_kp) THEN
1281 band_edges_scf = bs_env%band_edges_scf
1282 band_edges_scf_guess = band_edges_scf
1283 ELSE
1284
1285 IF (bs_env%n_spin == 1) THEN
1286 homo = bs_env%n_occ(1)
1287 band_edges_scf_guess%VBM = bs_env%eigenval_scf_Gamma(homo, 1)
1288 band_edges_scf_guess%CBM = bs_env%eigenval_scf_Gamma(homo + 1, 1)
1289 ELSE
1290 homo_1 = bs_env%n_occ(1)
1291 homo_2 = bs_env%n_occ(2)
1292 band_edges_scf_guess%VBM = max(bs_env%eigenval_scf_Gamma(homo_1, 1), &
1293 bs_env%eigenval_scf_Gamma(homo_2, 2))
1294 band_edges_scf_guess%CBM = min(bs_env%eigenval_scf_Gamma(homo_1 + 1, 1), &
1295 bs_env%eigenval_scf_Gamma(homo_2 + 1, 2))
1296 END IF
1297
1298 ! initialization
1299 band_edges_scf%VBM = -1000.0_dp
1300 band_edges_scf%CBM = 1000.0_dp
1301 band_edges_scf%DBG = 1000.0_dp
1302 END IF
1303
1304 e_min = band_edges_scf_guess%VBM - 0.5_dp*energy_window_dos
1305 e_max = band_edges_scf_guess%CBM + 0.5_dp*energy_window_dos
1306
1307 IF (bs_env%do_gw) THEN
1308 band_edges_g0w0 = bs_env%band_edges_GW
1309 e_min_g0w0 = band_edges_g0w0%VBM - 0.5_dp*energy_window_dos
1310 e_max_g0w0 = band_edges_g0w0%CBM + 0.5_dp*energy_window_dos
1311 e_min = min(e_min, e_min_g0w0)
1312 e_max = max(e_max, e_max_g0w0)
1313 END IF
1314
1315 e_total_window = e_max - e_min
1316
1317 n_e = int(e_total_window/energy_step_dos)
1318
1319 CALL get_qs_env(qs_env, nkind=nkind)
1320
1321 ALLOCATE (proj_mo_on_kind(n_ao, nkind))
1322 proj_mo_on_kind(:, :) = 0.0_dp
1323
1324 ALLOCATE (eigenval(n_ao))
1325 ALLOCATE (eigenval_spinor(2*n_ao))
1326 ALLOCATE (eigenval_spinor_no_soc(2*n_ao))
1327 ALLOCATE (eigenval_spinor_g0w0(2*n_ao))
1328
1329 ! collect the SOC spinor eigenvalues of every DOS k-point (band-structure-path
1330 ! k-points are included but carry wkp=0, see get_fermi_level_1c) to later obtain
1331 ! the SCF+SOC (and, if requested, G0W0+SOC) Fermi level from get_fermi_level_1c
1332 IF (bs_env%do_soc) THEN
1333 ALLOCATE (eigenval_spinor_all_scf(2*n_ao, bs_env%nkp_bs_and_DOS))
1334 IF (bs_env%do_gw) THEN
1335 ALLOCATE (eigenval_spinor_all_gw(2*n_ao, bs_env%nkp_bs_and_DOS))
1336 END IF
1337 END IF
1338
1339 IF (bs_env%do_dos_pdos) THEN
1340
1341 ALLOCATE (dos_scf(n_e))
1342 dos_scf(:) = 0.0_dp
1343 ALLOCATE (pdos_scf(n_e, nkind))
1344 pdos_scf(:, :) = 0.0_dp
1345
1346 IF (bs_env%do_soc) THEN
1347
1348 ALLOCATE (dos_scf_soc(n_e))
1349 dos_scf_soc(:) = 0.0_dp
1350 ALLOCATE (pdos_scf_soc(n_e, nkind))
1351 pdos_scf_soc(:, :) = 0.0_dp
1352
1353 END IF
1354
1355 IF (bs_env%do_gw) THEN
1356
1357 ALLOCATE (dos_g0w0(n_e))
1358 dos_g0w0(:) = 0.0_dp
1359 ALLOCATE (pdos_g0w0(n_e, nkind))
1360 pdos_g0w0(:, :) = 0.0_dp
1361
1362 IF (bs_env%do_soc) THEN
1363
1364 ALLOCATE (dos_g0w0_soc(n_e))
1365 dos_g0w0_soc(:) = 0.0_dp
1366 ALLOCATE (pdos_g0w0_soc(n_e, nkind))
1367 pdos_g0w0_soc(:, :) = 0.0_dp
1368
1369 END IF
1370 END IF
1371 END IF
1372
1373 CALL cp_cfm_create(cfm_mos_ikp(1), bs_env%fm_ks_Gamma(1)%matrix_struct)
1374 CALL cp_cfm_create(cfm_mos_ikp(2), bs_env%fm_ks_Gamma(1)%matrix_struct)
1375 CALL cp_cfm_create(cfm_work_ikp, bs_env%fm_ks_Gamma(1)%matrix_struct)
1376 CALL cp_cfm_create(cfm_s_ikp_copy, bs_env%fm_ks_Gamma(1)%matrix_struct)
1377
1378 IF (bs_env%do_soc) THEN
1379
1380 CALL cp_cfm_create(cfm_mos_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1381 CALL cp_cfm_create(cfm_work_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1382 CALL cp_cfm_create(cfm_s_ikp_spinor_copy, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1383 CALL cp_cfm_create(cfm_ks_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1384 CALL cp_cfm_create(cfm_soc_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1385 CALL cp_cfm_create(cfm_s_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1386 CALL cp_cfm_create(cfm_spinor_wf_ikp, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1387
1388 homo_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin)
1389
1390 SELECT CASE (bs_env%n_spin)
1391 CASE (1)
1392 n_el_spinor = real(bs_env%n_electron(1), kind=dp)
1393 CASE (2)
1394 n_el_spinor = real(bs_env%n_electron(1) + bs_env%n_electron(2), kind=dp)
1395 END SELECT
1396
1397 band_edges_scf_soc%VBM = -1000.0_dp
1398 band_edges_scf_soc%CBM = 1000.0_dp
1399 band_edges_scf_soc%DBG = 1000.0_dp
1400
1401 IF (bs_env%do_gw) THEN
1402 band_edges_g0w0_soc%VBM = -1000.0_dp
1403 band_edges_g0w0_soc%CBM = 1000.0_dp
1404 band_edges_g0w0_soc%DBG = 1000.0_dp
1405 END IF
1406
1407 IF (bs_env%unit_nr > 0) THEN
1408 WRITE (bs_env%unit_nr, '(A)') ''
1409 IF (bs_env%soc_window_occ > 0.0_dp) THEN
1410 WRITE (bs_env%unit_nr, '(T2,A,T71,F10.2)') 'SOC requested, SOC energy window occ (eV):', &
1411 bs_env%soc_window_occ*evolt
1412 ELSE
1413 WRITE (bs_env%unit_nr, '(T2,A,T71,A10)') 'SOC requested, SOC energy window occ (eV):', &
1414 ' no window'
1415 END IF
1416 IF (bs_env%soc_window_virt > 0.0_dp) THEN
1417 WRITE (bs_env%unit_nr, '(T2,A,T71,F10.2)') 'SOC requested, SOC energy window virt (eV):', &
1418 bs_env%soc_window_virt*evolt
1419 ELSE
1420 WRITE (bs_env%unit_nr, '(T2,A,T71,A10)') 'SOC requested, SOC energy window virt (eV):', &
1421 ' no window'
1422 END IF
1423 IF (bs_env%soc_window_occ > 0.0_dp .OR. bs_env%soc_window_virt > 0.0_dp) THEN
1424 WRITE (bs_env%unit_nr, '(T2,A,T71,F10.2)') 'SOC requested, SOC window smearing (eV):', &
1425 bs_env%soc_window_smearing*evolt
1426 END IF
1427 END IF
1428 END IF
1429
1430 IF (bs_env%do_ldos) THEN
1431 cpassert(bs_env%int_ldos_xyz == int_ldos_z)
1432 END IF
1433
1434 IF (bs_env%gw_implementation == tensor_small_cell_full_kp) THEN
1435 CALL cp_cfm_create(cfm_ks_ikp, bs_env%cfm_ks_kp(1, 1)%matrix_struct)
1436 CALL cp_cfm_create(cfm_s_ikp, bs_env%cfm_ks_kp(1, 1)%matrix_struct)
1437 END IF
1438
1439 DO ikp = 1, bs_env%nkp_bs_and_DOS
1440
1441 t1 = m_walltime()
1442
1443 DO ispin = 1, bs_env%n_spin
1444
1445 SELECT CASE (bs_env%gw_implementation)
1446 CASE (tensor_large_cell_gamma, ri_rs_large_cell_gamma, ri_rs_non_periodic, tensor_static_cohsex)
1447
1448 ! 1. get H^KS_µν(k_i) from H^KS_µν(k=0)
1449 CALL cfm_ikp_from_fm_gamma(cfm_ks_ikp, bs_env%fm_ks_Gamma(ispin), &
1450 ikp, qs_env, bs_env%kpoints_DOS, "ORB")
1451
1452 ! 2. get S_µν(k_i) from S_µν(k=0)
1453 CALL cfm_ikp_from_fm_gamma(cfm_s_ikp, bs_env%fm_s_Gamma, &
1454 ikp, qs_env, bs_env%kpoints_DOS, "ORB")
1455 CALL cp_cfm_to_cfm(cfm_s_ikp, cfm_s_ikp_copy)
1456
1457 ! 3. Diagonalize (Roothaan-Hall): H_KS(k_i)*C(k_i) = S(k_i)*C(k_i)*ϵ(k_i)
1458 CALL cp_cfm_geeig(cfm_ks_ikp, cfm_s_ikp_copy, cfm_mos_ikp(ispin), &
1459 eigenval, cfm_work_ikp)
1460
1461 bs_env%eigenval_scf(:, ikp, ispin) = eigenval(:)
1462
1463 CASE (tensor_small_cell_full_kp)
1464
1465 ! 1. get H^KS_µν(k_i)
1466 CALL cp_cfm_to_cfm(bs_env%cfm_ks_kp(ikp, ispin), cfm_ks_ikp)
1467
1468 ! 2. get S_µν(k_i)
1469 CALL cp_cfm_to_cfm(bs_env%cfm_s_kp(ikp), cfm_s_ikp)
1470
1471 ! 3. get C_µn(k_i) and ϵ_n(k_i)
1472 CALL cp_cfm_to_cfm(bs_env%cfm_mo_coeff_kp(ikp, ispin), cfm_mos_ikp(ispin))
1473 eigenval(:) = bs_env%eigenval_scf(:, ikp, ispin)
1474
1475 END SELECT
1476
1477 ! 4. Projection p_nk^A of MO ψ_nk(r) on atom type A (inspired by Mulliken charge)
1478 ! p_nk^A = sum_µ^A,ν C*_µ^A,n(k) S_µ^A,ν(k) C_ν,n(k)
1479 CALL compute_proj_mo_on_kind(proj_mo_on_kind, qs_env, cfm_mos_ikp(ispin), cfm_s_ikp)
1480
1481 ! 5. DOS and PDOS
1482 IF (bs_env%do_dos_pdos) THEN
1483 CALL add_to_dos_pdos(dos_scf, pdos_scf, eigenval, ikp, bs_env, n_e, e_min, &
1484 proj_mo_on_kind)
1485
1486 IF (bs_env%do_gw) THEN
1487 CALL add_to_dos_pdos(dos_g0w0, pdos_g0w0, bs_env%eigenval_GW(:, ikp, ispin), &
1488 ikp, bs_env, n_e, e_min, proj_mo_on_kind)
1489 END IF
1490 END IF
1491
1492 IF (bs_env%do_ldos) THEN
1493 CALL add_to_ldos_2d(ldos_scf_2d, qs_env, ikp, bs_env, cfm_mos_ikp(ispin), &
1494 eigenval(:), band_edges_scf_guess)
1495
1496 IF (bs_env%do_gw) THEN
1497 CALL add_to_ldos_2d(ldos_g0w0_2d, qs_env, ikp, bs_env, cfm_mos_ikp(ispin), &
1498 bs_env%eigenval_GW(:, ikp, 1), band_edges_g0w0)
1499 END IF
1500
1501 END IF
1502
1503 homo = bs_env%n_occ(ispin)
1504
1505 band_edges_scf%VBM = max(band_edges_scf%VBM, eigenval(homo))
1506 band_edges_scf%CBM = min(band_edges_scf%CBM, eigenval(homo + 1))
1507 band_edges_scf%DBG = min(band_edges_scf%DBG, eigenval(homo + 1) - eigenval(homo))
1508
1509 END DO ! spin
1510
1511 ! now the same with spin-orbit coupling
1512 IF (bs_env%do_soc) THEN
1513
1514 ! only print eigenvalues of DOS k-points in case no bandstructure path has been given
1515 print_dos_kpoints = (bs_env%nkp_only_bs <= 0)
1516 ! in kpoints_DOS, the last nkp_only_bs are bandstructure k-points
1517 is_bandstruc_kpoint = (ikp > bs_env%nkp_only_DOS)
1518 print_ikp = print_dos_kpoints .OR. is_bandstruc_kpoint
1519
1520 IF (print_dos_kpoints) THEN
1521 nkp = bs_env%nkp_only_DOS
1522 ikp_for_file = ikp
1523 ELSE
1524 nkp = bs_env%nkp_only_bs
1525 ikp_for_file = ikp - bs_env%nkp_only_DOS
1526 END IF
1527
1528 ! compute DFT+SOC eigenvalues; based on these, compute band edges, DOS and LDOS
1529 CALL soc_ev(bs_env, qs_env, ikp, bs_env%eigenval_scf, &
1530 e_min, cfm_mos_ikp, dos_scf_soc, pdos_scf_soc, &
1531 band_edges_scf_soc, eigenval_spinor, cfm_spinor_wf_ikp)
1532 eigenval_spinor_all_scf(:, ikp) = eigenval_spinor
1533 IF (.NOT. bs_env%do_gw .AND. print_ikp) THEN
1534 CALL write_soc_eigenvalues(eigenval_spinor, ikp_for_file, ikp, bs_env)
1535 END IF
1536
1537 IF (bs_env%do_ldos) THEN
1538 CALL add_to_ldos_2d(ldos_scf_2d_soc, qs_env, ikp, bs_env, cfm_spinor_wf_ikp, &
1539 eigenval_spinor, band_edges_scf_guess, .true., cfm_work_ikp)
1540 END IF
1541
1542 IF (bs_env%do_gw) THEN
1543
1544 ! compute G0W0+SOC eigenvalues; based on these, compute band edges, DOS and LDOS
1545 CALL soc_ev(bs_env, qs_env, ikp, bs_env%eigenval_GW, &
1546 e_min, cfm_mos_ikp, dos_g0w0_soc, pdos_g0w0_soc, &
1547 band_edges_g0w0_soc, eigenval_spinor_g0w0, cfm_spinor_wf_ikp)
1548 eigenval_spinor_all_gw(:, ikp) = eigenval_spinor_g0w0
1549 IF (print_ikp) THEN
1550 ! write SCF+SOC and G0W0+SOC eigenvalues to file
1551 ! SCF_and_G0W0_band_structure_for_kpoint_<ikp>_+_SOC
1552 CALL write_soc_eigenvalues(eigenval_spinor, ikp_for_file, ikp, bs_env, &
1553 eigenval_spinor_g0w0)
1554 END IF
1555
1556 END IF ! do_gw
1557
1558 END IF ! do_soc
1559
1560 IF (bs_env%unit_nr > 0 .AND. m_walltime() - t1 > 20.0_dp) THEN
1561 WRITE (bs_env%unit_nr, '(T2,A,T43,I5,A,I3,A,F7.1,A)') &
1562 'Compute DOS, LDOS for k-point ', ikp, ' /', bs_env%nkp_bs_and_DOS, &
1563 ', Execution time', m_walltime() - t1, ' s'
1564 END IF
1565
1566 END DO ! ikp_DOS
1567
1568 band_edges_scf%IDBG = band_edges_scf%CBM - band_edges_scf%VBM
1569
1570 SELECT CASE (bs_env%n_spin)
1571 CASE (1)
1572 CALL get_fermi_level_1c(bs_env, bs_env%eigenval_scf(:, :, 1), &
1573 REAL(bs_env%n_electron(1), kind=dp), &
1574 bs_env%spin_degeneracy, band_edges_scf%VBM, &
1575 band_edges_scf%CBM, band_edges_scf%E_fermi)
1576 CASE (2)
1577 CALL get_fermi_level_2c(bs_env, bs_env%eigenval_scf(:, :, 1:2), &
1578 REAL(bs_env%n_electron(1) + bs_env%n_electron(2), kind=dp), &
1579 band_edges_scf%VBM, band_edges_scf%CBM, band_edges_scf%E_fermi)
1580 END SELECT
1581
1582 IF (bs_env%do_soc) THEN
1583 band_edges_scf_soc%IDBG = band_edges_scf_soc%CBM - band_edges_scf_soc%VBM
1584 ! SOC merges both spins into 2*n_ao spinor states, each holding at most one
1585 ! electron (maxocc=1); homo_spinor = n_occ(1)+n_occ(n_spin) is the total spinor
1586 ! electron count, same quantity SOC_ev uses to locate the occupied/unoccupied edge
1587 CALL get_fermi_level_1c(bs_env, eigenval_spinor_all_scf, &
1588 n_el_spinor, 1.0_dp, &
1589 band_edges_scf_soc%VBM, band_edges_scf_soc%CBM, band_edges_scf_soc%E_fermi)
1590
1591 IF (bs_env%do_gw) THEN
1592 band_edges_g0w0_soc%IDBG = band_edges_g0w0_soc%CBM - band_edges_g0w0_soc%VBM
1593 CALL get_fermi_level_1c(bs_env, eigenval_spinor_all_gw, &
1594 n_el_spinor, 1.0_dp, &
1595 band_edges_g0w0_soc%VBM, band_edges_g0w0_soc%CBM, band_edges_g0w0_soc%E_fermi)
1596
1597 END IF
1598 END IF
1599
1600 CALL write_band_edges(band_edges_scf, "SCF", bs_env)
1601 IF (bs_env%do_dos_pdos) THEN
1602 CALL write_dos_pdos(dos_scf, pdos_scf, bs_env, qs_env, "SCF", e_min, band_edges_scf%VBM)
1603 END IF
1604 IF (bs_env%do_ldos) THEN
1605 CALL print_ldos_main(ldos_scf_2d, bs_env, band_edges_scf, "SCF")
1606 END IF
1607
1608 IF (bs_env%do_soc) THEN
1609 CALL write_band_edges(band_edges_scf_soc, "SCF+SOC", bs_env)
1610 IF (bs_env%do_dos_pdos) THEN
1611 CALL write_dos_pdos(dos_scf_soc, pdos_scf_soc, bs_env, qs_env, "SCF_SOC", &
1612 e_min, band_edges_scf_soc%VBM)
1613 END IF
1614 IF (bs_env%do_ldos) THEN
1615 ! argument band_edges_scf is actually correct because the non-SOC band edges
1616 ! have been used as reference in add_to_LDOS_2d
1617 CALL print_ldos_main(ldos_scf_2d_soc, bs_env, band_edges_scf, &
1618 "SCF_SOC")
1619 END IF
1620 END IF
1621
1622 ! the printed band edges carry the GW flavour that was actually run; the labels handed to
1623 ! write_dos_pdos / print_LDOS_main below become file names and stay as they are
1624 gw_label = gw_flavour_label(bs_env)
1625 gw_soc_label = trim(gw_label)//"+SOC"
1626
1627 IF (bs_env%do_gw) THEN
1628 CALL write_band_edges(band_edges_g0w0, trim(gw_label), bs_env)
1629 CALL write_band_edges(bs_env%band_edges_HF, "Hartree-Fock with SCF orbitals", bs_env)
1630 IF (bs_env%do_dos_pdos) THEN
1631 CALL write_dos_pdos(dos_g0w0, pdos_g0w0, bs_env, qs_env, "G0W0", e_min, &
1632 band_edges_g0w0%VBM)
1633 END IF
1634 IF (bs_env%do_ldos) THEN
1635 CALL print_ldos_main(ldos_g0w0_2d, bs_env, band_edges_g0w0, "G0W0")
1636 END IF
1637 END IF
1638
1639 IF (bs_env%do_soc .AND. bs_env%do_gw) THEN
1640 CALL write_band_edges(band_edges_g0w0_soc, trim(gw_soc_label), bs_env)
1641 IF (bs_env%do_dos_pdos) THEN
1642 CALL write_dos_pdos(dos_g0w0_soc, pdos_g0w0_soc, bs_env, qs_env, "G0W0_SOC", e_min, &
1643 band_edges_g0w0_soc%VBM)
1644 END IF
1645 END IF
1646
1647 CALL cp_cfm_release(cfm_s_ikp)
1648 CALL cp_cfm_release(cfm_ks_ikp)
1649 CALL cp_cfm_release(cfm_mos_ikp(1))
1650 CALL cp_cfm_release(cfm_mos_ikp(2))
1651 CALL cp_cfm_release(cfm_work_ikp)
1652 CALL cp_cfm_release(cfm_s_ikp_copy)
1653
1654 CALL cp_cfm_release(cfm_s_ikp_spinor)
1655 CALL cp_cfm_release(cfm_ks_ikp_spinor)
1656 CALL cp_cfm_release(cfm_soc_ikp_spinor)
1657 CALL cp_cfm_release(cfm_mos_ikp_spinor)
1658 CALL cp_cfm_release(cfm_work_ikp_spinor)
1659 CALL cp_cfm_release(cfm_s_ikp_spinor_copy)
1660 CALL cp_cfm_release(cfm_spinor_wf_ikp)
1661
1662 CALL timestop(handle)
1663
1664 END SUBROUTINE eval_bandstructure_properties
1665
1666! **************************************************************************************************
1667!> \brief ...
1668!> \param LDOS_2d ...
1669!> \param bs_env ...
1670!> \param band_edges ...
1671!> \param scf_gw_soc ...
1672! **************************************************************************************************
1673 SUBROUTINE print_ldos_main(LDOS_2d, bs_env, band_edges, scf_gw_soc)
1674 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ldos_2d
1675 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1676 TYPE(band_edges_type) :: band_edges
1677 CHARACTER(LEN=*) :: scf_gw_soc
1678
1679 CHARACTER(LEN=*), PARAMETER :: routinen = 'print_LDOS_main'
1680
1681 INTEGER :: handle, i_x, i_x_bin, i_x_end, i_x_end_bin, i_x_end_glob, i_x_start, &
1682 i_x_start_bin, i_x_start_glob, i_y, i_y_bin, i_y_end, i_y_end_bin, i_y_end_glob, &
1683 i_y_start, i_y_start_bin, i_y_start_glob, n_e
1684 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: n_sum_for_bins
1685 INTEGER, DIMENSION(2) :: bin_mesh
1686 LOGICAL :: do_xy_bins
1687 REAL(kind=dp) :: e_min, energy_step, energy_window
1688 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ldos_2d_bins
1689
1690 CALL timeset(routinen, handle)
1691
1692 n_e = SIZE(ldos_2d, 3)
1693
1694 energy_window = bs_env%energy_window_DOS
1695 energy_step = bs_env%energy_step_DOS
1696 e_min = band_edges%VBM - 0.5_dp*energy_window
1697
1698 bin_mesh(1:2) = bs_env%bin_mesh(1:2)
1699 do_xy_bins = (bin_mesh(1) > 0 .AND. bin_mesh(2) > 0)
1700
1701 i_x_start = lbound(ldos_2d, 1)
1702 i_x_end = ubound(ldos_2d, 1)
1703 i_y_start = lbound(ldos_2d, 2)
1704 i_y_end = ubound(ldos_2d, 2)
1705
1706 IF (do_xy_bins) THEN
1707 i_x_start_bin = 1
1708 i_x_end_bin = bin_mesh(1)
1709 i_y_start_bin = 1
1710 i_y_end_bin = bin_mesh(2)
1711 ELSE
1712 i_x_start_bin = i_x_start
1713 i_x_end_bin = i_x_end
1714 i_y_start_bin = i_y_start
1715 i_y_end_bin = i_y_end
1716 END IF
1717
1718 ALLOCATE (ldos_2d_bins(i_x_start_bin:i_x_end_bin, i_y_start_bin:i_y_end_bin, n_e))
1719 ldos_2d_bins(:, :, :) = 0.0_dp
1720
1721 IF (do_xy_bins) THEN
1722
1723 i_x_start_glob = i_x_start
1724 i_x_end_glob = i_x_end
1725 i_y_start_glob = i_y_start
1726 i_y_end_glob = i_y_end
1727
1728 CALL bs_env%para_env%min(i_x_start_glob)
1729 CALL bs_env%para_env%max(i_x_end_glob)
1730 CALL bs_env%para_env%min(i_y_start_glob)
1731 CALL bs_env%para_env%max(i_y_end_glob)
1732
1733 ALLOCATE (n_sum_for_bins(bin_mesh(1), bin_mesh(2)), source=0)
1734
1735 ! transform interval [i_x_start, i_x_end] to [1, bin_mesh(1)] (and same for y)
1736 DO i_y = i_y_start, i_y_end
1737 DO i_x = i_x_start, i_x_end
1738 i_x_bin = bin_mesh(1)*(i_x - i_x_start_glob)/(i_x_end_glob - i_x_start_glob + 1) + 1
1739 i_y_bin = bin_mesh(2)*(i_y - i_y_start_glob)/(i_y_end_glob - i_y_start_glob + 1) + 1
1740 ldos_2d_bins(i_x_bin, i_y_bin, :) = ldos_2d_bins(i_x_bin, i_y_bin, :) + &
1741 ldos_2d(i_x, i_y, :)
1742 n_sum_for_bins(i_x_bin, i_y_bin) = n_sum_for_bins(i_x_bin, i_y_bin) + 1
1743 END DO
1744 END DO
1745
1746 CALL bs_env%para_env%sum(ldos_2d_bins)
1747 CALL bs_env%para_env%sum(n_sum_for_bins)
1748
1749 ! divide by number of terms in the sum so we have the average LDOS(x,y,E)
1750 DO i_y_bin = 1, bin_mesh(2)
1751 DO i_x_bin = 1, bin_mesh(1)
1752 ldos_2d_bins(i_x_bin, i_y_bin, :) = ldos_2d_bins(i_x_bin, i_y_bin, :)/ &
1753 REAL(n_sum_for_bins(i_x_bin, i_y_bin), kind=dp)
1754 END DO
1755 END DO
1756
1757 ELSE
1758
1759 ldos_2d_bins(:, :, :) = ldos_2d(:, :, :)
1760
1761 END IF
1762
1763 IF (bin_mesh(1)*bin_mesh(2) < bs_env%n_bins_max_for_printing) THEN
1764 CALL print_ldos_2d_bins(ldos_2d_bins, bs_env, e_min, scf_gw_soc)
1765 ELSE
1766 cpwarn("The number of bins for the LDOS is too large. Decrease BIN_MESH.")
1767 END IF
1768
1769 CALL timestop(handle)
1770
1771 END SUBROUTINE print_ldos_main
1772
1773! **************************************************************************************************
1774!> \brief ...
1775!> \param LDOS_2d_bins ...
1776!> \param bs_env ...
1777!> \param E_min ...
1778!> \param scf_gw_soc ...
1779! **************************************************************************************************
1780 SUBROUTINE print_ldos_2d_bins(LDOS_2d_bins, bs_env, E_min, scf_gw_soc)
1781 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ldos_2d_bins
1782 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1783 REAL(kind=dp) :: e_min
1784 CHARACTER(LEN=*) :: scf_gw_soc
1785
1786 CHARACTER(LEN=*), PARAMETER :: routinen = 'print_LDOS_2d_bins'
1787
1788 CHARACTER(LEN=18) :: print_format
1789 CHARACTER(LEN=4) :: print_format_1, print_format_2
1790 CHARACTER(len=default_string_length) :: fname
1791 INTEGER :: handle, i_e, i_x, i_x_end, i_x_start, &
1792 i_y, i_y_end, i_y_start, iunit, n_e, &
1793 n_x, n_y
1794 REAL(kind=dp) :: energy
1795 REAL(kind=dp), DIMENSION(3) :: coord, idx
1796
1797 CALL timeset(routinen, handle)
1798
1799 i_x_start = lbound(ldos_2d_bins, 1)
1800 i_x_end = ubound(ldos_2d_bins, 1)
1801 i_y_start = lbound(ldos_2d_bins, 2)
1802 i_y_end = ubound(ldos_2d_bins, 2)
1803 n_e = SIZE(ldos_2d_bins, 3)
1804
1805 n_x = i_x_end - i_x_start + 1
1806 n_y = i_y_end - i_y_start + 1
1807
1808 IF (bs_env%para_env%is_source()) THEN
1809
1810 DO i_y = i_y_start, i_y_end
1811 DO i_x = i_x_start, i_x_end
1812
1813 idx(1) = (real(i_x, kind=dp) - 0.5_dp)/real(n_x, kind=dp)
1814 idx(2) = (real(i_y, kind=dp) - 0.5_dp)/real(n_y, kind=dp)
1815 idx(3) = 0.0_dp
1816 coord(1:3) = matmul(bs_env%hmat, idx)
1817
1818 CALL get_print_format(coord(1), print_format_1)
1819 CALL get_print_format(coord(2), print_format_2)
1820
1821 print_format = "(3A,"//print_format_1//",A,"//print_format_2//",A)"
1822
1823 WRITE (fname, print_format) "LDOS_", scf_gw_soc, &
1824 "_at_x_", coord(1)*angstrom, '_A_and_y_', coord(2)*angstrom, '_A'
1825
1826 CALL open_file(trim(fname), unit_number=iunit, file_status="REPLACE", &
1827 file_action="WRITE")
1828
1829 WRITE (iunit, "(2A)") " Energy E (eV) average LDOS(x,y,E) (1/(eV*Å^2), ", &
1830 "integrated over z, averaged inside bin)"
1831
1832 DO i_e = 1, n_e
1833 energy = e_min + i_e*bs_env%energy_step_DOS
1834 WRITE (iunit, "(2F17.3)") energy*evolt, &
1835 ldos_2d_bins(i_x, i_y, i_e)* &
1836 bs_env%unit_ldos_int_z_inv_Ang2_eV
1837 END DO
1838
1839 CALL close_file(iunit)
1840
1841 END DO
1842 END DO
1843
1844 END IF
1845
1846 CALL timestop(handle)
1847
1848 END SUBROUTINE print_ldos_2d_bins
1849
1850! **************************************************************************************************
1851!> \brief ...
1852!> \param coord ...
1853!> \param print_format ...
1854! **************************************************************************************************
1855 SUBROUTINE get_print_format(coord, print_format)
1856 REAL(kind=dp) :: coord
1857 CHARACTER(LEN=4) :: print_format
1858
1859 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_print_format'
1860
1861 INTEGER :: handle
1862
1863 CALL timeset(routinen, handle)
1864
1865 IF (coord < -10000/angstrom) THEN
1866 print_format = "F9.2"
1867 ELSE IF (coord < -1000/angstrom) THEN
1868 print_format = "F8.2"
1869 ELSE IF (coord < -100/angstrom) THEN
1870 print_format = "F7.2"
1871 ELSE IF (coord < -10/angstrom) THEN
1872 print_format = "F6.2"
1873 ELSE IF (coord < -1/angstrom) THEN
1874 print_format = "F5.2"
1875 ELSE IF (coord < 10/angstrom) THEN
1876 print_format = "F4.2"
1877 ELSE IF (coord < 100/angstrom) THEN
1878 print_format = "F5.2"
1879 ELSE IF (coord < 1000/angstrom) THEN
1880 print_format = "F6.2"
1881 ELSE IF (coord < 10000/angstrom) THEN
1882 print_format = "F7.2"
1883 ELSE
1884 print_format = "F8.2"
1885 END IF
1886
1887 CALL timestop(handle)
1888
1889 END SUBROUTINE get_print_format
1890
1891! **************************************************************************************************
1892!> \brief ...
1893!> \param bs_env ...
1894!> \param qs_env ...
1895!> \param ikp ...
1896!> \param eigenval_no_SOC non-SOC band energies (n_ao, n_kpoints, n_spin)
1897!> \param E_min ...
1898!> \param cfm_mos_ikp ...
1899!> \param DOS ...
1900!> \param PDOS ...
1901!> \param band_edges ...
1902!> \param eigenval_spinor ...
1903!> \param cfm_spinor_wf_ikp ...
1904! **************************************************************************************************
1905 SUBROUTINE soc_ev(bs_env, qs_env, ikp, eigenval_no_SOC, E_min, cfm_mos_ikp, &
1906 DOS, PDOS, band_edges, eigenval_spinor, cfm_spinor_wf_ikp)
1907
1908 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1909 TYPE(qs_environment_type), POINTER :: qs_env
1910 INTEGER :: ikp
1911 REAL(kind=dp), DIMENSION(:, :, :) :: eigenval_no_soc
1912 REAL(kind=dp) :: e_min
1913 TYPE(cp_cfm_type), DIMENSION(2) :: cfm_mos_ikp
1914 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dos
1915 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: pdos
1916 TYPE(band_edges_type) :: band_edges
1917 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_spinor
1918 TYPE(cp_cfm_type) :: cfm_spinor_wf_ikp
1919
1920 CHARACTER(LEN=*), PARAMETER :: routinen = 'SOC_ev'
1921
1922 INTEGER :: handle
1923 LOGICAL :: do_truncate
1924
1925 CALL timeset(routinen, handle)
1926
1927 !for the window truncation part
1928 do_truncate = (bs_env%soc_nbands_occ > 0 .AND. bs_env%soc_nbands_virt > 0)
1929
1930 IF (do_truncate) THEN
1931 CALL soc_ev_truncate_mos(bs_env, qs_env, ikp, eigenval_no_soc, cfm_mos_ikp, &
1932 band_edges, eigenval_spinor, cfm_spinor_wf_ikp)
1933 ELSE
1934 CALL soc_ev_all_mos(bs_env, qs_env, ikp, eigenval_no_soc, e_min, cfm_mos_ikp, &
1935 dos, pdos, band_edges, eigenval_spinor, cfm_spinor_wf_ikp)
1936 END IF
1937
1938 CALL timestop(handle)
1939
1940 END SUBROUTINE soc_ev
1941
1942! **************************************************************************************************
1943!> \brief SOC second-variational diagonalization over the full 2*n_ao spinor space (no truncation
1944!> window) -- the SOC_NBANDS_OCC/VIRT-truncated path lives in SOC_ev_truncate_mos.
1945!> \param bs_env ...
1946!> \param qs_env ...
1947!> \param ikp ...
1948!> \param eigenval_no_SOC non-SOC band energies (n_ao, n_kpoints, n_spin)
1949!> \param E_min ...
1950!> \param cfm_mos_ikp ...
1951!> \param DOS ...
1952!> \param PDOS ...
1953!> \param band_edges ...
1954!> \param eigenval_spinor ...
1955!> \param cfm_spinor_wf_ikp ...
1956! **************************************************************************************************
1957 SUBROUTINE soc_ev_all_mos(bs_env, qs_env, ikp, eigenval_no_SOC, E_min, cfm_mos_ikp, &
1958 DOS, PDOS, band_edges, eigenval_spinor, cfm_spinor_wf_ikp)
1959
1960 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1961 TYPE(qs_environment_type), POINTER :: qs_env
1962 INTEGER :: ikp
1963 REAL(kind=dp), DIMENSION(:, :, :) :: eigenval_no_soc
1964 REAL(kind=dp) :: e_min
1965 TYPE(cp_cfm_type), DIMENSION(2) :: cfm_mos_ikp
1966 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dos
1967 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: pdos
1968 TYPE(band_edges_type) :: band_edges
1969 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_spinor
1970 TYPE(cp_cfm_type) :: cfm_spinor_wf_ikp
1971
1972 CHARACTER(LEN=*), PARAMETER :: routinen = 'SOC_ev_all_mos'
1973
1974 INTEGER :: handle, homo_spinor, n_ao, n_e, nkind
1975 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_spinor_no_soc
1976 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: proj_mo_on_kind_spinor
1977 TYPE(cp_cfm_type) :: cfm_eigenvec_ikp_spinor, &
1978 cfm_ks_ikp_spinor, cfm_mos_ikp_spinor, &
1979 cfm_soc_ikp_spinor, cfm_work_ikp_spinor
1980
1981 CALL timeset(routinen, handle)
1982
1983 n_ao = bs_env%n_ao
1984 homo_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin)
1985 CALL get_qs_env(qs_env, nkind=nkind)
1986
1987 CALL cp_cfm_create(cfm_ks_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1988 CALL cp_cfm_create(cfm_soc_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1989 CALL cp_cfm_create(cfm_mos_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1990 CALL cp_cfm_create(cfm_work_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1991 CALL cp_cfm_create(cfm_eigenvec_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1992
1993 ALLOCATE (eigenval_spinor_no_soc(2*n_ao))
1994 ALLOCATE (proj_mo_on_kind_spinor(2*n_ao, nkind))
1995 ! PDOS not yet implemented -> projection is just zero -> PDOS is zero
1996 proj_mo_on_kind_spinor(:, :) = 0.0_dp
1997
1998 ! 1. get V^SOC_µν,σσ'(k_i)
1999 SELECT CASE (bs_env%gw_implementation)
2000 CASE (tensor_large_cell_gamma, ri_rs_large_cell_gamma, ri_rs_non_periodic, tensor_static_cohsex)
2001
2002 ! 1. get V^SOC_µν,σσ'(k_i) from V^SOC_µν,σσ'(k=0)
2003 CALL cfm_ikp_from_cfm_spinor_gamma(cfm_soc_ikp_spinor, &
2004 bs_env%cfm_SOC_spinor_ao(1), &
2005 bs_env%fm_s_Gamma%matrix_struct, &
2006 ikp, qs_env, bs_env%kpoints_DOS, "ORB")
2007
2008 CASE (tensor_small_cell_full_kp)
2009
2010 ! 1. V^SOC_µν,σσ'(k_i) already there
2011 CALL cp_cfm_to_cfm(bs_env%cfm_SOC_spinor_ao(ikp), cfm_soc_ikp_spinor)
2012
2013 END SELECT
2014
2015 ! 2. V^SOC_nn',σσ'(k_i) = sum_µν C^*_µn,σ(k_i) V^SOC_µν,σσ'(k_i) C_νn'(k_i),
2016 ! C_µn,σ(k_i): MO coefficiencts from diagonalizing KS-matrix h^KS_nn',σσ'(k_i)
2017
2018 ! 2.1 build matrix C_µn,σ(k_i)
2019 CALL cp_cfm_set_all(cfm_mos_ikp_spinor, z_zero)
2020 CALL add_cfm_submat(cfm_mos_ikp_spinor, cfm_mos_ikp(1), 1, 1)
2021 CALL add_cfm_submat(cfm_mos_ikp_spinor, cfm_mos_ikp(bs_env%n_spin), n_ao + 1, n_ao + 1)
2022
2023 ! 2.2 work_nν,σσ' = sum_µ C^*_µn,σ(k_i) V^SOC_µν,σσ'(k_i)
2024 CALL parallel_gemm('C', 'N', 2*n_ao, 2*n_ao, 2*n_ao, z_one, &
2025 cfm_mos_ikp_spinor, cfm_soc_ikp_spinor, &
2026 z_zero, cfm_work_ikp_spinor)
2027
2028 ! 2.3 V^SOC_nn',σσ'(k_i) = sum_ν work_nν,σσ' C_νn'(k_i)
2029 CALL parallel_gemm('N', 'N', 2*n_ao, 2*n_ao, 2*n_ao, z_one, &
2030 cfm_work_ikp_spinor, cfm_mos_ikp_spinor, &
2031 z_zero, cfm_ks_ikp_spinor)
2032
2033 ! 3. remove SOC outside of energy window (otherwise, numerical problems arise
2034 ! because energetically low semicore states and energetically very high
2035 ! unbound states couple to the states around the Fermi level)
2036 eigenval_spinor_no_soc(1:n_ao) = eigenval_no_soc(1:n_ao, ikp, 1)
2037 eigenval_spinor_no_soc(n_ao + 1:) = eigenval_no_soc(1:n_ao, ikp, bs_env%n_spin)
2038 IF (bs_env%soc_window_occ > 0.0_dp .OR. bs_env%soc_window_virt > 0.0_dp) THEN
2039 CALL remove_soc_outside_energy_window_mo(cfm_ks_ikp_spinor, &
2040 bs_env%soc_window_virt, &
2041 bs_env%soc_window_smearing, &
2042 eigenval_spinor_no_soc, &
2043 bs_env%e_fermi(1))
2044
2045 END IF
2046
2047 ! 4. h^G0W0+SOC_nn',σσ'(k_i) = ε_nσ^G0W0(k_i) δ_nn' δ_σσ' + V^SOC_nn',σσ'(k_i)
2048 CALL cfm_add_on_diag(cfm_ks_ikp_spinor, eigenval_spinor_no_soc)
2049
2050 ! 5. diagonalize h^G0W0+SOC_nn',σσ'(k_i) to get eigenvalues
2051 CALL cp_cfm_heevd(cfm_ks_ikp_spinor, cfm_eigenvec_ikp_spinor, eigenval_spinor)
2052
2053 ! 6. DOS from spinors, no PDOS
2054 IF (bs_env%do_dos_pdos) THEN
2055 n_e = SIZE(dos)
2056 CALL add_to_dos_pdos(dos, pdos, eigenval_spinor, &
2057 ikp, bs_env, n_e, e_min, proj_mo_on_kind_spinor)
2058 END IF
2059
2060 ! 7. valence band max. (VBM), conduction band min. (CBM) and direct bandgap (DBG)
2061 band_edges%VBM = max(band_edges%VBM, eigenval_spinor(homo_spinor))
2062 band_edges%CBM = min(band_edges%CBM, eigenval_spinor(homo_spinor + 1))
2063 band_edges%DBG = min(band_edges%DBG, eigenval_spinor(homo_spinor + 1) &
2064 - eigenval_spinor(homo_spinor))
2065
2066 ! 8. spinor wavefunctions:
2067 CALL cp_cfm_set_all(cfm_spinor_wf_ikp, z_zero)
2068 CALL parallel_gemm('N', 'N', 2*n_ao, 2*n_ao, 2*n_ao, z_one, &
2069 cfm_mos_ikp_spinor, cfm_eigenvec_ikp_spinor, &
2070 z_zero, cfm_spinor_wf_ikp)
2071
2072 CALL cp_cfm_release(cfm_ks_ikp_spinor)
2073 CALL cp_cfm_release(cfm_soc_ikp_spinor)
2074 CALL cp_cfm_release(cfm_work_ikp_spinor)
2075 CALL cp_cfm_release(cfm_eigenvec_ikp_spinor)
2076 CALL cp_cfm_release(cfm_mos_ikp_spinor)
2077
2078 CALL timestop(handle)
2079
2080 END SUBROUTINE soc_ev_all_mos
2081
2082! **************************************************************************************************
2083!> \brief SOC second-variational diagonalization truncated to a window of SOC_NBANDS_OCC occupied
2084!> and SOC_NBANDS_VIRT virtual states around the Fermi level -- the untruncated, full
2085!> 2*n_ao spinor-space path lives in SOC_ev_all_mos.
2086!> \param bs_env ...
2087!> \param qs_env ...
2088!> \param ikp ...
2089!> \param eigenval_no_SOC non-SOC band energies (n_ao, n_kpoints, n_spin)
2090!> \param cfm_mos_ikp ...
2091!> \param band_edges ...
2092!> \param eigenval_spinor ...
2093!> \param cfm_spinor_wf_ikp ...
2094! **************************************************************************************************
2095 SUBROUTINE soc_ev_truncate_mos(bs_env, qs_env, ikp, eigenval_no_SOC, cfm_mos_ikp, &
2096 band_edges, eigenval_spinor, cfm_spinor_wf_ikp)
2097
2098 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2099 TYPE(qs_environment_type), POINTER :: qs_env
2100 INTEGER :: ikp
2101 REAL(kind=dp), DIMENSION(:, :, :) :: eigenval_no_soc
2102 TYPE(cp_cfm_type), DIMENSION(2) :: cfm_mos_ikp
2103 TYPE(band_edges_type) :: band_edges
2104 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_spinor
2105 TYPE(cp_cfm_type) :: cfm_spinor_wf_ikp
2106
2107 CHARACTER(LEN=*), PARAMETER :: routinen = 'SOC_ev_truncate_mos'
2108
2109 INTEGER :: handle, homo_spinor, homo_used, &
2110 mo_index_start_1, mo_index_start_2, &
2111 n_ao, n_window, n_window2tot, nkind, &
2112 nocc, nvirt
2113 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eig_win_sorted, eigenval_spinor_no_soc
2114 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: proj_mo_on_kind_spinor
2115 TYPE(cp_cfm_type) :: cfm_eigenvec_ikp_spinor, &
2116 cfm_ks_ikp_spinor, cfm_mos_ikp_spinor, &
2117 cfm_soc_ikp_spinor, cfm_work_ikp_spinor
2118 TYPE(cp_fm_struct_type), POINTER :: fst_mos, fst_wa, fst_ww
2119 TYPE(cp_fm_type) :: fm_src_im, fm_src_re, fm_tgt_im, &
2120 fm_tgt_re
2121
2122 CALL timeset(routinen, handle)
2123
2124 n_ao = bs_env%n_ao
2125 homo_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin)
2126 CALL get_qs_env(qs_env, nkind=nkind)
2127
2128 nocc = bs_env%soc_nbands_occ
2129 nvirt = bs_env%soc_nbands_virt
2130 cpassert(nocc <= bs_env%n_occ_nosmear(1))
2131 cpassert(nocc <= bs_env%n_occ_nosmear(bs_env%n_spin))
2132 cpassert(nvirt <= n_ao - bs_env%n_occ_nosmear(1))
2133 cpassert(nvirt <= n_ao - bs_env%n_occ_nosmear(bs_env%n_spin))
2134 cpassert(.NOT. bs_env%do_dos_pdos)
2135 mo_index_start_1 = bs_env%n_occ_nosmear(1) - nocc + 1
2136 mo_index_start_2 = bs_env%n_occ_nosmear(bs_env%n_spin) - nocc + 1
2137 n_window = nocc + nvirt
2138 n_window2tot = 2*n_window
2139 homo_used = 2*nocc
2140
2141 NULLIFY (fst_ww, fst_wa, fst_mos)
2142 CALL cp_fm_struct_create(fst_ww, nrow_global=n_window2tot, ncol_global=n_window2tot, &
2143 template_fmstruct=bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
2144 CALL cp_fm_struct_create(fst_wa, nrow_global=n_window2tot, ncol_global=2*n_ao, &
2145 template_fmstruct=bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
2146 CALL cp_fm_struct_create(fst_mos, nrow_global=2*n_ao, ncol_global=n_window2tot, &
2147 template_fmstruct=bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
2148 CALL cp_cfm_create(cfm_ks_ikp_spinor, fst_ww)
2149 CALL cp_cfm_create(cfm_eigenvec_ikp_spinor, fst_ww)
2150 CALL cp_cfm_create(cfm_work_ikp_spinor, fst_wa)
2151 CALL cp_cfm_create(cfm_soc_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
2152 CALL cp_cfm_create(cfm_mos_ikp_spinor, fst_mos)
2153
2154 ALLOCATE (eigenval_spinor_no_soc(2*n_ao))
2155 ALLOCATE (proj_mo_on_kind_spinor(2*n_ao, nkind))
2156 ! PDOS not yet implemented -> projection is just zero -> PDOS is zero
2157 proj_mo_on_kind_spinor(:, :) = 0.0_dp
2158
2159 ! 1. get V^SOC_µν,σσ'(k_i)
2160 SELECT CASE (bs_env%gw_implementation)
2161 CASE (tensor_large_cell_gamma, ri_rs_large_cell_gamma, ri_rs_non_periodic, tensor_static_cohsex)
2162
2163 ! 1. get V^SOC_µν,σσ'(k_i) from V^SOC_µν,σσ'(k=0)
2164 CALL cfm_ikp_from_cfm_spinor_gamma(cfm_soc_ikp_spinor, &
2165 bs_env%cfm_SOC_spinor_ao(1), &
2166 bs_env%fm_s_Gamma%matrix_struct, &
2167 ikp, qs_env, bs_env%kpoints_DOS, "ORB")
2168
2169 CASE (tensor_small_cell_full_kp)
2170
2171 ! 1. V^SOC_µν,σσ'(k_i) already there
2172 CALL cp_cfm_to_cfm(bs_env%cfm_SOC_spinor_ao(ikp), cfm_soc_ikp_spinor)
2173
2174 END SELECT
2175
2176 ! 2. V^SOC_nn',σσ'(k_i) = sum_µν C^*_µn,σ(k_i) V^SOC_µν,σσ'(k_i) C_νn'(k_i),
2177 ! C_µn,σ(k_i): MO coefficiencts from diagonalizing KS-matrix h^KS_nn',σσ'(k_i)
2178
2179 ! 2.1 build matrix C_µn,σ(k_i)
2180 CALL cp_cfm_set_all(cfm_mos_ikp_spinor, z_zero)
2181 ! get_cfm_submat/add_cfm_submat only copy SQUARE blocks, so the rectangular
2182 ! (n_ao x n_window) column windows are copied with cp_fm_to_fm_submat (separate
2183 ! nrow/ncol) on the real and imaginary parts separately, then recombined.
2184 CALL cp_fm_create(fm_src_re, cfm_mos_ikp(1)%matrix_struct)
2185 CALL cp_fm_create(fm_src_im, cfm_mos_ikp(1)%matrix_struct)
2186 CALL cp_fm_create(fm_tgt_re, fst_mos)
2187 CALL cp_fm_create(fm_tgt_im, fst_mos)
2188 CALL cp_fm_set_all(fm_tgt_re, 0.0_dp)
2189 CALL cp_fm_set_all(fm_tgt_im, 0.0_dp)
2190 ! spin 1: MO columns mo_index_start_1..mo_index_start_1+n_window-1 -> rows 1..n_ao, cols 1..n_window
2191 CALL cp_cfm_to_fm(cfm_mos_ikp(1), fm_src_re, fm_src_im)
2192 CALL cp_fm_to_fm_submat(msource=fm_src_re, mtarget=fm_tgt_re, nrow=n_ao, ncol=n_window, &
2193 s_firstrow=1, s_firstcol=mo_index_start_1, t_firstrow=1, t_firstcol=1)
2194 CALL cp_fm_to_fm_submat(msource=fm_src_im, mtarget=fm_tgt_im, nrow=n_ao, ncol=n_window, &
2195 s_firstrow=1, s_firstcol=mo_index_start_1, t_firstrow=1, t_firstcol=1)
2196
2197 ! spin 2: MO columns mo_index_start_2..mo_index_start_2+n_window-1 -> rows n_ao+1..2*n_ao, cols n_window+1..2*n_window
2198 CALL cp_cfm_to_fm(cfm_mos_ikp(bs_env%n_spin), fm_src_re, fm_src_im)
2199 CALL cp_fm_to_fm_submat(msource=fm_src_re, mtarget=fm_tgt_re, nrow=n_ao, ncol=n_window, &
2200 s_firstrow=1, s_firstcol=mo_index_start_2, t_firstrow=n_ao + 1, t_firstcol=n_window + 1)
2201 CALL cp_fm_to_fm_submat(msource=fm_src_im, mtarget=fm_tgt_im, nrow=n_ao, ncol=n_window, &
2202 s_firstrow=1, s_firstcol=mo_index_start_2, t_firstrow=n_ao + 1, t_firstcol=n_window + 1)
2203
2204 CALL cp_fm_to_cfm(fm_tgt_re, fm_tgt_im, cfm_mos_ikp_spinor)
2205
2206 CALL cp_fm_release(fm_src_re)
2207 CALL cp_fm_release(fm_src_im)
2208 CALL cp_fm_release(fm_tgt_re)
2209 CALL cp_fm_release(fm_tgt_im)
2210
2211 ! 2.2 work_nν,σσ' = sum_µ C^*_µn,σ(k_i) V^SOC_µν,σσ'(k_i)
2212 CALL parallel_gemm('C', 'N', n_window2tot, 2*n_ao, 2*n_ao, z_one, &
2213 cfm_mos_ikp_spinor, cfm_soc_ikp_spinor, &
2214 z_zero, cfm_work_ikp_spinor)
2215 ! 2.3 V^SOC_nn',σσ'(k_i) = sum_ν work_nν,σσ' C_νn'(k_i)
2216 CALL parallel_gemm('N', 'N', n_window2tot, n_window2tot, 2*n_ao, z_one, &
2217 cfm_work_ikp_spinor, cfm_mos_ikp_spinor, &
2218 z_zero, cfm_ks_ikp_spinor)
2219
2220 ! 3. remove SOC outside of energy window: not applicable here -- the truncation window
2221 ! itself already restricts SOC to the states around the Fermi level
2222 eigenval_spinor_no_soc(1:n_ao) = eigenval_no_soc(1:n_ao, ikp, 1)
2223 eigenval_spinor_no_soc(n_ao + 1:) = eigenval_no_soc(1:n_ao, ikp, bs_env%n_spin)
2224
2225 ! 4. h^G0W0+SOC_nn',σσ'(k_i) = ε_nσ^G0W0(k_i) δ_nn' δ_σσ' + V^SOC_nn',σσ'(k_i)
2226 CALL cfm_add_on_diag(cfm_ks_ikp_spinor, &
2227 [eigenval_no_soc(mo_index_start_1:mo_index_start_1 + n_window - 1, ikp, 1), &
2228 eigenval_no_soc(mo_index_start_2:mo_index_start_2 + n_window - 1, ikp, bs_env%n_spin)])
2229
2230 ! 5. diagonalize h^G0W0+SOC_nn',σσ'(k_i) to get eigenvalues
2231 ALLOCATE (eig_win_sorted(n_window2tot))
2232 CALL cp_cfm_heevd(cfm_ks_ikp_spinor, cfm_eigenvec_ikp_spinor, eig_win_sorted)
2233
2234 ! 6. DOS from spinors: not supported for the truncated window (CPASSERT'd above)
2235
2236 ! 7. valence band max. (VBM), conduction band min. (CBM) and direct bandgap (DBG)
2237 band_edges%VBM = max(band_edges%VBM, eig_win_sorted(homo_used))
2238 band_edges%CBM = min(band_edges%CBM, eig_win_sorted(homo_used + 1))
2239 band_edges%DBG = min(band_edges%DBG, eig_win_sorted(homo_used + 1) - eig_win_sorted(homo_used))
2240
2241 CALL fill_truncated_eigenval_spinor(eigenval_no_soc, eig_win_sorted, ikp, bs_env%n_spin, &
2242 n_ao, mo_index_start_1, mo_index_start_2, n_window, &
2243 n_window2tot, eigenval_spinor)
2244
2245 ! 8. spinor wavefunctions:
2246 CALL cp_cfm_set_all(cfm_spinor_wf_ikp, z_zero)
2247 cpassert(.NOT. bs_env%do_ldos)
2248 CALL parallel_gemm('N', 'N', 2*n_ao, n_window2tot, n_window2tot, z_one, &
2249 cfm_mos_ikp_spinor, cfm_eigenvec_ikp_spinor, &
2250 z_zero, cfm_spinor_wf_ikp)
2251
2252 CALL cp_cfm_release(cfm_ks_ikp_spinor)
2253 CALL cp_cfm_release(cfm_soc_ikp_spinor)
2254 CALL cp_cfm_release(cfm_work_ikp_spinor)
2255 CALL cp_cfm_release(cfm_eigenvec_ikp_spinor)
2256 CALL cp_cfm_release(cfm_mos_ikp_spinor)
2257
2258 CALL cp_fm_struct_release(fst_ww)
2259 CALL cp_fm_struct_release(fst_wa)
2260 CALL cp_fm_struct_release(fst_mos)
2261
2262 CALL timestop(handle)
2263
2264 END SUBROUTINE soc_ev_truncate_mos
2265! **************************************************************************************************
2266!> \brief ...
2267!> \param eigenval_no_SOC non-SOC band energies (n_ao, n_kpoints, n_spin)
2268!> \param eig_win_sorted SOC-corrected eigenvalues of the window, ascending
2269!> \param ikp k-point index into eigenval_no_SOC
2270!> \param n_spin number of spin channels
2271!> \param n_ao number of AO basis functions (per spin)
2272!> \param mo_index_start_1 first retained state index for spin 1
2273!> \param mo_index_start_2 first retained state index for spin n_spin
2274!> \param n_window per-spin window size
2275!> \param n_window2tot total window size
2276!> \param eigenval_spinor output: full repacked spinor spectrum (2*n_ao)
2277! **************************************************************************************************
2278 SUBROUTINE fill_truncated_eigenval_spinor(eigenval_no_SOC, eig_win_sorted, ikp, n_spin, n_ao, &
2279 mo_index_start_1, mo_index_start_2, n_window, &
2280 n_window2tot, eigenval_spinor)
2281 REAL(kind=dp), DIMENSION(:, :, :) :: eigenval_no_soc
2282 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eig_win_sorted
2283 INTEGER :: ikp, n_spin, n_ao, mo_index_start_1, &
2284 mo_index_start_2, n_window, &
2285 n_window2tot
2286 REAL(kind=dp), DIMENSION(:) :: eigenval_spinor
2287
2288 CHARACTER(LEN=*), PARAMETER :: routinen = 'fill_truncated_eigenval_spinor'
2289
2290 INTEGER :: handle, i, ia, ib, idx_end, idx_start, &
2291 mo_index_start_1_m1, &
2292 mo_index_start_2_m1
2293
2294 CALL timeset(routinen, handle)
2295
2296 ! verify the window was wide enough that SOC did not push a retained boundary state
2297 ! past an excluded neighbor on either side -- this would invalidate the merge below
2298 IF (mo_index_start_1 > 1) THEN
2299 IF (eig_win_sorted(1) < eigenval_no_soc(mo_index_start_1 - 1, ikp, 1)) THEN
2300 cpabort("SOC_NBANDS_OCC too small.")
2301 END IF
2302 END IF
2303 IF (mo_index_start_2 > 1) THEN
2304 IF (eig_win_sorted(1) < eigenval_no_soc(mo_index_start_2 - 1, ikp, n_spin)) THEN
2305 cpabort("SOC_NBANDS_OCC too small.")
2306 END IF
2307 END IF
2308 IF (mo_index_start_1 + n_window <= n_ao) THEN
2309 IF (eig_win_sorted(n_window2tot) > eigenval_no_soc(mo_index_start_1 + n_window, ikp, 1)) THEN
2310 cpabort("SOC_NBANDS_VIRT too small.")
2311 END IF
2312 END IF
2313 IF (mo_index_start_2 + n_window <= n_ao) THEN
2314 IF (eig_win_sorted(n_window2tot) > eigenval_no_soc(mo_index_start_2 + n_window, ikp, n_spin)) THEN
2315 cpabort("SOC_NBANDS_VIRT too small.")
2316 END IF
2317 END IF
2318
2319 mo_index_start_1_m1 = mo_index_start_1 - 1
2320 mo_index_start_2_m1 = mo_index_start_2 - 1
2321
2322 ia = 1
2323 ib = 1
2324 DO i = 1, mo_index_start_1_m1 + mo_index_start_2_m1
2325 IF (ia <= mo_index_start_1_m1 .AND. ib <= mo_index_start_2_m1) THEN
2326 IF (eigenval_no_soc(ia, ikp, 1) <= eigenval_no_soc(ib, ikp, n_spin)) THEN
2327 eigenval_spinor(i) = eigenval_no_soc(ia, ikp, 1)
2328 ia = ia + 1
2329 ELSE
2330 eigenval_spinor(i) = eigenval_no_soc(ib, ikp, n_spin)
2331 ib = ib + 1
2332 END IF
2333 ELSE IF (ia <= mo_index_start_1_m1) THEN
2334 eigenval_spinor(i) = eigenval_no_soc(ia, ikp, 1)
2335 ia = ia + 1
2336 ELSE
2337 eigenval_spinor(i) = eigenval_no_soc(ib, ikp, n_spin)
2338 ib = ib + 1
2339 END IF
2340 END DO
2341
2342 idx_start = mo_index_start_1_m1 + mo_index_start_2_m1 + 1
2343 idx_end = mo_index_start_1_m1 + mo_index_start_2_m1 + n_window2tot
2344
2345 eigenval_spinor(idx_start:idx_end) = eig_win_sorted(:)
2346
2347 ia = mo_index_start_1 + n_window
2348 ib = mo_index_start_2 + n_window
2349 DO i = mo_index_start_1_m1 + mo_index_start_2_m1 + n_window2tot + 1, 2*n_ao
2350 IF (ia <= n_ao .AND. ib <= n_ao) THEN
2351 IF (eigenval_no_soc(ia, ikp, 1) <= eigenval_no_soc(ib, ikp, n_spin)) THEN
2352 eigenval_spinor(i) = eigenval_no_soc(ia, ikp, 1)
2353 ia = ia + 1
2354 ELSE
2355 eigenval_spinor(i) = eigenval_no_soc(ib, ikp, n_spin)
2356 ib = ib + 1
2357 END IF
2358 ELSE IF (ia <= n_ao) THEN
2359 eigenval_spinor(i) = eigenval_no_soc(ia, ikp, 1)
2360 ia = ia + 1
2361 ELSE
2362 eigenval_spinor(i) = eigenval_no_soc(ib, ikp, n_spin)
2363 ib = ib + 1
2364 END IF
2365 END DO
2366
2367 DEALLOCATE (eig_win_sorted)
2368
2369 CALL timestop(handle)
2370
2371 END SUBROUTINE fill_truncated_eigenval_spinor
2372! **************************************************************************************************
2373!> \brief ...
2374!> \param DOS ...
2375!> \param PDOS ...
2376!> \param eigenval ...
2377!> \param ikp ...
2378!> \param bs_env ...
2379!> \param n_E ...
2380!> \param E_min ...
2381!> \param proj_mo_on_kind ...
2382! **************************************************************************************************
2383 SUBROUTINE add_to_dos_pdos(DOS, PDOS, eigenval, ikp, bs_env, n_E, E_min, proj_mo_on_kind)
2384
2385 REAL(kind=dp), DIMENSION(:) :: dos
2386 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: pdos
2387 REAL(kind=dp), DIMENSION(:) :: eigenval
2388 INTEGER :: ikp
2389 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2390 INTEGER :: n_e
2391 REAL(kind=dp) :: e_min
2392 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: proj_mo_on_kind
2393
2394 CHARACTER(LEN=*), PARAMETER :: routinen = 'add_to_DOS_PDOS'
2395
2396 INTEGER :: handle, i_e, i_kind, i_mo, n_mo, nkind
2397 REAL(kind=dp) :: broadening, energy, energy_step_dos, wkp
2398
2399 CALL timeset(routinen, handle)
2400
2401 energy_step_dos = bs_env%energy_step_DOS
2402 broadening = bs_env%broadening_DOS
2403
2404 n_mo = SIZE(eigenval)
2405 nkind = SIZE(proj_mo_on_kind, 2)
2406
2407 ! normalize to closed-shell / open-shell
2408 wkp = bs_env%kpoints_DOS%wkp(ikp)*bs_env%spin_degeneracy
2409 DO i_e = 1, n_e
2410 energy = e_min + i_e*energy_step_dos
2411 DO i_mo = 1, n_mo
2412 ! DOS
2413 dos(i_e) = dos(i_e) + wkp*gaussian(energy - eigenval(i_mo), broadening)
2414
2415 ! PDOS
2416 DO i_kind = 1, nkind
2417 IF (proj_mo_on_kind(i_mo, i_kind) > 0.0_dp) THEN
2418 pdos(i_e, i_kind) = pdos(i_e, i_kind) + &
2419 proj_mo_on_kind(i_mo, i_kind)*wkp* &
2420 gaussian(energy - eigenval(i_mo), broadening)
2421 END IF
2422 END DO
2423 END DO
2424 END DO
2425
2426 CALL timestop(handle)
2427
2428 END SUBROUTINE add_to_dos_pdos
2429
2430! **************************************************************************************************
2431!> \brief ...
2432!> \param LDOS_2d ...
2433!> \param qs_env ...
2434!> \param ikp ...
2435!> \param bs_env ...
2436!> \param cfm_mos_ikp ...
2437!> \param eigenval ...
2438!> \param band_edges ...
2439!> \param do_spinor ...
2440!> \param cfm_non_spinor ...
2441! **************************************************************************************************
2442 SUBROUTINE add_to_ldos_2d(LDOS_2d, qs_env, ikp, bs_env, cfm_mos_ikp, eigenval, &
2443 band_edges, do_spinor, cfm_non_spinor)
2444 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ldos_2d
2445 TYPE(qs_environment_type), POINTER :: qs_env
2446 INTEGER :: ikp
2447 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2448 TYPE(cp_cfm_type) :: cfm_mos_ikp
2449 REAL(kind=dp), DIMENSION(:) :: eigenval
2450 TYPE(band_edges_type) :: band_edges
2451 LOGICAL, OPTIONAL :: do_spinor
2452 TYPE(cp_cfm_type), OPTIONAL :: cfm_non_spinor
2453
2454 CHARACTER(LEN=*), PARAMETER :: routinen = 'add_to_LDOS_2d'
2455
2456 INTEGER :: handle, i_e, i_x_end, i_x_start, i_y_end, i_y_start, i_z, i_z_end, i_z_start, &
2457 j_col, j_mo, n_e, n_mo, n_z, ncol_local, nimages, z_end_global, z_start_global
2458 INTEGER, DIMENSION(:), POINTER :: col_indices
2459 LOGICAL :: is_any_weight_non_zero, my_do_spinor
2460 REAL(kind=dp) :: broadening, e_max, e_min, &
2461 e_total_window, energy, energy_step, &
2462 energy_window, spin_degeneracy, weight
2463 TYPE(cp_cfm_type) :: cfm_weighted_dm_ikp, cfm_work
2464 TYPE(cp_fm_type) :: fm_non_spinor, fm_weighted_dm_mic
2465 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: weighted_dm_mic
2466 TYPE(dft_control_type), POINTER :: dft_control
2467 TYPE(pw_c1d_gs_type) :: rho_g
2468 TYPE(pw_env_type), POINTER :: pw_env
2469 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2470 TYPE(pw_r3d_rs_type) :: ldos_3d
2471 TYPE(qs_ks_env_type), POINTER :: ks_env
2472
2473 CALL timeset(routinen, handle)
2474
2475 my_do_spinor = .false.
2476 IF (PRESENT(do_spinor)) my_do_spinor = do_spinor
2477
2478 CALL get_qs_env(qs_env, ks_env=ks_env, pw_env=pw_env, dft_control=dft_control)
2479
2480 ! previously, dft_control%nimages set to # neighbor cells, revert for Γ-only KS matrix
2481 nimages = dft_control%nimages
2482 dft_control%nimages = bs_env%nimages_scf
2483
2484 energy_window = bs_env%energy_window_DOS
2485 energy_step = bs_env%energy_step_DOS
2486 broadening = bs_env%broadening_DOS
2487
2488 e_min = band_edges%VBM - 0.5_dp*energy_window
2489 e_max = band_edges%CBM + 0.5_dp*energy_window
2490 e_total_window = e_max - e_min
2491
2492 n_e = int(e_total_window/energy_step)
2493
2494 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
2495
2496 CALL auxbas_pw_pool%create_pw(ldos_3d)
2497 CALL auxbas_pw_pool%create_pw(rho_g)
2498
2499 i_x_start = lbound(ldos_3d%array, 1)
2500 i_x_end = ubound(ldos_3d%array, 1)
2501 i_y_start = lbound(ldos_3d%array, 2)
2502 i_y_end = ubound(ldos_3d%array, 2)
2503 i_z_start = lbound(ldos_3d%array, 3)
2504 i_z_end = ubound(ldos_3d%array, 3)
2505
2506 z_start_global = i_z_start
2507 z_end_global = i_z_end
2508
2509 CALL bs_env%para_env%min(z_start_global)
2510 CALL bs_env%para_env%max(z_end_global)
2511 n_z = z_end_global - z_start_global + 1
2512
2513 IF (any(abs(bs_env%hmat(1:2, 3)) > 1.0e-6_dp) .OR. any(abs(bs_env%hmat(3, 1:2)) > 1.0e-6_dp)) THEN
2514 cpabort("Please choose a cell that has 90° angles to the z-direction.")
2515 END IF
2516 ! for integration, we need the dz and the conversion from H -> eV and a_Bohr -> Å
2517 bs_env%unit_ldos_int_z_inv_Ang2_eV = bs_env%hmat(3, 3)/real(n_z, kind=dp)/evolt/angstrom**2
2518
2519 IF (ikp == 1) THEN
2520 ALLOCATE (ldos_2d(i_x_start:i_x_end, i_y_start:i_y_end, n_e))
2521 ldos_2d(:, :, :) = 0.0_dp
2522 END IF
2523
2524 CALL cp_cfm_create(cfm_work, cfm_mos_ikp%matrix_struct)
2525 CALL cp_cfm_create(cfm_weighted_dm_ikp, cfm_mos_ikp%matrix_struct)
2526 CALL cp_fm_create(fm_weighted_dm_mic, cfm_mos_ikp%matrix_struct)
2527 IF (my_do_spinor) THEN
2528 CALL cp_fm_create(fm_non_spinor, cfm_non_spinor%matrix_struct)
2529 END IF
2530
2531 CALL cp_cfm_get_info(matrix=cfm_mos_ikp, &
2532 ncol_global=n_mo, &
2533 ncol_local=ncol_local, &
2534 col_indices=col_indices)
2535
2536 NULLIFY (weighted_dm_mic)
2537 CALL dbcsr_allocate_matrix_set(weighted_dm_mic, 1)
2538 ALLOCATE (weighted_dm_mic(1)%matrix)
2539 CALL dbcsr_create(weighted_dm_mic(1)%matrix, template=bs_env%mat_ao_ao%matrix, &
2540 matrix_type=dbcsr_type_symmetric)
2541
2542 DO i_e = 1, n_e
2543
2544 energy = e_min + i_e*energy_step
2545
2546 is_any_weight_non_zero = .false.
2547
2548 DO j_col = 1, ncol_local
2549
2550 j_mo = col_indices(j_col)
2551
2552 IF (my_do_spinor) THEN
2553 spin_degeneracy = 1.0_dp
2554 ELSE
2555 spin_degeneracy = bs_env%spin_degeneracy
2556 END IF
2557
2558 weight = gaussian(energy - eigenval(j_mo), broadening)*spin_degeneracy
2559
2560 cfm_work%local_data(:, j_col) = cfm_mos_ikp%local_data(:, j_col)*weight
2561
2562 IF (weight > 1.0e-5_dp) is_any_weight_non_zero = .true.
2563
2564 END DO
2565
2566 CALL bs_env%para_env%sync()
2567 CALL bs_env%para_env%sum(is_any_weight_non_zero)
2568 CALL bs_env%para_env%sync()
2569
2570 ! cycle if there are no states at the energy i_E
2571 IF (is_any_weight_non_zero) THEN
2572
2573 CALL parallel_gemm('N', 'C', n_mo, n_mo, n_mo, z_one, &
2574 cfm_mos_ikp, cfm_work, z_zero, cfm_weighted_dm_ikp)
2575
2576 IF (my_do_spinor) THEN
2577
2578 ! contribution from up,up to fm_non_spinor
2579 CALL get_cfm_submat(cfm_non_spinor, cfm_weighted_dm_ikp, 1, 1)
2580 CALL cp_fm_set_all(fm_non_spinor, 0.0_dp)
2581 CALL mic_contribution_from_ikp(bs_env, qs_env, fm_non_spinor, &
2582 cfm_non_spinor, ikp, bs_env%kpoints_DOS, &
2583 "ORB", bs_env%kpoints_DOS%wkp(ikp))
2584
2585 ! add contribution from down,down to fm_non_spinor
2586 CALL get_cfm_submat(cfm_non_spinor, cfm_weighted_dm_ikp, n_mo/2, n_mo/2)
2587 CALL mic_contribution_from_ikp(bs_env, qs_env, fm_non_spinor, &
2588 cfm_non_spinor, ikp, bs_env%kpoints_DOS, &
2589 "ORB", bs_env%kpoints_DOS%wkp(ikp))
2590 CALL copy_fm_to_dbcsr(fm_non_spinor, weighted_dm_mic(1)%matrix, &
2591 keep_sparsity=.false.)
2592 ELSE
2593 CALL cp_fm_set_all(fm_weighted_dm_mic, 0.0_dp)
2594 CALL mic_contribution_from_ikp(bs_env, qs_env, fm_weighted_dm_mic, &
2595 cfm_weighted_dm_ikp, ikp, bs_env%kpoints_DOS, &
2596 "ORB", bs_env%kpoints_DOS%wkp(ikp))
2597 CALL copy_fm_to_dbcsr(fm_weighted_dm_mic, weighted_dm_mic(1)%matrix, &
2598 keep_sparsity=.false.)
2599 END IF
2600
2601 ldos_3d%array(:, :, :) = 0.0_dp
2602
2603 CALL calculate_rho_elec(matrix_p_kp=weighted_dm_mic, &
2604 rho=ldos_3d, &
2605 rho_gspace=rho_g, &
2606 ks_env=ks_env)
2607
2608 DO i_z = i_z_start, i_z_end
2609 ldos_2d(:, :, i_e) = ldos_2d(:, :, i_e) + ldos_3d%array(:, :, i_z)
2610 END DO
2611
2612 END IF
2613
2614 END DO
2615
2616 ! set back nimages
2617 dft_control%nimages = nimages
2618
2619 CALL auxbas_pw_pool%give_back_pw(ldos_3d)
2620 CALL auxbas_pw_pool%give_back_pw(rho_g)
2621
2622 CALL cp_cfm_release(cfm_work)
2623 CALL cp_cfm_release(cfm_weighted_dm_ikp)
2624
2625 CALL cp_fm_release(fm_weighted_dm_mic)
2626
2627 CALL dbcsr_deallocate_matrix_set(weighted_dm_mic)
2628
2629 IF (my_do_spinor) THEN
2630 CALL cp_fm_release(fm_non_spinor)
2631 END IF
2632
2633 CALL timestop(handle)
2634
2635 END SUBROUTINE add_to_ldos_2d
2636
2637! **************************************************************************************************
2638!> \brief ...
2639!> \param eigenval_spinor ...
2640!> \param ikp_for_file ...
2641!> \param ikp ...
2642!> \param bs_env ...
2643!> \param eigenval_spinor_G0W0 ...
2644! **************************************************************************************************
2645 SUBROUTINE write_soc_eigenvalues(eigenval_spinor, ikp_for_file, ikp, bs_env, eigenval_spinor_G0W0)
2646
2647 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_spinor
2648 INTEGER :: ikp_for_file, ikp
2649 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2650 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: eigenval_spinor_g0w0
2651
2652 CHARACTER(LEN=*), PARAMETER :: routinen = 'write_SOC_eigenvalues'
2653
2654 CHARACTER(len=3) :: occ_vir
2655 CHARACTER(LEN=default_string_length) :: fname
2656 INTEGER :: handle, i_mo, iunit, n_occ_spinor
2657
2658 CALL timeset(routinen, handle)
2659
2660 fname = "bandstructure_SCF_and_G0W0_plus_SOC"
2661
2662 IF (bs_env%para_env%is_source()) THEN
2663
2664 IF (ikp_for_file == 1) THEN
2665 CALL open_file(trim(fname), unit_number=iunit, file_status="REPLACE", &
2666 file_action="WRITE")
2667 ELSE
2668 CALL open_file(trim(fname), unit_number=iunit, file_status="OLD", &
2669 file_action="WRITE", file_position="APPEND")
2670 END IF
2671
2672 WRITE (iunit, "(A)") " "
2673 WRITE (iunit, "(A10,I7,A25,3F10.4)") "kpoint: ", ikp_for_file, "coordinate: ", &
2674 bs_env%kpoints_DOS%xkp(:, ikp)
2675 WRITE (iunit, "(A)") " "
2676
2677 IF (PRESENT(eigenval_spinor_g0w0)) THEN
2678 ! SCF+SOC and G0W0+SOC eigenvalues
2679 WRITE (iunit, "(A5,A12,2A22)") "n", "k", "ϵ_nk^DFT+SOC (eV)", "ϵ_nk^G0W0+SOC (eV)"
2680 ELSE
2681 ! SCF+SOC eigenvalues only
2682 WRITE (iunit, "(A5,A12,A22)") "n", "k", "ϵ_nk^DFT+SOC (eV)"
2683 END IF
2684
2685 n_occ_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin)
2686
2687 DO i_mo = 1, SIZE(eigenval_spinor)
2688 IF (i_mo <= n_occ_spinor) occ_vir = 'occ'
2689 IF (i_mo > n_occ_spinor) occ_vir = 'vir'
2690 IF (PRESENT(eigenval_spinor_g0w0)) THEN
2691 ! SCF+SOC and G0W0+SOC eigenvalues
2692 WRITE (iunit, "(I5,3A,I5,4F16.3,2F17.3)") i_mo, ' (', occ_vir, ') ', &
2693 ikp_for_file, eigenval_spinor(i_mo)*evolt, eigenval_spinor_g0w0(i_mo)*evolt
2694 ELSE
2695 ! SCF+SOC eigenvalues only
2696 WRITE (iunit, "(I5,3A,I5,4F16.3,F17.3)") i_mo, ' (', occ_vir, ') ', &
2697 ikp_for_file, eigenval_spinor(i_mo)*evolt
2698 END IF
2699 END DO
2700
2701 CALL close_file(iunit)
2702
2703 END IF
2704
2705 CALL timestop(handle)
2706
2707 END SUBROUTINE write_soc_eigenvalues
2708
2709! **************************************************************************************************
2710!> \brief ...
2711!> \param int_number ...
2712!> \return ...
2713! **************************************************************************************************
2714 PURE FUNCTION count_digits(int_number)
2715
2716 INTEGER, INTENT(IN) :: int_number
2717 INTEGER :: count_digits
2718
2719 INTEGER :: digitcount, tempint
2720
2721 digitcount = 0
2722
2723 tempint = int_number
2724
2725 DO WHILE (tempint /= 0)
2726 tempint = tempint/10
2727 digitcount = digitcount + 1
2728 END DO
2729
2730 count_digits = digitcount
2731
2732 END FUNCTION count_digits
2733
2734! **************************************************************************************************
2735!> \brief ...
2736!> \param band_edges ...
2737!> \param scf_gw_soc ...
2738!> \param bs_env ...
2739! **************************************************************************************************
2740 SUBROUTINE write_band_edges(band_edges, scf_gw_soc, bs_env)
2741
2742 TYPE(band_edges_type) :: band_edges
2743 CHARACTER(LEN=*) :: scf_gw_soc
2744 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2745
2746 CHARACTER(LEN=*), PARAMETER :: routinen = 'write_band_edges'
2747
2748 CHARACTER(LEN=17) :: print_format
2749 INTEGER :: handle, u
2750
2751 CALL timeset(routinen, handle)
2752
2753 ! print format
2754 print_format = "(T2,2A,T61,F20.3)"
2755
2756 u = bs_env%unit_nr
2757 IF (u > 0) THEN
2758 WRITE (u, '(T2,A)') ''
2759 WRITE (u, print_format) scf_gw_soc, ' valence band maximum (eV):', band_edges%VBM*evolt
2760 WRITE (u, print_format) scf_gw_soc, ' conduction band minimum (eV):', band_edges%CBM*evolt
2761 WRITE (u, print_format) scf_gw_soc, ' indirect band gap (eV):', band_edges%IDBG*evolt
2762 WRITE (u, print_format) scf_gw_soc, ' direct band gap (eV):', band_edges%DBG*evolt
2763 END IF
2764
2765 CALL timestop(handle)
2766
2767 END SUBROUTINE write_band_edges
2768
2769! **************************************************************************************************
2770!> \brief ...
2771!> \param DOS ...
2772!> \param PDOS ...
2773!> \param bs_env ...
2774!> \param qs_env ...
2775!> \param scf_gw_soc ...
2776!> \param E_min ...
2777!> \param E_VBM ...
2778! **************************************************************************************************
2779 SUBROUTINE write_dos_pdos(DOS, PDOS, bs_env, qs_env, scf_gw_soc, E_min, E_VBM)
2780 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dos
2781 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: pdos
2782 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2783 TYPE(qs_environment_type), POINTER :: qs_env
2784 CHARACTER(LEN=*) :: scf_gw_soc
2785 REAL(kind=dp) :: e_min, e_vbm
2786
2787 CHARACTER(LEN=*), PARAMETER :: routinen = 'write_dos_pdos'
2788
2789 CHARACTER(LEN=3), DIMENSION(100) :: elements
2790 CHARACTER(LEN=default_string_length) :: atom_name, fname, output_string
2791 INTEGER :: handle, i_e, i_kind, iatom, iunit, n_a, &
2792 n_e, nkind
2793 REAL(kind=dp) :: energy
2794 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2795
2796 CALL timeset(routinen, handle)
2797
2798 WRITE (fname, "(3A)") "DOS_PDOS_", scf_gw_soc, ".out"
2799
2800 n_e = SIZE(pdos, 1)
2801 nkind = SIZE(pdos, 2)
2802 CALL get_qs_env(qs_env, particle_set=particle_set)
2803
2804 IF (bs_env%para_env%is_source()) THEN
2805
2806 CALL open_file(trim(fname), unit_number=iunit, file_status="REPLACE", file_action="WRITE")
2807
2808 n_a = 2 + nkind
2809
2810 DO iatom = 1, bs_env%n_atom
2811 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
2812 kind_number=i_kind, name=atom_name)
2813 elements(i_kind) = atom_name(1:3)
2814 END DO
2815
2816 WRITE (output_string, "(A,I1,A)") "(", n_a, "A)"
2817
2818 WRITE (iunit, trim(output_string)) "Energy-E_F (eV) DOS (1/eV) PDOS (1/eV) ", &
2819 " of atom type ", elements(1:nkind)
2820
2821 WRITE (output_string, "(A,I1,A)") "(", n_a, "F13.5)"
2822
2823 DO i_e = 1, n_e
2824 ! energy is relative to valence band maximum => - E_VBM
2825 energy = e_min + i_e*bs_env%energy_step_DOS - e_vbm
2826 WRITE (iunit, trim(output_string)) energy*evolt, dos(i_e)/evolt, pdos(i_e, :)/evolt
2827 END DO
2828
2829 CALL close_file(iunit)
2830
2831 END IF
2832
2833 CALL timestop(handle)
2834
2835 END SUBROUTINE write_dos_pdos
2836
2837! **************************************************************************************************
2838!> \brief ...
2839!> \param energy ...
2840!> \param broadening ...
2841!> \return ...
2842! **************************************************************************************************
2843 PURE FUNCTION gaussian(energy, broadening)
2844
2845 REAL(kind=dp), INTENT(IN) :: energy, broadening
2846 REAL(kind=dp) :: gaussian
2847
2848 IF (abs(energy) < 5*broadening) THEN
2849 gaussian = 1.0_dp/broadening/sqrt(twopi)*exp(-0.5_dp*energy**2/broadening**2)
2850 ELSE
2851 gaussian = 0.0_dp
2852 END IF
2853
2854 END FUNCTION gaussian
2855
2856! **************************************************************************************************
2857!> \brief ...
2858!> \param proj_mo_on_kind ...
2859!> \param qs_env ...
2860!> \param cfm_mos ...
2861!> \param cfm_s ...
2862! **************************************************************************************************
2863 SUBROUTINE compute_proj_mo_on_kind(proj_mo_on_kind, qs_env, cfm_mos, cfm_s)
2864 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: proj_mo_on_kind
2865 TYPE(qs_environment_type), POINTER :: qs_env
2866 TYPE(cp_cfm_type) :: cfm_mos, cfm_s
2867
2868 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_proj_mo_on_kind'
2869
2870 INTEGER :: handle, i_atom, i_global, i_kind, i_row, &
2871 j_col, n_ao, n_mo, ncol_local, nkind, &
2872 nrow_local
2873 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_bf, kind_of
2874 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2875 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2876 TYPE(cp_cfm_type) :: cfm_proj, cfm_s_i_kind, cfm_work
2877 TYPE(cp_fm_type) :: fm_proj_im, fm_proj_re
2878
2879 CALL timeset(routinen, handle)
2880
2881 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, nkind=nkind)
2882 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
2883
2884 CALL cp_cfm_get_info(matrix=cfm_mos, &
2885 nrow_global=n_mo, &
2886 nrow_local=nrow_local, &
2887 ncol_local=ncol_local, &
2888 row_indices=row_indices, &
2889 col_indices=col_indices)
2890
2891 n_ao = qs_env%bs_env%n_ao
2892
2893 ALLOCATE (atom_from_bf(n_ao))
2894 CALL get_atom_index_from_basis_function_index(qs_env, atom_from_bf, n_ao, "ORB")
2895
2896 proj_mo_on_kind(:, :) = 0.0_dp
2897
2898 CALL cp_cfm_create(cfm_s_i_kind, cfm_s%matrix_struct)
2899 CALL cp_cfm_create(cfm_work, cfm_s%matrix_struct)
2900 CALL cp_cfm_create(cfm_proj, cfm_s%matrix_struct)
2901 CALL cp_fm_create(fm_proj_re, cfm_s%matrix_struct)
2902 CALL cp_fm_create(fm_proj_im, cfm_s%matrix_struct)
2903
2904 DO i_kind = 1, nkind
2905
2906 CALL cp_cfm_to_cfm(cfm_s, cfm_s_i_kind)
2907
2908 ! set entries in overlap matrix to zero which do not belong to atoms of i_kind
2909 DO j_col = 1, ncol_local
2910 DO i_row = 1, nrow_local
2911
2912 i_global = row_indices(i_row)
2913
2914 IF (i_global <= n_ao) THEN
2915 i_atom = atom_from_bf(i_global)
2916 ELSE IF (i_global <= 2*n_ao) THEN
2917 i_atom = atom_from_bf(i_global - n_ao)
2918 ELSE
2919 cpabort("Wrong indices.")
2920 END IF
2921
2922 IF (i_kind /= kind_of(i_atom)) THEN
2923 cfm_s_i_kind%local_data(i_row, j_col) = z_zero
2924 END IF
2925
2926 END DO
2927 END DO
2928
2929 CALL parallel_gemm('N', 'N', n_mo, n_mo, n_mo, z_one, &
2930 cfm_s_i_kind, cfm_mos, z_zero, cfm_work)
2931 CALL parallel_gemm('C', 'N', n_mo, n_mo, n_mo, z_one, &
2932 cfm_mos, cfm_work, z_zero, cfm_proj)
2933
2934 CALL cp_cfm_to_fm(cfm_proj, fm_proj_re, fm_proj_im)
2935
2936 CALL cp_fm_get_diag(fm_proj_im, proj_mo_on_kind(:, i_kind))
2937 CALL cp_fm_get_diag(fm_proj_re, proj_mo_on_kind(:, i_kind))
2938
2939 END DO ! i_kind
2940
2941 CALL cp_cfm_release(cfm_s_i_kind)
2942 CALL cp_cfm_release(cfm_work)
2943 CALL cp_cfm_release(cfm_proj)
2944 CALL cp_fm_release(fm_proj_re)
2945 CALL cp_fm_release(fm_proj_im)
2946
2947 CALL timestop(handle)
2948
2949 END SUBROUTINE compute_proj_mo_on_kind
2950
2951! **************************************************************************************************
2952!> \brief ...
2953!> \param cfm_spinor_ikp ...
2954!> \param cfm_spinor_Gamma ...
2955!> \param fm_struct_non_spinor ...
2956!> \param ikp ...
2957!> \param qs_env ...
2958!> \param kpoints ...
2959!> \param basis_type ...
2960! **************************************************************************************************
2961 SUBROUTINE cfm_ikp_from_cfm_spinor_gamma(cfm_spinor_ikp, cfm_spinor_Gamma, fm_struct_non_spinor, &
2962 ikp, qs_env, kpoints, basis_type)
2963 TYPE(cp_cfm_type) :: cfm_spinor_ikp, cfm_spinor_gamma
2964 TYPE(cp_fm_struct_type), POINTER :: fm_struct_non_spinor
2965 INTEGER :: ikp
2966 TYPE(qs_environment_type), POINTER :: qs_env
2967 TYPE(kpoint_type), POINTER :: kpoints
2968 CHARACTER(LEN=*) :: basis_type
2969
2970 CHARACTER(LEN=*), PARAMETER :: routinen = 'cfm_ikp_from_cfm_spinor_Gamma'
2971
2972 INTEGER :: handle, i_block, i_offset, j_block, &
2973 j_offset, n_ao
2974 TYPE(cp_cfm_type) :: cfm_non_spinor_gamma, cfm_non_spinor_ikp
2975 TYPE(cp_fm_type) :: fm_non_spinor_gamma_im, &
2976 fm_non_spinor_gamma_re
2977
2978 CALL timeset(routinen, handle)
2979
2980 CALL cp_cfm_create(cfm_non_spinor_gamma, fm_struct_non_spinor)
2981 CALL cp_cfm_create(cfm_non_spinor_ikp, fm_struct_non_spinor)
2982 CALL cp_fm_create(fm_non_spinor_gamma_re, fm_struct_non_spinor)
2983 CALL cp_fm_create(fm_non_spinor_gamma_im, fm_struct_non_spinor)
2984
2985 CALL cp_cfm_get_info(cfm_non_spinor_gamma, nrow_global=n_ao)
2986
2987 CALL cp_cfm_set_all(cfm_spinor_ikp, z_zero)
2988
2989 DO i_block = 0, 1
2990 DO j_block = 0, 1
2991 i_offset = i_block*n_ao + 1
2992 j_offset = j_block*n_ao + 1
2993 CALL get_cfm_submat(cfm_non_spinor_gamma, cfm_spinor_gamma, i_offset, j_offset)
2994 CALL cp_cfm_to_fm(cfm_non_spinor_gamma, fm_non_spinor_gamma_re, fm_non_spinor_gamma_im)
2995
2996 ! transform real part of Gamma-point matrix to ikp
2997 CALL cfm_ikp_from_fm_gamma(cfm_non_spinor_ikp, fm_non_spinor_gamma_re, &
2998 ikp, qs_env, kpoints, basis_type)
2999 CALL add_cfm_submat(cfm_spinor_ikp, cfm_non_spinor_ikp, i_offset, j_offset)
3000
3001 ! transform imag part of Gamma-point matrix to ikp
3002 CALL cfm_ikp_from_fm_gamma(cfm_non_spinor_ikp, fm_non_spinor_gamma_im, &
3003 ikp, qs_env, kpoints, basis_type)
3004 CALL add_cfm_submat(cfm_spinor_ikp, cfm_non_spinor_ikp, i_offset, j_offset, gaussi)
3005
3006 END DO
3007 END DO
3008
3009 CALL cp_cfm_release(cfm_non_spinor_gamma)
3010 CALL cp_cfm_release(cfm_non_spinor_ikp)
3011 CALL cp_fm_release(fm_non_spinor_gamma_re)
3012 CALL cp_fm_release(fm_non_spinor_gamma_im)
3013
3014 CALL timestop(handle)
3015
3016 END SUBROUTINE cfm_ikp_from_cfm_spinor_gamma
3017
3018! **************************************************************************************************
3019!> \brief Reconstruct a k-point operator from a real Gamma matrix using one MIC phase.
3020!> \param cfm_ikp ...
3021!> \param fm_Gamma ...
3022!> \param ikp ...
3023!> \param qs_env ...
3024!> \param kpoints ...
3025!> \param basis_type ...
3026! **************************************************************************************************
3027 SUBROUTINE cfm_ikp_from_fm_gamma(cfm_ikp, fm_Gamma, ikp, qs_env, kpoints, basis_type)
3028 TYPE(cp_cfm_type) :: cfm_ikp
3029 TYPE(cp_fm_type) :: fm_gamma
3030 INTEGER :: ikp
3031 TYPE(qs_environment_type), POINTER :: qs_env
3032 TYPE(kpoint_type), POINTER :: kpoints
3033 CHARACTER(LEN=*) :: basis_type
3034
3035 CHARACTER(LEN=*), PARAMETER :: routinen = 'cfm_ikp_from_fm_Gamma'
3036
3037 INTEGER :: handle
3038 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_bf
3039 TYPE(cell_type), POINTER :: cell
3040 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3041
3042 CALL timeset(routinen, handle)
3043
3044 cpassert(basis_type == "ORB" .OR. basis_type == "RI_AUX")
3045 ALLOCATE (atom_from_bf(fm_gamma%matrix_struct%nrow_global))
3046 CALL get_atom_index_from_basis_function_index(qs_env, atom_from_bf, SIZE(atom_from_bf), basis_type)
3047 NULLIFY (cell, particle_set)
3048 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
3049 CALL fm_gamma_to_kp(cfm_ikp, fm_gamma, cell, particle_set, atom_from_bf, &
3050 kpoints%index_to_cell, kpoints%xkp(:, ikp))
3051
3052 CALL timestop(handle)
3053
3054 END SUBROUTINE cfm_ikp_from_fm_gamma
3055
3056! **************************************************************************************************
3057!> \brief ...
3058!> \param bs_env ...
3059!> \param qs_env ...
3060!> \param fm_W_MIC_freq_j ...
3061!> \param cfm_W_ikp_freq_j ...
3062!> \param ikp ...
3063!> \param kpoints ...
3064!> \param basis_type ...
3065!> \param wkp_ext ...
3066! **************************************************************************************************
3067 SUBROUTINE mic_contribution_from_ikp(bs_env, qs_env, fm_W_MIC_freq_j, &
3068 cfm_W_ikp_freq_j, ikp, kpoints, basis_type, wkp_ext)
3069 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3070 TYPE(qs_environment_type), POINTER :: qs_env
3071 TYPE(cp_fm_type) :: fm_w_mic_freq_j
3072 TYPE(cp_cfm_type) :: cfm_w_ikp_freq_j
3073 INTEGER, INTENT(IN) :: ikp
3074 TYPE(kpoint_type), POINTER :: kpoints
3075 CHARACTER(LEN=*) :: basis_type
3076 REAL(kind=dp), OPTIONAL :: wkp_ext
3077
3078 CHARACTER(LEN=*), PARAMETER :: routinen = 'MIC_contribution_from_ikp'
3079
3080 INTEGER :: handle, i_bf, iatom, iatom_old, irow, &
3081 j_bf, jatom, jatom_old, jcol, n_bf, &
3082 ncol_local, nrow_local, num_cells
3083 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_bf_index
3084 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
3085 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
3086 REAL(kind=dp) :: contribution, weight_im, weight_re, &
3087 wkp_of_ikp
3088 REAL(kind=dp), DIMENSION(3, 3) :: hmat
3089 REAL(kind=dp), DIMENSION(:), POINTER :: wkp
3090 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
3091 TYPE(cell_type), POINTER :: cell
3092 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3093
3094 CALL timeset(routinen, handle)
3095
3096 ! get number of basis functions (bf) for different basis sets
3097 IF (basis_type == "ORB") THEN
3098 n_bf = qs_env%bs_env%n_ao
3099 ELSE IF (basis_type == "RI_AUX") THEN
3100 n_bf = qs_env%bs_env%n_RI
3101 ELSE
3102 cpabort("Only ORB and RI_AUX basis implemented.")
3103 END IF
3104
3105 ALLOCATE (atom_from_bf_index(n_bf))
3106 CALL get_atom_index_from_basis_function_index(qs_env, atom_from_bf_index, n_bf, basis_type)
3107
3108 NULLIFY (cell, particle_set)
3109 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
3110 CALL get_cell(cell=cell, h=hmat)
3111
3112 CALL cp_cfm_get_info(matrix=cfm_w_ikp_freq_j, &
3113 nrow_local=nrow_local, &
3114 ncol_local=ncol_local, &
3115 row_indices=row_indices, &
3116 col_indices=col_indices)
3117
3118 CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp)
3119 index_to_cell => kpoints%index_to_cell
3120 num_cells = SIZE(index_to_cell, 2)
3121
3122 iatom_old = 0
3123 jatom_old = 0
3124
3125 DO jcol = 1, ncol_local
3126 DO irow = 1, nrow_local
3127
3128 i_bf = row_indices(irow)
3129 j_bf = col_indices(jcol)
3130
3131 iatom = atom_from_bf_index(i_bf)
3132 jatom = atom_from_bf_index(j_bf)
3133
3134 IF (PRESENT(wkp_ext)) THEN
3135 wkp_of_ikp = wkp_ext
3136 ELSE
3137 SELECT CASE (bs_env%l_RI(i_bf) + bs_env%l_RI(j_bf))
3138 CASE (0)
3139 ! both RI functions are s-functions, k-extrapolation for 2D and 3D
3140 wkp_of_ikp = wkp(ikp)
3141 CASE (1)
3142 ! one function is an s-function, the other a p-function, k-extrapolation for 3D
3143 wkp_of_ikp = bs_env%wkp_s_p(ikp)
3144 CASE DEFAULT
3145 ! for any other matrix element of W, there is no need for extrapolation
3146 wkp_of_ikp = bs_env%wkp_no_extra(ikp)
3147 END SELECT
3148 END IF
3149
3150 IF (iatom /= iatom_old .OR. jatom /= jatom_old) THEN
3151
3152 CALL compute_weight_re_im(weight_re, weight_im, &
3153 num_cells, iatom, jatom, xkp(1:3, ikp), wkp_of_ikp, &
3154 cell, index_to_cell, hmat, particle_set)
3155
3156 iatom_old = iatom
3157 jatom_old = jatom
3158
3159 END IF
3160
3161 contribution = weight_re*real(cfm_w_ikp_freq_j%local_data(irow, jcol)) + &
3162 weight_im*aimag(cfm_w_ikp_freq_j%local_data(irow, jcol))
3163
3164 fm_w_mic_freq_j%local_data(irow, jcol) = fm_w_mic_freq_j%local_data(irow, jcol) &
3165 + contribution
3166
3167 END DO
3168 END DO
3169
3170 CALL timestop(handle)
3171
3172 END SUBROUTINE mic_contribution_from_ikp
3173
3174! **************************************************************************************************
3175!> \brief ...
3176!> \param xkp ...
3177!> \param ikp_start ...
3178!> \param ikp_end ...
3179!> \param grid ...
3180! **************************************************************************************************
3181 SUBROUTINE compute_xkp(xkp, ikp_start, ikp_end, grid)
3182
3183 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
3184 INTEGER :: ikp_start, ikp_end
3185 INTEGER, DIMENSION(3) :: grid
3186
3187 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_xkp'
3188
3189 INTEGER :: handle, i, ix, iy, iz
3190
3191 CALL timeset(routinen, handle)
3192
3193 i = ikp_start
3194 DO ix = 1, grid(1)
3195 DO iy = 1, grid(2)
3196 DO iz = 1, grid(3)
3197
3198 IF (i > ikp_end) cycle
3199
3200 xkp(1, i) = real(2*ix - grid(1) - 1, kind=dp)/(2._dp*real(grid(1), kind=dp))
3201 xkp(2, i) = real(2*iy - grid(2) - 1, kind=dp)/(2._dp*real(grid(2), kind=dp))
3202 xkp(3, i) = real(2*iz - grid(3) - 1, kind=dp)/(2._dp*real(grid(3), kind=dp))
3203 i = i + 1
3204
3205 END DO
3206 END DO
3207 END DO
3208
3209 CALL timestop(handle)
3210
3211 END SUBROUTINE compute_xkp
3212
3213! **************************************************************************************************
3214!> \brief ...
3215!> \param kpoints ...
3216!> \param qs_env ...
3217! **************************************************************************************************
3218 SUBROUTINE kpoint_init_cell_index_simple(kpoints, qs_env)
3219
3220 TYPE(kpoint_type), POINTER :: kpoints
3221 TYPE(qs_environment_type), POINTER :: qs_env
3222
3223 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_init_cell_index_simple'
3224
3225 INTEGER :: handle, nimages
3226 TYPE(mp_para_env_type), POINTER :: para_env
3227 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3228 POINTER :: sab_orb
3229
3230 CALL timeset(routinen, handle)
3231
3232 NULLIFY (para_env, sab_orb)
3233 CALL get_qs_env(qs_env=qs_env, para_env=para_env, sab_orb=sab_orb)
3234 CALL kpoint_init_cell_index(kpoints, sab_orb, para_env, nimages)
3235
3236 CALL timestop(handle)
3237
3238 END SUBROUTINE kpoint_init_cell_index_simple
3239
3240! **************************************************************************************************
3241!> \brief ...
3242!> \param qs_env ...
3243!> \param bs_env ...
3244! **************************************************************************************************
3245 SUBROUTINE soc(qs_env, bs_env)
3246 TYPE(qs_environment_type), POINTER :: qs_env
3247 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3248
3249 CHARACTER(LEN=*), PARAMETER :: routinen = 'soc'
3250
3251 INTEGER :: handle
3252
3253 CALL timeset(routinen, handle)
3254
3255 ! V^SOC_µν^(α),R = ħ/2 < ϕ_µ cell O | sum_ℓ ΔV_ℓ^SO(r,r') L^(α) | ϕ_ν cell R>, α = x,y,z
3256 ! see Hartwigsen, Goedecker, Hutter, Eq.(18), (19) (doi.org/10.1103/PhysRevB.58.3641)
3257 CALL v_soc_xyz_from_pseudopotential(qs_env, bs_env%mat_V_SOC_xyz)
3258
3259 ! Calculate H^SOC_µν,σσ'(k) = sum_α V^SOC_µν^(α)(k)*Pauli-matrix^(α)_σσ'
3260 ! see Hartwigsen, Goedecker, Hutter, Eq.(18) (doi.org/10.1103/PhysRevB.58.3641)
3261 SELECT CASE (bs_env%gw_implementation)
3262 CASE (tensor_large_cell_gamma, ri_rs_large_cell_gamma, ri_rs_non_periodic, tensor_static_cohsex)
3263
3264 ! H^SOC_µν,σσ' = sum_α V^SOC_µν^(α)*Pauli-matrix^(α)_σσ'
3265 CALL h_ks_spinor_gamma(bs_env)
3266
3267 CASE (tensor_small_cell_full_kp)
3268
3269 ! V^SOC_µν^(α),R -> V^SOC_µν^(α)(k); then calculate spinor H^SOC_µν,σσ'(k) (see above)
3270 CALL h_ks_spinor_kp(qs_env, bs_env)
3271
3272 END SELECT
3273
3274 CALL timestop(handle)
3275
3276 END SUBROUTINE soc
3277
3278! **************************************************************************************************
3279!> \brief ...
3280!> \param bs_env ...
3281! **************************************************************************************************
3282 SUBROUTINE h_ks_spinor_gamma(bs_env)
3283
3284 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3285
3286 CHARACTER(LEN=*), PARAMETER :: routinen = 'H_KS_spinor_Gamma'
3287
3288 INTEGER :: handle, nao, s
3289 TYPE(cp_fm_struct_type), POINTER :: str
3290
3291 CALL timeset(routinen, handle)
3292
3293 CALL cp_fm_get_info(bs_env%fm_ks_Gamma(1), nrow_global=nao)
3294
3295 ALLOCATE (bs_env%cfm_SOC_spinor_ao(1))
3296 CALL create_cfm_double(bs_env%cfm_SOC_spinor_ao(1), fm_orig=bs_env%fm_ks_Gamma(1))
3297 CALL cp_cfm_set_all(bs_env%cfm_SOC_spinor_ao(1), z_zero)
3298
3299 str => bs_env%fm_ks_Gamma(1)%matrix_struct
3300
3301 s = nao + 1
3302
3303 ! careful: inside add_dbcsr_submat, mat_V_SOC_xyz is multiplied by i because the real matrix
3304 ! mat_V_SOC_xyz is antisymmetric as V_SOC matrix is purely imaginary and Hermitian
3305 ! V_x * sigma_x: sigma_x = ((0,1),(1,0))
3306 ! ud block (1,s): +i*V_x
3307 CALL add_dbcsr_submat(bs_env%cfm_SOC_spinor_ao(1), bs_env%mat_V_SOC_xyz(1, 1)%matrix, &
3308 str, 1, s, z_one, .false.)
3309 ! du block (s,1): +i*V_x
3310 CALL add_dbcsr_submat(bs_env%cfm_SOC_spinor_ao(1), bs_env%mat_V_SOC_xyz(1, 1)%matrix, &
3311 str, s, 1, z_one, .false.)
3312
3313 ! V_y * sigma_y: sigma_y = ((0,-i),(i,0))
3314 ! ud block (1,s): i*(i*V_y) = -V_y (extra gaussi factor)
3315 CALL add_dbcsr_submat(bs_env%cfm_SOC_spinor_ao(1), bs_env%mat_V_SOC_xyz(2, 1)%matrix, &
3316 str, 1, s, gaussi, .false.)
3317 ! du block (s,1): -i*(i*V_y) = +V_y (extra -gaussi factor)
3318 CALL add_dbcsr_submat(bs_env%cfm_SOC_spinor_ao(1), bs_env%mat_V_SOC_xyz(2, 1)%matrix, &
3319 str, s, 1, -gaussi, .false.)
3320
3321 ! V_z * sigma_z: sigma_z = ((1,0),(0,-1))
3322 ! uu block (1,1): +i*V_z
3323 CALL add_dbcsr_submat(bs_env%cfm_SOC_spinor_ao(1), bs_env%mat_V_SOC_xyz(3, 1)%matrix, &
3324 str, 1, 1, z_one, .false.)
3325 ! dd block (s,s): -i*V_z
3326 CALL add_dbcsr_submat(bs_env%cfm_SOC_spinor_ao(1), bs_env%mat_V_SOC_xyz(3, 1)%matrix, &
3327 str, s, s, -z_one, .false.)
3328
3329 CALL timestop(handle)
3330
3331 END SUBROUTINE h_ks_spinor_gamma
3332
3333! **************************************************************************************************
3334!> \brief ...
3335!> \param qs_env ...
3336!> \param bs_env ...
3337! **************************************************************************************************
3338 SUBROUTINE h_ks_spinor_kp(qs_env, bs_env)
3339 TYPE(qs_environment_type), POINTER :: qs_env
3340 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3341
3342 CHARACTER(LEN=*), PARAMETER :: routinen = 'H_KS_spinor_kp'
3343
3344 INTEGER :: handle, i_dim, ikp, n_spin, &
3345 nkp_bs_and_dos, s
3346 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_scf
3347 REAL(kind=dp), DIMENSION(3) :: xkp
3348 TYPE(cp_cfm_type) :: cfm_v_soc_xyz_ikp
3349 TYPE(cp_fm_struct_type), POINTER :: str
3350 TYPE(kpoint_type), POINTER :: kpoints_scf
3351 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3352 POINTER :: sab_nl
3353
3354 CALL timeset(routinen, handle)
3355
3356 nkp_bs_and_dos = bs_env%nkp_bs_and_DOS
3357 n_spin = bs_env%n_spin
3358 s = bs_env%n_ao + 1
3359 str => bs_env%cfm_ks_kp(1, 1)%matrix_struct
3360
3361 CALL cp_cfm_create(cfm_v_soc_xyz_ikp, bs_env%cfm_work_mo%matrix_struct)
3362
3363 CALL alloc_cfm_double_array_1d(bs_env%cfm_SOC_spinor_ao, bs_env%cfm_ks_kp(1, 1), nkp_bs_and_dos)
3364
3365 CALL get_qs_env(qs_env, kpoints=kpoints_scf)
3366
3367 NULLIFY (sab_nl)
3368 CALL get_kpoint_info(kpoints_scf, sab_nl=sab_nl, cell_to_index=cell_to_index_scf)
3369
3370 DO i_dim = 1, 3
3371
3372 DO ikp = 1, nkp_bs_and_dos
3373
3374 xkp(1:3) = bs_env%kpoints_DOS%xkp(1:3, ikp)
3375
3376 CALL cp_cfm_set_all(cfm_v_soc_xyz_ikp, z_zero)
3377
3378 CALL rsmat_to_kp(bs_env%mat_V_SOC_xyz, i_dim, xkp, cell_to_index_scf, &
3379 sab_nl, bs_env, cfm_v_soc_xyz_ikp, imag_rs_mat=.true.)
3380
3381 ! multiply V_SOC with i because bs_env%mat_V_SOC_xyz stores imag. part (real part = 0)
3382 CALL cp_cfm_scale(gaussi, cfm_v_soc_xyz_ikp)
3383
3384 SELECT CASE (i_dim)
3385 CASE (1)
3386 ! add V^SOC_x * σ_x for σ_x = ( (0,1) (1,0) )
3387 CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_v_soc_xyz_ikp, 1, s)
3388 CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_v_soc_xyz_ikp, s, 1)
3389 CASE (2)
3390 ! add V^SOC_y * σ_y for σ_y = ( (0,-i) (i,0) )
3391 CALL cp_cfm_scale(gaussi, cfm_v_soc_xyz_ikp)
3392 CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_v_soc_xyz_ikp, 1, s)
3393 CALL cp_cfm_scale(-z_one, cfm_v_soc_xyz_ikp)
3394 CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_v_soc_xyz_ikp, s, 1)
3395 CASE (3)
3396 ! add V^SOC_z * σ_z for σ_z = ( (1,0) (0,1) )
3397 CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_v_soc_xyz_ikp, 1, 1)
3398 CALL cp_cfm_scale(-z_one, cfm_v_soc_xyz_ikp)
3399 CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_v_soc_xyz_ikp, s, s)
3400 END SELECT
3401
3402 END DO
3403
3404 END DO ! ikp
3405
3406 CALL cp_cfm_release(cfm_v_soc_xyz_ikp)
3407
3408 CALL timestop(handle)
3409
3410 END SUBROUTINE h_ks_spinor_kp
3411
3412! **************************************************************************************************
3413!> \brief ...
3414!> \param cfm_array ...
3415!> \param cfm_template ...
3416!> \param n ...
3417! **************************************************************************************************
3418 SUBROUTINE alloc_cfm_double_array_1d(cfm_array, cfm_template, n)
3419 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_array
3420 TYPE(cp_cfm_type) :: cfm_template
3421 INTEGER :: n
3422
3423 CHARACTER(LEN=*), PARAMETER :: routinen = 'alloc_cfm_double_array_1d'
3424
3425 INTEGER :: handle, i
3426
3427 CALL timeset(routinen, handle)
3428
3429 ALLOCATE (cfm_array(n))
3430 DO i = 1, n
3431 CALL create_cfm_double(cfm_array(i), cfm_orig=cfm_template)
3432 CALL cp_cfm_set_all(cfm_array(i), z_zero)
3433 END DO
3434
3435 CALL timestop(handle)
3436
3437 END SUBROUTINE alloc_cfm_double_array_1d
3438
3439! **************************************************************************************************
3440!> \brief ...
3441!> \param bs_env ...
3442! **************************************************************************************************
3443 SUBROUTINE get_all_vbm_cbm_bandgaps(bs_env)
3444
3445 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3446
3447 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_all_VBM_CBM_bandgaps'
3448
3449 INTEGER :: handle
3450
3451 CALL timeset(routinen, handle)
3452
3453 CALL get_vbm_cbm_bandgaps(bs_env%band_edges_scf, bs_env%eigenval_scf, bs_env)
3454 CALL get_vbm_cbm_bandgaps(bs_env%band_edges_GW, bs_env%eigenval_GW, bs_env)
3455 CALL get_vbm_cbm_bandgaps(bs_env%band_edges_HF, bs_env%eigenval_HF, bs_env)
3456
3457 CALL check_qp_gap_sanity(bs_env)
3458
3459 CALL check_scf_gw_level_ordering(bs_env)
3460
3461 CALL timestop(handle)
3462
3463 END SUBROUTINE get_all_vbm_cbm_bandgaps
3464
3465! **************************************************************************************************
3466!> \brief Warn if the GW fundamental band gap is inverted or implausibly large, i.e. if the
3467!> quasiparticle solve has produced a spectrum that cannot be physical.
3468!> \param bs_env ...
3469! **************************************************************************************************
3470 SUBROUTINE check_qp_gap_sanity(bs_env)
3471
3472 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3473
3474 CHARACTER(LEN=default_string_length) :: gw_label
3475 REAL(kind=dp) :: gap, gap_scf
3476
3477 gap = bs_env%band_edges_GW%IDBG
3478 gap_scf = bs_env%band_edges_scf%IDBG
3479 gw_label = gw_flavour_label(bs_env)
3480
3481 ! requiring a healthy SCF gap keeps the inversion test from firing on a genuine metal
3482 IF (gap < -eps_qp_gap .AND. gap_scf > eps_qp_gap) THEN
3483 CALL cp_warn(__location__, &
3484 trim(gw_label)//" band gap is negative ("// &
3485 trim(adjustl(cp_to_string(gap*evolt, '(F12.3)')))//" eV): the quasiparticle "// &
3486 "spectrum is inverted. Check numerical parameters.")
3487 ELSE IF (abs(gap) > max_qp_gap) THEN
3488 CALL cp_warn(__location__, &
3489 trim(gw_label)//" band gap is implausibly large ("// &
3490 trim(adjustl(cp_to_string(gap*evolt, '(F12.3)')))//" eV): the quasiparticle "// &
3491 "solve has likely diverged. Check numerical parameters.")
3492 END IF
3493
3494 END SUBROUTINE check_qp_gap_sanity
3495
3496! **************************************************************************************************
3497!> \brief Check whether the GW correction reorders the frontier levels with respect to the SCF
3498!> eigenvalues, i.e. whether the SCF HOMO (LUMO) is a different level than the GW HOMO
3499!> (LUMO). If so, the SCF and GW band edges printed afterwards belong to different
3500!> orbitals, which is easily overlooked => warn.
3501!> Levels closer than eps_degeneracy in the SCF spectrum count as degenerate, i.e. their
3502!> reordering by GW is not reported.
3503!> \param bs_env ...
3504! **************************************************************************************************
3505 SUBROUTINE check_scf_gw_level_ordering(bs_env)
3506
3507 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3508
3509 CHARACTER(LEN=*), PARAMETER :: routinen = 'check_scf_gw_level_ordering'
3510 REAL(kind=dp), PARAMETER :: eps_degeneracy = 1.0e-05_dp
3511
3512 CHARACTER(LEN=default_string_length) :: gw_label, level_gw, level_scf, &
3513 spin_string
3514 CHARACTER(LEN=max_line_length) :: msg
3515 INTEGER :: handle, homo, ispin, n_mo
3516 INTEGER, DIMENSION(2) :: loc_gw, loc_scf
3517
3518 CALL timeset(routinen, handle)
3519
3520 n_mo = bs_env%n_ao
3521 gw_label = gw_flavour_label(bs_env)
3522
3523 DO ispin = 1, bs_env%n_spin
3524
3525 homo = bs_env%n_occ(ispin)
3526
3527 IF (bs_env%n_spin == 2) THEN
3528 WRITE (spin_string, '(A,I0,A)') " (spin ", ispin, ")"
3529 ELSE
3530 spin_string = ""
3531 END IF
3532
3533 ! occupied levels: level carrying the SCF VBM versus level carrying the GW VBM
3534 loc_scf(:) = maxloc(bs_env%eigenval_scf(1:homo, :, ispin))
3535 loc_gw(:) = maxloc(bs_env%eigenval_GW(1:homo, :, ispin))
3536
3537 ! only report a genuine reordering, not a permutation within a degenerate SCF shell
3538 IF (bs_env%eigenval_scf(loc_scf(1), loc_scf(2), ispin) - &
3539 bs_env%eigenval_scf(loc_gw(1), loc_gw(2), ispin) > eps_degeneracy) THEN
3540
3541 CALL level_string(level_scf, loc_scf, bs_env)
3542 CALL level_string(level_gw, loc_gw, bs_env)
3543
3544 msg = trim(gw_label)//" reorders the occupied levels"//trim(spin_string)// &
3545 ": the SCF valence band maximum (SCF HOMO) is "//trim(level_scf)// &
3546 ", whereas the "//trim(gw_label)//" valence band maximum ("// &
3547 trim(gw_label)//" HOMO) is "//trim(level_gw)//". The SCF and "// &
3548 trim(gw_label)//" band edges printed below therefore belong to different "// &
3549 "orbitals; see the eigenvalues in the file bandstructure_SCF_and_G0W0."
3550 cpwarn(trim(msg))
3551
3552 END IF
3553
3554 ! empty levels: level carrying the SCF CBM versus level carrying the GW CBM
3555 loc_scf(:) = minloc(bs_env%eigenval_scf(homo + 1:n_mo, :, ispin))
3556 loc_gw(:) = minloc(bs_env%eigenval_GW(homo + 1:n_mo, :, ispin))
3557 loc_scf(1) = loc_scf(1) + homo
3558 loc_gw(1) = loc_gw(1) + homo
3559
3560 IF (bs_env%eigenval_scf(loc_gw(1), loc_gw(2), ispin) - &
3561 bs_env%eigenval_scf(loc_scf(1), loc_scf(2), ispin) > eps_degeneracy) THEN
3562
3563 CALL level_string(level_scf, loc_scf, bs_env)
3564 CALL level_string(level_gw, loc_gw, bs_env)
3565
3566 msg = trim(gw_label)//" reorders the empty levels"//trim(spin_string)// &
3567 ": the SCF conduction band minimum (SCF LUMO) is "//trim(level_scf)// &
3568 ", whereas the "//trim(gw_label)//" conduction band minimum ("// &
3569 trim(gw_label)//" LUMO) is "//trim(level_gw)//". The SCF and "// &
3570 trim(gw_label)//" band edges printed below therefore belong to different "// &
3571 "orbitals; see the eigenvalues in the file bandstructure_SCF_and_G0W0."
3572 cpwarn(trim(msg))
3573
3574 END IF
3575
3576 END DO
3577
3578 CALL timestop(handle)
3579
3580 END SUBROUTINE check_scf_gw_level_ordering
3581
3582! **************************************************************************************************
3583!> \brief Human readable name of the level (i_mo, ikp) given as loc = [i_mo, ikp]
3584!> \param string ...
3585!> \param loc ...
3586!> \param bs_env ...
3587! **************************************************************************************************
3588 SUBROUTINE level_string(string, loc, bs_env)
3589
3590 CHARACTER(LEN=*), INTENT(OUT) :: string
3591 INTEGER, DIMENSION(2), INTENT(IN) :: loc
3592 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3593
3594 IF (bs_env%nkp_bs_and_DOS > 1) THEN
3595 WRITE (string, '(A,I0,A,I0)') "level ", loc(1), " at k-point ", loc(2)
3596 ELSE
3597 WRITE (string, '(A,I0)') "level ", loc(1)
3598 END IF
3599
3600 END SUBROUTINE level_string
3601
3602! **************************************************************************************************
3603!> \brief ...
3604!> \param band_edges ...
3605!> \param ev ...
3606!> \param bs_env ...
3607! **************************************************************************************************
3608 SUBROUTINE get_vbm_cbm_bandgaps(band_edges, ev, bs_env)
3609 TYPE(band_edges_type) :: band_edges
3610 REAL(kind=dp), DIMENSION(:, :, :) :: ev
3611 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3612
3613 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_VBM_CBM_bandgaps'
3614
3615 INTEGER :: handle, homo, homo_1, homo_2, ikp, &
3616 ispin, lumo, lumo_1, lumo_2, n_mo
3617 REAL(kind=dp) :: e_dbg_at_ikp
3618
3619 CALL timeset(routinen, handle)
3620
3621 n_mo = bs_env%n_ao
3622
3623 band_edges%DBG = 1000.0_dp
3624
3625 SELECT CASE (bs_env%n_spin)
3626 CASE (1)
3627 homo = bs_env%n_occ(1)
3628 lumo = homo + 1
3629 band_edges%VBM = maxval(ev(1:homo, :, 1))
3630 band_edges%CBM = minval(ev(homo + 1:n_mo, :, 1))
3631 CASE (2)
3632 homo_1 = bs_env%n_occ(1)
3633 lumo_1 = homo_1 + 1
3634 homo_2 = bs_env%n_occ(2)
3635 lumo_2 = homo_2 + 1
3636 band_edges%VBM = max(maxval(ev(1:homo_1, :, 1)), maxval(ev(1:homo_2, :, 2)))
3637 band_edges%CBM = min(minval(ev(homo_1 + 1:n_mo, :, 1)), minval(ev(homo_2 + 1:n_mo, :, 2)))
3638 CASE DEFAULT
3639 cpabort("Error with number of spins.")
3640 END SELECT
3641
3642 band_edges%IDBG = band_edges%CBM - band_edges%VBM
3643
3644 DO ispin = 1, bs_env%n_spin
3645
3646 homo = bs_env%n_occ(ispin)
3647
3648 DO ikp = 1, bs_env%nkp_bs_and_DOS
3649
3650 e_dbg_at_ikp = -maxval(ev(1:homo, ikp, ispin)) + minval(ev(homo + 1:n_mo, ikp, ispin))
3651
3652 IF (e_dbg_at_ikp < band_edges%DBG) band_edges%DBG = e_dbg_at_ikp
3653
3654 END DO
3655
3656 END DO
3657 SELECT CASE (bs_env%n_spin)
3658 CASE (1)
3659 CALL get_fermi_level_1c(bs_env, ev(:, :, 1), &
3660 REAL(bs_env%n_electron(1), kind=dp), &
3661 bs_env%spin_degeneracy, band_edges%VBM, band_edges%CBM, band_edges%E_fermi)
3662 CASE (2)
3663 CALL get_fermi_level_2c(bs_env, ev(:, :, 1:2), &
3664 REAL(bs_env%n_electron(1) + bs_env%n_electron(2), kind=dp), &
3665 band_edges%VBM, band_edges%CBM, band_edges%E_fermi)
3666 END SELECT
3667
3668 CALL timestop(handle)
3669
3670 END SUBROUTINE get_vbm_cbm_bandgaps
3671
3672! **************************************************************************************************
3673!> \brief Fermi level (chemical potential) for one set of band energies on the (weighted)
3674!> DOS k-point mesh bs_env%kpoints_DOS (band-structure-path k-points carry wkp=0 and
3675!> therefore do not contribute). Used both for a single spin-degenerate channel
3676!> (n_spin=1, maxocc=bs_env%spin_degeneracy=2) and for the spin-orbit-coupled spinor
3677!> space, where the two spins are already merged into one set of 2*n_ao states, each
3678!> holding at most one electron (maxocc=1).
3679!> IMPORTANT: Smearkp's T=0 bisection does *not* converge to the gap midpoint for a
3680!> real (flat) gap: N(mu) is exactly n_el for every mu in [VBM,CBM), so the bisection
3681!> exits on its very first probe, wherever that lands inside the gap, not at the
3682!> midpoint. Therefore Smearkp is only invoked when bs_env%do_smear is active (a
3683!> genuine metal, smeared with the ground-state SCF's &SMEAR settings); otherwise
3684!> E_fermi = VBM (highest occupied eigenvalue), matching the convention ABINIT uses
3685!> for its fixed-occupation-number case (occopt<3), rather than the gap midpoint.
3686!> \param bs_env ...
3687!> \param eigenval band energies (n_states, n_kpoints)
3688!> \param n_el target (weighted) number of occupied states
3689!> \param maxocc maximum occupation per state
3690!> \param VBM valence band maximum, used as E_fermi when bs_env%do_smear = .FALSE.
3691!> \param CBM conduction band minimum, used only to detect a wide gap
3692!> \param e_fermi ...
3693! **************************************************************************************************
3694 SUBROUTINE get_fermi_level_1c(bs_env, eigenval, n_el, maxocc, VBM, CBM, e_fermi)
3695 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3696 REAL(kind=dp), DIMENSION(:, :) :: eigenval
3697 REAL(kind=dp) :: n_el, maxocc, vbm, cbm, e_fermi
3698
3699 REAL(kind=dp), PARAMETER :: gap_smear_ratio = 20.0_dp
3700
3701 REAL(kind=dp) :: kts
3702 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: f_occ
3703
3704 IF (bs_env%do_smear .AND. .NOT. (cbm > vbm .AND. &
3705 cbm - vbm > gap_smear_ratio*bs_env%smear_width)) THEN
3706 ALLOCATE (f_occ(SIZE(eigenval, 1), SIZE(eigenval, 2)))
3707 CALL smearkp(f_occ, e_fermi, kts, eigenval, n_el, &
3708 bs_env%kpoints_DOS%wkp(1:SIZE(eigenval, 2)), &
3709 bs_env%smear_width, maxocc, bs_env%smear_method)
3710 DEALLOCATE (f_occ)
3711 ELSE
3712 e_fermi = vbm
3713 END IF
3714
3715 END SUBROUTINE get_fermi_level_1c
3716
3717! **************************************************************************************************
3718!> \brief Fermi level (chemical potential) shared across two spin channels with independent
3719!> band energies but the same chemical potential (standard spin-polarized DFT); see
3720!> get_fermi_level_1c for the k-point weighting and for why Smearkp2 is only used
3721!> when bs_env%do_smear is active (E_fermi = VBM otherwise).
3722!> \param bs_env ...
3723!> \param eigenval band energies (n_states, n_kpoints, 2)
3724!> \param n_el target (weighted) total number of occupied states, summed over both spins
3725!> \param VBM valence band maximum, used as E_fermi when bs_env%do_smear = .FALSE.
3726!> \param CBM conduction band minimum, used only to detect a wide gap
3727!> \param e_fermi ...
3728! **************************************************************************************************
3729 SUBROUTINE get_fermi_level_2c(bs_env, eigenval, n_el, VBM, CBM, e_fermi)
3730 TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3731 REAL(kind=dp), DIMENSION(:, :, :) :: eigenval
3732 REAL(kind=dp) :: n_el, vbm, cbm, e_fermi
3733
3734 REAL(kind=dp), PARAMETER :: gap_smear_ratio = 20.0_dp
3735
3736 REAL(kind=dp) :: kts
3737 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: f_occ
3738
3739 IF (bs_env%do_smear .AND. .NOT. (cbm > vbm .AND. &
3740 cbm - vbm > gap_smear_ratio*bs_env%smear_width)) THEN
3741 ALLOCATE (f_occ(SIZE(eigenval, 1), SIZE(eigenval, 2), 2))
3742 CALL smearkp2(f_occ, e_fermi, kts, eigenval, n_el, &
3743 bs_env%kpoints_DOS%wkp(1:SIZE(eigenval, 2)), &
3744 bs_env%smear_width, bs_env%smear_method)
3745 DEALLOCATE (f_occ)
3746 ELSE
3747 e_fermi = vbm
3748 END IF
3749
3750 END SUBROUTINE get_fermi_level_2c
3751END MODULE post_scf_bandstructure_utils
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:126
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition grid_common.h:162
cp_cfm_basic_linalg::cp_cfm_scale
Definition cp_cfm_basic_linalg.F:64
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:57
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:95
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:103
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:89
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:84
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
bse_input
Reads the &BSE input section into the BSE environment; called by read_mp2_section with the section un...
Definition bse_input.F:15
bse_input::read_bse_section
subroutine, public read_bse_section(bse_section, bse_env)
Reads every keyword of the &BSE section into the settings of the BSE environment.
Definition bse_input.F:38
bse_main
Main routines for GW + Bethe-Salpeter for computing electronic excitations: a GW path prepares the BS...
Definition bse_main.F:16
bse_main::bse_fulldiag_memory_check
subroutine, public bse_fulldiag_memory_check(bse_env, window, homo, para_env, unit_nr)
Early memory check for the full diagonalization, before any GW work: the window on the DFT axis gives...
Definition bse_main.F:506
bse_types
The BSE environment: the settings of the &BSE section and the state a GW path prepares for the solver...
Definition bse_types.F:14
bse_util
Auxiliary routines for GW + Bethe-Salpeter for computing electronic excitations.
Definition bse_util.F:13
bse_util::bse_window_from_dft
subroutine, public bse_window_from_dft(eigenval_scf, n_mo, homo, cutoff_occ, cutoff_empty, window, warn)
The active MO window on the DFT axis: one absolute MO range for every spin; when the spin windows dif...
Definition bse_util.F:576
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:234
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_cholesky
various cholesky decomposition related routines
Definition cp_cfm_cholesky.F:13
cp_cfm_cholesky::cp_cfm_cholesky_decompose
subroutine, public cp_cfm_cholesky_decompose(matrix, n, info_out)
Used to replace a symmetric positive definite matrix M with its Cholesky decomposition U: M = U^T * U...
Definition cp_cfm_cholesky.F:53
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_geeig_canon
subroutine, public cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig, nmo_retained)
General Eigenvalue Problem AX = BXE Use canonical orthogonalization.
Definition cp_cfm_diag.F:686
cp_cfm_diag::cp_cfm_heevd
subroutine, public cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
Perform a diagonalisation of a complex matrix.
Definition cp_cfm_diag.F:94
cp_cfm_diag::cp_cfm_geeig
subroutine, public cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work, lowest_subset)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition cp_cfm_diag.F:439
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:183
cp_cfm_types::cp_fm_to_cfm
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
Definition cp_cfm_types.F:845
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:126
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:632
cp_cfm_types::cp_cfm_set_all
subroutine, public cp_cfm_set_all(matrix, alpha, beta)
Set all elements of the full matrix to alpha. Besides, set all diagonal matrix elements to beta (if g...
Definition cp_cfm_types.F:203
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition cp_cfm_types.F:793
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_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:525
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
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, plan, check_plan)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:233
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:129
cp_files
Utility routines to open and close files.
Definition cp_files.F:22
cp_files::open_file
subroutine, public open_file(file_name, unit_number, file_status, file_form, file_action, file_position, file_pad, file_access, debug)
Opens a file and returns its logical unit number.
Definition cp_files.F:197
cp_files::close_file
subroutine, public close_file(unit_number, file_status)
Closes an open file given by its logical unit number.
Definition cp_files.F:64
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::cp_fm_geeig_canon
subroutine, public cp_fm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig, nmo_retained)
General Eigenvalue Problem AX = BXE Use canonical diagonalization : U*s**(-1/2).
Definition cp_fm_diag.F:1718
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
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:1045
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag, n_active)
returns the diagonal elements of a fm
Definition cp_fm_types.F:566
cp_fm_types::cp_fm_to_fm_submat
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
Definition cp_fm_types.F:1500
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:166
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_parser_methods
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition cp_parser_methods.F:20
cp_parser_methods::read_float_object
elemental subroutine, public read_float_object(string, object, error_message)
Returns a floating point number read from a string including fraction like z1/z2.
Definition cp_parser_methods.F:1265
dft_plus_u
Add the DFT+U contribution to the Hamiltonian matrix.
Definition dft_plus_u.F:18
dft_plus_u::lowdin_kp_prepare
subroutine, public lowdin_kp_prepare(kp, orbitals, use_real_wfn, smat)
Prepare selected Lowdin projectors and their reusable product workspace.
Definition dft_plus_u.F:1015
dft_plus_u::lowdin_kp_apply
subroutine, public lowdin_kp_apply(kp, v, rmat, cmat)
Add B diag(v) B**H before any solver, DIIS residual, or history storage uses H(k).
Definition dft_plus_u.F:1267
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::smear_fermi_dirac
integer, parameter, public smear_fermi_dirac
Definition input_constants.F:951
input_constants::tensor_small_cell_full_kp
integer, parameter, public tensor_small_cell_full_kp
Definition input_constants.F:1151
input_constants::ri_rs_non_periodic
integer, parameter, public ri_rs_non_periodic
Definition input_constants.F:1151
input_constants::ri_rs_large_cell_gamma
integer, parameter, public ri_rs_large_cell_gamma
Definition input_constants.F:1151
input_constants::tensor_large_cell_gamma
integer, parameter, public tensor_large_cell_gamma
Definition input_constants.F:1151
input_constants::bse_fulldiag
integer, parameter, public bse_fulldiag
Definition input_constants.F:1393
input_constants::smear_gaussian
integer, parameter, public smear_gaussian
Definition input_constants.F:951
input_constants::smear_mv
integer, parameter, public smear_mv
Definition input_constants.F:951
input_constants::g0w0
integer, parameter, public g0w0
Definition input_constants.F:1151
input_constants::int_ldos_z
integer, parameter, public int_ldos_z
Definition input_constants.F:1151
input_constants::evgw0
integer, parameter, public evgw0
Definition input_constants.F:1151
input_constants::gaussian
integer, parameter, public gaussian
Definition input_constants.F:164
input_constants::tensor_static_cohsex
integer, parameter, public tensor_static_cohsex
Definition input_constants.F:1151
input_constants::cohsex_off
integer, parameter, public cohsex_off
Definition input_constants.F:1181
input_constants::smear_mp
integer, parameter, public smear_mp
Definition input_constants.F:951
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_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:735
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:708
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:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::max_line_length
integer, parameter, public max_line_length
Definition kinds.F:59
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kpoint_k_r_trafo_simple
Implements simple k-space/R-space transformations and Gamma-point MIC reconstruction.
Definition kpoint_k_r_trafo_simple.F:16
kpoint_k_r_trafo_simple::fm_gamma_to_kp
subroutine, public fm_gamma_to_kp(cfm_kp, fm_gamma, cell, particle_set, atom_from_bf, index_to_cell, xkp)
Reconstruct a complex k-point matrix from a real Gamma matrix using one MIC phase....
Definition kpoint_k_r_trafo_simple.F:283
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:1138
kpoint_methods::kpoint_init_cell_index
subroutine, public kpoint_init_cell_index(kpoint, sab_nl, para_env, nimages)
Generates the mapping of cell indices and linear RS index CELL (0,0,0) is always mapped to index 1.
Definition kpoint_methods.F:1004
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, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:585
kpoint_types::kpoint_create
subroutine, public kpoint_create(kpoint)
Create a kpoint environment.
Definition kpoint_types.F:265
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:144
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:144
mathconstants::gaussi
complex(kind=dp), parameter, public gaussi
Definition mathconstants.F:143
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:144
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
mo_window
Common selection and union operations for contiguous molecular-orbital windows.
Definition mo_window.F:11
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
post_scf_bandstructure_types
Definition post_scf_bandstructure_types.F:13
post_scf_bandstructure_types::max_qp_gap
real(kind=dp), parameter, public max_qp_gap
Definition post_scf_bandstructure_types.F:392
post_scf_bandstructure_types::eps_qp_gap
real(kind=dp), parameter, public eps_qp_gap
Definition post_scf_bandstructure_types.F:392
post_scf_bandstructure_utils
Definition post_scf_bandstructure_utils.F:13
post_scf_bandstructure_utils::eval_bandstructure_properties
subroutine, public eval_bandstructure_properties(qs_env, bs_env)
...
Definition post_scf_bandstructure_utils.F:1243
post_scf_bandstructure_utils::rsmat_to_kp
subroutine, public rsmat_to_kp(mat_rs, ispin, xkp, cell_to_index_scf, sab_nl, bs_env, cfm_kp, imag_rs_mat)
...
Definition post_scf_bandstructure_utils.F:863
post_scf_bandstructure_utils::kpoint_init_cell_index_simple
subroutine, public kpoint_init_cell_index_simple(kpoints, qs_env)
...
Definition post_scf_bandstructure_utils.F:3219
post_scf_bandstructure_utils::allocate_gw_eigenvalues
subroutine, public allocate_gw_eigenvalues(bs_env)
Allocate the arrays holding the GW quasiparticle energies.
Definition post_scf_bandstructure_utils.F:144
post_scf_bandstructure_utils::gw_flavour_label
character(len=default_string_length) function, public gw_flavour_label(bs_env)
Name of the GW flavour that was requested, for printing.
Definition post_scf_bandstructure_utils.F:184
post_scf_bandstructure_utils::cfm_ikp_from_fm_gamma
subroutine, public cfm_ikp_from_fm_gamma(cfm_ikp, fm_gamma, ikp, qs_env, kpoints, basis_type)
Reconstruct a k-point operator from a real Gamma matrix using one MIC phase.
Definition post_scf_bandstructure_utils.F:3028
post_scf_bandstructure_utils::get_all_vbm_cbm_bandgaps
subroutine, public get_all_vbm_cbm_bandgaps(bs_env)
...
Definition post_scf_bandstructure_utils.F:3444
post_scf_bandstructure_utils::soc
subroutine, public soc(qs_env, bs_env)
...
Definition post_scf_bandstructure_utils.F:3246
post_scf_bandstructure_utils::mic_contribution_from_ikp
subroutine, public mic_contribution_from_ikp(bs_env, qs_env, fm_w_mic_freq_j, cfm_w_ikp_freq_j, ikp, kpoints, basis_type, wkp_ext)
...
Definition post_scf_bandstructure_utils.F:3069
post_scf_bandstructure_utils::compute_xkp
subroutine, public compute_xkp(xkp, ikp_start, ikp_end, grid)
...
Definition post_scf_bandstructure_utils.F:3182
post_scf_bandstructure_utils::create_and_init_bs_env
subroutine, public create_and_init_bs_env(qs_env, bs_env, post_scf_bandstructure_section)
...
Definition post_scf_bandstructure_utils.F:204
post_scf_bandstructure_utils::get_vbm_cbm_bandgaps
subroutine, public get_vbm_cbm_bandgaps(band_edges, ev, bs_env)
...
Definition post_scf_bandstructure_utils.F:3609
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition qs_collocate_density.F:38
qs_collocate_density::calculate_rho_elec
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
Definition qs_collocate_density.F:1713
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, mimic, 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_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, 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, matrix_vhxc, 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, xcint_weights, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, harris_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, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:544
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, cmo_coeff)
Get the components of a MO set data structure.
Definition qs_mo_types.F:445
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
rpa_gw_im_time_util
Utility routines for GW with imaginary time.
Definition rpa_gw_im_time_util.F:13
rpa_gw_im_time_util::compute_weight_re_im
subroutine, public compute_weight_re_im(weight_re, weight_im, num_cells, iatom, jatom, xkp, wkp_w, cell, index_to_cell, hmat, particle_set)
...
Definition rpa_gw_im_time_util.F:634
rpa_gw_im_time_util::get_atom_index_from_basis_function_index
subroutine, public get_atom_index_from_basis_function_index(qs_env, atom_from_basis_index, basis_size, basis_type, first_bf_from_atom)
...
Definition rpa_gw_im_time_util.F:579
scf_control_types
parameters that control an scf iteration
Definition scf_control_types.F:17
smearing_utils
Unified smearing module supporting four methods: smear_fermi_dirac — Fermi-Dirac distribution smear_g...
Definition smearing_utils.F:24
smearing_utils::smearkp
subroutine, public smearkp(f, mu, kts, e, nel, wk, sigma, maxocc, method)
Bisection search for mu given a target electron count (k-point case, single spin channel or spin-dege...
Definition smearing_utils.F:532
smearing_utils::smearkp2
subroutine, public smearkp2(f, mu, kts, e, nel, wk, sigma, method)
Bisection search for mu (k-point, spin-polarised with a shared chemical potential across both spin ch...
Definition smearing_utils.F:676
soc_pseudopotential_methods
Definition soc_pseudopotential_methods.F:8
soc_pseudopotential_methods::v_soc_xyz_from_pseudopotential
subroutine, public v_soc_xyz_from_pseudopotential(qs_env, mat_v_soc_xyz)
V^SOC_µν^(α),R = ħ/2 < ϕ_µ cell O | sum_ℓ ΔV_ℓ^SO(r,r') L^(α) | ϕ_ν cell R>, α = x,...
Definition soc_pseudopotential_methods.F:63
soc_pseudopotential_methods::remove_soc_outside_energy_window_mo
subroutine, public remove_soc_outside_energy_window_mo(cfm_ks_spinor, e_win_cbm, temp_smear, eigenval, e_fermi)
...
Definition soc_pseudopotential_methods.F:194
soc_pseudopotential_utils
Definition soc_pseudopotential_utils.F:8
soc_pseudopotential_utils::create_cfm_double
subroutine, public create_cfm_double(cfm_double, fm_orig, cfm_orig)
...
Definition soc_pseudopotential_utils.F:325
soc_pseudopotential_utils::add_dbcsr_submat
subroutine, public add_dbcsr_submat(cfm_mat_target, mat_source, fm_struct_source, nstart_row, nstart_col, factor, add_also_herm_conj)
...
Definition soc_pseudopotential_utils.F:60
soc_pseudopotential_utils::add_cfm_submat
subroutine, public add_cfm_submat(cfm_mat_target, cfm_mat_source, nstart_row, nstart_col, factor)
...
Definition soc_pseudopotential_utils.F:209
soc_pseudopotential_utils::get_cfm_submat
subroutine, public get_cfm_submat(cfm_mat_target, cfm_mat_source, nstart_row, nstart_col)
...
Definition soc_pseudopotential_utils.F:273
soc_pseudopotential_utils::cfm_add_on_diag
subroutine, public cfm_add_on_diag(cfm, alpha)
...
Definition soc_pseudopotential_utils.F:122
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3362
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
bse_types::bse_env_type
Settings of the &BSE section (read_bse_section, re-read before every BSE run, normalised in place by ...
Definition bse_types.F:44
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:70
cp_control_types::dft_control_type
Definition cp_control_types.F:786
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
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:115
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_env_type
Keeps information about a specific k-point.
Definition kpoint_types.F:94
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:196
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
mo_window::mo_window_type
Definition mo_window.F:21
particle_types::particle_type
Definition particle_types.F:35
post_scf_bandstructure_types::band_edges_type
Definition post_scf_bandstructure_types.F:57
post_scf_bandstructure_types::post_scf_bandstructure_type
Definition post_scf_bandstructure_types.F:65
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:223
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:137
qs_mo_types::mo_set_type
Definition qs_mo_types.F:53
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
scf_control_types::scf_control_type
Definition scf_control_types.F:135