Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief
10 : !> \author Jan Wilhelm
11 : !> \date 07.2023
12 : ! **************************************************************************************************
13 : MODULE 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 :
137 : CONTAINS
138 :
139 : ! **************************************************************************************************
140 : !> \brief Allocate the arrays holding the GW quasiparticle energies
141 : !> \param bs_env ...
142 : ! **************************************************************************************************
143 156 : 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 156 : CALL timeset(routineN, handle)
151 :
152 156 : n_ao = bs_env%n_ao
153 156 : nkp = bs_env%nkp_bs_and_DOS
154 156 : n_spin = bs_env%n_spin
155 :
156 156 : IF (.NOT. ALLOCATED(bs_env%eigenval_GW)) THEN
157 700 : ALLOCATE (bs_env%eigenval_GW(n_ao, nkp, n_spin))
158 : END IF
159 156 : IF (.NOT. ALLOCATED(bs_env%eigenval_HF)) THEN
160 700 : ALLOCATE (bs_env%eigenval_HF(n_ao, nkp, n_spin))
161 : END IF
162 :
163 312 : SELECT CASE (bs_env%gw_flavour)
164 : CASE (G0W0, evGW0)
165 156 : IF (.NOT. ALLOCATED(bs_env%eigenval_G0W0)) THEN
166 700 : ALLOCATE (bs_env%eigenval_G0W0(n_ao, nkp, n_spin))
167 : END IF
168 : END SELECT
169 :
170 156 : IF (bs_env%gw_flavour == evGW0 .AND. .NOT. ALLOCATED(bs_env%eigenval_evGW0)) THEN
171 20 : ALLOCATE (bs_env%eigenval_evGW0(n_ao, nkp, n_spin))
172 : END IF
173 :
174 156 : CALL timestop(handle)
175 :
176 156 : 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 689 : 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 744 : SELECT CASE (bs_env%gw_flavour)
188 : CASE (evGW0)
189 55 : label = "evGW0"
190 : CASE DEFAULT
191 689 : label = "G0W0"
192 : END SELECT
193 689 : IF (bs_env%do_static_cohsex) label = "static COHSEX"
194 :
195 689 : 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 140 : 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 140 : CALL timeset(routineN, handle)
213 :
214 13720 : ALLOCATE (bs_env)
215 :
216 140 : CALL print_header(bs_env)
217 :
218 140 : CALL read_bandstructure_input_parameters(bs_env, post_scf_bandstructure_section, qs_env)
219 :
220 140 : CALL get_parameters_from_qs_env(qs_env, bs_env)
221 :
222 140 : CALL set_heuristic_parameters(bs_env)
223 :
224 262 : 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 122 : CALL setup_kpoints_DOS_large_cell_Gamma(qs_env, bs_env, bs_env%kpoints_DOS)
228 :
229 122 : CALL allocate_and_fill_fm_ks_fm_s(qs_env, bs_env)
230 :
231 122 : CALL diagonalize_ks_matrix(bs_env)
232 :
233 122 : CALL check_positive_definite_overlap_mat(bs_env, qs_env)
234 :
235 : CASE (tensor_small_cell_full_kp)
236 :
237 18 : CALL setup_kpoints_scf_desymm(qs_env, bs_env, bs_env%kpoints_scf_desymm, .TRUE.)
238 :
239 18 : CALL setup_kpoints_DOS_small_cell_full_kp(bs_env, bs_env%kpoints_DOS)
240 :
241 18 : CALL allocate_and_fill_fm_ks_fm_s(qs_env, bs_env)
242 :
243 158 : CALL compute_cfm_mo_coeff_kp_and_eigenval_scf_kp(qs_env, bs_env)
244 :
245 : END SELECT
246 :
247 140 : CALL read_bse_input(qs_env, bs_env)
248 :
249 140 : CALL timestop(handle)
250 :
251 140 : 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 140 : 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 560 : TYPE(bse_env_type) :: bse_env
269 : TYPE(mo_window_type) :: window
270 : TYPE(section_vals_type), POINTER :: bse_sec, input
271 :
272 140 : CALL timeset(routineN, handle)
273 :
274 140 : NULLIFY (bse_sec, input)
275 140 : CALL get_qs_env(qs_env, input=input)
276 140 : bse_sec => section_vals_get_subs_vals(input, "PROPERTIES%BSE")
277 140 : CALL section_vals_get(bse_sec, explicit=do_bse_props)
278 :
279 140 : IF (do_bse_props) THEN
280 : ! the full diagonalization's footprint is known here, before any GW work
281 6 : CALL read_bse_section(bse_sec, bse_env)
282 6 : IF (bse_env%bse_diag_method == bse_fulldiag) THEN
283 4 : n_spin = bs_env%n_spin
284 : ! the canonical orthogonalisation parks the linearly dependent directions above n_mo_retained
285 12 : 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 4 : 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 4 : bs_env%unit_nr)
292 : END IF
293 : END IF
294 140 : bs_env%keep_W_freq_zero = do_bse_props
295 140 : bs_env%keep_3c_nl = do_bse_props
296 :
297 140 : CALL timestop(handle)
298 :
299 140 : END SUBROUTINE read_bse_input
300 :
301 : ! **************************************************************************************************
302 : !> \brief ...
303 : !> \param bs_env ...
304 : !> \param bs_sec ...
305 : !> \param qs_env ...
306 : ! **************************************************************************************************
307 140 : 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 140 : 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 140 : CALL timeset(routineN, handle)
326 140 : NULLIFY (cell)
327 140 : CALL get_qs_env(qs_env=qs_env, cell=cell)
328 1820 : cart_hmat(:, :) = cell%hmat(:, :)
329 140 : IF (cell%input_cell_canonicalized) cart_hmat(:, :) = cell%input_hmat(:, :)
330 :
331 140 : NULLIFY (gw_sec)
332 140 : gw_sec => section_vals_get_subs_vals(bs_sec, "GW")
333 140 : CALL section_vals_get(gw_sec, explicit=bs_env%do_gw)
334 140 : CALL section_vals_val_get(gw_sec, "COHSEX", i_val=bs_env%cohsex_flavor)
335 140 : bs_env%do_static_cohsex = (bs_env%cohsex_flavor /= cohsex_off)
336 140 : CALL section_vals_val_get(gw_sec, "RI_RS%_SECTION_PARAMETERS_", l_val=bs_env%do_gw_ri_rs)
337 :
338 140 : NULLIFY (soc_sec)
339 140 : soc_sec => section_vals_get_subs_vals(bs_sec, "SOC")
340 140 : CALL section_vals_get(soc_sec, explicit=bs_env%do_soc)
341 :
342 140 : CALL section_vals_val_get(soc_sec, "SOC_WINDOW_OCC", r_val=bs_env%soc_window_occ)
343 140 : CALL section_vals_val_get(soc_sec, "SOC_WINDOW_VIRT", r_val=bs_env%soc_window_virt)
344 140 : CALL section_vals_val_get(soc_sec, "SOC_WINDOW_SMEARING", r_val=bs_env%soc_window_smearing)
345 :
346 140 : CALL section_vals_val_get(soc_sec, "SOC_NBANDS_OCC", i_val=bs_env%soc_nbands_occ)
347 140 : CALL section_vals_val_get(soc_sec, "SOC_NBANDS_VIRT", i_val=bs_env%soc_nbands_virt)
348 :
349 140 : NULLIFY (dos_pdos_sec)
350 140 : dos_pdos_sec => section_vals_get_subs_vals(bs_sec, "DOS")
351 140 : CALL section_vals_get(dos_pdos_sec, explicit=bs_env%do_dos_pdos)
352 :
353 140 : CALL section_vals_val_get(bs_sec, "DOS%KPOINTS", i_vals=bs_env%nkp_grid_DOS_input)
354 140 : CALL section_vals_val_get(bs_sec, "DOS%ENERGY_WINDOW", r_val=bs_env%energy_window_DOS)
355 140 : CALL section_vals_val_get(bs_sec, "DOS%ENERGY_STEP", r_val=bs_env%energy_step_DOS)
356 140 : CALL section_vals_val_get(bs_sec, "DOS%BROADENING", r_val=bs_env%broadening_DOS)
357 :
358 140 : NULLIFY (ldos_sec)
359 140 : ldos_sec => section_vals_get_subs_vals(bs_sec, "DOS%LDOS")
360 140 : CALL section_vals_get(ldos_sec, explicit=bs_env%do_ldos)
361 :
362 140 : CALL section_vals_val_get(ldos_sec, "INTEGRATION", i_val=bs_env%int_ldos_xyz)
363 140 : CALL section_vals_val_get(ldos_sec, "BIN_MESH", i_vals=bs_env%bin_mesh)
364 :
365 140 : NULLIFY (kp_bs_sec)
366 140 : kp_bs_sec => section_vals_get_subs_vals(bs_sec, "BANDSTRUCTURE_PATH")
367 140 : CALL section_vals_val_get(kp_bs_sec, "NPOINTS", i_val=bs_env%input_kp_bs_npoints)
368 140 : CALL section_vals_val_get(kp_bs_sec, "UNITS", c_val=ustr)
369 140 : CALL uppercase(ustr)
370 140 : CALL section_vals_val_get(kp_bs_sec, "SPECIAL_POINT", n_rep_val=bs_env%input_kp_bs_n_sp_pts)
371 :
372 140 : NULLIFY (floquet_sec)
373 140 : floquet_sec => section_vals_get_subs_vals(bs_sec, "FLOQUET")
374 140 : CALL section_vals_get(floquet_sec, explicit=bs_env%do_floquet)
375 140 : CALL section_vals_val_get(floquet_sec, "AMPLITUDE", r_val=bs_env%floquet_amplitude)
376 140 : CALL section_vals_val_get(floquet_sec, "FREQUENCY", r_val=bs_env%floquet_omega)
377 140 : CALL section_vals_val_get(floquet_sec, "POLARISATION", r_vals=bs_env%floquet_polarisation)
378 140 : CALL section_vals_val_get(floquet_sec, "PHASE_OFFSETS", r_vals=bs_env%floquet_phi)
379 140 : CALL section_vals_val_get(floquet_sec, "MAX_FLOQUET_INDEX", i_val=bs_env%max_floquet_index)
380 140 : CALL section_vals_val_get(floquet_sec, "MAX_NMO", i_val=bs_env%max_nmo_floquet)
381 140 : CALL section_vals_val_get(floquet_sec, "EPS_FLOQUET", r_val=bs_env%eps_floquet)
382 140 : CALL section_vals_val_get(floquet_sec, "ENERGY_WINDOW", r_val=bs_env%energy_window_floquet)
383 140 : CALL section_vals_val_get(floquet_sec, "ENERGY_STEP", r_val=bs_env%energy_step_floquet)
384 140 : CALL section_vals_val_get(floquet_sec, "BROADENING", r_val=bs_env%broadening_floquet)
385 140 : CALL section_vals_val_get(floquet_sec, "MEM_FILL_FRACTION", r_val=bs_env%floquet_mem_fill_fraction)
386 140 : CALL section_vals_val_get(floquet_sec, "TEMPERATURE", r_val=bs_env%floquet_temperature)
387 140 : CALL section_vals_val_get(floquet_sec, "FLOQUET_DOS_FILE_NAME", c_val=bs_env%floquet_dos_file)
388 140 : CALL section_vals_val_get(floquet_sec, "QUASI_ENERGIES_FILE_NAME", c_val=bs_env%floquet_qe_file)
389 140 : 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 284 : ALLOCATE (bs_env%xkp_special(3, bs_env%input_kp_bs_n_sp_pts))
393 150 : DO ikp = 1, bs_env%input_kp_bs_n_sp_pts
394 10 : CALL section_vals_val_get(kp_bs_sec, "SPECIAL_POINT", i_rep_val=ikp, c_vals=string_ptr)
395 10 : CPASSERT(SIZE(string_ptr(:), 1) == 4)
396 40 : DO i = 1, 3
397 30 : CALL read_float_object(string_ptr(i + 1), kpptr(i), error_msg)
398 40 : IF (LEN_TRIM(error_msg) > 0) CPABORT(TRIM(error_msg))
399 : END DO
400 140 : SELECT CASE (ustr)
401 : CASE ("B_VECTOR")
402 40 : 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 0 : 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 0 : kpptr(3)*cart_hmat(3, 1:3))/twopi
411 : CASE DEFAULT
412 10 : CPABORT("Unknown unit <"//TRIM(ustr)//"> specified for k-point definition")
413 : END SELECT
414 : END DO
415 :
416 140 : CALL timestop(handle)
417 :
418 140 : END SUBROUTINE read_bandstructure_input_parameters
419 :
420 : ! **************************************************************************************************
421 : !> \brief ...
422 : !> \param bs_env ...
423 : ! **************************************************************************************************
424 140 : 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 140 : CALL timeset(routineN, handle)
433 :
434 140 : bs_env%unit_nr = cp_logger_get_default_io_unit()
435 :
436 140 : u = bs_env%unit_nr
437 :
438 140 : IF (u > 0) THEN
439 70 : WRITE (u, '(T2,A)') ' '
440 70 : WRITE (u, '(T2,A)') REPEAT('-', 79)
441 70 : WRITE (u, '(T2,A,A78)') '-', '-'
442 70 : WRITE (u, '(T2,A,A51,A27)') '-', 'BANDSTRUCTURE CALCULATION', '-'
443 70 : WRITE (u, '(T2,A,A78)') '-', '-'
444 70 : WRITE (u, '(T2,A)') REPEAT('-', 79)
445 70 : WRITE (u, '(T2,A)') ' '
446 : END IF
447 :
448 140 : CALL timestop(handle)
449 :
450 140 : END SUBROUTINE print_header
451 :
452 : ! **************************************************************************************************
453 : !> \brief ...
454 : !> \param qs_env ...
455 : !> \param bs_env ...
456 : !> \param kpoints ...
457 : ! **************************************************************************************************
458 122 : 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 122 : CALL timeset(routineN, handle)
473 :
474 : ! routine adapted from mp2_integrals.F
475 122 : NULLIFY (kpoints)
476 122 : CALL kpoint_create(kpoints)
477 :
478 122 : kpoints%kp_scheme = "GENERAL"
479 :
480 122 : n_special_kp = bs_env%input_kp_bs_n_sp_pts
481 122 : n_kp_in_line = bs_env%input_kp_bs_npoints
482 :
483 488 : periodic(1:3) = bs_env%periodic(1:3)
484 :
485 488 : DO i_dim = 1, 3
486 :
487 366 : CPASSERT(periodic(i_dim) == 0 .OR. periodic(i_dim) == 1)
488 :
489 488 : IF (bs_env%nkp_grid_DOS_input(i_dim) < 0) THEN
490 330 : IF (periodic(i_dim) == 1) nkp_grid(i_dim) = 2
491 330 : IF (periodic(i_dim) == 0) nkp_grid(i_dim) = 1
492 : ELSE
493 36 : 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 122 : IF (nkp_grid(1) > 1) THEN
500 8 : nkp_only_DOS = (nkp_grid(1) + 1)/2*nkp_grid(2)*nkp_grid(3)
501 114 : ELSE IF (nkp_grid(2) > 1) THEN
502 6 : nkp_only_DOS = nkp_grid(1)*(nkp_grid(2) + 1)/2*nkp_grid(3)
503 108 : ELSE IF (nkp_grid(3) > 1) THEN
504 2 : nkp_only_DOS = nkp_grid(1)*nkp_grid(2)*(nkp_grid(3) + 1)/2
505 : ELSE
506 106 : 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 122 : IF (n_special_kp > 0) THEN
512 0 : nkp_only_bs = n_kp_in_line*(n_special_kp - 1) + 1
513 : ELSE
514 : nkp_only_bs = 0
515 : END IF
516 :
517 122 : nkp = nkp_only_DOS + nkp_only_bs
518 :
519 488 : kpoints%nkp_grid(1:3) = nkp_grid(1:3)
520 122 : kpoints%nkp = nkp
521 :
522 122 : bs_env%nkp_bs_and_DOS = nkp
523 122 : bs_env%nkp_only_bs = nkp_only_bs
524 122 : bs_env%nkp_only_DOS = nkp_only_DOS
525 :
526 610 : ALLOCATE (kpoints%xkp(3, nkp), kpoints%wkp(nkp))
527 262 : kpoints%wkp(:) = 0.0_dp
528 262 : kpoints%wkp(1:nkp_only_DOS) = 1.0_dp/REAL(nkp_only_DOS, KIND=dp)
529 :
530 122 : CALL compute_xkp(kpoints%xkp, 1, nkp_only_DOS, nkp_grid)
531 :
532 122 : IF (n_special_kp > 0) THEN
533 0 : kpoints%xkp(1:3, nkp_only_DOS + 1) = bs_env%xkp_special(1:3, 1)
534 0 : ikk = nkp_only_DOS + 1
535 0 : DO i_special_kp = 2, n_special_kp
536 0 : DO i_kp_in_line = 1, n_kp_in_line
537 0 : 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 0 : bs_env%xkp_special(1:3, i_special_kp - 1))
542 0 : kpoints%wkp(ikk) = 0.0_dp
543 : END DO
544 : END DO
545 : END IF
546 :
547 122 : CALL kpoint_init_cell_index_simple(kpoints, qs_env)
548 :
549 122 : u = bs_env%unit_nr
550 :
551 122 : IF (u > 0 .AND. bs_env%do_periodic) THEN
552 8 : IF (nkp_only_bs > 0) THEN
553 : WRITE (u, FMT="(T2,1A,T77,I4)") &
554 0 : "Number of special k-points for the bandstructure", n_special_kp
555 0 : WRITE (u, FMT="(T2,1A,T77,I4)") "Number of k-points for the bandstructure", nkp
556 : WRITE (u, FMT="(T2,1A,T69,3I4)") &
557 0 : "K-point mesh for the density of states (DOS)", nkp_grid(1:3)
558 : ELSE
559 : WRITE (u, FMT="(T2,1A,T69,3I4)") &
560 8 : "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 122 : CALL timestop(handle)
565 :
566 122 : 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 18 : 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 18 : CALL timeset(routineN, handle)
590 :
591 18 : NULLIFY (kpoints)
592 18 : CALL kpoint_create(kpoints)
593 :
594 18 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints_scf)
595 :
596 72 : nkp_grid(1:3) = kpoints_scf%nkp_grid(1:3)
597 18 : 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 72 : DO i_dim = 1, 3
601 72 : IF (bs_env%periodic(i_dim) == 1) THEN
602 36 : CPASSERT(nkp_grid(i_dim) >= 4)
603 : END IF
604 : END DO
605 :
606 18 : kpoints%kp_scheme = "GENERAL"
607 72 : kpoints%nkp_grid(1:3) = nkp_grid(1:3)
608 18 : kpoints%nkp = nkp
609 18 : bs_env%nkp_scf_desymm = nkp
610 :
611 54 : ALLOCATE (kpoints%xkp(1:3, nkp))
612 18 : CALL compute_xkp(kpoints%xkp, 1, nkp, nkp_grid)
613 :
614 54 : ALLOCATE (kpoints%wkp(nkp))
615 306 : 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 72 : 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 72 : cixd(1:3) = cell_grid(1:3)/2
623 :
624 18 : nimages = cell_grid(1)*cell_grid(2)*cell_grid(3)
625 :
626 18 : bs_env%nimages_scf_desymm = nimages
627 72 : bs_env%cell_grid_scf_desymm(1:3) = cell_grid(1:3)
628 :
629 18 : IF (ASSOCIATED(kpoints%index_to_cell)) DEALLOCATE (kpoints%index_to_cell)
630 18 : IF (ASSOCIATED(kpoints%cell_to_index)) DEALLOCATE (kpoints%cell_to_index)
631 :
632 90 : ALLOCATE (kpoints%cell_to_index(-cixd(1):cixd(1), -cixd(2):cixd(2), -cixd(3):cixd(3)))
633 54 : ALLOCATE (kpoints%index_to_cell(3, nimages))
634 :
635 18 : img = 0
636 44 : DO i_cell_x = -cixd(1), cixd(1)
637 122 : DO j_cell_y = -cixd(2), cixd(2)
638 266 : DO k_cell_z = -cixd(3), cixd(3)
639 162 : img = img + 1
640 162 : kpoints%cell_to_index(i_cell_x, j_cell_y, k_cell_z) = img
641 726 : 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 18 : u = bs_env%unit_nr
647 18 : IF (u > 0 .AND. do_print) THEN
648 9 : WRITE (u, FMT="(T2,A,I49)") "Number of cells for G, χ, W, Σ", nimages
649 : END IF
650 :
651 18 : CALL timestop(handle)
652 :
653 18 : END SUBROUTINE setup_kpoints_scf_desymm
654 :
655 : ! **************************************************************************************************
656 : !> \brief ...
657 : !> \param bs_env ...
658 : !> \param kpoints ...
659 : ! **************************************************************************************************
660 18 : 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 18 : CALL timeset(routineN, handle)
672 :
673 : ! routine adapted from mp2_integrals.F
674 18 : NULLIFY (kpoints)
675 18 : CALL kpoint_create(kpoints)
676 :
677 18 : n_special_kp = bs_env%input_kp_bs_n_sp_pts
678 18 : n_kp_in_line = bs_env%input_kp_bs_npoints
679 18 : 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 18 : IF (n_special_kp > 0) THEN
684 4 : nkp_only_bs = n_kp_in_line*(n_special_kp - 1) + 1
685 : ELSE
686 : nkp_only_bs = 0
687 : END IF
688 18 : nkp = nkp_only_bs + nkp_scf_desymm
689 :
690 54 : ALLOCATE (kpoints%xkp(3, nkp))
691 54 : ALLOCATE (kpoints%wkp(nkp))
692 :
693 18 : kpoints%nkp = nkp
694 :
695 18 : bs_env%nkp_bs_and_DOS = nkp
696 18 : bs_env%nkp_only_bs = nkp_only_bs
697 18 : bs_env%nkp_only_DOS = nkp_scf_desymm
698 :
699 2340 : kpoints%xkp(1:3, 1:nkp_scf_desymm) = bs_env%kpoints_scf_desymm%xkp(1:3, 1:nkp_scf_desymm)
700 306 : kpoints%wkp(1:nkp_scf_desymm) = 1.0_dp/REAL(nkp_scf_desymm, KIND=dp)
701 :
702 18 : IF (n_special_kp > 0) THEN
703 32 : kpoints%xkp(1:3, nkp_scf_desymm + 1) = bs_env%xkp_special(1:3, 1)
704 4 : ikk = nkp_scf_desymm + 1
705 10 : DO i_special_kp = 2, n_special_kp
706 70 : DO i_kp_in_line = 1, n_kp_in_line
707 60 : 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 480 : bs_env%xkp_special(1:3, i_special_kp - 1))
712 66 : kpoints%wkp(ikk) = 0.0_dp
713 : END DO
714 : END DO
715 : END IF
716 :
717 18 : IF (ASSOCIATED(kpoints%index_to_cell)) DEALLOCATE (kpoints%index_to_cell)
718 :
719 54 : ALLOCATE (kpoints%index_to_cell(3, bs_env%nimages_scf_desymm))
720 1332 : kpoints%index_to_cell(:, :) = bs_env%kpoints_scf_desymm%index_to_cell(:, :)
721 :
722 18 : u = bs_env%unit_nr
723 :
724 18 : IF (u > 0) THEN
725 9 : WRITE (u, FMT="(T2,1A,T77,I4)") "Number of special k-points for the bandstructure", &
726 18 : n_special_kp
727 9 : WRITE (u, FMT="(T2,1A,T77,I4)") "Number of k-points for the bandstructure", nkp
728 : END IF
729 :
730 18 : CALL timestop(handle)
731 :
732 18 : END SUBROUTINE setup_kpoints_DOS_small_cell_full_kp
733 :
734 : ! **************************************************************************************************
735 : !> \brief ...
736 : !> \param qs_env ...
737 : !> \param bs_env ...
738 : ! **************************************************************************************************
739 18 : 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 18 : 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 18 : 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 18 : POINTER :: sab_nl
756 :
757 18 : 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 18 : kpoints=kpoints_scf)
763 :
764 18 : NULLIFY (sab_nl, lowdin_kp)
765 18 : IF (ALLOCATED(kpoints_scf%lowdin_v)) THEN
766 0 : IF (ANY(kpoints_scf%lowdin_v /= 0.0_dp)) ALLOCATE (lowdin_kp)
767 : END IF
768 18 : CALL get_kpoint_info(kpoints_scf, sab_nl=sab_nl, cell_to_index=cell_to_index_scf)
769 :
770 18 : CALL cp_cfm_create(cfm_ks, bs_env%cfm_work_mo%matrix_struct)
771 18 : CALL cp_cfm_create(cfm_s, bs_env%cfm_work_mo%matrix_struct)
772 18 : 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 18 : nkp_bs_and_DOS = bs_env%nkp_bs_and_DOS
776 :
777 18 : CALL allocate_GW_eigenvalues(bs_env)
778 :
779 442 : ALLOCATE (bs_env%cfm_mo_coeff_kp(nkp_bs_and_DOS, bs_env%n_spin))
780 442 : ALLOCATE (bs_env%cfm_ks_kp(nkp_bs_and_DOS, bs_env%n_spin))
781 406 : ALLOCATE (bs_env%cfm_s_kp(nkp_bs_and_DOS))
782 370 : DO ikp = 1, nkp_bs_and_DOS
783 704 : DO ispin = 1, bs_env%n_spin
784 352 : CALL cp_cfm_create(bs_env%cfm_mo_coeff_kp(ikp, ispin), bs_env%cfm_work_mo%matrix_struct)
785 704 : CALL cp_cfm_create(bs_env%cfm_ks_kp(ikp, ispin), bs_env%cfm_work_mo%matrix_struct)
786 : END DO
787 370 : CALL cp_cfm_create(bs_env%cfm_s_kp(ikp), bs_env%cfm_work_mo%matrix_struct)
788 : END DO
789 :
790 370 : DO ikp = 1, nkp_bs_and_DOS
791 :
792 1408 : 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 352 : CALL rsmat_to_kp(matrix_s, 1, xkp, cell_to_index_scf, sab_nl, bs_env, cfm_s)
796 352 : CALL cp_cfm_to_cfm(cfm_s, bs_env%cfm_s_kp(ikp))
797 352 : IF (ASSOCIATED(lowdin_kp)) THEN
798 0 : CALL lowdin_kp_prepare(lowdin_kp, kpoints_scf%lowdin_orbitals, .FALSE., cfm_s)
799 : END IF
800 722 : DO ispin = 1, bs_env%n_spin
801 : ! h^KS^R -> h^KS(k)
802 352 : CALL rsmat_to_kp(matrix_ks, ispin, xkp, cell_to_index_scf, sab_nl, bs_env, cfm_ks)
803 352 : IF (ASSOCIATED(lowdin_kp)) THEN
804 0 : 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 352 : CALL cp_cfm_to_cfm(cfm_ks, bs_env%cfm_ks_kp(ikp, ispin))
810 352 : 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 352 : nmo_retained=nmo_retained)
818 352 : 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 704 : CALL cp_cfm_to_cfm(cfm_mos, bs_env%cfm_mo_coeff_kp(ikp, ispin))
823 :
824 : END DO
825 : END DO
826 :
827 36 : DO ispin = 1, bs_env%n_spin
828 370 : VBM = MAXVAL(bs_env%eigenval_scf(bs_env%n_occ(ispin), :, ispin))
829 370 : CBM = MINVAL(bs_env%eigenval_scf(bs_env%n_occ(ispin) + 1, :, ispin))
830 :
831 36 : bs_env%e_fermi(ispin) = 0.5_dp*(VBM + CBM)
832 :
833 : END DO
834 :
835 18 : CALL get_VBM_CBM_bandgaps(bs_env%band_edges_scf, bs_env%eigenval_scf, bs_env)
836 :
837 18 : CALL cp_cfm_release(cfm_ks)
838 18 : CALL cp_cfm_release(cfm_s)
839 18 : CALL cp_cfm_release(cfm_mos)
840 18 : IF (ASSOCIATED(lowdin_kp)) THEN
841 0 : CALL cp_cfm_release(lowdin_kp%lowdin_cb)
842 0 : CALL cp_cfm_release(lowdin_kp%lowdin_cwork)
843 0 : CALL lowdin_kp%lowdin_gemm%destroy()
844 0 : DEALLOCATE (lowdin_kp)
845 : END IF
846 :
847 18 : CALL timestop(handle)
848 :
849 36 : 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 1976 : 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 1976 : CALL timeset(routineN, handle)
880 :
881 1976 : ALLOCATE (rmat, cmat, nsmat)
882 :
883 1976 : imag_rs_mat_private = .FALSE.
884 1976 : IF (PRESENT(imag_rs_mat)) imag_rs_mat_private = imag_rs_mat
885 :
886 954 : IF (imag_rs_mat_private) THEN
887 954 : CALL dbcsr_create(rmat, template=mat_rs(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
888 954 : CALL dbcsr_create(cmat, template=mat_rs(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
889 : ELSE
890 1022 : CALL dbcsr_create(rmat, template=mat_rs(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
891 1022 : CALL dbcsr_create(cmat, template=mat_rs(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
892 : END IF
893 1976 : CALL dbcsr_create(nsmat, template=mat_rs(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
894 1976 : CALL cp_dbcsr_alloc_block_from_nbl(rmat, sab_nl)
895 1976 : CALL cp_dbcsr_alloc_block_from_nbl(cmat, sab_nl)
896 :
897 1976 : CALL dbcsr_set(rmat, 0.0_dp)
898 1976 : CALL dbcsr_set(cmat, 0.0_dp)
899 : CALL rskp_transform(rmatrix=rmat, cmatrix=cmat, rsmat=mat_rs, ispin=ispin, &
900 1976 : xkp=xkp, cell_to_index=cell_to_index_scf, sab_nl=sab_nl)
901 :
902 1976 : CALL dbcsr_desymmetrize(rmat, nsmat)
903 1976 : CALL copy_dbcsr_to_fm(nsmat, bs_env%fm_work_mo(1))
904 1976 : CALL dbcsr_desymmetrize(cmat, nsmat)
905 1976 : CALL copy_dbcsr_to_fm(nsmat, bs_env%fm_work_mo(2))
906 1976 : CALL cp_fm_to_cfm(bs_env%fm_work_mo(1), bs_env%fm_work_mo(2), cfm_kp)
907 :
908 1976 : CALL dbcsr_deallocate_matrix(rmat)
909 1976 : CALL dbcsr_deallocate_matrix(cmat)
910 1976 : CALL dbcsr_deallocate_matrix(nsmat)
911 :
912 1976 : CALL timestop(handle)
913 :
914 1976 : END SUBROUTINE rsmat_to_kp
915 :
916 : ! **************************************************************************************************
917 : !> \brief ...
918 : !> \param bs_env ...
919 : ! **************************************************************************************************
920 122 : 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 122 : CALL timeset(routineN, handle)
929 :
930 488 : ALLOCATE (bs_env%eigenval_scf_Gamma(bs_env%n_ao, bs_env%n_spin))
931 :
932 262 : DO ispin = 1, bs_env%n_spin
933 :
934 : ! use work matrices because the matrices are overwritten in cp_fm_geeig_canon
935 140 : CALL cp_fm_to_fm(bs_env%fm_ks_Gamma(ispin), bs_env%fm_work_mo(1))
936 140 : 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 140 : nmo_retained=nmo_retained)
947 140 : bs_env%n_mo_retained = MIN(bs_env%n_mo_retained, nmo_retained)
948 :
949 140 : VBM = bs_env%eigenval_scf_Gamma(bs_env%n_occ(ispin), ispin)
950 140 : CBM = bs_env%eigenval_scf_Gamma(bs_env%n_occ(ispin) + 1, ispin)
951 :
952 262 : bs_env%e_fermi(ispin) = 0.5_dp*(VBM + CBM)
953 :
954 : END DO
955 :
956 122 : CALL timestop(handle)
957 :
958 : ! Gamma-only path for molecules: eigenval_scf is filled here from the Gamma eigenvalues
959 262 : DO ispin = 1, bs_env%n_spin
960 424 : DO ikp = 1, bs_env%nkp_bs_and_DOS
961 2280 : bs_env%eigenval_scf(:, ikp, ispin) = bs_env%eigenval_scf_Gamma(:, ispin)
962 : END DO
963 : END DO
964 :
965 122 : END SUBROUTINE diagonalize_ks_matrix
966 :
967 : ! **************************************************************************************************
968 : !> \brief ...
969 : !> \param bs_env ...
970 : !> \param qs_env ...
971 : ! **************************************************************************************************
972 122 : 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 122 : CALL timeset(routineN, handle)
982 :
983 262 : 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 140 : ikp, qs_env, bs_env%kpoints_DOS, "ORB")
988 :
989 : ! check whether S_µν(k_i) is positive definite
990 140 : 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 262 : IF (info /= 0) THEN
995 0 : u = bs_env%unit_nr
996 :
997 0 : IF (u > 0) THEN
998 0 : WRITE (u, FMT="(T2,A)") ""
999 : WRITE (u, FMT="(T2,A)") "ERROR: The Cholesky decomposition "// &
1000 0 : "of the k-point overlap matrix failed. This is"
1001 : WRITE (u, FMT="(T2,A)") "because the algorithm is "// &
1002 0 : "only correct in the limit of large cells. The cell of "
1003 : WRITE (u, FMT="(T2,A)") "the calculation is too small. "// &
1004 0 : "Use MULTIPLE_UNIT_CELL to create a larger cell "
1005 0 : WRITE (u, FMT="(T2,A)") "and to prevent this error."
1006 : END IF
1007 :
1008 0 : CALL bs_env%para_env%sync()
1009 0 : CPABORT("Please see information on the error above.")
1010 :
1011 : END IF ! Cholesky decomposition failed
1012 :
1013 : END DO ! ikp
1014 :
1015 122 : CALL cp_cfm_release(cfm_s_ikp)
1016 :
1017 122 : CALL timestop(handle)
1018 :
1019 122 : END SUBROUTINE check_positive_definite_overlap_mat
1020 :
1021 : ! **************************************************************************************************
1022 : !> \brief ...
1023 : !> \param qs_env ...
1024 : !> \param bs_env ...
1025 : ! **************************************************************************************************
1026 280 : 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 140 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1040 : TYPE(mp_para_env_type), POINTER :: para_env
1041 140 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1042 : TYPE(scf_control_type), POINTER :: scf_control
1043 : TYPE(section_vals_type), POINTER :: input
1044 :
1045 140 : CALL timeset(routineN, handle)
1046 :
1047 : CALL get_qs_env(qs_env, &
1048 : dft_control=dft_control, &
1049 : scf_control=scf_control, &
1050 140 : mos=mos)
1051 :
1052 140 : bs_env%n_spin = dft_control%nspins
1053 140 : IF (bs_env%n_spin == 1) bs_env%spin_degeneracy = 2.0_dp
1054 140 : IF (bs_env%n_spin == 2) bs_env%spin_degeneracy = 1.0_dp
1055 :
1056 140 : CALL get_mo_set(mo_set=mos(1), nao=n_ao, homo=homo, nelectron=nelectron, maxocc=maxocc)
1057 140 : bs_env%n_ao = n_ao
1058 420 : bs_env%n_occ_nosmear(1:2) = NINT(nelectron/maxocc) !fix for moire large-cell Gamma
1059 420 : bs_env%n_electron(1:2) = nelectron
1060 : ! lowered by the diagonalization below if the canonical orthogonalization removes modes
1061 140 : bs_env%n_mo_retained = n_ao
1062 420 : bs_env%n_occ(1:2) = homo
1063 420 : bs_env%n_vir(1:2) = n_ao - homo
1064 140 : IF (bs_env%n_spin == 2) THEN
1065 18 : CALL get_mo_set(mo_set=mos(2), homo=homo, nelectron=nelectron, maxocc=maxocc)
1066 18 : bs_env%n_occ_nosmear(2) = NINT(nelectron/maxocc)
1067 18 : bs_env%n_occ(2) = homo
1068 18 : bs_env%n_vir(2) = n_ao - homo
1069 18 : bs_env%n_electron(2) = nelectron
1070 : END IF
1071 140 : 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 140 : bs_env%do_smear = .FALSE.
1076 140 : bs_env%smear_method = smear_fermi_dirac
1077 140 : bs_env%smear_width = 0.0_dp
1078 140 : IF (scf_control%smear%do_smear) THEN
1079 4 : SELECT CASE (scf_control%smear%method)
1080 : CASE (smear_fermi_dirac)
1081 2 : bs_env%do_smear = .TRUE.
1082 2 : bs_env%smear_method = smear_fermi_dirac
1083 2 : bs_env%smear_width = scf_control%smear%electronic_temperature
1084 : CASE (smear_gaussian, smear_mp, smear_mv)
1085 0 : bs_env%do_smear = .TRUE.
1086 0 : bs_env%smear_method = scf_control%smear%method
1087 0 : bs_env%smear_width = scf_control%smear%smearing_width
1088 : CASE DEFAULT
1089 2 : 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 0 : 'falling back to a T=0 bisection.'
1093 : END IF
1094 : END SELECT
1095 : END IF
1096 140 : 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 140 : CALL get_qs_env(qs_env, para_env=para_env)
1099 140 : color_sub = 0
1100 140 : ALLOCATE (bs_env%para_env)
1101 140 : CALL bs_env%para_env%from_split(para_env, color_sub)
1102 :
1103 140 : CALL get_qs_env(qs_env, particle_set=particle_set)
1104 :
1105 140 : n_atom = SIZE(particle_set)
1106 140 : bs_env%n_atom = n_atom
1107 :
1108 140 : CALL get_qs_env(qs_env=qs_env, cell=cell)
1109 140 : CALL get_cell(cell=cell, periodic=periodic, h=hmat)
1110 560 : bs_env%periodic(1:3) = periodic(1:3)
1111 480 : bs_env%do_periodic = ANY(bs_env%periodic /= 0)
1112 1820 : bs_env%hmat(1:3, 1:3) = hmat
1113 140 : bs_env%nimages_scf = dft_control%nimages
1114 140 : IF (dft_control%nimages == 1) THEN
1115 122 : IF (bs_env%do_gw_ri_rs) THEN
1116 40 : IF (bs_env%do_periodic) THEN
1117 0 : bs_env%gw_implementation = ri_rs_large_cell_Gamma
1118 : ELSE
1119 40 : bs_env%gw_implementation = ri_rs_non_periodic
1120 : END IF
1121 82 : ELSE IF (bs_env%do_static_cohsex) THEN
1122 8 : bs_env%gw_implementation = tensor_static_cohsex
1123 : ELSE
1124 74 : bs_env%gw_implementation = tensor_large_cell_Gamma
1125 : END IF
1126 18 : ELSE IF (dft_control%nimages > 1) THEN
1127 18 : IF (bs_env%do_gw_ri_rs) THEN
1128 0 : CPABORT("RI-RS Not Implemented for K-point Calculations")
1129 18 : ELSE IF (bs_env%do_static_cohsex) THEN
1130 0 : CPABORT("COHSEX Not Implemented for K-point Calculations")
1131 : ELSE
1132 18 : bs_env%gw_implementation = tensor_small_cell_full_kp
1133 : END IF
1134 : ELSE
1135 0 : CPABORT("Wrong number of cells from DFT calculation.")
1136 : END IF
1137 :
1138 140 : u = bs_env%unit_nr
1139 :
1140 : ! Marek : Get and save the rtp method
1141 140 : CALL get_qs_env(qs_env=qs_env, input=input)
1142 140 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%_SECTION_PARAMETERS_", i_val=bs_env%rtp_method)
1143 :
1144 140 : IF (u > 0) THEN
1145 70 : WRITE (u, FMT="(T2,2A,T73,I8)") "Number of occupied molecular orbitals (MOs) ", &
1146 140 : "= Number of occupied bands", homo
1147 70 : WRITE (u, FMT="(T2,2A,T73,I8)") "Number of unoccupied (= virtual) MOs ", &
1148 140 : "= Number of unoccupied bands", n_ao - homo
1149 70 : WRITE (u, FMT="(T2,A,T73,I8)") "Number of Gaussian basis functions for MOs", n_ao
1150 70 : IF (bs_env%gw_implementation == tensor_small_cell_full_kp) THEN
1151 9 : WRITE (u, FMT="(T2,2A,T73,I8)") "Number of cells considered in the DFT ", &
1152 18 : "calculation", bs_env%nimages_scf
1153 : END IF
1154 : END IF
1155 :
1156 140 : CALL timestop(handle)
1157 :
1158 140 : END SUBROUTINE get_parameters_from_qs_env
1159 :
1160 : ! **************************************************************************************************
1161 : !> \brief ...
1162 : !> \param bs_env ...
1163 : ! **************************************************************************************************
1164 140 : 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 140 : CALL timeset(routineN, handle)
1172 :
1173 140 : bs_env%n_bins_max_for_printing = 5000
1174 :
1175 140 : CALL timestop(handle)
1176 :
1177 140 : END SUBROUTINE set_heuristic_parameters
1178 :
1179 : ! **************************************************************************************************
1180 : !> \brief ...
1181 : !> \param qs_env ...
1182 : !> \param bs_env ...
1183 : ! **************************************************************************************************
1184 140 : 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 140 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
1194 : TYPE(mp_para_env_type), POINTER :: para_env
1195 :
1196 140 : 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 140 : matrix_s_kp=matrix_s)
1203 :
1204 140 : NULLIFY (fm_struct)
1205 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=bs_env%n_ao, &
1206 140 : ncol_global=bs_env%n_ao, para_env=para_env)
1207 :
1208 700 : DO i_work = 1, SIZE(bs_env%fm_work_mo)
1209 700 : CALL cp_fm_create(bs_env%fm_work_mo(i_work), fm_struct)
1210 : END DO
1211 :
1212 140 : CALL cp_cfm_create(bs_env%cfm_work_mo, fm_struct)
1213 140 : CALL cp_cfm_create(bs_env%cfm_work_mo_2, fm_struct)
1214 :
1215 140 : CALL cp_fm_create(bs_env%fm_s_Gamma, fm_struct)
1216 140 : CALL copy_dbcsr_to_fm(matrix_s(1, 1)%matrix, bs_env%fm_s_Gamma)
1217 :
1218 298 : DO ispin = 1, bs_env%n_spin
1219 158 : CALL cp_fm_create(bs_env%fm_ks_Gamma(ispin), fm_struct)
1220 158 : CALL copy_dbcsr_to_fm(matrix_ks(ispin, 1)%matrix, bs_env%fm_ks_Gamma(ispin))
1221 298 : CALL cp_fm_create(bs_env%fm_mo_coeff_Gamma(ispin), fm_struct)
1222 : END DO
1223 :
1224 140 : CALL cp_fm_struct_release(fm_struct)
1225 :
1226 140 : NULLIFY (bs_env%mat_ao_ao%matrix)
1227 140 : ALLOCATE (bs_env%mat_ao_ao%matrix)
1228 : CALL dbcsr_create(bs_env%mat_ao_ao%matrix, template=matrix_s(1, 1)%matrix, &
1229 140 : matrix_type=dbcsr_type_no_symmetry)
1230 :
1231 700 : ALLOCATE (bs_env%eigenval_scf(bs_env%n_ao, bs_env%nkp_bs_and_DOS, bs_env%n_spin))
1232 :
1233 140 : CALL timestop(handle)
1234 :
1235 140 : END SUBROUTINE allocate_and_fill_fm_ks_fm_s
1236 :
1237 : ! **************************************************************************************************
1238 : !> \brief ...
1239 : !> \param qs_env ...
1240 : !> \param bs_env ...
1241 : ! **************************************************************************************************
1242 140 : 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 140 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: DOS_G0W0, DOS_G0W0_SOC, DOS_scf, DOS_scf_SOC, &
1257 140 : eigenval, eigenval_spinor, eigenval_spinor_G0W0, eigenval_spinor_no_SOC
1258 140 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenval_spinor_all_gw, &
1259 140 : eigenval_spinor_all_scf, PDOS_G0W0, PDOS_G0W0_SOC, PDOS_scf, PDOS_scf_SOC, proj_mo_on_kind
1260 140 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: LDOS_G0W0_2d, LDOS_scf_2d, &
1261 140 : 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 420 : TYPE(cp_cfm_type), DIMENSION(2) :: cfm_mos_ikp
1269 :
1270 140 : CALL timeset(routineN, handle)
1271 :
1272 140 : n_ao = bs_env%n_ao
1273 :
1274 140 : energy_window_DOS = bs_env%energy_window_DOS
1275 140 : energy_step_DOS = bs_env%energy_step_DOS
1276 140 : broadening = bs_env%broadening_DOS
1277 :
1278 : ! if we have done GW or a full kpoint SCF, we already have the band edges
1279 140 : IF (bs_env%do_gw .OR. &
1280 : bs_env%gw_implementation == tensor_small_cell_full_kp) THEN
1281 140 : band_edges_scf = bs_env%band_edges_scf
1282 140 : band_edges_scf_guess = band_edges_scf
1283 : ELSE
1284 :
1285 0 : IF (bs_env%n_spin == 1) THEN
1286 0 : homo = bs_env%n_occ(1)
1287 0 : band_edges_scf_guess%VBM = bs_env%eigenval_scf_Gamma(homo, 1)
1288 0 : band_edges_scf_guess%CBM = bs_env%eigenval_scf_Gamma(homo + 1, 1)
1289 : ELSE
1290 0 : homo_1 = bs_env%n_occ(1)
1291 0 : homo_2 = bs_env%n_occ(2)
1292 : band_edges_scf_guess%VBM = MAX(bs_env%eigenval_scf_Gamma(homo_1, 1), &
1293 0 : 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 0 : bs_env%eigenval_scf_Gamma(homo_2 + 1, 2))
1296 : END IF
1297 :
1298 : ! initialization
1299 0 : band_edges_scf%VBM = -1000.0_dp
1300 0 : band_edges_scf%CBM = 1000.0_dp
1301 0 : band_edges_scf%DBG = 1000.0_dp
1302 : END IF
1303 :
1304 140 : E_min = band_edges_scf_guess%VBM - 0.5_dp*energy_window_DOS
1305 140 : E_max = band_edges_scf_guess%CBM + 0.5_dp*energy_window_DOS
1306 :
1307 140 : IF (bs_env%do_gw) THEN
1308 138 : band_edges_G0W0 = bs_env%band_edges_GW
1309 138 : E_min_G0W0 = band_edges_G0W0%VBM - 0.5_dp*energy_window_DOS
1310 138 : E_max_G0W0 = band_edges_G0W0%CBM + 0.5_dp*energy_window_DOS
1311 138 : E_min = MIN(E_min, E_min_G0W0)
1312 138 : E_max = MAX(E_max, E_max_G0W0)
1313 : END IF
1314 :
1315 140 : E_total_window = E_max - E_min
1316 :
1317 140 : n_E = INT(E_total_window/energy_step_DOS)
1318 :
1319 140 : CALL get_qs_env(qs_env, nkind=nkind)
1320 :
1321 560 : ALLOCATE (proj_mo_on_kind(n_ao, nkind))
1322 140 : proj_mo_on_kind(:, :) = 0.0_dp
1323 :
1324 420 : ALLOCATE (eigenval(n_ao))
1325 420 : ALLOCATE (eigenval_spinor(2*n_ao))
1326 280 : ALLOCATE (eigenval_spinor_no_SOC(2*n_ao))
1327 280 : 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 140 : IF (bs_env%do_soc) THEN
1333 88 : ALLOCATE (eigenval_spinor_all_scf(2*n_ao, bs_env%nkp_bs_and_DOS))
1334 22 : IF (bs_env%do_gw) THEN
1335 66 : ALLOCATE (eigenval_spinor_all_gw(2*n_ao, bs_env%nkp_bs_and_DOS))
1336 : END IF
1337 : END IF
1338 :
1339 140 : IF (bs_env%do_dos_pdos) THEN
1340 :
1341 78 : ALLOCATE (DOS_scf(n_E))
1342 26 : DOS_scf(:) = 0.0_dp
1343 104 : ALLOCATE (PDOS_scf(n_E, nkind))
1344 26 : PDOS_scf(:, :) = 0.0_dp
1345 :
1346 26 : IF (bs_env%do_soc) THEN
1347 :
1348 32 : ALLOCATE (DOS_scf_SOC(n_E))
1349 16 : DOS_scf_SOC(:) = 0.0_dp
1350 48 : ALLOCATE (PDOS_scf_SOC(n_E, nkind))
1351 16 : PDOS_scf_SOC(:, :) = 0.0_dp
1352 :
1353 : END IF
1354 :
1355 26 : IF (bs_env%do_gw) THEN
1356 :
1357 52 : ALLOCATE (DOS_G0W0(n_E))
1358 26 : DOS_G0W0(:) = 0.0_dp
1359 78 : ALLOCATE (PDOS_G0W0(n_E, nkind))
1360 26 : PDOS_G0W0(:, :) = 0.0_dp
1361 :
1362 26 : IF (bs_env%do_soc) THEN
1363 :
1364 32 : ALLOCATE (DOS_G0W0_SOC(n_E))
1365 16 : DOS_G0W0_SOC(:) = 0.0_dp
1366 48 : ALLOCATE (PDOS_G0W0_SOC(n_E, nkind))
1367 16 : PDOS_G0W0_SOC(:, :) = 0.0_dp
1368 :
1369 : END IF
1370 : END IF
1371 : END IF
1372 :
1373 140 : CALL cp_cfm_create(cfm_mos_ikp(1), bs_env%fm_ks_Gamma(1)%matrix_struct)
1374 140 : CALL cp_cfm_create(cfm_mos_ikp(2), bs_env%fm_ks_Gamma(1)%matrix_struct)
1375 140 : CALL cp_cfm_create(cfm_work_ikp, bs_env%fm_ks_Gamma(1)%matrix_struct)
1376 140 : CALL cp_cfm_create(cfm_s_ikp_copy, bs_env%fm_ks_Gamma(1)%matrix_struct)
1377 :
1378 140 : IF (bs_env%do_soc) THEN
1379 :
1380 22 : CALL cp_cfm_create(cfm_mos_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1381 22 : CALL cp_cfm_create(cfm_work_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1382 22 : CALL cp_cfm_create(cfm_s_ikp_spinor_copy, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1383 22 : CALL cp_cfm_create(cfm_ks_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1384 22 : CALL cp_cfm_create(cfm_SOC_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1385 22 : CALL cp_cfm_create(cfm_s_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1386 22 : CALL cp_cfm_create(cfm_spinor_wf_ikp, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1387 :
1388 22 : homo_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin)
1389 :
1390 18 : SELECT CASE (bs_env%n_spin)
1391 : CASE (1)
1392 18 : n_el_spinor = REAL(bs_env%n_electron(1), KIND=dp)
1393 : CASE (2)
1394 22 : n_el_spinor = REAL(bs_env%n_electron(1) + bs_env%n_electron(2), KIND=dp)
1395 : END SELECT
1396 :
1397 22 : band_edges_scf_SOC%VBM = -1000.0_dp
1398 22 : band_edges_scf_SOC%CBM = 1000.0_dp
1399 22 : band_edges_scf_SOC%DBG = 1000.0_dp
1400 :
1401 22 : IF (bs_env%do_gw) THEN
1402 22 : band_edges_G0W0_SOC%VBM = -1000.0_dp
1403 22 : band_edges_G0W0_SOC%CBM = 1000.0_dp
1404 22 : band_edges_G0W0_SOC%DBG = 1000.0_dp
1405 : END IF
1406 :
1407 22 : IF (bs_env%unit_nr > 0) THEN
1408 11 : WRITE (bs_env%unit_nr, '(A)') ''
1409 11 : IF (bs_env%soc_window_occ > 0.0_dp) THEN
1410 3 : WRITE (bs_env%unit_nr, '(T2,A,T71,F10.2)') 'SOC requested, SOC energy window occ (eV):', &
1411 6 : bs_env%soc_window_occ*evolt
1412 : ELSE
1413 8 : WRITE (bs_env%unit_nr, '(T2,A,T71,A10)') 'SOC requested, SOC energy window occ (eV):', &
1414 16 : ' no window'
1415 : END IF
1416 11 : IF (bs_env%soc_window_virt > 0.0_dp) THEN
1417 3 : WRITE (bs_env%unit_nr, '(T2,A,T71,F10.2)') 'SOC requested, SOC energy window virt (eV):', &
1418 6 : bs_env%soc_window_virt*evolt
1419 : ELSE
1420 8 : WRITE (bs_env%unit_nr, '(T2,A,T71,A10)') 'SOC requested, SOC energy window virt (eV):', &
1421 16 : ' no window'
1422 : END IF
1423 11 : IF (bs_env%soc_window_occ > 0.0_dp .OR. bs_env%soc_window_virt > 0.0_dp) THEN
1424 4 : WRITE (bs_env%unit_nr, '(T2,A,T71,F10.2)') 'SOC requested, SOC window smearing (eV):', &
1425 8 : bs_env%soc_window_smearing*evolt
1426 : END IF
1427 : END IF
1428 : END IF
1429 :
1430 140 : IF (bs_env%do_ldos) THEN
1431 2 : CPASSERT(bs_env%int_ldos_xyz == int_ldos_z)
1432 : END IF
1433 :
1434 140 : IF (bs_env%gw_implementation == tensor_small_cell_full_kp) THEN
1435 18 : CALL cp_cfm_create(cfm_ks_ikp, bs_env%cfm_ks_kp(1, 1)%matrix_struct)
1436 18 : CALL cp_cfm_create(cfm_s_ikp, bs_env%cfm_ks_kp(1, 1)%matrix_struct)
1437 : END IF
1438 :
1439 632 : DO ikp = 1, bs_env%nkp_bs_and_DOS
1440 :
1441 492 : t1 = m_walltime()
1442 :
1443 1006 : DO ispin = 1, bs_env%n_spin
1444 :
1445 676 : 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 162 : 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 162 : ikp, qs_env, bs_env%kpoints_DOS, "ORB")
1455 162 : 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 162 : eigenval, cfm_work_ikp)
1460 :
1461 2140 : bs_env%eigenval_scf(:, ikp, ispin) = eigenval(:)
1462 :
1463 : CASE (tensor_small_cell_full_kp)
1464 :
1465 : ! 1. get H^KS_µν(k_i)
1466 352 : CALL cp_cfm_to_cfm(bs_env%cfm_ks_kp(ikp, ispin), cfm_ks_ikp)
1467 :
1468 : ! 2. get S_µν(k_i)
1469 352 : 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 352 : CALL cp_cfm_to_cfm(bs_env%cfm_mo_coeff_kp(ikp, ispin), cfm_mos_ikp(ispin))
1473 5120 : 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 514 : 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 514 : 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 250 : proj_mo_on_kind)
1485 :
1486 250 : IF (bs_env%do_gw) THEN
1487 : CALL add_to_DOS_PDOS(DOS_G0W0, PDOS_G0W0, bs_env%eigenval_GW(:, ikp, ispin), &
1488 250 : ikp, bs_env, n_E, E_min, proj_mo_on_kind)
1489 : END IF
1490 : END IF
1491 :
1492 514 : 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 2 : eigenval(:), band_edges_scf_guess)
1495 :
1496 2 : 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 2 : bs_env%eigenval_GW(:, ikp, 1), band_edges_G0W0)
1499 : END IF
1500 :
1501 : END IF
1502 :
1503 514 : homo = bs_env%n_occ(ispin)
1504 :
1505 514 : band_edges_scf%VBM = MAX(band_edges_scf%VBM, eigenval(homo))
1506 514 : band_edges_scf%CBM = MIN(band_edges_scf%CBM, eigenval(homo + 1))
1507 1006 : 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 492 : IF (bs_env%do_soc) THEN
1513 :
1514 : ! only print eigenvalues of DOS k-points in case no bandstructure path has been given
1515 328 : print_DOS_kpoints = (bs_env%nkp_only_bs <= 0)
1516 : ! in kpoints_DOS, the last nkp_only_bs are bandstructure k-points
1517 328 : is_bandstruc_kpoint = (ikp > bs_env%nkp_only_DOS)
1518 328 : print_ikp = print_DOS_kpoints .OR. is_bandstruc_kpoint
1519 :
1520 328 : IF (print_DOS_kpoints) THEN
1521 234 : nkp = bs_env%nkp_only_DOS
1522 234 : ikp_for_file = ikp
1523 : ELSE
1524 94 : nkp = bs_env%nkp_only_bs
1525 94 : 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 328 : band_edges_scf_SOC, eigenval_spinor, cfm_spinor_wf_ikp)
1532 8472 : eigenval_spinor_all_scf(:, ikp) = eigenval_spinor
1533 328 : IF (.NOT. bs_env%do_gw .AND. print_ikp) THEN
1534 0 : CALL write_SOC_eigenvalues(eigenval_spinor, ikp_for_file, ikp, bs_env)
1535 : END IF
1536 :
1537 328 : 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 2 : eigenval_spinor, band_edges_scf_guess, .TRUE., cfm_work_ikp)
1540 : END IF
1541 :
1542 328 : 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 328 : band_edges_G0W0_SOC, eigenval_spinor_G0W0, cfm_spinor_wf_ikp)
1548 8472 : eigenval_spinor_all_gw(:, ikp) = eigenval_spinor_G0W0
1549 328 : 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 296 : eigenval_spinor_G0W0)
1554 : END IF
1555 :
1556 : END IF ! do_gw
1557 :
1558 : END IF ! do_soc
1559 :
1560 632 : 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 0 : 'Compute DOS, LDOS for k-point ', ikp, ' /', bs_env%nkp_bs_and_DOS, &
1563 0 : ', Execution time', m_walltime() - t1, ' s'
1564 : END IF
1565 :
1566 : END DO ! ikp_DOS
1567 :
1568 140 : band_edges_scf%IDBG = band_edges_scf%CBM - band_edges_scf%VBM
1569 :
1570 262 : 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 122 : 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 140 : band_edges_scf%VBM, band_edges_scf%CBM, band_edges_scf%E_fermi)
1580 : END SELECT
1581 :
1582 140 : IF (bs_env%do_soc) THEN
1583 22 : 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 22 : band_edges_scf_SOC%VBM, band_edges_scf_SOC%CBM, band_edges_scf_SOC%E_fermi)
1590 :
1591 22 : IF (bs_env%do_gw) THEN
1592 22 : 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 22 : 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 140 : CALL write_band_edges(band_edges_scf, "SCF", bs_env)
1601 140 : IF (bs_env%do_dos_pdos) THEN
1602 26 : CALL write_dos_pdos(DOS_scf, PDOS_scf, bs_env, qs_env, "SCF", E_min, band_edges_scf%VBM)
1603 : END IF
1604 140 : IF (bs_env%do_ldos) THEN
1605 2 : CALL print_LDOS_main(LDOS_scf_2d, bs_env, band_edges_scf, "SCF")
1606 : END IF
1607 :
1608 140 : IF (bs_env%do_soc) THEN
1609 22 : CALL write_band_edges(band_edges_scf_SOC, "SCF+SOC", bs_env)
1610 22 : 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 16 : E_min, band_edges_scf_SOC%VBM)
1613 : END IF
1614 22 : 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 2 : "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 140 : gw_label = gw_flavour_label(bs_env)
1625 140 : gw_soc_label = TRIM(gw_label)//"+SOC"
1626 :
1627 140 : IF (bs_env%do_gw) THEN
1628 138 : CALL write_band_edges(band_edges_G0W0, TRIM(gw_label), bs_env)
1629 138 : CALL write_band_edges(bs_env%band_edges_HF, "Hartree-Fock with SCF orbitals", bs_env)
1630 138 : IF (bs_env%do_dos_pdos) THEN
1631 : CALL write_dos_pdos(DOS_G0W0, PDOS_G0W0, bs_env, qs_env, "G0W0", E_min, &
1632 26 : band_edges_G0W0%VBM)
1633 : END IF
1634 138 : IF (bs_env%do_ldos) THEN
1635 2 : CALL print_LDOS_main(LDOS_G0W0_2d, bs_env, band_edges_G0W0, "G0W0")
1636 : END IF
1637 : END IF
1638 :
1639 140 : IF (bs_env%do_soc .AND. bs_env%do_gw) THEN
1640 22 : CALL write_band_edges(band_edges_G0W0_SOC, TRIM(gw_soc_label), bs_env)
1641 22 : 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 16 : band_edges_G0W0_SOC%VBM)
1644 : END IF
1645 : END IF
1646 :
1647 140 : CALL cp_cfm_release(cfm_s_ikp)
1648 140 : CALL cp_cfm_release(cfm_ks_ikp)
1649 140 : CALL cp_cfm_release(cfm_mos_ikp(1))
1650 140 : CALL cp_cfm_release(cfm_mos_ikp(2))
1651 140 : CALL cp_cfm_release(cfm_work_ikp)
1652 140 : CALL cp_cfm_release(cfm_s_ikp_copy)
1653 :
1654 140 : CALL cp_cfm_release(cfm_s_ikp_spinor)
1655 140 : CALL cp_cfm_release(cfm_ks_ikp_spinor)
1656 140 : CALL cp_cfm_release(cfm_SOC_ikp_spinor)
1657 140 : CALL cp_cfm_release(cfm_mos_ikp_spinor)
1658 140 : CALL cp_cfm_release(cfm_work_ikp_spinor)
1659 140 : CALL cp_cfm_release(cfm_s_ikp_spinor_copy)
1660 140 : CALL cp_cfm_release(cfm_spinor_wf_ikp)
1661 :
1662 140 : CALL timestop(handle)
1663 :
1664 560 : 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 6 : 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 6 : 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 6 : CALL timeset(routineN, handle)
1691 :
1692 6 : n_E = SIZE(LDOS_2d, 3)
1693 :
1694 6 : energy_window = bs_env%energy_window_DOS
1695 6 : energy_step = bs_env%energy_step_DOS
1696 6 : E_min = band_edges%VBM - 0.5_dp*energy_window
1697 :
1698 18 : bin_mesh(1:2) = bs_env%bin_mesh(1:2)
1699 6 : do_xy_bins = (bin_mesh(1) > 0 .AND. bin_mesh(2) > 0)
1700 :
1701 6 : i_x_start = LBOUND(LDOS_2d, 1)
1702 6 : i_x_end = UBOUND(LDOS_2d, 1)
1703 6 : i_y_start = LBOUND(LDOS_2d, 2)
1704 6 : i_y_end = UBOUND(LDOS_2d, 2)
1705 :
1706 6 : IF (do_xy_bins) THEN
1707 6 : i_x_start_bin = 1
1708 6 : i_x_end_bin = bin_mesh(1)
1709 6 : i_y_start_bin = 1
1710 6 : 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 30 : ALLOCATE (LDOS_2d_bins(i_x_start_bin:i_x_end_bin, i_y_start_bin:i_y_end_bin, n_E))
1719 6 : LDOS_2d_bins(:, :, :) = 0.0_dp
1720 :
1721 6 : IF (do_xy_bins) THEN
1722 :
1723 6 : i_x_start_glob = i_x_start
1724 6 : i_x_end_glob = i_x_end
1725 6 : i_y_start_glob = i_y_start
1726 6 : i_y_end_glob = i_y_end
1727 :
1728 6 : CALL bs_env%para_env%min(i_x_start_glob)
1729 6 : CALL bs_env%para_env%max(i_x_end_glob)
1730 6 : CALL bs_env%para_env%min(i_y_start_glob)
1731 6 : CALL bs_env%para_env%max(i_y_end_glob)
1732 :
1733 24 : 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 390 : DO i_y = i_y_start, i_y_end
1737 4230 : DO i_x = i_x_start, i_x_end
1738 3840 : i_x_bin = bin_mesh(1)*(i_x - i_x_start_glob)/(i_x_end_glob - i_x_start_glob + 1) + 1
1739 3840 : 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 1073920 : LDOS_2d(i_x, i_y, :)
1742 4224 : 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 6 : CALL bs_env%para_env%sum(LDOS_2d_bins)
1747 6 : 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 30 : DO i_y_bin = 1, bin_mesh(2)
1751 126 : 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 26872 : REAL(n_sum_for_bins(i_x_bin, i_y_bin), KIND=dp)
1754 : END DO
1755 : END DO
1756 :
1757 : ELSE
1758 :
1759 0 : LDOS_2d_bins(:, :, :) = LDOS_2d(:, :, :)
1760 :
1761 : END IF
1762 :
1763 6 : IF (bin_mesh(1)*bin_mesh(2) < bs_env%n_bins_max_for_printing) THEN
1764 6 : CALL print_LDOS_2d_bins(LDOS_2d_bins, bs_env, E_min, scf_gw_soc)
1765 : ELSE
1766 0 : CPWARN("The number of bins for the LDOS is too large. Decrease BIN_MESH.")
1767 : END IF
1768 :
1769 6 : CALL timestop(handle)
1770 :
1771 12 : 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 6 : 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 6 : CALL timeset(routineN, handle)
1798 :
1799 6 : i_x_start = LBOUND(LDOS_2d_bins, 1)
1800 6 : i_x_end = UBOUND(LDOS_2d_bins, 1)
1801 6 : i_y_start = LBOUND(LDOS_2d_bins, 2)
1802 6 : i_y_end = UBOUND(LDOS_2d_bins, 2)
1803 6 : n_E = SIZE(LDOS_2d_bins, 3)
1804 :
1805 6 : n_x = i_x_end - i_x_start + 1
1806 6 : n_y = i_y_end - i_y_start + 1
1807 :
1808 6 : IF (bs_env%para_env%is_source()) THEN
1809 :
1810 15 : DO i_y = i_y_start, i_y_end
1811 63 : DO i_x = i_x_start, i_x_end
1812 :
1813 48 : idx(1) = (REAL(i_x, KIND=dp) - 0.5_dp)/REAL(n_x, KIND=dp)
1814 48 : idx(2) = (REAL(i_y, KIND=dp) - 0.5_dp)/REAL(n_y, KIND=dp)
1815 48 : idx(3) = 0.0_dp
1816 624 : coord(1:3) = MATMUL(bs_env%hmat, idx)
1817 :
1818 48 : CALL get_print_format(coord(1), print_format_1)
1819 48 : CALL get_print_format(coord(2), print_format_2)
1820 :
1821 48 : print_format = "(3A,"//print_format_1//",A,"//print_format_2//",A)"
1822 :
1823 48 : WRITE (fname, print_format) "LDOS_", scf_gw_soc, &
1824 96 : "_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 48 : file_action="WRITE")
1828 :
1829 48 : WRITE (iunit, "(2A)") " Energy E (eV) average LDOS(x,y,E) (1/(eV*Å^2), ", &
1830 96 : "integrated over z, averaged inside bin)"
1831 :
1832 13424 : DO i_E = 1, n_E
1833 13376 : energy = E_min + i_E*bs_env%energy_step_DOS
1834 13376 : WRITE (iunit, "(2F17.3)") energy*evolt, &
1835 : LDOS_2d_bins(i_x, i_y, i_E)* &
1836 26800 : bs_env%unit_ldos_int_z_inv_Ang2_eV
1837 : END DO
1838 :
1839 60 : CALL close_file(iunit)
1840 :
1841 : END DO
1842 : END DO
1843 :
1844 : END IF
1845 :
1846 6 : CALL timestop(handle)
1847 :
1848 6 : END SUBROUTINE print_LDOS_2d_bins
1849 :
1850 : ! **************************************************************************************************
1851 : !> \brief ...
1852 : !> \param coord ...
1853 : !> \param print_format ...
1854 : ! **************************************************************************************************
1855 96 : 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 96 : CALL timeset(routineN, handle)
1864 :
1865 96 : IF (coord < -10000/angstrom) THEN
1866 0 : print_format = "F9.2"
1867 96 : ELSE IF (coord < -1000/angstrom) THEN
1868 0 : print_format = "F8.2"
1869 96 : ELSE IF (coord < -100/angstrom) THEN
1870 0 : print_format = "F7.2"
1871 96 : ELSE IF (coord < -10/angstrom) THEN
1872 0 : print_format = "F6.2"
1873 96 : ELSE IF (coord < -1/angstrom) THEN
1874 0 : print_format = "F5.2"
1875 96 : ELSE IF (coord < 10/angstrom) THEN
1876 96 : print_format = "F4.2"
1877 0 : ELSE IF (coord < 100/angstrom) THEN
1878 0 : print_format = "F5.2"
1879 0 : ELSE IF (coord < 1000/angstrom) THEN
1880 0 : print_format = "F6.2"
1881 0 : ELSE IF (coord < 10000/angstrom) THEN
1882 0 : print_format = "F7.2"
1883 : ELSE
1884 0 : print_format = "F8.2"
1885 : END IF
1886 :
1887 96 : CALL timestop(handle)
1888 :
1889 96 : 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 656 : 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 656 : CALL timeset(routineN, handle)
1926 :
1927 : !for the window truncation part
1928 656 : 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 0 : 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 656 : DOS, PDOS, band_edges, eigenval_spinor, cfm_spinor_wf_ikp)
1936 : END IF
1937 :
1938 656 : CALL timestop(handle)
1939 :
1940 656 : 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 656 : 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 656 : CALL timeset(routineN, handle)
1982 :
1983 656 : n_ao = bs_env%n_ao
1984 656 : homo_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin)
1985 656 : CALL get_qs_env(qs_env, nkind=nkind)
1986 :
1987 656 : CALL cp_cfm_create(cfm_ks_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1988 656 : CALL cp_cfm_create(cfm_SOC_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1989 656 : CALL cp_cfm_create(cfm_mos_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1990 656 : CALL cp_cfm_create(cfm_work_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1991 656 : CALL cp_cfm_create(cfm_eigenvec_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1992 :
1993 1968 : ALLOCATE (eigenval_spinor_no_SOC(2*n_ao))
1994 2624 : ALLOCATE (proj_mo_on_kind_spinor(2*n_ao, nkind))
1995 : ! PDOS not yet implemented -> projection is just zero -> PDOS is zero
1996 656 : proj_mo_on_kind_spinor(:, :) = 0.0_dp
1997 :
1998 : ! 1. get V^SOC_µν,σσ'(k_i)
1999 676 : 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 20 : 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 656 : 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 656 : CALL cp_cfm_set_all(cfm_mos_ikp_spinor, z_zero)
2020 656 : CALL add_cfm_submat(cfm_mos_ikp_spinor, cfm_mos_ikp(1), 1, 1)
2021 656 : 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 656 : 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 656 : 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 8800 : eigenval_spinor_no_SOC(1:n_ao) = eigenval_no_SOC(1:n_ao, ikp, 1)
2037 8800 : eigenval_spinor_no_SOC(n_ao + 1:) = eigenval_no_SOC(1:n_ao, ikp, bs_env%n_spin)
2038 656 : 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 256 : 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 656 : 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 656 : CALL cp_cfm_heevd(cfm_ks_ikp_spinor, cfm_eigenvec_ikp_spinor, eigenval_spinor)
2052 :
2053 : ! 6. DOS from spinors, no PDOS
2054 656 : IF (bs_env%do_dos_pdos) THEN
2055 452 : n_E = SIZE(DOS)
2056 : CALL add_to_DOS_PDOS(DOS, PDOS, eigenval_spinor, &
2057 452 : 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 656 : band_edges%VBM = MAX(band_edges%VBM, eigenval_spinor(homo_spinor))
2062 656 : 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 656 : - eigenval_spinor(homo_spinor))
2065 :
2066 : ! 8. spinor wavefunctions:
2067 656 : 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 656 : z_zero, cfm_spinor_wf_ikp)
2071 :
2072 656 : CALL cp_cfm_release(cfm_ks_ikp_spinor)
2073 656 : CALL cp_cfm_release(cfm_SOC_ikp_spinor)
2074 656 : CALL cp_cfm_release(cfm_work_ikp_spinor)
2075 656 : CALL cp_cfm_release(cfm_eigenvec_ikp_spinor)
2076 656 : CALL cp_cfm_release(cfm_mos_ikp_spinor)
2077 :
2078 656 : CALL timestop(handle)
2079 :
2080 1968 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig_win_sorted, eigenval_spinor_no_SOC
2114 0 : 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 0 : CALL timeset(routineN, handle)
2123 :
2124 0 : n_ao = bs_env%n_ao
2125 0 : homo_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin)
2126 0 : CALL get_qs_env(qs_env, nkind=nkind)
2127 :
2128 0 : nocc = bs_env%soc_nbands_occ
2129 0 : nvirt = bs_env%soc_nbands_virt
2130 0 : CPASSERT(nocc <= bs_env%n_occ_nosmear(1))
2131 0 : CPASSERT(nocc <= bs_env%n_occ_nosmear(bs_env%n_spin))
2132 0 : CPASSERT(nvirt <= n_ao - bs_env%n_occ_nosmear(1))
2133 0 : CPASSERT(nvirt <= n_ao - bs_env%n_occ_nosmear(bs_env%n_spin))
2134 0 : CPASSERT(.NOT. bs_env%do_dos_pdos)
2135 0 : mo_index_start_1 = bs_env%n_occ_nosmear(1) - nocc + 1
2136 0 : mo_index_start_2 = bs_env%n_occ_nosmear(bs_env%n_spin) - nocc + 1
2137 0 : n_window = nocc + nvirt
2138 0 : n_window2tot = 2*n_window
2139 0 : homo_used = 2*nocc
2140 :
2141 0 : NULLIFY (fst_ww, fst_wa, fst_mos)
2142 : CALL cp_fm_struct_create(fst_ww, nrow_global=n_window2tot, ncol_global=n_window2tot, &
2143 0 : 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 0 : 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 0 : template_fmstruct=bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
2148 0 : CALL cp_cfm_create(cfm_ks_ikp_spinor, fst_ww)
2149 0 : CALL cp_cfm_create(cfm_eigenvec_ikp_spinor, fst_ww)
2150 0 : CALL cp_cfm_create(cfm_work_ikp_spinor, fst_wa)
2151 0 : CALL cp_cfm_create(cfm_SOC_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
2152 0 : CALL cp_cfm_create(cfm_mos_ikp_spinor, fst_mos)
2153 :
2154 0 : ALLOCATE (eigenval_spinor_no_SOC(2*n_ao))
2155 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : CALL cp_fm_create(fm_src_re, cfm_mos_ikp(1)%matrix_struct)
2185 0 : CALL cp_fm_create(fm_src_im, cfm_mos_ikp(1)%matrix_struct)
2186 0 : CALL cp_fm_create(fm_tgt_re, fst_mos)
2187 0 : CALL cp_fm_create(fm_tgt_im, fst_mos)
2188 0 : CALL cp_fm_set_all(fm_tgt_re, 0.0_dp)
2189 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : s_firstrow=1, s_firstcol=mo_index_start_2, t_firstrow=n_ao + 1, t_firstcol=n_window + 1)
2203 :
2204 0 : CALL cp_fm_to_cfm(fm_tgt_re, fm_tgt_im, cfm_mos_ikp_spinor)
2205 :
2206 0 : CALL cp_fm_release(fm_src_re)
2207 0 : CALL cp_fm_release(fm_src_im)
2208 0 : CALL cp_fm_release(fm_tgt_re)
2209 0 : 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 0 : 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 0 : 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 0 : eigenval_spinor_no_SOC(1:n_ao) = eigenval_no_SOC(1:n_ao, ikp, 1)
2223 0 : 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 0 : 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 0 : ALLOCATE (eig_win_sorted(n_window2tot))
2232 0 : 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 0 : band_edges%VBM = MAX(band_edges%VBM, eig_win_sorted(homo_used))
2238 0 : band_edges%CBM = MIN(band_edges%CBM, eig_win_sorted(homo_used + 1))
2239 0 : 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 0 : n_window2tot, eigenval_spinor)
2244 :
2245 : ! 8. spinor wavefunctions:
2246 0 : CALL cp_cfm_set_all(cfm_spinor_wf_ikp, z_zero)
2247 0 : 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 0 : z_zero, cfm_spinor_wf_ikp)
2251 :
2252 0 : CALL cp_cfm_release(cfm_ks_ikp_spinor)
2253 0 : CALL cp_cfm_release(cfm_SOC_ikp_spinor)
2254 0 : CALL cp_cfm_release(cfm_work_ikp_spinor)
2255 0 : CALL cp_cfm_release(cfm_eigenvec_ikp_spinor)
2256 0 : CALL cp_cfm_release(cfm_mos_ikp_spinor)
2257 :
2258 0 : CALL cp_fm_struct_release(fst_ww)
2259 0 : CALL cp_fm_struct_release(fst_wa)
2260 0 : CALL cp_fm_struct_release(fst_mos)
2261 :
2262 0 : CALL timestop(handle)
2263 :
2264 0 : 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 0 : 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 0 : 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 0 : 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 0 : IF (mo_index_start_1 > 1) THEN
2299 0 : IF (eig_win_sorted(1) < eigenval_no_SOC(mo_index_start_1 - 1, ikp, 1)) THEN
2300 0 : CPABORT("SOC_NBANDS_OCC too small.")
2301 : END IF
2302 : END IF
2303 0 : IF (mo_index_start_2 > 1) THEN
2304 0 : IF (eig_win_sorted(1) < eigenval_no_SOC(mo_index_start_2 - 1, ikp, n_spin)) THEN
2305 0 : CPABORT("SOC_NBANDS_OCC too small.")
2306 : END IF
2307 : END IF
2308 0 : IF (mo_index_start_1 + n_window <= n_ao) THEN
2309 0 : IF (eig_win_sorted(n_window2tot) > eigenval_no_SOC(mo_index_start_1 + n_window, ikp, 1)) THEN
2310 0 : CPABORT("SOC_NBANDS_VIRT too small.")
2311 : END IF
2312 : END IF
2313 0 : IF (mo_index_start_2 + n_window <= n_ao) THEN
2314 0 : IF (eig_win_sorted(n_window2tot) > eigenval_no_SOC(mo_index_start_2 + n_window, ikp, n_spin)) THEN
2315 0 : CPABORT("SOC_NBANDS_VIRT too small.")
2316 : END IF
2317 : END IF
2318 :
2319 0 : mo_index_start_1_m1 = mo_index_start_1 - 1
2320 0 : mo_index_start_2_m1 = mo_index_start_2 - 1
2321 :
2322 0 : ia = 1
2323 0 : ib = 1
2324 0 : DO i = 1, mo_index_start_1_m1 + mo_index_start_2_m1
2325 0 : IF (ia <= mo_index_start_1_m1 .AND. ib <= mo_index_start_2_m1) THEN
2326 0 : IF (eigenval_no_SOC(ia, ikp, 1) <= eigenval_no_SOC(ib, ikp, n_spin)) THEN
2327 0 : eigenval_spinor(i) = eigenval_no_SOC(ia, ikp, 1)
2328 0 : ia = ia + 1
2329 : ELSE
2330 0 : eigenval_spinor(i) = eigenval_no_SOC(ib, ikp, n_spin)
2331 0 : ib = ib + 1
2332 : END IF
2333 0 : ELSE IF (ia <= mo_index_start_1_m1) THEN
2334 0 : eigenval_spinor(i) = eigenval_no_SOC(ia, ikp, 1)
2335 0 : ia = ia + 1
2336 : ELSE
2337 0 : eigenval_spinor(i) = eigenval_no_SOC(ib, ikp, n_spin)
2338 0 : ib = ib + 1
2339 : END IF
2340 : END DO
2341 :
2342 0 : idx_start = mo_index_start_1_m1 + mo_index_start_2_m1 + 1
2343 0 : idx_end = mo_index_start_1_m1 + mo_index_start_2_m1 + n_window2tot
2344 :
2345 0 : eigenval_spinor(idx_start:idx_end) = eig_win_sorted(:)
2346 :
2347 0 : ia = mo_index_start_1 + n_window
2348 0 : ib = mo_index_start_2 + n_window
2349 0 : DO i = mo_index_start_1_m1 + mo_index_start_2_m1 + n_window2tot + 1, 2*n_ao
2350 0 : IF (ia <= n_ao .AND. ib <= n_ao) THEN
2351 0 : IF (eigenval_no_SOC(ia, ikp, 1) <= eigenval_no_SOC(ib, ikp, n_spin)) THEN
2352 0 : eigenval_spinor(i) = eigenval_no_SOC(ia, ikp, 1)
2353 0 : ia = ia + 1
2354 : ELSE
2355 0 : eigenval_spinor(i) = eigenval_no_SOC(ib, ikp, n_spin)
2356 0 : ib = ib + 1
2357 : END IF
2358 0 : ELSE IF (ia <= n_ao) THEN
2359 0 : eigenval_spinor(i) = eigenval_no_SOC(ia, ikp, 1)
2360 0 : ia = ia + 1
2361 : ELSE
2362 0 : eigenval_spinor(i) = eigenval_no_SOC(ib, ikp, n_spin)
2363 0 : ib = ib + 1
2364 : END IF
2365 : END DO
2366 :
2367 0 : DEALLOCATE (eig_win_sorted)
2368 :
2369 0 : CALL timestop(handle)
2370 :
2371 0 : 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 952 : 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 952 : CALL timeset(routineN, handle)
2400 :
2401 952 : energy_step_DOS = bs_env%energy_step_DOS
2402 952 : broadening = bs_env%broadening_DOS
2403 :
2404 952 : n_mo = SIZE(eigenval)
2405 952 : nkind = SIZE(proj_mo_on_kind, 2)
2406 :
2407 : ! normalize to closed-shell / open-shell
2408 952 : wkp = bs_env%kpoints_DOS%wkp(ikp)*bs_env%spin_degeneracy
2409 2147288 : DO i_E = 1, n_E
2410 2146336 : energy = E_min + i_E*energy_step_DOS
2411 47466732 : DO i_mo = 1, n_mo
2412 : ! DOS
2413 45319444 : DOS(i_E) = DOS(i_E) + wkp*Gaussian(energy - eigenval(i_mo), broadening)
2414 :
2415 : ! PDOS
2416 138104668 : DO i_kind = 1, nkind
2417 135958332 : 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 30345484 : Gaussian(energy - eigenval(i_mo), broadening)
2421 : END IF
2422 : END DO
2423 : END DO
2424 : END DO
2425 :
2426 952 : CALL timestop(handle)
2427 :
2428 952 : 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 6 : 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 6 : 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 6 : 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 6 : CALL timeset(routineN, handle)
2474 :
2475 6 : my_do_spinor = .FALSE.
2476 6 : IF (PRESENT(do_spinor)) my_do_spinor = do_spinor
2477 :
2478 6 : 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 6 : nimages = dft_control%nimages
2482 6 : dft_control%nimages = bs_env%nimages_scf
2483 :
2484 6 : energy_window = bs_env%energy_window_DOS
2485 6 : energy_step = bs_env%energy_step_DOS
2486 6 : broadening = bs_env%broadening_DOS
2487 :
2488 6 : E_min = band_edges%VBM - 0.5_dp*energy_window
2489 6 : E_max = band_edges%CBM + 0.5_dp*energy_window
2490 6 : E_total_window = E_max - E_min
2491 :
2492 6 : n_E = INT(E_total_window/energy_step)
2493 :
2494 6 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
2495 :
2496 6 : CALL auxbas_pw_pool%create_pw(LDOS_3d)
2497 6 : CALL auxbas_pw_pool%create_pw(rho_g)
2498 :
2499 6 : i_x_start = LBOUND(LDOS_3d%array, 1)
2500 6 : i_x_end = UBOUND(LDOS_3d%array, 1)
2501 6 : i_y_start = LBOUND(LDOS_3d%array, 2)
2502 6 : i_y_end = UBOUND(LDOS_3d%array, 2)
2503 6 : i_z_start = LBOUND(LDOS_3d%array, 3)
2504 6 : i_z_end = UBOUND(LDOS_3d%array, 3)
2505 :
2506 6 : z_start_global = i_z_start
2507 6 : z_end_global = i_z_end
2508 :
2509 6 : CALL bs_env%para_env%min(z_start_global)
2510 6 : CALL bs_env%para_env%max(z_end_global)
2511 6 : n_z = z_end_global - z_start_global + 1
2512 :
2513 36 : 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 0 : 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 6 : bs_env%unit_ldos_int_z_inv_Ang2_eV = bs_env%hmat(3, 3)/REAL(n_z, KIND=dp)/evolt/angstrom**2
2518 :
2519 6 : IF (ikp == 1) THEN
2520 30 : ALLOCATE (LDOS_2d(i_x_start:i_x_end, i_y_start:i_y_end, n_E))
2521 6 : LDOS_2d(:, :, :) = 0.0_dp
2522 : END IF
2523 :
2524 6 : CALL cp_cfm_create(cfm_work, cfm_mos_ikp%matrix_struct)
2525 6 : CALL cp_cfm_create(cfm_weighted_dm_ikp, cfm_mos_ikp%matrix_struct)
2526 6 : CALL cp_fm_create(fm_weighted_dm_MIC, cfm_mos_ikp%matrix_struct)
2527 6 : IF (my_do_spinor) THEN
2528 2 : 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 6 : col_indices=col_indices)
2535 :
2536 6 : NULLIFY (weighted_dm_MIC)
2537 6 : CALL dbcsr_allocate_matrix_set(weighted_dm_MIC, 1)
2538 6 : ALLOCATE (weighted_dm_MIC(1)%matrix)
2539 : CALL dbcsr_create(weighted_dm_MIC(1)%matrix, template=bs_env%mat_ao_ao%matrix, &
2540 6 : matrix_type=dbcsr_type_symmetric)
2541 :
2542 1678 : DO i_E = 1, n_E
2543 :
2544 1672 : energy = E_min + i_E*energy_step
2545 :
2546 1672 : is_any_weight_non_zero = .FALSE.
2547 :
2548 20950 : DO j_col = 1, ncol_local
2549 :
2550 19278 : j_mo = col_indices(j_col)
2551 :
2552 19278 : IF (my_do_spinor) THEN
2553 : spin_degeneracy = 1.0_dp
2554 : ELSE
2555 10818 : spin_degeneracy = bs_env%spin_degeneracy
2556 : END IF
2557 :
2558 19278 : weight = Gaussian(energy - eigenval(j_mo), broadening)*spin_degeneracy
2559 :
2560 144099 : cfm_work%local_data(:, j_col) = cfm_mos_ikp%local_data(:, j_col)*weight
2561 :
2562 20950 : IF (weight > 1.0E-5_dp) is_any_weight_non_zero = .TRUE.
2563 :
2564 : END DO
2565 :
2566 1672 : CALL bs_env%para_env%sync()
2567 1672 : CALL bs_env%para_env%sum(is_any_weight_non_zero)
2568 1672 : CALL bs_env%para_env%sync()
2569 :
2570 : ! cycle if there are no states at the energy i_E
2571 1678 : IF (is_any_weight_non_zero) THEN
2572 :
2573 : CALL parallel_gemm('N', 'C', n_mo, n_mo, n_mo, z_one, &
2574 24 : cfm_mos_ikp, cfm_work, z_zero, cfm_weighted_dm_ikp)
2575 :
2576 24 : IF (my_do_spinor) THEN
2577 :
2578 : ! contribution from up,up to fm_non_spinor
2579 8 : CALL get_cfm_submat(cfm_non_spinor, cfm_weighted_dm_ikp, 1, 1)
2580 8 : 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 8 : "ORB", bs_env%kpoints_DOS%wkp(ikp))
2584 :
2585 : ! add contribution from down,down to fm_non_spinor
2586 8 : 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 8 : "ORB", bs_env%kpoints_DOS%wkp(ikp))
2590 : CALL copy_fm_to_dbcsr(fm_non_spinor, weighted_dm_MIC(1)%matrix, &
2591 8 : keep_sparsity=.FALSE.)
2592 : ELSE
2593 16 : 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 16 : "ORB", bs_env%kpoints_DOS%wkp(ikp))
2597 : CALL copy_fm_to_dbcsr(fm_weighted_dm_MIC, weighted_dm_MIC(1)%matrix, &
2598 16 : keep_sparsity=.FALSE.)
2599 : END IF
2600 :
2601 338424 : 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 24 : ks_env=ks_env)
2607 :
2608 504 : DO i_z = i_z_start, i_z_end
2609 338424 : 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 6 : dft_control%nimages = nimages
2618 :
2619 6 : CALL auxbas_pw_pool%give_back_pw(LDOS_3d)
2620 6 : CALL auxbas_pw_pool%give_back_pw(rho_g)
2621 :
2622 6 : CALL cp_cfm_release(cfm_work)
2623 6 : CALL cp_cfm_release(cfm_weighted_dm_ikp)
2624 :
2625 6 : CALL cp_fm_release(fm_weighted_dm_MIC)
2626 :
2627 6 : CALL dbcsr_deallocate_matrix_set(weighted_dm_MIC)
2628 :
2629 6 : IF (my_do_spinor) THEN
2630 2 : CALL cp_fm_release(fm_non_spinor)
2631 : END IF
2632 :
2633 6 : CALL timestop(handle)
2634 :
2635 6 : 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 296 : 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 296 : CALL timeset(routineN, handle)
2659 :
2660 296 : fname = "bandstructure_SCF_and_G0W0_plus_SOC"
2661 :
2662 296 : IF (bs_env%para_env%is_source()) THEN
2663 :
2664 148 : IF (ikp_for_file == 1) THEN
2665 : CALL open_file(TRIM(fname), unit_number=iunit, file_status="REPLACE", &
2666 11 : file_action="WRITE")
2667 : ELSE
2668 : CALL open_file(TRIM(fname), unit_number=iunit, file_status="OLD", &
2669 137 : file_action="WRITE", file_position="APPEND")
2670 : END IF
2671 :
2672 148 : WRITE (iunit, "(A)") " "
2673 148 : WRITE (iunit, "(A10,I7,A25,3F10.4)") "kpoint: ", ikp_for_file, "coordinate: ", &
2674 740 : bs_env%kpoints_DOS%xkp(:, ikp)
2675 148 : WRITE (iunit, "(A)") " "
2676 :
2677 148 : IF (PRESENT(eigenval_spinor_G0W0)) THEN
2678 : ! SCF+SOC and G0W0+SOC eigenvalues
2679 148 : WRITE (iunit, "(A5,A12,2A22)") "n", "k", "ϵ_nk^DFT+SOC (eV)", "ϵ_nk^G0W0+SOC (eV)"
2680 : ELSE
2681 : ! SCF+SOC eigenvalues only
2682 0 : WRITE (iunit, "(A5,A12,A22)") "n", "k", "ϵ_nk^DFT+SOC (eV)"
2683 : END IF
2684 :
2685 148 : n_occ_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin)
2686 :
2687 3932 : DO i_mo = 1, SIZE(eigenval_spinor)
2688 3784 : IF (i_mo <= n_occ_spinor) occ_vir = 'occ'
2689 3784 : IF (i_mo > n_occ_spinor) occ_vir = 'vir'
2690 3932 : IF (PRESENT(eigenval_spinor_G0W0)) THEN
2691 : ! SCF+SOC and G0W0+SOC eigenvalues
2692 3784 : WRITE (iunit, "(I5,3A,I5,4F16.3,2F17.3)") i_mo, ' (', occ_vir, ') ', &
2693 7568 : ikp_for_file, eigenval_spinor(i_mo)*evolt, eigenval_spinor_G0W0(i_mo)*evolt
2694 : ELSE
2695 : ! SCF+SOC eigenvalues only
2696 0 : WRITE (iunit, "(I5,3A,I5,4F16.3,F17.3)") i_mo, ' (', occ_vir, ') ', &
2697 0 : ikp_for_file, eigenval_spinor(i_mo)*evolt
2698 : END IF
2699 : END DO
2700 :
2701 148 : CALL close_file(iunit)
2702 :
2703 : END IF
2704 :
2705 296 : CALL timestop(handle)
2706 :
2707 296 : END SUBROUTINE write_SOC_eigenvalues
2708 :
2709 : ! **************************************************************************************************
2710 : !> \brief ...
2711 : !> \param int_number ...
2712 : !> \return ...
2713 : ! **************************************************************************************************
2714 0 : PURE FUNCTION count_digits(int_number)
2715 :
2716 : INTEGER, INTENT(IN) :: int_number
2717 : INTEGER :: count_digits
2718 :
2719 : INTEGER :: digitCount, tempInt
2720 :
2721 0 : digitCount = 0
2722 :
2723 0 : tempInt = int_number
2724 :
2725 0 : DO WHILE (tempInt /= 0)
2726 0 : tempInt = tempInt/10
2727 0 : digitCount = digitCount + 1
2728 : END DO
2729 :
2730 0 : count_digits = digitCount
2731 :
2732 0 : END FUNCTION count_digits
2733 :
2734 : ! **************************************************************************************************
2735 : !> \brief ...
2736 : !> \param band_edges ...
2737 : !> \param scf_gw_soc ...
2738 : !> \param bs_env ...
2739 : ! **************************************************************************************************
2740 460 : 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 460 : CALL timeset(routineN, handle)
2752 :
2753 : ! print format
2754 460 : print_format = "(T2,2A,T61,F20.3)"
2755 :
2756 460 : u = bs_env%unit_nr
2757 460 : IF (u > 0) THEN
2758 230 : WRITE (u, '(T2,A)') ''
2759 230 : WRITE (u, print_format) scf_gw_soc, ' valence band maximum (eV):', band_edges%VBM*evolt
2760 230 : WRITE (u, print_format) scf_gw_soc, ' conduction band minimum (eV):', band_edges%CBM*evolt
2761 230 : WRITE (u, print_format) scf_gw_soc, ' indirect band gap (eV):', band_edges%IDBG*evolt
2762 230 : WRITE (u, print_format) scf_gw_soc, ' direct band gap (eV):', band_edges%DBG*evolt
2763 : END IF
2764 :
2765 460 : CALL timestop(handle)
2766 :
2767 460 : 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 84 : 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 84 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2795 :
2796 84 : CALL timeset(routineN, handle)
2797 :
2798 84 : WRITE (fname, "(3A)") "DOS_PDOS_", scf_gw_soc, ".out"
2799 :
2800 84 : n_E = SIZE(PDOS, 1)
2801 84 : nkind = SIZE(PDOS, 2)
2802 84 : CALL get_qs_env(qs_env, particle_set=particle_set)
2803 :
2804 84 : IF (bs_env%para_env%is_source()) THEN
2805 :
2806 42 : CALL open_file(TRIM(fname), unit_number=iunit, file_status="REPLACE", file_action="WRITE")
2807 :
2808 42 : n_A = 2 + nkind
2809 :
2810 162 : DO iatom = 1, bs_env%n_atom
2811 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
2812 120 : kind_number=i_kind, name=atom_name)
2813 162 : elements(i_kind) = atom_name(1:3)
2814 : END DO
2815 :
2816 42 : WRITE (output_string, "(A,I1,A)") "(", n_A, "A)"
2817 :
2818 42 : WRITE (iunit, TRIM(output_string)) "Energy-E_F (eV) DOS (1/eV) PDOS (1/eV) ", &
2819 84 : " of atom type ", elements(1:nkind)
2820 :
2821 42 : WRITE (output_string, "(A,I1,A)") "(", n_A, "F13.5)"
2822 :
2823 87510 : DO i_E = 1, n_E
2824 : ! energy is relative to valence band maximum => - E_VBM
2825 87468 : energy = E_min + i_E*bs_env%energy_step_DOS - E_VBM
2826 262446 : WRITE (iunit, TRIM(output_string)) energy*evolt, DOS(i_E)/evolt, PDOS(i_E, :)/evolt
2827 : END DO
2828 :
2829 42 : CALL close_file(iunit)
2830 :
2831 : END IF
2832 :
2833 84 : CALL timestop(handle)
2834 :
2835 84 : END SUBROUTINE write_dos_pdos
2836 :
2837 : ! **************************************************************************************************
2838 : !> \brief ...
2839 : !> \param energy ...
2840 : !> \param broadening ...
2841 : !> \return ...
2842 : ! **************************************************************************************************
2843 75684206 : PURE FUNCTION Gaussian(energy, broadening)
2844 :
2845 : REAL(KIND=dp), INTENT(IN) :: energy, broadening
2846 : REAL(KIND=dp) :: Gaussian
2847 :
2848 75684206 : IF (ABS(energy) < 5*broadening) THEN
2849 116164 : 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 75684206 : 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 514 : 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 514 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_bf, kind_of
2874 514 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2875 514 : 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 514 : CALL timeset(routineN, handle)
2880 :
2881 514 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, nkind=nkind)
2882 514 : 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 514 : col_indices=col_indices)
2890 :
2891 514 : n_ao = qs_env%bs_env%n_ao
2892 :
2893 1542 : ALLOCATE (atom_from_bf(n_ao))
2894 514 : CALL get_atom_index_from_basis_function_index(qs_env, atom_from_bf, n_ao, "ORB")
2895 :
2896 514 : proj_mo_on_kind(:, :) = 0.0_dp
2897 :
2898 514 : CALL cp_cfm_create(cfm_s_i_kind, cfm_s%matrix_struct)
2899 514 : CALL cp_cfm_create(cfm_work, cfm_s%matrix_struct)
2900 514 : CALL cp_cfm_create(cfm_proj, cfm_s%matrix_struct)
2901 514 : CALL cp_fm_create(fm_proj_re, cfm_s%matrix_struct)
2902 514 : CALL cp_fm_create(fm_proj_im, cfm_s%matrix_struct)
2903 :
2904 1420 : DO i_kind = 1, nkind
2905 :
2906 906 : 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 12198 : DO j_col = 1, ncol_local
2910 90134 : DO i_row = 1, nrow_local
2911 :
2912 77936 : i_global = row_indices(i_row)
2913 :
2914 77936 : IF (i_global <= n_ao) THEN
2915 77936 : i_atom = atom_from_bf(i_global)
2916 0 : ELSE IF (i_global <= 2*n_ao) THEN
2917 0 : i_atom = atom_from_bf(i_global - n_ao)
2918 : ELSE
2919 0 : CPABORT("Wrong indices.")
2920 : END IF
2921 :
2922 89228 : IF (i_kind /= kind_of(i_atom)) THEN
2923 35834 : 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 906 : 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 906 : cfm_mos, cfm_work, z_zero, cfm_proj)
2933 :
2934 906 : CALL cp_cfm_to_fm(cfm_proj, fm_proj_re, fm_proj_im)
2935 :
2936 906 : CALL cp_fm_get_diag(fm_proj_im, proj_mo_on_kind(:, i_kind))
2937 1420 : CALL cp_fm_get_diag(fm_proj_re, proj_mo_on_kind(:, i_kind))
2938 :
2939 : END DO ! i_kind
2940 :
2941 514 : CALL cp_cfm_release(cfm_s_i_kind)
2942 514 : CALL cp_cfm_release(cfm_work)
2943 514 : CALL cp_cfm_release(cfm_proj)
2944 514 : CALL cp_fm_release(fm_proj_re)
2945 514 : CALL cp_fm_release(fm_proj_im)
2946 :
2947 514 : CALL timestop(handle)
2948 :
2949 2056 : 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 120 : 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 20 : CALL timeset(routineN, handle)
2979 :
2980 20 : CALL cp_cfm_create(cfm_non_spinor_Gamma, fm_struct_non_spinor)
2981 20 : CALL cp_cfm_create(cfm_non_spinor_ikp, fm_struct_non_spinor)
2982 20 : CALL cp_fm_create(fm_non_spinor_Gamma_re, fm_struct_non_spinor)
2983 20 : CALL cp_fm_create(fm_non_spinor_Gamma_im, fm_struct_non_spinor)
2984 :
2985 20 : CALL cp_cfm_get_info(cfm_non_spinor_Gamma, nrow_global=n_ao)
2986 :
2987 20 : CALL cp_cfm_set_all(cfm_spinor_ikp, z_zero)
2988 :
2989 60 : DO i_block = 0, 1
2990 140 : DO j_block = 0, 1
2991 80 : i_offset = i_block*n_ao + 1
2992 80 : j_offset = j_block*n_ao + 1
2993 80 : CALL get_cfm_submat(cfm_non_spinor_Gamma, cfm_spinor_Gamma, i_offset, j_offset)
2994 80 : 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 80 : ikp, qs_env, kpoints, basis_type)
2999 80 : 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 80 : ikp, qs_env, kpoints, basis_type)
3004 120 : 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 20 : CALL cp_cfm_release(cfm_non_spinor_Gamma)
3010 20 : CALL cp_cfm_release(cfm_non_spinor_ikp)
3011 20 : CALL cp_fm_release(fm_non_spinor_Gamma_re)
3012 20 : CALL cp_fm_release(fm_non_spinor_Gamma_im)
3013 :
3014 20 : CALL timestop(handle)
3015 :
3016 20 : 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 7788 : 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 7788 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_bf
3039 : TYPE(cell_type), POINTER :: cell
3040 7788 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3041 :
3042 7788 : CALL timeset(routineN, handle)
3043 :
3044 7788 : CPASSERT(basis_type == "ORB" .OR. basis_type == "RI_AUX")
3045 23364 : ALLOCATE (atom_from_bf(fm_Gamma%matrix_struct%nrow_global))
3046 7788 : CALL get_atom_index_from_basis_function_index(qs_env, atom_from_bf, SIZE(atom_from_bf), basis_type)
3047 7788 : NULLIFY (cell, particle_set)
3048 7788 : 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 7788 : kpoints%index_to_cell, kpoints%xkp(:, ikp))
3051 :
3052 7788 : CALL timestop(handle)
3053 :
3054 15576 : 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 3410 : 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 3410 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_bf_index
3084 3410 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
3085 3410 : 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 3410 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
3090 3410 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
3091 : TYPE(cell_type), POINTER :: cell
3092 3410 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3093 :
3094 3410 : CALL timeset(routineN, handle)
3095 :
3096 : ! get number of basis functions (bf) for different basis sets
3097 3410 : IF (basis_type == "ORB") THEN
3098 32 : n_bf = qs_env%bs_env%n_ao
3099 3378 : ELSE IF (basis_type == "RI_AUX") THEN
3100 3378 : n_bf = qs_env%bs_env%n_RI
3101 : ELSE
3102 0 : CPABORT("Only ORB and RI_AUX basis implemented.")
3103 : END IF
3104 :
3105 10230 : ALLOCATE (atom_from_bf_index(n_bf))
3106 3410 : CALL get_atom_index_from_basis_function_index(qs_env, atom_from_bf_index, n_bf, basis_type)
3107 :
3108 3410 : NULLIFY (cell, particle_set)
3109 3410 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
3110 3410 : 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 3410 : col_indices=col_indices)
3117 :
3118 3410 : CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp)
3119 3410 : index_to_cell => kpoints%index_to_cell
3120 3410 : num_cells = SIZE(index_to_cell, 2)
3121 :
3122 3410 : iatom_old = 0
3123 3410 : jatom_old = 0
3124 :
3125 73304 : DO jcol = 1, ncol_local
3126 1529527 : DO irow = 1, nrow_local
3127 :
3128 1456223 : i_bf = row_indices(irow)
3129 1456223 : j_bf = col_indices(jcol)
3130 :
3131 1456223 : iatom = atom_from_bf_index(i_bf)
3132 1456223 : jatom = atom_from_bf_index(j_bf)
3133 :
3134 1456223 : IF (PRESENT(wkp_ext)) THEN
3135 3496 : wkp_of_ikp = wkp_ext
3136 : ELSE
3137 1497526 : 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 44799 : 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 193398 : 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 1452727 : wkp_of_ikp = bs_env%wkp_no_extra(ikp)
3147 : END SELECT
3148 : END IF
3149 :
3150 1456223 : 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 138140 : cell, index_to_cell, hmat, particle_set)
3155 :
3156 138140 : iatom_old = iatom
3157 138140 : jatom_old = jatom
3158 :
3159 : END IF
3160 :
3161 : contribution = weight_re*REAL(cfm_W_ikp_freq_j%local_data(irow, jcol)) + &
3162 1456223 : 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 1526117 : + contribution
3166 :
3167 : END DO
3168 : END DO
3169 :
3170 3410 : CALL timestop(handle)
3171 :
3172 10230 : 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 140 : 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 140 : CALL timeset(routineN, handle)
3192 :
3193 140 : i = ikp_start
3194 300 : DO ix = 1, grid(1)
3195 572 : DO iy = 1, grid(2)
3196 894 : DO iz = 1, grid(3)
3197 :
3198 462 : IF (i > ikp_end) CYCLE
3199 :
3200 428 : xkp(1, i) = REAL(2*ix - grid(1) - 1, KIND=dp)/(2._dp*REAL(grid(1), KIND=dp))
3201 428 : xkp(2, i) = REAL(2*iy - grid(2) - 1, KIND=dp)/(2._dp*REAL(grid(2), KIND=dp))
3202 428 : xkp(3, i) = REAL(2*iz - grid(3) - 1, KIND=dp)/(2._dp*REAL(grid(3), KIND=dp))
3203 734 : i = i + 1
3204 :
3205 : END DO
3206 : END DO
3207 : END DO
3208 :
3209 140 : CALL timestop(handle)
3210 :
3211 140 : END SUBROUTINE compute_xkp
3212 :
3213 : ! **************************************************************************************************
3214 : !> \brief ...
3215 : !> \param kpoints ...
3216 : !> \param qs_env ...
3217 : ! **************************************************************************************************
3218 244 : 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 122 : POINTER :: sab_orb
3229 :
3230 122 : CALL timeset(routineN, handle)
3231 :
3232 122 : NULLIFY (para_env, sab_orb)
3233 122 : CALL get_qs_env(qs_env=qs_env, para_env=para_env, sab_orb=sab_orb)
3234 122 : CALL kpoint_init_cell_index(kpoints, sab_orb, para_env, nimages)
3235 :
3236 122 : CALL timestop(handle)
3237 :
3238 122 : END SUBROUTINE kpoint_init_cell_index_simple
3239 :
3240 : ! **************************************************************************************************
3241 : !> \brief ...
3242 : !> \param qs_env ...
3243 : !> \param bs_env ...
3244 : ! **************************************************************************************************
3245 22 : 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 22 : 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 22 : 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 28 : 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 6 : 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 22 : CALL H_KS_spinor_kp(qs_env, bs_env)
3271 :
3272 : END SELECT
3273 :
3274 22 : CALL timestop(handle)
3275 :
3276 22 : END SUBROUTINE soc
3277 :
3278 : ! **************************************************************************************************
3279 : !> \brief ...
3280 : !> \param bs_env ...
3281 : ! **************************************************************************************************
3282 6 : 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 6 : CALL timeset(routineN, handle)
3292 :
3293 6 : CALL cp_fm_get_info(bs_env%fm_ks_Gamma(1), nrow_global=nao)
3294 :
3295 12 : ALLOCATE (bs_env%cfm_SOC_spinor_ao(1))
3296 6 : CALL create_cfm_double(bs_env%cfm_SOC_spinor_ao(1), fm_orig=bs_env%fm_ks_Gamma(1))
3297 6 : CALL cp_cfm_set_all(bs_env%cfm_SOC_spinor_ao(1), z_zero)
3298 :
3299 6 : str => bs_env%fm_ks_Gamma(1)%matrix_struct
3300 :
3301 6 : 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 6 : 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 6 : 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 6 : 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 6 : 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 6 : 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 6 : str, s, s, -z_one, .FALSE.)
3328 :
3329 6 : CALL timestop(handle)
3330 :
3331 6 : END SUBROUTINE H_KS_spinor_Gamma
3332 :
3333 : ! **************************************************************************************************
3334 : !> \brief ...
3335 : !> \param qs_env ...
3336 : !> \param bs_env ...
3337 : ! **************************************************************************************************
3338 32 : 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 16 : 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 16 : POINTER :: sab_nl
3353 :
3354 16 : CALL timeset(routineN, handle)
3355 :
3356 16 : nkp_bs_and_DOS = bs_env%nkp_bs_and_DOS
3357 16 : n_spin = bs_env%n_spin
3358 16 : s = bs_env%n_ao + 1
3359 16 : str => bs_env%cfm_ks_kp(1, 1)%matrix_struct
3360 :
3361 16 : CALL cp_cfm_create(cfm_V_SOC_xyz_ikp, bs_env%cfm_work_mo%matrix_struct)
3362 :
3363 16 : 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 16 : CALL get_qs_env(qs_env, kpoints=kpoints_scf)
3366 :
3367 16 : NULLIFY (sab_nl)
3368 16 : CALL get_kpoint_info(kpoints_scf, sab_nl=sab_nl, cell_to_index=cell_to_index_scf)
3369 :
3370 64 : DO i_dim = 1, 3
3371 :
3372 1018 : DO ikp = 1, nkp_bs_and_DOS
3373 :
3374 3816 : xkp(1:3) = bs_env%kpoints_DOS%xkp(1:3, ikp)
3375 :
3376 954 : 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 954 : 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 954 : CALL cp_cfm_scale(gaussi, cfm_V_SOC_xyz_ikp)
3383 :
3384 48 : SELECT CASE (i_dim)
3385 : CASE (1)
3386 : ! add V^SOC_x * σ_x for σ_x = ( (0,1) (1,0) )
3387 318 : CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_V_SOC_xyz_ikp, 1, s)
3388 318 : 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 318 : CALL cp_cfm_scale(gaussi, cfm_V_SOC_xyz_ikp)
3392 318 : CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_V_SOC_xyz_ikp, 1, s)
3393 318 : CALL cp_cfm_scale(-z_one, cfm_V_SOC_xyz_ikp)
3394 318 : 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 318 : CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_V_SOC_xyz_ikp, 1, 1)
3398 318 : CALL cp_cfm_scale(-z_one, cfm_V_SOC_xyz_ikp)
3399 1272 : 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 16 : CALL cp_cfm_release(cfm_V_SOC_xyz_ikp)
3407 :
3408 16 : CALL timestop(handle)
3409 :
3410 16 : END SUBROUTINE H_KS_spinor_kp
3411 :
3412 : ! **************************************************************************************************
3413 : !> \brief ...
3414 : !> \param cfm_array ...
3415 : !> \param cfm_template ...
3416 : !> \param n ...
3417 : ! **************************************************************************************************
3418 16 : 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 16 : CALL timeset(routineN, handle)
3428 :
3429 366 : ALLOCATE (cfm_array(n))
3430 334 : DO i = 1, n
3431 318 : CALL create_cfm_double(cfm_array(i), cfm_orig=cfm_template)
3432 334 : CALL cp_cfm_set_all(cfm_array(i), z_zero)
3433 : END DO
3434 :
3435 16 : CALL timestop(handle)
3436 :
3437 16 : END SUBROUTINE alloc_cfm_double_array_1d
3438 :
3439 : ! **************************************************************************************************
3440 : !> \brief ...
3441 : !> \param bs_env ...
3442 : ! **************************************************************************************************
3443 152 : 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 152 : CALL timeset(routineN, handle)
3452 :
3453 152 : CALL get_VBM_CBM_bandgaps(bs_env%band_edges_scf, bs_env%eigenval_scf, bs_env)
3454 152 : CALL get_VBM_CBM_bandgaps(bs_env%band_edges_GW, bs_env%eigenval_GW, bs_env)
3455 152 : CALL get_VBM_CBM_bandgaps(bs_env%band_edges_HF, bs_env%eigenval_HF, bs_env)
3456 :
3457 152 : CALL check_qp_gap_sanity(bs_env)
3458 :
3459 152 : CALL check_scf_gw_level_ordering(bs_env)
3460 :
3461 152 : CALL timestop(handle)
3462 :
3463 152 : 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 152 : 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 152 : gap = bs_env%band_edges_GW%IDBG
3478 152 : gap_scf = bs_env%band_edges_scf%IDBG
3479 152 : 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 152 : 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 0 : "spectrum is inverted. Check numerical parameters.")
3487 152 : 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 0 : "solve has likely diverged. Check numerical parameters.")
3492 : END IF
3493 :
3494 152 : 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 152 : 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 152 : CALL timeset(routineN, handle)
3519 :
3520 152 : n_mo = bs_env%n_ao
3521 152 : gw_label = gw_flavour_label(bs_env)
3522 :
3523 328 : DO ispin = 1, bs_env%n_spin
3524 :
3525 176 : homo = bs_env%n_occ(ispin)
3526 :
3527 176 : IF (bs_env%n_spin == 2) THEN
3528 48 : WRITE (spin_string, '(A,I0,A)') " (spin ", ispin, ")"
3529 : ELSE
3530 128 : spin_string = ""
3531 : END IF
3532 :
3533 : ! occupied levels: level carrying the SCF VBM versus level carrying the GW VBM
3534 2876 : loc_scf(:) = MAXLOC(bs_env%eigenval_scf(1:homo, :, ispin))
3535 2876 : 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 176 : 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 2 : CALL level_string(level_scf, loc_scf, bs_env)
3542 2 : 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 2 : "orbitals; see the eigenvalues in the file bandstructure_SCF_and_G0W0."
3550 2 : CPWARN(TRIM(msg))
3551 :
3552 : END IF
3553 :
3554 : ! empty levels: level carrying the SCF CBM versus level carrying the GW CBM
3555 5316 : loc_scf(:) = MINLOC(bs_env%eigenval_scf(homo + 1:n_mo, :, ispin))
3556 5316 : loc_gw(:) = MINLOC(bs_env%eigenval_GW(homo + 1:n_mo, :, ispin))
3557 176 : loc_scf(1) = loc_scf(1) + homo
3558 176 : loc_gw(1) = loc_gw(1) + homo
3559 :
3560 176 : IF (bs_env%eigenval_scf(loc_gw(1), loc_gw(2), ispin) - &
3561 152 : bs_env%eigenval_scf(loc_scf(1), loc_scf(2), ispin) > eps_degeneracy) THEN
3562 :
3563 0 : CALL level_string(level_scf, loc_scf, bs_env)
3564 0 : 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 0 : "orbitals; see the eigenvalues in the file bandstructure_SCF_and_G0W0."
3572 0 : CPWARN(TRIM(msg))
3573 :
3574 : END IF
3575 :
3576 : END DO
3577 :
3578 152 : CALL timestop(handle)
3579 :
3580 152 : 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 4 : 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 4 : IF (bs_env%nkp_bs_and_DOS > 1) THEN
3595 4 : WRITE (string, '(A,I0,A,I0)') "level ", loc(1), " at k-point ", loc(2)
3596 : ELSE
3597 0 : WRITE (string, '(A,I0)') "level ", loc(1)
3598 : END IF
3599 :
3600 4 : END SUBROUTINE level_string
3601 :
3602 : ! **************************************************************************************************
3603 : !> \brief ...
3604 : !> \param band_edges ...
3605 : !> \param ev ...
3606 : !> \param bs_env ...
3607 : ! **************************************************************************************************
3608 474 : 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 474 : CALL timeset(routineN, handle)
3620 :
3621 474 : n_mo = bs_env%n_ao
3622 :
3623 474 : band_edges%DBG = 1000.0_dp
3624 :
3625 876 : SELECT CASE (bs_env%n_spin)
3626 : CASE (1)
3627 402 : homo = bs_env%n_occ(1)
3628 402 : lumo = homo + 1
3629 8294 : band_edges%VBM = MAXVAL(ev(1:homo, :, 1))
3630 16656 : band_edges%CBM = MINVAL(ev(homo + 1:n_mo, :, 1))
3631 : CASE (2)
3632 72 : homo_1 = bs_env%n_occ(1)
3633 72 : lumo_1 = homo_1 + 1
3634 72 : homo_2 = bs_env%n_occ(2)
3635 72 : lumo_2 = homo_2 + 1
3636 984 : band_edges%VBM = MAX(MAXVAL(ev(1:homo_1, :, 1)), MAXVAL(ev(1:homo_2, :, 2)))
3637 1380 : band_edges%CBM = MIN(MINVAL(ev(homo_1 + 1:n_mo, :, 1)), MINVAL(ev(homo_2 + 1:n_mo, :, 2)))
3638 : CASE DEFAULT
3639 474 : CPABORT("Error with number of spins.")
3640 : END SELECT
3641 :
3642 474 : band_edges%IDBG = band_edges%CBM - band_edges%VBM
3643 :
3644 1020 : DO ispin = 1, bs_env%n_spin
3645 :
3646 546 : homo = bs_env%n_occ(ispin)
3647 :
3648 2872 : DO ikp = 1, bs_env%nkp_bs_and_DOS
3649 :
3650 24514 : E_DBG_at_ikp = -MAXVAL(ev(1:homo, ikp, ispin)) + MINVAL(ev(homo + 1:n_mo, ikp, ispin))
3651 :
3652 2398 : IF (E_DBG_at_ikp < band_edges%DBG) band_edges%DBG = E_DBG_at_ikp
3653 :
3654 : END DO
3655 :
3656 : END DO
3657 402 : 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 402 : 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 474 : band_edges%VBM, band_edges%CBM, band_edges%E_fermi)
3666 : END SELECT
3667 :
3668 474 : CALL timestop(handle)
3669 :
3670 474 : 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 568 : 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 568 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: f_occ
3703 :
3704 568 : IF (bs_env%do_smear .AND. .NOT. (CBM > VBM .AND. &
3705 : CBM - VBM > gap_smear_ratio*bs_env%smear_width)) THEN
3706 16 : 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 4 : bs_env%smear_width, maxocc, bs_env%smear_method)
3710 4 : DEALLOCATE (f_occ)
3711 : ELSE
3712 564 : e_fermi = VBM
3713 : END IF
3714 :
3715 568 : 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 90 : 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 90 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: f_occ
3738 :
3739 90 : IF (bs_env%do_smear .AND. .NOT. (CBM > VBM .AND. &
3740 : CBM - VBM > gap_smear_ratio*bs_env%smear_width)) THEN
3741 0 : 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 0 : bs_env%smear_width, bs_env%smear_method)
3745 0 : DEALLOCATE (f_occ)
3746 : ELSE
3747 90 : e_fermi = VBM
3748 : END IF
3749 :
3750 90 : END SUBROUTINE get_fermi_level_2c
3751 : END MODULE post_scf_bandstructure_utils
|
