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 cell_types, ONLY: cell_type,&
18 : get_cell,&
19 : pbc
20 : USE cp_blacs_env, ONLY: cp_blacs_env_type
21 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale
22 : USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose
23 : USE cp_cfm_diag, ONLY: cp_cfm_geeig,&
24 : cp_cfm_geeig_canon,&
25 : cp_cfm_heevd
26 : USE cp_cfm_types, ONLY: cp_cfm_create,&
27 : cp_cfm_get_info,&
28 : cp_cfm_release,&
29 : cp_cfm_set_all,&
30 : cp_cfm_to_cfm,&
31 : cp_cfm_to_fm,&
32 : cp_cfm_type,&
33 : cp_fm_to_cfm
34 : USE cp_control_types, ONLY: dft_control_type
35 : USE cp_dbcsr_api, ONLY: &
36 : dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_p_type, dbcsr_set, &
37 : dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
38 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
39 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
40 : copy_fm_to_dbcsr,&
41 : dbcsr_allocate_matrix_set,&
42 : dbcsr_deallocate_matrix_set
43 : USE cp_files, ONLY: close_file,&
44 : open_file
45 : USE cp_fm_diag, ONLY: cp_fm_geeig_canon
46 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
47 : cp_fm_struct_release,&
48 : cp_fm_struct_type
49 : USE cp_fm_types, ONLY: cp_fm_create,&
50 : cp_fm_get_diag,&
51 : cp_fm_get_info,&
52 : cp_fm_release,&
53 : cp_fm_set_all,&
54 : cp_fm_to_fm,&
55 : cp_fm_type
56 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit,&
57 : cp_to_string
58 : USE cp_parser_methods, ONLY: read_float_object
59 : USE input_constants, ONLY: int_ldos_z,&
60 : large_cell_Gamma,&
61 : large_cell_Gamma_ri_rs,&
62 : non_periodic_ri_rs,&
63 : small_cell_full_kp
64 : USE input_section_types, ONLY: section_vals_get,&
65 : section_vals_get_subs_vals,&
66 : section_vals_type,&
67 : section_vals_val_get
68 : USE kinds, ONLY: default_string_length,&
69 : dp,&
70 : max_line_length
71 : USE kpoint_methods, ONLY: kpoint_init_cell_index,&
72 : rskp_transform
73 : USE kpoint_types, ONLY: get_kpoint_info,&
74 : kpoint_create,&
75 : kpoint_type
76 : USE machine, ONLY: m_walltime
77 : USE mathconstants, ONLY: gaussi,&
78 : twopi,&
79 : z_one,&
80 : z_zero
81 : USE message_passing, ONLY: mp_para_env_type
82 : USE parallel_gemm_api, ONLY: parallel_gemm
83 : USE particle_types, ONLY: particle_type
84 : USE physcon, ONLY: angstrom,&
85 : evolt
86 : USE post_scf_bandstructure_types, ONLY: band_edges_type,&
87 : eps_qp_gap,&
88 : max_qp_gap,&
89 : post_scf_bandstructure_type
90 : USE pw_env_types, ONLY: pw_env_get,&
91 : pw_env_type
92 : USE pw_pool_types, ONLY: pw_pool_type
93 : USE pw_types, ONLY: pw_c1d_gs_type,&
94 : pw_r3d_rs_type
95 : USE qs_collocate_density, ONLY: calculate_rho_elec
96 : USE qs_environment_types, ONLY: get_qs_env,&
97 : qs_environment_type
98 : USE qs_ks_types, ONLY: qs_ks_env_type
99 : USE qs_mo_types, ONLY: get_mo_set,&
100 : mo_set_type
101 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
102 : USE rpa_gw_im_time_util, ONLY: compute_weight_re_im,&
103 : get_atom_index_from_basis_function_index
104 : USE scf_control_types, ONLY: scf_control_type
105 : USE soc_pseudopotential_methods, ONLY: V_SOC_xyz_from_pseudopotential,&
106 : remove_soc_outside_energy_window_mo
107 : USE soc_pseudopotential_utils, ONLY: add_cfm_submat,&
108 : add_dbcsr_submat,&
109 : cfm_add_on_diag,&
110 : create_cfm_double,&
111 : get_cfm_submat
112 : USE string_utilities, ONLY: uppercase
113 : #include "base/base_uses.f90"
114 :
115 : IMPLICIT NONE
116 :
117 : PRIVATE
118 :
119 : PUBLIC :: create_and_init_bs_env, &
120 : eval_bandstructure_properties, cfm_ikp_from_fm_Gamma, &
121 : MIC_contribution_from_ikp, compute_xkp, kpoint_init_cell_index_simple, &
122 : rsmat_to_kp, soc, get_VBM_CBM_bandgaps, get_all_VBM_CBM_bandgaps
123 :
124 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'post_scf_bandstructure_utils'
125 :
126 : CONTAINS
127 :
128 : ! **************************************************************************************************
129 : !> \brief ...
130 : !> \param qs_env ...
131 : !> \param bs_env ...
132 : !> \param post_scf_bandstructure_section ...
133 : ! **************************************************************************************************
134 106 : SUBROUTINE create_and_init_bs_env(qs_env, bs_env, post_scf_bandstructure_section)
135 : TYPE(qs_environment_type), POINTER :: qs_env
136 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
137 : TYPE(section_vals_type), POINTER :: post_scf_bandstructure_section
138 :
139 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_and_init_bs_env'
140 :
141 : INTEGER :: handle
142 :
143 106 : CALL timeset(routineN, handle)
144 :
145 10388 : ALLOCATE (bs_env)
146 :
147 106 : CALL print_header(bs_env)
148 :
149 106 : CALL read_bandstructure_input_parameters(bs_env, post_scf_bandstructure_section, qs_env)
150 :
151 106 : CALL get_parameters_from_qs_env(qs_env, bs_env)
152 :
153 106 : CALL set_heuristic_parameters(bs_env)
154 :
155 194 : SELECT CASE (bs_env%small_cell_full_kp_or_large_cell_Gamma)
156 : CASE (large_cell_Gamma, large_cell_Gamma_ri_rs, non_periodic_ri_rs)
157 :
158 88 : CALL setup_kpoints_DOS_large_cell_Gamma(qs_env, bs_env, bs_env%kpoints_DOS)
159 :
160 88 : CALL allocate_and_fill_fm_ks_fm_s(qs_env, bs_env)
161 :
162 88 : CALL diagonalize_ks_matrix(bs_env)
163 :
164 88 : CALL check_positive_definite_overlap_mat(bs_env, qs_env)
165 :
166 : CASE (small_cell_full_kp)
167 :
168 18 : CALL setup_kpoints_scf_desymm(qs_env, bs_env, bs_env%kpoints_scf_desymm, .TRUE.)
169 18 : CALL setup_kpoints_scf_desymm(qs_env, bs_env, bs_env%kpoints_scf_desymm_2, .FALSE.)
170 :
171 18 : CALL setup_kpoints_DOS_small_cell_full_kp(bs_env, bs_env%kpoints_DOS)
172 :
173 18 : CALL allocate_and_fill_fm_ks_fm_s(qs_env, bs_env)
174 :
175 124 : CALL compute_cfm_mo_coeff_kp_and_eigenval_scf_kp(qs_env, bs_env)
176 :
177 : END SELECT
178 :
179 106 : CALL timestop(handle)
180 :
181 106 : END SUBROUTINE create_and_init_bs_env
182 :
183 : ! **************************************************************************************************
184 : !> \brief ...
185 : !> \param bs_env ...
186 : !> \param bs_sec ...
187 : !> \param qs_env ...
188 : ! **************************************************************************************************
189 106 : SUBROUTINE read_bandstructure_input_parameters(bs_env, bs_sec, qs_env)
190 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
191 : TYPE(section_vals_type), POINTER :: bs_sec
192 : TYPE(qs_environment_type), POINTER :: qs_env
193 :
194 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_bandstructure_input_parameters'
195 :
196 : CHARACTER(LEN=default_string_length) :: ustr
197 : CHARACTER(LEN=default_string_length), &
198 106 : DIMENSION(:), POINTER :: string_ptr
199 : CHARACTER(LEN=max_line_length) :: error_msg
200 : INTEGER :: handle, i, ikp
201 : REAL(KIND=dp), DIMENSION(3) :: kpptr
202 : REAL(KIND=dp), DIMENSION(3, 3) :: cart_hmat
203 : TYPE(cell_type), POINTER :: cell
204 : TYPE(section_vals_type), POINTER :: dos_pdos_sec, floquet_sec, gw_sec, &
205 : kp_bs_sec, ldos_sec, soc_sec
206 :
207 106 : CALL timeset(routineN, handle)
208 106 : NULLIFY (cell)
209 106 : CALL get_qs_env(qs_env=qs_env, cell=cell)
210 1378 : cart_hmat(:, :) = cell%hmat(:, :)
211 106 : IF (cell%input_cell_canonicalized) cart_hmat(:, :) = cell%input_hmat(:, :)
212 :
213 106 : NULLIFY (gw_sec)
214 106 : gw_sec => section_vals_get_subs_vals(bs_sec, "GW")
215 106 : CALL section_vals_get(gw_sec, explicit=bs_env%do_gw)
216 106 : CALL section_vals_val_get(gw_sec, "RI_RS", l_val=bs_env%do_gw_ri_rs)
217 :
218 106 : NULLIFY (soc_sec)
219 106 : soc_sec => section_vals_get_subs_vals(bs_sec, "SOC")
220 106 : CALL section_vals_get(soc_sec, explicit=bs_env%do_soc)
221 :
222 106 : CALL section_vals_val_get(soc_sec, "SOC_WINDOW_OCC", r_val=bs_env%soc_window_occ)
223 106 : CALL section_vals_val_get(soc_sec, "SOC_WINDOW_VIRT", r_val=bs_env%soc_window_virt)
224 106 : CALL section_vals_val_get(soc_sec, "SOC_WINDOW_SMEARING", r_val=bs_env%soc_window_smearing)
225 :
226 106 : NULLIFY (dos_pdos_sec)
227 106 : dos_pdos_sec => section_vals_get_subs_vals(bs_sec, "DOS")
228 106 : CALL section_vals_get(dos_pdos_sec, explicit=bs_env%do_dos_pdos)
229 :
230 106 : CALL section_vals_val_get(bs_sec, "DOS%KPOINTS", i_vals=bs_env%nkp_grid_DOS_input)
231 106 : CALL section_vals_val_get(bs_sec, "DOS%ENERGY_WINDOW", r_val=bs_env%energy_window_DOS)
232 106 : CALL section_vals_val_get(bs_sec, "DOS%ENERGY_STEP", r_val=bs_env%energy_step_DOS)
233 106 : CALL section_vals_val_get(bs_sec, "DOS%BROADENING", r_val=bs_env%broadening_DOS)
234 :
235 106 : NULLIFY (ldos_sec)
236 106 : ldos_sec => section_vals_get_subs_vals(bs_sec, "DOS%LDOS")
237 106 : CALL section_vals_get(ldos_sec, explicit=bs_env%do_ldos)
238 :
239 106 : CALL section_vals_val_get(ldos_sec, "INTEGRATION", i_val=bs_env%int_ldos_xyz)
240 106 : CALL section_vals_val_get(ldos_sec, "BIN_MESH", i_vals=bs_env%bin_mesh)
241 :
242 106 : NULLIFY (kp_bs_sec)
243 106 : kp_bs_sec => section_vals_get_subs_vals(bs_sec, "BANDSTRUCTURE_PATH")
244 106 : CALL section_vals_val_get(kp_bs_sec, "NPOINTS", i_val=bs_env%input_kp_bs_npoints)
245 106 : CALL section_vals_val_get(kp_bs_sec, "UNITS", c_val=ustr)
246 106 : CALL uppercase(ustr)
247 106 : CALL section_vals_val_get(kp_bs_sec, "SPECIAL_POINT", n_rep_val=bs_env%input_kp_bs_n_sp_pts)
248 :
249 106 : NULLIFY (floquet_sec)
250 106 : floquet_sec => section_vals_get_subs_vals(bs_sec, "FLOQUET")
251 106 : CALL section_vals_get(floquet_sec, explicit=bs_env%do_floquet)
252 106 : CALL section_vals_val_get(floquet_sec, "AMPLITUDE", r_val=bs_env%floquet_amplitude)
253 106 : CALL section_vals_val_get(floquet_sec, "FREQUENCY", r_val=bs_env%floquet_omega)
254 106 : CALL section_vals_val_get(floquet_sec, "POLARISATION", r_vals=bs_env%floquet_polarisation)
255 106 : CALL section_vals_val_get(floquet_sec, "PHASE_OFFSETS", r_vals=bs_env%floquet_phi)
256 106 : CALL section_vals_val_get(floquet_sec, "MAX_FLOQUET_INDEX", i_val=bs_env%max_floquet_index)
257 106 : CALL section_vals_val_get(floquet_sec, "EPS_FLOQUET", r_val=bs_env%eps_floquet)
258 106 : CALL section_vals_val_get(floquet_sec, "ENERGY_WINDOW", r_val=bs_env%energy_window_floquet)
259 106 : CALL section_vals_val_get(floquet_sec, "ENERGY_STEP", r_val=bs_env%energy_step_floquet)
260 106 : CALL section_vals_val_get(floquet_sec, "BROADENING", r_val=bs_env%broadening_floquet)
261 106 : CALL section_vals_val_get(floquet_sec, "MEM_FILL_FRACTION", r_val=bs_env%floquet_mem_fill_fraction)
262 106 : CALL section_vals_val_get(floquet_sec, "TEMPERATURE", r_val=bs_env%floquet_temperature)
263 106 : CALL section_vals_val_get(floquet_sec, "FLOQUET_DOS_FILE_NAME", c_val=bs_env%floquet_dos_file)
264 106 : CALL section_vals_val_get(floquet_sec, "QUASI_ENERGIES_FILE_NAME", c_val=bs_env%floquet_qe_file)
265 106 : CALL section_vals_val_get(floquet_sec, "FLOQUET_BS_FILE_NAME", c_val=bs_env%floquet_bs_file)
266 :
267 : ! read special points for band structure
268 216 : ALLOCATE (bs_env%xkp_special(3, bs_env%input_kp_bs_n_sp_pts))
269 116 : DO ikp = 1, bs_env%input_kp_bs_n_sp_pts
270 10 : CALL section_vals_val_get(kp_bs_sec, "SPECIAL_POINT", i_rep_val=ikp, c_vals=string_ptr)
271 10 : CPASSERT(SIZE(string_ptr(:), 1) == 4)
272 40 : DO i = 1, 3
273 30 : CALL read_float_object(string_ptr(i + 1), kpptr(i), error_msg)
274 40 : IF (LEN_TRIM(error_msg) > 0) CPABORT(TRIM(error_msg))
275 : END DO
276 106 : SELECT CASE (ustr)
277 : CASE ("B_VECTOR")
278 40 : bs_env%xkp_special(1:3, ikp) = kpptr(1:3)
279 : CASE ("CART_ANGSTROM")
280 : bs_env%xkp_special(1:3, ikp) = (kpptr(1)*cart_hmat(1, 1:3) + &
281 : kpptr(2)*cart_hmat(2, 1:3) + &
282 0 : kpptr(3)*cart_hmat(3, 1:3))/twopi*angstrom
283 : CASE ("CART_BOHR")
284 : bs_env%xkp_special(1:3, ikp) = (kpptr(1)*cart_hmat(1, 1:3) + &
285 : kpptr(2)*cart_hmat(2, 1:3) + &
286 0 : kpptr(3)*cart_hmat(3, 1:3))/twopi
287 : CASE DEFAULT
288 10 : CPABORT("Unknown unit <"//TRIM(ustr)//"> specified for k-point definition")
289 : END SELECT
290 : END DO
291 :
292 106 : CALL timestop(handle)
293 :
294 106 : END SUBROUTINE read_bandstructure_input_parameters
295 :
296 : ! **************************************************************************************************
297 : !> \brief ...
298 : !> \param bs_env ...
299 : ! **************************************************************************************************
300 106 : SUBROUTINE print_header(bs_env)
301 :
302 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
303 :
304 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_header'
305 :
306 : INTEGER :: handle, u
307 :
308 106 : CALL timeset(routineN, handle)
309 :
310 106 : bs_env%unit_nr = cp_logger_get_default_io_unit()
311 :
312 106 : u = bs_env%unit_nr
313 :
314 106 : IF (u > 0) THEN
315 53 : WRITE (u, '(T2,A)') ' '
316 53 : WRITE (u, '(T2,A)') REPEAT('-', 79)
317 53 : WRITE (u, '(T2,A,A78)') '-', '-'
318 53 : WRITE (u, '(T2,A,A51,A27)') '-', 'BANDSTRUCTURE CALCULATION', '-'
319 53 : WRITE (u, '(T2,A,A78)') '-', '-'
320 53 : WRITE (u, '(T2,A)') REPEAT('-', 79)
321 53 : WRITE (u, '(T2,A)') ' '
322 : END IF
323 :
324 106 : CALL timestop(handle)
325 :
326 106 : END SUBROUTINE print_header
327 :
328 : ! **************************************************************************************************
329 : !> \brief ...
330 : !> \param qs_env ...
331 : !> \param bs_env ...
332 : !> \param kpoints ...
333 : ! **************************************************************************************************
334 88 : SUBROUTINE setup_kpoints_DOS_large_cell_Gamma(qs_env, bs_env, kpoints)
335 :
336 : TYPE(qs_environment_type), POINTER :: qs_env
337 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
338 : TYPE(kpoint_type), POINTER :: kpoints
339 :
340 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_kpoints_DOS_large_cell_Gamma'
341 :
342 : INTEGER :: handle, i_dim, i_kp_in_line, &
343 : i_special_kp, ikk, n_kp_in_line, &
344 : n_special_kp, nkp, nkp_only_bs, &
345 : nkp_only_DOS, u
346 : INTEGER, DIMENSION(3) :: nkp_grid, periodic
347 :
348 88 : CALL timeset(routineN, handle)
349 :
350 : ! routine adapted from mp2_integrals.F
351 88 : NULLIFY (kpoints)
352 88 : CALL kpoint_create(kpoints)
353 :
354 88 : kpoints%kp_scheme = "GENERAL"
355 :
356 88 : n_special_kp = bs_env%input_kp_bs_n_sp_pts
357 88 : n_kp_in_line = bs_env%input_kp_bs_npoints
358 :
359 352 : periodic(1:3) = bs_env%periodic(1:3)
360 :
361 352 : DO i_dim = 1, 3
362 :
363 264 : CPASSERT(periodic(i_dim) == 0 .OR. periodic(i_dim) == 1)
364 :
365 352 : IF (bs_env%nkp_grid_DOS_input(i_dim) < 0) THEN
366 246 : IF (periodic(i_dim) == 1) nkp_grid(i_dim) = 2
367 246 : IF (periodic(i_dim) == 0) nkp_grid(i_dim) = 1
368 : ELSE
369 18 : nkp_grid(i_dim) = bs_env%nkp_grid_DOS_input(i_dim)
370 : END IF
371 :
372 : END DO
373 :
374 : ! use the k <-> -k symmetry to reduce the number of kpoints
375 88 : IF (nkp_grid(1) > 1) THEN
376 4 : nkp_only_DOS = (nkp_grid(1) + 1)/2*nkp_grid(2)*nkp_grid(3)
377 84 : ELSE IF (nkp_grid(2) > 1) THEN
378 4 : nkp_only_DOS = nkp_grid(1)*(nkp_grid(2) + 1)/2*nkp_grid(3)
379 80 : ELSE IF (nkp_grid(3) > 1) THEN
380 2 : nkp_only_DOS = nkp_grid(1)*nkp_grid(2)*(nkp_grid(3) + 1)/2
381 : ELSE
382 78 : nkp_only_DOS = 1
383 : END IF
384 :
385 : ! we will compute the GW QP levels for all k's in the bandstructure path but also
386 : ! for all k-points from the SCF (e.g. for DOS or for self-consistent GW)
387 88 : IF (n_special_kp > 0) THEN
388 0 : nkp_only_bs = n_kp_in_line*(n_special_kp - 1) + 1
389 : ELSE
390 : nkp_only_bs = 0
391 : END IF
392 :
393 88 : nkp = nkp_only_DOS + nkp_only_bs
394 :
395 352 : kpoints%nkp_grid(1:3) = nkp_grid(1:3)
396 88 : kpoints%nkp = nkp
397 :
398 88 : bs_env%nkp_bs_and_DOS = nkp
399 88 : bs_env%nkp_only_bs = nkp_only_bs
400 88 : bs_env%nkp_only_DOS = nkp_only_DOS
401 :
402 440 : ALLOCATE (kpoints%xkp(3, nkp), kpoints%wkp(nkp))
403 184 : kpoints%wkp(1:nkp_only_DOS) = 1.0_dp/REAL(nkp_only_DOS, KIND=dp)
404 :
405 88 : CALL compute_xkp(kpoints%xkp, 1, nkp_only_DOS, nkp_grid)
406 :
407 88 : IF (n_special_kp > 0) THEN
408 0 : kpoints%xkp(1:3, nkp_only_DOS + 1) = bs_env%xkp_special(1:3, 1)
409 0 : ikk = nkp_only_DOS + 1
410 0 : DO i_special_kp = 2, n_special_kp
411 0 : DO i_kp_in_line = 1, n_kp_in_line
412 0 : ikk = ikk + 1
413 : kpoints%xkp(1:3, ikk) = bs_env%xkp_special(1:3, i_special_kp - 1) + &
414 : REAL(i_kp_in_line, KIND=dp)/REAL(n_kp_in_line, KIND=dp)* &
415 : (bs_env%xkp_special(1:3, i_special_kp) - &
416 0 : bs_env%xkp_special(1:3, i_special_kp - 1))
417 0 : kpoints%wkp(ikk) = 0.0_dp
418 : END DO
419 : END DO
420 : END IF
421 :
422 88 : CALL kpoint_init_cell_index_simple(kpoints, qs_env)
423 :
424 88 : u = bs_env%unit_nr
425 :
426 88 : IF (u > 0) THEN
427 44 : IF (nkp_only_bs > 0) THEN
428 : WRITE (u, FMT="(T2,1A,T77,I4)") &
429 0 : "Number of special k-points for the bandstructure", n_special_kp
430 0 : WRITE (u, FMT="(T2,1A,T77,I4)") "Number of k-points for the bandstructure", nkp
431 : WRITE (u, FMT="(T2,1A,T69,3I4)") &
432 0 : "K-point mesh for the density of states (DOS)", nkp_grid(1:3)
433 : ELSE
434 : WRITE (u, FMT="(T2,1A,T69,3I4)") &
435 44 : "K-point mesh for the density of states (DOS) and the self-energy", nkp_grid(1:3)
436 : END IF
437 : END IF
438 :
439 88 : CALL timestop(handle)
440 :
441 88 : END SUBROUTINE setup_kpoints_DOS_large_cell_Gamma
442 :
443 : ! **************************************************************************************************
444 : !> \brief ...
445 : !> \param qs_env ...
446 : !> \param bs_env ...
447 : !> \param kpoints ...
448 : !> \param do_print ...
449 : ! **************************************************************************************************
450 36 : SUBROUTINE setup_kpoints_scf_desymm(qs_env, bs_env, kpoints, do_print)
451 : TYPE(qs_environment_type), POINTER :: qs_env
452 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
453 : TYPE(kpoint_type), POINTER :: kpoints
454 :
455 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_kpoints_scf_desymm'
456 :
457 : INTEGER :: handle, i_cell_x, i_dim, img, j_cell_y, &
458 : k_cell_z, nimages, nkp, u
459 : INTEGER, DIMENSION(3) :: cell_grid, cixd, nkp_grid
460 : TYPE(kpoint_type), POINTER :: kpoints_scf
461 :
462 : LOGICAL:: do_print
463 :
464 36 : CALL timeset(routineN, handle)
465 :
466 36 : NULLIFY (kpoints)
467 36 : CALL kpoint_create(kpoints)
468 :
469 36 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints_scf)
470 :
471 144 : nkp_grid(1:3) = kpoints_scf%nkp_grid(1:3)
472 36 : nkp = nkp_grid(1)*nkp_grid(2)*nkp_grid(3)
473 :
474 : ! we need in periodic directions at least 4 k-points in the SCF
475 144 : DO i_dim = 1, 3
476 144 : IF (bs_env%periodic(i_dim) == 1) THEN
477 72 : CPASSERT(nkp_grid(i_dim) >= 4)
478 : END IF
479 : END DO
480 :
481 36 : kpoints%kp_scheme = "GENERAL"
482 144 : kpoints%nkp_grid(1:3) = nkp_grid(1:3)
483 36 : kpoints%nkp = nkp
484 36 : bs_env%nkp_scf_desymm = nkp
485 :
486 108 : ALLOCATE (kpoints%xkp(1:3, nkp))
487 36 : CALL compute_xkp(kpoints%xkp, 1, nkp, nkp_grid)
488 :
489 108 : ALLOCATE (kpoints%wkp(nkp))
490 612 : kpoints%wkp(:) = 1.0_dp/REAL(nkp, KIND=dp)
491 :
492 : ! for example 4x3x6 kpoint grid -> 3x3x5 cell grid because we need the same number of
493 : ! neighbor cells on both sides of the unit cell
494 144 : cell_grid(1:3) = nkp_grid(1:3) - MODULO(nkp_grid(1:3) + 1, 2)
495 :
496 : ! cell index: for example for x: from -n_x/2 to +n_x/2, n_x: number of cells in x direction
497 144 : cixd(1:3) = cell_grid(1:3)/2
498 :
499 36 : nimages = cell_grid(1)*cell_grid(2)*cell_grid(3)
500 :
501 36 : bs_env%nimages_scf_desymm = nimages
502 144 : bs_env%cell_grid_scf_desymm(1:3) = cell_grid(1:3)
503 :
504 36 : IF (ASSOCIATED(kpoints%index_to_cell)) DEALLOCATE (kpoints%index_to_cell)
505 36 : IF (ASSOCIATED(kpoints%cell_to_index)) DEALLOCATE (kpoints%cell_to_index)
506 :
507 180 : ALLOCATE (kpoints%cell_to_index(-cixd(1):cixd(1), -cixd(2):cixd(2), -cixd(3):cixd(3)))
508 108 : ALLOCATE (kpoints%index_to_cell(3, nimages))
509 :
510 36 : img = 0
511 88 : DO i_cell_x = -cixd(1), cixd(1)
512 244 : DO j_cell_y = -cixd(2), cixd(2)
513 532 : DO k_cell_z = -cixd(3), cixd(3)
514 324 : img = img + 1
515 324 : kpoints%cell_to_index(i_cell_x, j_cell_y, k_cell_z) = img
516 1452 : kpoints%index_to_cell(1:3, img) = [i_cell_x, j_cell_y, k_cell_z]
517 : END DO
518 : END DO
519 : END DO
520 :
521 36 : u = bs_env%unit_nr
522 36 : IF (u > 0 .AND. do_print) THEN
523 9 : WRITE (u, FMT="(T2,A,I49)") "Number of cells for G, χ, W, Σ", nimages
524 : END IF
525 :
526 36 : CALL timestop(handle)
527 :
528 36 : END SUBROUTINE setup_kpoints_scf_desymm
529 :
530 : ! **************************************************************************************************
531 : !> \brief ...
532 : !> \param bs_env ...
533 : !> \param kpoints ...
534 : ! **************************************************************************************************
535 18 : SUBROUTINE setup_kpoints_DOS_small_cell_full_kp(bs_env, kpoints)
536 :
537 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
538 : TYPE(kpoint_type), POINTER :: kpoints
539 :
540 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_kpoints_DOS_small_cell_full_kp'
541 :
542 : INTEGER :: handle, i_kp_in_line, i_special_kp, ikk, &
543 : n_kp_in_line, n_special_kp, nkp, &
544 : nkp_only_bs, nkp_scf_desymm, u
545 :
546 18 : CALL timeset(routineN, handle)
547 :
548 : ! routine adapted from mp2_integrals.F
549 18 : NULLIFY (kpoints)
550 18 : CALL kpoint_create(kpoints)
551 :
552 18 : n_special_kp = bs_env%input_kp_bs_n_sp_pts
553 18 : n_kp_in_line = bs_env%input_kp_bs_npoints
554 18 : nkp_scf_desymm = bs_env%nkp_scf_desymm
555 :
556 : ! we will compute the GW QP levels for all k's in the bandstructure path but also
557 : ! for all k-points from the SCF (e.g. for DOS or for self-consistent GW)
558 18 : IF (n_special_kp > 0) THEN
559 4 : nkp_only_bs = n_kp_in_line*(n_special_kp - 1) + 1
560 : ELSE
561 : nkp_only_bs = 0
562 : END IF
563 18 : nkp = nkp_only_bs + nkp_scf_desymm
564 :
565 54 : ALLOCATE (kpoints%xkp(3, nkp))
566 54 : ALLOCATE (kpoints%wkp(nkp))
567 :
568 18 : kpoints%nkp = nkp
569 :
570 18 : bs_env%nkp_bs_and_DOS = nkp
571 18 : bs_env%nkp_only_bs = nkp_only_bs
572 18 : bs_env%nkp_only_DOS = nkp_scf_desymm
573 :
574 2340 : kpoints%xkp(1:3, 1:nkp_scf_desymm) = bs_env%kpoints_scf_desymm%xkp(1:3, 1:nkp_scf_desymm)
575 306 : kpoints%wkp(1:nkp_scf_desymm) = 1.0_dp/REAL(nkp_scf_desymm, KIND=dp)
576 :
577 18 : IF (n_special_kp > 0) THEN
578 32 : kpoints%xkp(1:3, nkp_scf_desymm + 1) = bs_env%xkp_special(1:3, 1)
579 4 : ikk = nkp_scf_desymm + 1
580 10 : DO i_special_kp = 2, n_special_kp
581 70 : DO i_kp_in_line = 1, n_kp_in_line
582 60 : ikk = ikk + 1
583 : kpoints%xkp(1:3, ikk) = bs_env%xkp_special(1:3, i_special_kp - 1) + &
584 : REAL(i_kp_in_line, KIND=dp)/REAL(n_kp_in_line, KIND=dp)* &
585 : (bs_env%xkp_special(1:3, i_special_kp) - &
586 480 : bs_env%xkp_special(1:3, i_special_kp - 1))
587 66 : kpoints%wkp(ikk) = 0.0_dp
588 : END DO
589 : END DO
590 : END IF
591 :
592 18 : IF (ASSOCIATED(kpoints%index_to_cell)) DEALLOCATE (kpoints%index_to_cell)
593 :
594 54 : ALLOCATE (kpoints%index_to_cell(3, bs_env%nimages_scf_desymm))
595 1332 : kpoints%index_to_cell(:, :) = bs_env%kpoints_scf_desymm%index_to_cell(:, :)
596 :
597 18 : u = bs_env%unit_nr
598 :
599 18 : IF (u > 0) THEN
600 9 : WRITE (u, FMT="(T2,1A,T77,I4)") "Number of special k-points for the bandstructure", &
601 18 : n_special_kp
602 9 : WRITE (u, FMT="(T2,1A,T77,I4)") "Number of k-points for the bandstructure", nkp
603 : END IF
604 :
605 18 : CALL timestop(handle)
606 :
607 18 : END SUBROUTINE setup_kpoints_DOS_small_cell_full_kp
608 :
609 : ! **************************************************************************************************
610 : !> \brief ...
611 : !> \param qs_env ...
612 : !> \param bs_env ...
613 : ! **************************************************************************************************
614 18 : SUBROUTINE compute_cfm_mo_coeff_kp_and_eigenval_scf_kp(qs_env, bs_env)
615 : TYPE(qs_environment_type), POINTER :: qs_env
616 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
617 :
618 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_cfm_mo_coeff_kp_and_eigenval_scf_kp'
619 :
620 : INTEGER :: handle, ikp, ispin, nkp_bs_and_DOS
621 18 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_scf
622 : REAL(KIND=dp) :: CBM, VBM
623 : REAL(KIND=dp), DIMENSION(3) :: xkp
624 : TYPE(cp_cfm_type) :: cfm_ks, cfm_mos, cfm_s
625 18 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
626 : TYPE(kpoint_type), POINTER :: kpoints_scf
627 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
628 18 : POINTER :: sab_nl
629 :
630 18 : CALL timeset(routineN, handle)
631 :
632 : CALL get_qs_env(qs_env, &
633 : matrix_ks_kp=matrix_ks, &
634 : matrix_s_kp=matrix_s, &
635 18 : kpoints=kpoints_scf)
636 :
637 18 : NULLIFY (sab_nl)
638 18 : CALL get_kpoint_info(kpoints_scf, sab_nl=sab_nl, cell_to_index=cell_to_index_scf)
639 :
640 18 : CALL cp_cfm_create(cfm_ks, bs_env%cfm_work_mo%matrix_struct)
641 18 : CALL cp_cfm_create(cfm_s, bs_env%cfm_work_mo%matrix_struct)
642 18 : CALL cp_cfm_create(cfm_mos, bs_env%cfm_work_mo%matrix_struct)
643 :
644 : ! nkp_bs_and_DOS contains desymmetrized k-point mesh from SCF and k-points from GW bandstructure
645 18 : nkp_bs_and_DOS = bs_env%nkp_bs_and_DOS
646 :
647 90 : ALLOCATE (bs_env%eigenval_G0W0(bs_env%n_ao, nkp_bs_and_DOS, bs_env%n_spin))
648 90 : ALLOCATE (bs_env%eigenval_HF(bs_env%n_ao, nkp_bs_and_DOS, bs_env%n_spin))
649 442 : ALLOCATE (bs_env%cfm_mo_coeff_kp(nkp_bs_and_DOS, bs_env%n_spin))
650 442 : ALLOCATE (bs_env%cfm_ks_kp(nkp_bs_and_DOS, bs_env%n_spin))
651 406 : ALLOCATE (bs_env%cfm_s_kp(nkp_bs_and_DOS))
652 370 : DO ikp = 1, nkp_bs_and_DOS
653 704 : DO ispin = 1, bs_env%n_spin
654 352 : CALL cp_cfm_create(bs_env%cfm_mo_coeff_kp(ikp, ispin), bs_env%cfm_work_mo%matrix_struct)
655 704 : CALL cp_cfm_create(bs_env%cfm_ks_kp(ikp, ispin), bs_env%cfm_work_mo%matrix_struct)
656 : END DO
657 370 : CALL cp_cfm_create(bs_env%cfm_s_kp(ikp), bs_env%cfm_work_mo%matrix_struct)
658 : END DO
659 :
660 36 : DO ispin = 1, bs_env%n_spin
661 370 : DO ikp = 1, nkp_bs_and_DOS
662 :
663 1408 : xkp(1:3) = bs_env%kpoints_DOS%xkp(1:3, ikp)
664 :
665 : ! h^KS^R -> h^KS(k)
666 352 : CALL rsmat_to_kp(matrix_ks, ispin, xkp, cell_to_index_scf, sab_nl, bs_env, cfm_ks)
667 :
668 : ! S^R -> S(k)
669 352 : CALL rsmat_to_kp(matrix_s, 1, xkp, cell_to_index_scf, sab_nl, bs_env, cfm_s)
670 :
671 : ! we store the complex KS matrix as fm matrix because the infrastructure for fm is
672 : ! much nicer compared to cfm
673 352 : CALL cp_cfm_to_cfm(cfm_ks, bs_env%cfm_ks_kp(ikp, ispin))
674 352 : CALL cp_cfm_to_cfm(cfm_s, bs_env%cfm_s_kp(ikp))
675 :
676 : ! Diagonalize KS-matrix via Rothaan-Hall equation:
677 : ! H^KS(k) C(k) = S(k) C(k) ε(k)
678 : CALL cp_cfm_geeig_canon(cfm_ks, cfm_s, cfm_mos, &
679 : bs_env%eigenval_scf(:, ikp, ispin), &
680 352 : bs_env%cfm_work_mo, bs_env%eps_eigval_mat_s)
681 :
682 : ! we store the complex MO coeff as fm matrix because the infrastructure for fm is
683 : ! much nicer compared to cfm
684 370 : CALL cp_cfm_to_cfm(cfm_mos, bs_env%cfm_mo_coeff_kp(ikp, ispin))
685 :
686 : END DO
687 :
688 370 : VBM = MAXVAL(bs_env%eigenval_scf(bs_env%n_occ(ispin), :, ispin))
689 370 : CBM = MINVAL(bs_env%eigenval_scf(bs_env%n_occ(ispin) + 1, :, ispin))
690 :
691 36 : bs_env%e_fermi(ispin) = 0.5_dp*(VBM + CBM)
692 :
693 : END DO
694 :
695 18 : CALL get_VBM_CBM_bandgaps(bs_env%band_edges_scf, bs_env%eigenval_scf, bs_env)
696 :
697 18 : CALL cp_cfm_release(cfm_ks)
698 18 : CALL cp_cfm_release(cfm_s)
699 18 : CALL cp_cfm_release(cfm_mos)
700 :
701 18 : CALL timestop(handle)
702 :
703 36 : END SUBROUTINE compute_cfm_mo_coeff_kp_and_eigenval_scf_kp
704 :
705 : ! **************************************************************************************************
706 : !> \brief ...
707 : !> \param mat_rs ...
708 : !> \param ispin ...
709 : !> \param xkp ...
710 : !> \param cell_to_index_scf ...
711 : !> \param sab_nl ...
712 : !> \param bs_env ...
713 : !> \param cfm_kp ...
714 : !> \param imag_rs_mat ...
715 : ! **************************************************************************************************
716 1976 : SUBROUTINE rsmat_to_kp(mat_rs, ispin, xkp, cell_to_index_scf, sab_nl, bs_env, cfm_kp, imag_rs_mat)
717 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_rs
718 : INTEGER :: ispin
719 : REAL(KIND=dp), DIMENSION(3) :: xkp
720 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_scf
721 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
722 : POINTER :: sab_nl
723 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
724 : TYPE(cp_cfm_type) :: cfm_kp
725 : LOGICAL, OPTIONAL :: imag_rs_mat
726 :
727 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rsmat_to_kp'
728 :
729 : INTEGER :: handle
730 : LOGICAL :: imag_rs_mat_private
731 : TYPE(dbcsr_type), POINTER :: cmat, nsmat, rmat
732 :
733 1976 : CALL timeset(routineN, handle)
734 :
735 1976 : ALLOCATE (rmat, cmat, nsmat)
736 :
737 1976 : imag_rs_mat_private = .FALSE.
738 1976 : IF (PRESENT(imag_rs_mat)) imag_rs_mat_private = imag_rs_mat
739 :
740 954 : IF (imag_rs_mat_private) THEN
741 954 : CALL dbcsr_create(rmat, template=mat_rs(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
742 954 : CALL dbcsr_create(cmat, template=mat_rs(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
743 : ELSE
744 1022 : CALL dbcsr_create(rmat, template=mat_rs(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
745 1022 : CALL dbcsr_create(cmat, template=mat_rs(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
746 : END IF
747 1976 : CALL dbcsr_create(nsmat, template=mat_rs(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
748 1976 : CALL cp_dbcsr_alloc_block_from_nbl(rmat, sab_nl)
749 1976 : CALL cp_dbcsr_alloc_block_from_nbl(cmat, sab_nl)
750 :
751 1976 : CALL dbcsr_set(rmat, 0.0_dp)
752 1976 : CALL dbcsr_set(cmat, 0.0_dp)
753 : CALL rskp_transform(rmatrix=rmat, cmatrix=cmat, rsmat=mat_rs, ispin=ispin, &
754 1976 : xkp=xkp, cell_to_index=cell_to_index_scf, sab_nl=sab_nl)
755 :
756 1976 : CALL dbcsr_desymmetrize(rmat, nsmat)
757 1976 : CALL copy_dbcsr_to_fm(nsmat, bs_env%fm_work_mo(1))
758 1976 : CALL dbcsr_desymmetrize(cmat, nsmat)
759 1976 : CALL copy_dbcsr_to_fm(nsmat, bs_env%fm_work_mo(2))
760 1976 : CALL cp_fm_to_cfm(bs_env%fm_work_mo(1), bs_env%fm_work_mo(2), cfm_kp)
761 :
762 1976 : CALL dbcsr_deallocate_matrix(rmat)
763 1976 : CALL dbcsr_deallocate_matrix(cmat)
764 1976 : CALL dbcsr_deallocate_matrix(nsmat)
765 :
766 1976 : CALL timestop(handle)
767 :
768 1976 : END SUBROUTINE rsmat_to_kp
769 :
770 : ! **************************************************************************************************
771 : !> \brief ...
772 : !> \param bs_env ...
773 : ! **************************************************************************************************
774 88 : SUBROUTINE diagonalize_ks_matrix(bs_env)
775 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
776 :
777 : CHARACTER(LEN=*), PARAMETER :: routineN = 'diagonalize_ks_matrix'
778 :
779 : INTEGER :: handle, ikp, ispin
780 : REAL(KIND=dp) :: CBM, VBM
781 :
782 88 : CALL timeset(routineN, handle)
783 :
784 352 : ALLOCATE (bs_env%eigenval_scf_Gamma(bs_env%n_ao, bs_env%n_spin))
785 :
786 190 : DO ispin = 1, bs_env%n_spin
787 :
788 : ! use work matrices because the matrices are overwritten in cp_fm_geeig_canon
789 102 : CALL cp_fm_to_fm(bs_env%fm_ks_Gamma(ispin), bs_env%fm_work_mo(1))
790 102 : CALL cp_fm_to_fm(bs_env%fm_s_Gamma, bs_env%fm_work_mo(2))
791 :
792 : ! diagonalize the Kohn-Sham matrix to obtain MO coefficients and SCF eigenvalues
793 : ! (at the Gamma-point)
794 : CALL cp_fm_geeig_canon(bs_env%fm_work_mo(1), &
795 : bs_env%fm_work_mo(2), &
796 : bs_env%fm_mo_coeff_Gamma(ispin), &
797 : bs_env%eigenval_scf_Gamma(:, ispin), &
798 : bs_env%fm_work_mo(3), &
799 102 : bs_env%eps_eigval_mat_s)
800 :
801 102 : VBM = bs_env%eigenval_scf_Gamma(bs_env%n_occ(ispin), ispin)
802 102 : CBM = bs_env%eigenval_scf_Gamma(bs_env%n_occ(ispin) + 1, ispin)
803 :
804 102 : bs_env%band_edges_scf_Gamma(ispin)%VBM = VBM
805 102 : bs_env%band_edges_scf_Gamma(ispin)%CBM = CBM
806 190 : bs_env%e_fermi(ispin) = 0.5_dp*(VBM + CBM)
807 :
808 : END DO
809 :
810 88 : CALL timestop(handle)
811 :
812 : ! Gamma-only path for molecules: eigenval_scf is filled here from the Gamma eigenvalues
813 190 : DO ispin = 1, bs_env%n_spin
814 304 : DO ikp = 1, bs_env%nkp_bs_and_DOS
815 1420 : bs_env%eigenval_scf(:, ikp, ispin) = bs_env%eigenval_scf_Gamma(:, ispin)
816 : END DO
817 : END DO
818 :
819 88 : END SUBROUTINE diagonalize_ks_matrix
820 :
821 : ! **************************************************************************************************
822 : !> \brief ...
823 : !> \param bs_env ...
824 : !> \param qs_env ...
825 : ! **************************************************************************************************
826 88 : SUBROUTINE check_positive_definite_overlap_mat(bs_env, qs_env)
827 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
828 : TYPE(qs_environment_type), POINTER :: qs_env
829 :
830 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_positive_definite_overlap_mat'
831 :
832 : INTEGER :: handle, ikp, info, u
833 : TYPE(cp_cfm_type) :: cfm_s_ikp
834 :
835 88 : CALL timeset(routineN, handle)
836 :
837 184 : DO ikp = 1, bs_env%kpoints_DOS%nkp
838 :
839 : ! get S_µν(k_i) from S_µν(k=0)
840 : CALL cfm_ikp_from_fm_Gamma(cfm_s_ikp, bs_env%fm_s_Gamma, &
841 96 : ikp, qs_env, bs_env%kpoints_DOS, "ORB")
842 :
843 : ! check whether S_µν(k_i) is positive definite
844 96 : CALL cp_cfm_cholesky_decompose(matrix=cfm_s_ikp, n=bs_env%n_ao, info_out=info)
845 :
846 : ! check if Cholesky decomposition failed (Cholesky decomposition only works for
847 : ! positive definite matrices
848 184 : IF (info /= 0) THEN
849 0 : u = bs_env%unit_nr
850 :
851 0 : IF (u > 0) THEN
852 0 : WRITE (u, FMT="(T2,A)") ""
853 : WRITE (u, FMT="(T2,A)") "ERROR: The Cholesky decomposition "// &
854 0 : "of the k-point overlap matrix failed. This is"
855 : WRITE (u, FMT="(T2,A)") "because the algorithm is "// &
856 0 : "only correct in the limit of large cells. The cell of "
857 : WRITE (u, FMT="(T2,A)") "the calculation is too small. "// &
858 0 : "Use MULTIPLE_UNIT_CELL to create a larger cell "
859 0 : WRITE (u, FMT="(T2,A)") "and to prevent this error."
860 : END IF
861 :
862 0 : CALL bs_env%para_env%sync()
863 0 : CPABORT("Please see information on the error above.")
864 :
865 : END IF ! Cholesky decomposition failed
866 :
867 : END DO ! ikp
868 :
869 88 : CALL cp_cfm_release(cfm_s_ikp)
870 :
871 88 : CALL timestop(handle)
872 :
873 88 : END SUBROUTINE check_positive_definite_overlap_mat
874 :
875 : ! **************************************************************************************************
876 : !> \brief ...
877 : !> \param qs_env ...
878 : !> \param bs_env ...
879 : ! **************************************************************************************************
880 212 : SUBROUTINE get_parameters_from_qs_env(qs_env, bs_env)
881 : TYPE(qs_environment_type), POINTER :: qs_env
882 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
883 :
884 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_parameters_from_qs_env'
885 :
886 : INTEGER :: color_sub, handle, homo, n_ao, n_atom, u
887 : INTEGER, DIMENSION(3) :: periodic
888 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
889 : TYPE(cell_type), POINTER :: cell
890 : TYPE(dft_control_type), POINTER :: dft_control
891 106 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
892 : TYPE(mp_para_env_type), POINTER :: para_env
893 106 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
894 : TYPE(scf_control_type), POINTER :: scf_control
895 : TYPE(section_vals_type), POINTER :: input
896 :
897 106 : CALL timeset(routineN, handle)
898 :
899 : CALL get_qs_env(qs_env, &
900 : dft_control=dft_control, &
901 : scf_control=scf_control, &
902 106 : mos=mos)
903 :
904 106 : bs_env%n_spin = dft_control%nspins
905 106 : IF (bs_env%n_spin == 1) bs_env%spin_degeneracy = 2.0_dp
906 106 : IF (bs_env%n_spin == 2) bs_env%spin_degeneracy = 1.0_dp
907 :
908 106 : CALL get_mo_set(mo_set=mos(1), nao=n_ao, homo=homo)
909 106 : bs_env%n_ao = n_ao
910 318 : bs_env%n_occ(1:2) = homo
911 318 : bs_env%n_vir(1:2) = n_ao - homo
912 :
913 106 : IF (bs_env%n_spin == 2) THEN
914 14 : CALL get_mo_set(mo_set=mos(2), homo=homo)
915 14 : bs_env%n_occ(2) = homo
916 14 : bs_env%n_vir(2) = n_ao - homo
917 : END IF
918 :
919 106 : bs_env%eps_eigval_mat_s = scf_control%eps_eigval
920 :
921 : ! get para_env from qs_env (bs_env%para_env is identical to para_env in qs_env)
922 106 : CALL get_qs_env(qs_env, para_env=para_env)
923 106 : color_sub = 0
924 106 : ALLOCATE (bs_env%para_env)
925 106 : CALL bs_env%para_env%from_split(para_env, color_sub)
926 :
927 106 : CALL get_qs_env(qs_env, particle_set=particle_set)
928 :
929 106 : n_atom = SIZE(particle_set)
930 106 : bs_env%n_atom = n_atom
931 :
932 106 : CALL get_qs_env(qs_env=qs_env, cell=cell)
933 106 : CALL get_cell(cell=cell, periodic=periodic, h=hmat)
934 424 : bs_env%periodic(1:3) = periodic(1:3)
935 1378 : bs_env%hmat(1:3, 1:3) = hmat
936 106 : bs_env%nimages_scf = dft_control%nimages
937 106 : IF (dft_control%nimages == 1) THEN
938 88 : IF (bs_env%do_gw_ri_rs) THEN
939 80 : IF (ANY(periodic /= 0)) THEN
940 0 : bs_env%small_cell_full_kp_or_large_cell_Gamma = large_cell_Gamma_ri_rs
941 : ELSE
942 20 : bs_env%small_cell_full_kp_or_large_cell_Gamma = non_periodic_ri_rs
943 : END IF
944 : ELSE
945 68 : bs_env%small_cell_full_kp_or_large_cell_Gamma = large_cell_Gamma
946 : END IF
947 18 : ELSE IF (dft_control%nimages > 1) THEN
948 18 : IF (bs_env%do_gw_ri_rs) THEN
949 0 : CPABORT("RI-RS Not Implemented for K-point Calculations")
950 : ELSE
951 18 : bs_env%small_cell_full_kp_or_large_cell_Gamma = small_cell_full_kp
952 : END IF
953 : ELSE
954 0 : CPABORT("Wrong number of cells from DFT calculation.")
955 : END IF
956 :
957 106 : u = bs_env%unit_nr
958 :
959 : ! Marek : Get and save the rtp method
960 106 : CALL get_qs_env(qs_env=qs_env, input=input)
961 106 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%_SECTION_PARAMETERS_", i_val=bs_env%rtp_method)
962 :
963 106 : IF (u > 0) THEN
964 53 : WRITE (u, FMT="(T2,2A,T73,I8)") "Number of occupied molecular orbitals (MOs) ", &
965 106 : "= Number of occupied bands", homo
966 53 : WRITE (u, FMT="(T2,2A,T73,I8)") "Number of unoccupied (= virtual) MOs ", &
967 106 : "= Number of unoccupied bands", n_ao - homo
968 53 : WRITE (u, FMT="(T2,A,T73,I8)") "Number of Gaussian basis functions for MOs", n_ao
969 53 : IF (bs_env%small_cell_full_kp_or_large_cell_Gamma == small_cell_full_kp) THEN
970 9 : WRITE (u, FMT="(T2,2A,T73,I8)") "Number of cells considered in the DFT ", &
971 18 : "calculation", bs_env%nimages_scf
972 : END IF
973 : END IF
974 :
975 106 : CALL timestop(handle)
976 :
977 106 : END SUBROUTINE get_parameters_from_qs_env
978 :
979 : ! **************************************************************************************************
980 : !> \brief ...
981 : !> \param bs_env ...
982 : ! **************************************************************************************************
983 106 : SUBROUTINE set_heuristic_parameters(bs_env)
984 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
985 :
986 : CHARACTER(LEN=*), PARAMETER :: routineN = 'set_heuristic_parameters'
987 :
988 : INTEGER :: handle
989 :
990 106 : CALL timeset(routineN, handle)
991 :
992 106 : bs_env%n_bins_max_for_printing = 5000
993 :
994 106 : CALL timestop(handle)
995 :
996 106 : END SUBROUTINE set_heuristic_parameters
997 :
998 : ! **************************************************************************************************
999 : !> \brief ...
1000 : !> \param qs_env ...
1001 : !> \param bs_env ...
1002 : ! **************************************************************************************************
1003 106 : SUBROUTINE allocate_and_fill_fm_ks_fm_s(qs_env, bs_env)
1004 : TYPE(qs_environment_type), POINTER :: qs_env
1005 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1006 :
1007 : CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_and_fill_fm_ks_fm_s'
1008 :
1009 : INTEGER :: handle, i_work, ispin
1010 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1011 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1012 106 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
1013 : TYPE(mp_para_env_type), POINTER :: para_env
1014 :
1015 106 : CALL timeset(routineN, handle)
1016 :
1017 : CALL get_qs_env(qs_env, &
1018 : para_env=para_env, &
1019 : blacs_env=blacs_env, &
1020 : matrix_ks_kp=matrix_ks, &
1021 106 : matrix_s_kp=matrix_s)
1022 :
1023 106 : NULLIFY (fm_struct)
1024 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=bs_env%n_ao, &
1025 106 : ncol_global=bs_env%n_ao, para_env=para_env)
1026 :
1027 530 : DO i_work = 1, SIZE(bs_env%fm_work_mo)
1028 530 : CALL cp_fm_create(bs_env%fm_work_mo(i_work), fm_struct)
1029 : END DO
1030 :
1031 106 : CALL cp_cfm_create(bs_env%cfm_work_mo, fm_struct)
1032 106 : CALL cp_cfm_create(bs_env%cfm_work_mo_2, fm_struct)
1033 :
1034 106 : CALL cp_fm_create(bs_env%fm_s_Gamma, fm_struct)
1035 106 : CALL copy_dbcsr_to_fm(matrix_s(1, 1)%matrix, bs_env%fm_s_Gamma)
1036 :
1037 226 : DO ispin = 1, bs_env%n_spin
1038 120 : CALL cp_fm_create(bs_env%fm_ks_Gamma(ispin), fm_struct)
1039 120 : CALL copy_dbcsr_to_fm(matrix_ks(ispin, 1)%matrix, bs_env%fm_ks_Gamma(ispin))
1040 226 : CALL cp_fm_create(bs_env%fm_mo_coeff_Gamma(ispin), fm_struct)
1041 : END DO
1042 :
1043 106 : CALL cp_fm_struct_release(fm_struct)
1044 :
1045 106 : NULLIFY (bs_env%mat_ao_ao%matrix)
1046 106 : ALLOCATE (bs_env%mat_ao_ao%matrix)
1047 : CALL dbcsr_create(bs_env%mat_ao_ao%matrix, template=matrix_s(1, 1)%matrix, &
1048 106 : matrix_type=dbcsr_type_no_symmetry)
1049 :
1050 530 : ALLOCATE (bs_env%eigenval_scf(bs_env%n_ao, bs_env%nkp_bs_and_DOS, bs_env%n_spin))
1051 :
1052 106 : CALL timestop(handle)
1053 :
1054 106 : END SUBROUTINE allocate_and_fill_fm_ks_fm_s
1055 :
1056 : ! **************************************************************************************************
1057 : !> \brief ...
1058 : !> \param qs_env ...
1059 : !> \param bs_env ...
1060 : ! **************************************************************************************************
1061 106 : SUBROUTINE eval_bandstructure_properties(qs_env, bs_env)
1062 : TYPE(qs_environment_type), POINTER :: qs_env
1063 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1064 :
1065 : CHARACTER(LEN=*), PARAMETER :: routineN = 'eval_bandstructure_properties'
1066 :
1067 : INTEGER :: handle, homo, homo_1, homo_2, &
1068 : homo_spinor, ikp, ikp_for_file, ispin, &
1069 : n_ao, n_E, nkind, nkp
1070 : LOGICAL :: is_bandstruc_kpoint, print_DOS_kpoints, &
1071 : print_ikp
1072 : REAL(KIND=dp) :: broadening, E_max, E_max_G0W0, E_min, &
1073 : E_min_G0W0, E_total_window, &
1074 : energy_step_DOS, energy_window_DOS, t1
1075 106 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: DOS_G0W0, DOS_G0W0_SOC, DOS_scf, DOS_scf_SOC, &
1076 106 : eigenval, eigenval_spinor, eigenval_spinor_G0W0, eigenval_spinor_no_SOC
1077 106 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: PDOS_G0W0, PDOS_G0W0_SOC, PDOS_scf, &
1078 106 : PDOS_scf_SOC, proj_mo_on_kind
1079 106 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: LDOS_G0W0_2d, LDOS_scf_2d, &
1080 106 : LDOS_scf_2d_SOC
1081 : TYPE(band_edges_type) :: band_edges_G0W0, band_edges_G0W0_SOC, &
1082 : band_edges_scf, band_edges_scf_guess, &
1083 : band_edges_scf_SOC
1084 : TYPE(cp_cfm_type) :: cfm_ks_ikp, cfm_ks_ikp_spinor, cfm_mos_ikp_spinor, cfm_s_ikp, &
1085 : cfm_s_ikp_copy, cfm_s_ikp_spinor, cfm_s_ikp_spinor_copy, cfm_SOC_ikp_spinor, &
1086 : cfm_spinor_wf_ikp, cfm_work_ikp, cfm_work_ikp_spinor
1087 318 : TYPE(cp_cfm_type), DIMENSION(2) :: cfm_mos_ikp
1088 :
1089 106 : CALL timeset(routineN, handle)
1090 :
1091 106 : n_ao = bs_env%n_ao
1092 :
1093 106 : energy_window_DOS = bs_env%energy_window_DOS
1094 106 : energy_step_DOS = bs_env%energy_step_DOS
1095 106 : broadening = bs_env%broadening_DOS
1096 :
1097 : ! if we have done GW or a full kpoint SCF, we already have the band edges
1098 106 : IF (bs_env%do_gw .OR. &
1099 : bs_env%small_cell_full_kp_or_large_cell_Gamma == small_cell_full_kp) THEN
1100 106 : band_edges_scf = bs_env%band_edges_scf
1101 106 : band_edges_scf_guess = band_edges_scf
1102 : ELSE
1103 :
1104 0 : IF (bs_env%n_spin == 1) THEN
1105 0 : homo = bs_env%n_occ(1)
1106 0 : band_edges_scf_guess%VBM = bs_env%eigenval_scf_Gamma(homo, 1)
1107 0 : band_edges_scf_guess%CBM = bs_env%eigenval_scf_Gamma(homo + 1, 1)
1108 : ELSE
1109 0 : homo_1 = bs_env%n_occ(1)
1110 0 : homo_2 = bs_env%n_occ(2)
1111 : band_edges_scf_guess%VBM = MAX(bs_env%eigenval_scf_Gamma(homo_1, 1), &
1112 0 : bs_env%eigenval_scf_Gamma(homo_2, 2))
1113 : band_edges_scf_guess%CBM = MIN(bs_env%eigenval_scf_Gamma(homo_1 + 1, 1), &
1114 0 : bs_env%eigenval_scf_Gamma(homo_2 + 1, 2))
1115 : END IF
1116 :
1117 : ! initialization
1118 0 : band_edges_scf%VBM = -1000.0_dp
1119 0 : band_edges_scf%CBM = 1000.0_dp
1120 0 : band_edges_scf%DBG = 1000.0_dp
1121 : END IF
1122 :
1123 106 : E_min = band_edges_scf_guess%VBM - 0.5_dp*energy_window_DOS
1124 106 : E_max = band_edges_scf_guess%CBM + 0.5_dp*energy_window_DOS
1125 :
1126 106 : IF (bs_env%do_gw) THEN
1127 104 : band_edges_G0W0 = bs_env%band_edges_G0W0
1128 104 : E_min_G0W0 = band_edges_G0W0%VBM - 0.5_dp*energy_window_DOS
1129 104 : E_max_G0W0 = band_edges_G0W0%CBM + 0.5_dp*energy_window_DOS
1130 104 : E_min = MIN(E_min, E_min_G0W0)
1131 104 : E_max = MAX(E_max, E_max_G0W0)
1132 : END IF
1133 :
1134 106 : E_total_window = E_max - E_min
1135 :
1136 106 : n_E = INT(E_total_window/energy_step_DOS)
1137 :
1138 106 : CALL get_qs_env(qs_env, nkind=nkind)
1139 :
1140 424 : ALLOCATE (proj_mo_on_kind(n_ao, nkind))
1141 106 : proj_mo_on_kind(:, :) = 0.0_dp
1142 :
1143 318 : ALLOCATE (eigenval(n_ao))
1144 318 : ALLOCATE (eigenval_spinor(2*n_ao))
1145 212 : ALLOCATE (eigenval_spinor_no_SOC(2*n_ao))
1146 212 : ALLOCATE (eigenval_spinor_G0W0(2*n_ao))
1147 :
1148 106 : IF (bs_env%do_dos_pdos) THEN
1149 :
1150 60 : ALLOCATE (DOS_scf(n_E))
1151 20 : DOS_scf(:) = 0.0_dp
1152 80 : ALLOCATE (PDOS_scf(n_E, nkind))
1153 20 : PDOS_scf(:, :) = 0.0_dp
1154 :
1155 20 : IF (bs_env%do_soc) THEN
1156 :
1157 32 : ALLOCATE (DOS_scf_SOC(n_E))
1158 16 : DOS_scf_SOC(:) = 0.0_dp
1159 48 : ALLOCATE (PDOS_scf_SOC(n_E, nkind))
1160 16 : PDOS_scf_SOC(:, :) = 0.0_dp
1161 :
1162 : END IF
1163 :
1164 20 : IF (bs_env%do_gw) THEN
1165 :
1166 40 : ALLOCATE (DOS_G0W0(n_E))
1167 20 : DOS_G0W0(:) = 0.0_dp
1168 60 : ALLOCATE (PDOS_G0W0(n_E, nkind))
1169 20 : PDOS_G0W0(:, :) = 0.0_dp
1170 :
1171 20 : IF (bs_env%do_soc) THEN
1172 :
1173 32 : ALLOCATE (DOS_G0W0_SOC(n_E))
1174 16 : DOS_G0W0_SOC(:) = 0.0_dp
1175 48 : ALLOCATE (PDOS_G0W0_SOC(n_E, nkind))
1176 16 : PDOS_G0W0_SOC(:, :) = 0.0_dp
1177 :
1178 : END IF
1179 : END IF
1180 : END IF
1181 :
1182 106 : CALL cp_cfm_create(cfm_mos_ikp(1), bs_env%fm_ks_Gamma(1)%matrix_struct)
1183 106 : CALL cp_cfm_create(cfm_mos_ikp(2), bs_env%fm_ks_Gamma(1)%matrix_struct)
1184 106 : CALL cp_cfm_create(cfm_work_ikp, bs_env%fm_ks_Gamma(1)%matrix_struct)
1185 106 : CALL cp_cfm_create(cfm_s_ikp_copy, bs_env%fm_ks_Gamma(1)%matrix_struct)
1186 :
1187 106 : IF (bs_env%do_soc) THEN
1188 :
1189 22 : CALL cp_cfm_create(cfm_mos_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1190 22 : CALL cp_cfm_create(cfm_work_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1191 22 : CALL cp_cfm_create(cfm_s_ikp_spinor_copy, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1192 22 : CALL cp_cfm_create(cfm_ks_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1193 22 : CALL cp_cfm_create(cfm_SOC_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1194 22 : CALL cp_cfm_create(cfm_s_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1195 22 : CALL cp_cfm_create(cfm_spinor_wf_ikp, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1196 :
1197 22 : homo_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin)
1198 :
1199 22 : band_edges_scf_SOC%VBM = -1000.0_dp
1200 22 : band_edges_scf_SOC%CBM = 1000.0_dp
1201 22 : band_edges_scf_SOC%DBG = 1000.0_dp
1202 :
1203 22 : IF (bs_env%do_gw) THEN
1204 22 : band_edges_G0W0_SOC%VBM = -1000.0_dp
1205 22 : band_edges_G0W0_SOC%CBM = 1000.0_dp
1206 22 : band_edges_G0W0_SOC%DBG = 1000.0_dp
1207 : END IF
1208 :
1209 22 : IF (bs_env%unit_nr > 0) THEN
1210 11 : WRITE (bs_env%unit_nr, '(A)') ''
1211 11 : IF (bs_env%soc_window_occ > 0.0_dp) THEN
1212 3 : WRITE (bs_env%unit_nr, '(T2,A,T71,F10.2)') 'SOC requested, SOC energy window occ (eV):', &
1213 6 : bs_env%soc_window_occ*evolt
1214 : ELSE
1215 8 : WRITE (bs_env%unit_nr, '(T2,A,T71,A10)') 'SOC requested, SOC energy window occ (eV):', &
1216 16 : ' no window'
1217 : END IF
1218 11 : IF (bs_env%soc_window_virt > 0.0_dp) THEN
1219 3 : WRITE (bs_env%unit_nr, '(T2,A,T71,F10.2)') 'SOC requested, SOC energy window virt (eV):', &
1220 6 : bs_env%soc_window_virt*evolt
1221 : ELSE
1222 8 : WRITE (bs_env%unit_nr, '(T2,A,T71,A10)') 'SOC requested, SOC energy window virt (eV):', &
1223 16 : ' no window'
1224 : END IF
1225 11 : IF (bs_env%soc_window_occ > 0.0_dp .OR. bs_env%soc_window_virt > 0.0_dp) THEN
1226 4 : WRITE (bs_env%unit_nr, '(T2,A,T71,F10.2)') 'SOC requested, SOC window smearing (eV):', &
1227 8 : bs_env%soc_window_smearing*evolt
1228 : END IF
1229 : END IF
1230 : END IF
1231 :
1232 106 : IF (bs_env%do_ldos) THEN
1233 2 : CPASSERT(bs_env%int_ldos_xyz == int_ldos_z)
1234 : END IF
1235 :
1236 106 : IF (bs_env%unit_nr > 0) THEN
1237 53 : WRITE (bs_env%unit_nr, '(A)') ''
1238 : END IF
1239 :
1240 106 : IF (bs_env%small_cell_full_kp_or_large_cell_Gamma == small_cell_full_kp) THEN
1241 18 : CALL cp_cfm_create(cfm_ks_ikp, bs_env%cfm_ks_kp(1, 1)%matrix_struct)
1242 18 : CALL cp_cfm_create(cfm_s_ikp, bs_env%cfm_ks_kp(1, 1)%matrix_struct)
1243 : END IF
1244 :
1245 554 : DO ikp = 1, bs_env%nkp_bs_and_DOS
1246 :
1247 448 : t1 = m_walltime()
1248 :
1249 914 : DO ispin = 1, bs_env%n_spin
1250 :
1251 580 : SELECT CASE (bs_env%small_cell_full_kp_or_large_cell_Gamma)
1252 : CASE (large_cell_Gamma, large_cell_Gamma_ri_rs, non_periodic_ri_rs)
1253 :
1254 : ! 1. get H^KS_µν(k_i) from H^KS_µν(k=0)
1255 : CALL cfm_ikp_from_fm_Gamma(cfm_ks_ikp, bs_env%fm_ks_Gamma(ispin), &
1256 114 : ikp, qs_env, bs_env%kpoints_DOS, "ORB")
1257 :
1258 : ! 2. get S_µν(k_i) from S_µν(k=0)
1259 : CALL cfm_ikp_from_fm_Gamma(cfm_s_ikp, bs_env%fm_s_Gamma, &
1260 114 : ikp, qs_env, bs_env%kpoints_DOS, "ORB")
1261 114 : CALL cp_cfm_to_cfm(cfm_s_ikp, cfm_s_ikp_copy)
1262 :
1263 : ! 3. Diagonalize (Roothaan-Hall): H_KS(k_i)*C(k_i) = S(k_i)*C(k_i)*ϵ(k_i)
1264 : CALL cp_cfm_geeig(cfm_ks_ikp, cfm_s_ikp_copy, cfm_mos_ikp(ispin), &
1265 114 : eigenval, cfm_work_ikp)
1266 :
1267 : CASE (small_cell_full_kp)
1268 :
1269 : ! 1. get H^KS_µν(k_i)
1270 352 : CALL cp_cfm_to_cfm(bs_env%cfm_ks_kp(ikp, ispin), cfm_ks_ikp)
1271 :
1272 : ! 2. get S_µν(k_i)
1273 352 : CALL cp_cfm_to_cfm(bs_env%cfm_s_kp(ikp), cfm_s_ikp)
1274 :
1275 : ! 3. get C_µn(k_i) and ϵ_n(k_i)
1276 352 : CALL cp_cfm_to_cfm(bs_env%cfm_mo_coeff_kp(ikp, ispin), cfm_mos_ikp(ispin))
1277 5072 : eigenval(:) = bs_env%eigenval_scf(:, ikp, ispin)
1278 :
1279 : END SELECT
1280 :
1281 : ! 4. Projection p_nk^A of MO ψ_nk(r) on atom type A (inspired by Mulliken charge)
1282 : ! p_nk^A = sum_µ^A,ν C*_µ^A,n(k) S_µ^A,ν(k) C_ν,n(k)
1283 466 : CALL compute_proj_mo_on_kind(proj_mo_on_kind, qs_env, cfm_mos_ikp(ispin), cfm_s_ikp)
1284 :
1285 : ! 5. DOS and PDOS
1286 466 : IF (bs_env%do_dos_pdos) THEN
1287 : CALL add_to_DOS_PDOS(DOS_scf, PDOS_scf, eigenval, ikp, bs_env, n_E, E_min, &
1288 234 : proj_mo_on_kind)
1289 :
1290 234 : IF (bs_env%do_gw) THEN
1291 : CALL add_to_DOS_PDOS(DOS_G0W0, PDOS_G0W0, bs_env%eigenval_G0W0(:, ikp, ispin), &
1292 234 : ikp, bs_env, n_E, E_min, proj_mo_on_kind)
1293 : END IF
1294 : END IF
1295 :
1296 466 : IF (bs_env%do_ldos) THEN
1297 : CALL add_to_LDOS_2d(LDOS_scf_2d, qs_env, ikp, bs_env, cfm_mos_ikp(ispin), &
1298 2 : eigenval(:), band_edges_scf_guess)
1299 :
1300 2 : IF (bs_env%do_gw) THEN
1301 : CALL add_to_LDOS_2d(LDOS_G0W0_2d, qs_env, ikp, bs_env, cfm_mos_ikp(ispin), &
1302 2 : bs_env%eigenval_G0W0(:, ikp, 1), band_edges_G0W0)
1303 : END IF
1304 :
1305 : END IF
1306 :
1307 466 : homo = bs_env%n_occ(ispin)
1308 :
1309 466 : band_edges_scf%VBM = MAX(band_edges_scf%VBM, eigenval(homo))
1310 466 : band_edges_scf%CBM = MIN(band_edges_scf%CBM, eigenval(homo + 1))
1311 914 : band_edges_scf%DBG = MIN(band_edges_scf%DBG, eigenval(homo + 1) - eigenval(homo))
1312 :
1313 : END DO ! spin
1314 :
1315 : ! now the same with spin-orbit coupling
1316 448 : IF (bs_env%do_soc) THEN
1317 :
1318 : ! only print eigenvalues of DOS k-points in case no bandstructure path has been given
1319 328 : print_DOS_kpoints = (bs_env%nkp_only_bs <= 0)
1320 : ! in kpoints_DOS, the last nkp_only_bs are bandstructure k-points
1321 328 : is_bandstruc_kpoint = (ikp > bs_env%nkp_only_DOS)
1322 328 : print_ikp = print_DOS_kpoints .OR. is_bandstruc_kpoint
1323 :
1324 328 : IF (print_DOS_kpoints) THEN
1325 234 : nkp = bs_env%nkp_only_DOS
1326 234 : ikp_for_file = ikp
1327 : ELSE
1328 94 : nkp = bs_env%nkp_only_bs
1329 94 : ikp_for_file = ikp - bs_env%nkp_only_DOS
1330 : END IF
1331 :
1332 : ! compute DFT+SOC eigenvalues; based on these, compute band edges, DOS and LDOS
1333 : CALL SOC_ev(bs_env, qs_env, ikp, bs_env%eigenval_scf, &
1334 : E_min, cfm_mos_ikp, DOS_scf_SOC, PDOS_scf_SOC, &
1335 328 : band_edges_scf_SOC, eigenval_spinor, cfm_spinor_wf_ikp)
1336 :
1337 328 : IF (.NOT. bs_env%do_gw .AND. print_ikp) THEN
1338 0 : CALL write_SOC_eigenvalues(eigenval_spinor, ikp_for_file, ikp, bs_env)
1339 : END IF
1340 :
1341 328 : IF (bs_env%do_ldos) THEN
1342 : CALL add_to_LDOS_2d(LDOS_scf_2d_SOC, qs_env, ikp, bs_env, cfm_spinor_wf_ikp, &
1343 2 : eigenval_spinor, band_edges_scf_guess, .TRUE., cfm_work_ikp)
1344 : END IF
1345 :
1346 328 : IF (bs_env%do_gw) THEN
1347 :
1348 : ! compute G0W0+SOC eigenvalues; based on these, compute band edges, DOS and LDOS
1349 : CALL SOC_ev(bs_env, qs_env, ikp, bs_env%eigenval_G0W0, &
1350 : E_min, cfm_mos_ikp, DOS_G0W0_SOC, PDOS_G0W0_SOC, &
1351 328 : band_edges_G0W0_SOC, eigenval_spinor_G0W0, cfm_spinor_wf_ikp)
1352 :
1353 328 : IF (print_ikp) THEN
1354 : ! write SCF+SOC and G0W0+SOC eigenvalues to file
1355 : ! SCF_and_G0W0_band_structure_for_kpoint_<ikp>_+_SOC
1356 : CALL write_SOC_eigenvalues(eigenval_spinor, ikp_for_file, ikp, bs_env, &
1357 296 : eigenval_spinor_G0W0)
1358 : END IF
1359 :
1360 : END IF ! do_gw
1361 :
1362 : END IF ! do_soc
1363 :
1364 554 : IF (bs_env%unit_nr > 0 .AND. m_walltime() - t1 > 20.0_dp) THEN
1365 : WRITE (bs_env%unit_nr, '(T2,A,T43,I5,A,I3,A,F7.1,A)') &
1366 0 : 'Compute DOS, LDOS for k-point ', ikp, ' /', bs_env%nkp_bs_and_DOS, &
1367 0 : ', Execution time', m_walltime() - t1, ' s'
1368 : END IF
1369 :
1370 : END DO ! ikp_DOS
1371 :
1372 106 : band_edges_scf%IDBG = band_edges_scf%CBM - band_edges_scf%VBM
1373 106 : IF (bs_env%do_soc) THEN
1374 22 : band_edges_scf_SOC%IDBG = band_edges_scf_SOC%CBM - band_edges_scf_SOC%VBM
1375 22 : IF (bs_env%do_gw) THEN
1376 22 : band_edges_G0W0_SOC%IDBG = band_edges_G0W0_SOC%CBM - band_edges_G0W0_SOC%VBM
1377 : END IF
1378 : END IF
1379 :
1380 106 : CALL write_band_edges(band_edges_scf, "SCF", bs_env)
1381 106 : IF (bs_env%do_dos_pdos) THEN
1382 20 : CALL write_dos_pdos(DOS_scf, PDOS_scf, bs_env, qs_env, "SCF", E_min, band_edges_scf%VBM)
1383 : END IF
1384 106 : IF (bs_env%do_ldos) THEN
1385 2 : CALL print_LDOS_main(LDOS_scf_2d, bs_env, band_edges_scf, "SCF")
1386 : END IF
1387 :
1388 106 : IF (bs_env%do_soc) THEN
1389 22 : CALL write_band_edges(band_edges_scf_SOC, "SCF+SOC", bs_env)
1390 22 : IF (bs_env%do_dos_pdos) THEN
1391 : CALL write_dos_pdos(DOS_scf_SOC, PDOS_scf_SOC, bs_env, qs_env, "SCF_SOC", &
1392 16 : E_min, band_edges_scf_SOC%VBM)
1393 : END IF
1394 22 : IF (bs_env%do_ldos) THEN
1395 : ! argument band_edges_scf is actually correct because the non-SOC band edges
1396 : ! have been used as reference in add_to_LDOS_2d
1397 : CALL print_LDOS_main(LDOS_scf_2d_SOC, bs_env, band_edges_scf, &
1398 2 : "SCF_SOC")
1399 : END IF
1400 : END IF
1401 :
1402 106 : IF (bs_env%do_gw) THEN
1403 104 : CALL write_band_edges(band_edges_G0W0, "G0W0", bs_env)
1404 104 : CALL write_band_edges(bs_env%band_edges_HF, "Hartree-Fock with SCF orbitals", bs_env)
1405 104 : IF (bs_env%do_dos_pdos) THEN
1406 : CALL write_dos_pdos(DOS_G0W0, PDOS_G0W0, bs_env, qs_env, "G0W0", E_min, &
1407 20 : band_edges_G0W0%VBM)
1408 : END IF
1409 104 : IF (bs_env%do_ldos) THEN
1410 2 : CALL print_LDOS_main(LDOS_G0W0_2d, bs_env, band_edges_G0W0, "G0W0")
1411 : END IF
1412 : END IF
1413 :
1414 106 : IF (bs_env%do_soc .AND. bs_env%do_gw) THEN
1415 22 : CALL write_band_edges(band_edges_G0W0_SOC, "G0W0+SOC", bs_env)
1416 22 : IF (bs_env%do_dos_pdos) THEN
1417 : CALL write_dos_pdos(DOS_G0W0_SOC, PDOS_G0W0_SOC, bs_env, qs_env, "G0W0_SOC", E_min, &
1418 16 : band_edges_G0W0_SOC%VBM)
1419 : END IF
1420 : END IF
1421 :
1422 106 : CALL cp_cfm_release(cfm_s_ikp)
1423 106 : CALL cp_cfm_release(cfm_ks_ikp)
1424 106 : CALL cp_cfm_release(cfm_mos_ikp(1))
1425 106 : CALL cp_cfm_release(cfm_mos_ikp(2))
1426 106 : CALL cp_cfm_release(cfm_work_ikp)
1427 106 : CALL cp_cfm_release(cfm_s_ikp_copy)
1428 :
1429 106 : CALL cp_cfm_release(cfm_s_ikp_spinor)
1430 106 : CALL cp_cfm_release(cfm_ks_ikp_spinor)
1431 106 : CALL cp_cfm_release(cfm_SOC_ikp_spinor)
1432 106 : CALL cp_cfm_release(cfm_mos_ikp_spinor)
1433 106 : CALL cp_cfm_release(cfm_work_ikp_spinor)
1434 106 : CALL cp_cfm_release(cfm_s_ikp_spinor_copy)
1435 106 : CALL cp_cfm_release(cfm_spinor_wf_ikp)
1436 :
1437 106 : CALL timestop(handle)
1438 :
1439 424 : END SUBROUTINE eval_bandstructure_properties
1440 :
1441 : ! **************************************************************************************************
1442 : !> \brief ...
1443 : !> \param LDOS_2d ...
1444 : !> \param bs_env ...
1445 : !> \param band_edges ...
1446 : !> \param scf_gw_soc ...
1447 : ! **************************************************************************************************
1448 6 : SUBROUTINE print_LDOS_main(LDOS_2d, bs_env, band_edges, scf_gw_soc)
1449 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: LDOS_2d
1450 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1451 : TYPE(band_edges_type) :: band_edges
1452 : CHARACTER(LEN=*) :: scf_gw_soc
1453 :
1454 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_LDOS_main'
1455 :
1456 : INTEGER :: handle, i_x, i_x_bin, i_x_end, i_x_end_bin, i_x_end_glob, i_x_start, &
1457 : i_x_start_bin, i_x_start_glob, i_y, i_y_bin, i_y_end, i_y_end_bin, i_y_end_glob, &
1458 : i_y_start, i_y_start_bin, i_y_start_glob, n_E
1459 6 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: n_sum_for_bins
1460 : INTEGER, DIMENSION(2) :: bin_mesh
1461 : LOGICAL :: do_xy_bins
1462 : REAL(KIND=dp) :: E_min, energy_step, energy_window
1463 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: LDOS_2d_bins
1464 :
1465 6 : CALL timeset(routineN, handle)
1466 :
1467 6 : n_E = SIZE(LDOS_2d, 3)
1468 :
1469 6 : energy_window = bs_env%energy_window_DOS
1470 6 : energy_step = bs_env%energy_step_DOS
1471 6 : E_min = band_edges%VBM - 0.5_dp*energy_window
1472 :
1473 18 : bin_mesh(1:2) = bs_env%bin_mesh(1:2)
1474 6 : do_xy_bins = (bin_mesh(1) > 0 .AND. bin_mesh(2) > 0)
1475 :
1476 6 : i_x_start = LBOUND(LDOS_2d, 1)
1477 6 : i_x_end = UBOUND(LDOS_2d, 1)
1478 6 : i_y_start = LBOUND(LDOS_2d, 2)
1479 6 : i_y_end = UBOUND(LDOS_2d, 2)
1480 :
1481 6 : IF (do_xy_bins) THEN
1482 6 : i_x_start_bin = 1
1483 6 : i_x_end_bin = bin_mesh(1)
1484 6 : i_y_start_bin = 1
1485 6 : i_y_end_bin = bin_mesh(2)
1486 : ELSE
1487 : i_x_start_bin = i_x_start
1488 : i_x_end_bin = i_x_end
1489 : i_y_start_bin = i_y_start
1490 : i_y_end_bin = i_y_end
1491 : END IF
1492 :
1493 30 : ALLOCATE (LDOS_2d_bins(i_x_start_bin:i_x_end_bin, i_y_start_bin:i_y_end_bin, n_E))
1494 6 : LDOS_2d_bins(:, :, :) = 0.0_dp
1495 :
1496 6 : IF (do_xy_bins) THEN
1497 :
1498 6 : i_x_start_glob = i_x_start
1499 6 : i_x_end_glob = i_x_end
1500 6 : i_y_start_glob = i_y_start
1501 6 : i_y_end_glob = i_y_end
1502 :
1503 6 : CALL bs_env%para_env%min(i_x_start_glob)
1504 6 : CALL bs_env%para_env%max(i_x_end_glob)
1505 6 : CALL bs_env%para_env%min(i_y_start_glob)
1506 6 : CALL bs_env%para_env%max(i_y_end_glob)
1507 :
1508 24 : ALLOCATE (n_sum_for_bins(bin_mesh(1), bin_mesh(2)), SOURCE=0)
1509 :
1510 : ! transform interval [i_x_start, i_x_end] to [1, bin_mesh(1)] (and same for y)
1511 390 : DO i_y = i_y_start, i_y_end
1512 4230 : DO i_x = i_x_start, i_x_end
1513 3840 : i_x_bin = bin_mesh(1)*(i_x - i_x_start_glob)/(i_x_end_glob - i_x_start_glob + 1) + 1
1514 3840 : i_y_bin = bin_mesh(2)*(i_y - i_y_start_glob)/(i_y_end_glob - i_y_start_glob + 1) + 1
1515 : LDOS_2d_bins(i_x_bin, i_y_bin, :) = LDOS_2d_bins(i_x_bin, i_y_bin, :) + &
1516 1073920 : LDOS_2d(i_x, i_y, :)
1517 4224 : n_sum_for_bins(i_x_bin, i_y_bin) = n_sum_for_bins(i_x_bin, i_y_bin) + 1
1518 : END DO
1519 : END DO
1520 :
1521 6 : CALL bs_env%para_env%sum(LDOS_2d_bins)
1522 6 : CALL bs_env%para_env%sum(n_sum_for_bins)
1523 :
1524 : ! divide by number of terms in the sum so we have the average LDOS(x,y,E)
1525 30 : DO i_y_bin = 1, bin_mesh(2)
1526 126 : DO i_x_bin = 1, bin_mesh(1)
1527 : LDOS_2d_bins(i_x_bin, i_y_bin, :) = LDOS_2d_bins(i_x_bin, i_y_bin, :)/ &
1528 26872 : REAL(n_sum_for_bins(i_x_bin, i_y_bin), KIND=dp)
1529 : END DO
1530 : END DO
1531 :
1532 : ELSE
1533 :
1534 0 : LDOS_2d_bins(:, :, :) = LDOS_2d(:, :, :)
1535 :
1536 : END IF
1537 :
1538 6 : IF (bin_mesh(1)*bin_mesh(2) < bs_env%n_bins_max_for_printing) THEN
1539 6 : CALL print_LDOS_2d_bins(LDOS_2d_bins, bs_env, E_min, scf_gw_soc)
1540 : ELSE
1541 0 : CPWARN("The number of bins for the LDOS is too large. Decrease BIN_MESH.")
1542 : END IF
1543 :
1544 6 : CALL timestop(handle)
1545 :
1546 12 : END SUBROUTINE print_LDOS_main
1547 :
1548 : ! **************************************************************************************************
1549 : !> \brief ...
1550 : !> \param LDOS_2d_bins ...
1551 : !> \param bs_env ...
1552 : !> \param E_min ...
1553 : !> \param scf_gw_soc ...
1554 : ! **************************************************************************************************
1555 6 : SUBROUTINE print_LDOS_2d_bins(LDOS_2d_bins, bs_env, E_min, scf_gw_soc)
1556 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: LDOS_2d_bins
1557 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1558 : REAL(KIND=dp) :: E_min
1559 : CHARACTER(LEN=*) :: scf_gw_soc
1560 :
1561 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_LDOS_2d_bins'
1562 :
1563 : CHARACTER(LEN=18) :: print_format
1564 : CHARACTER(LEN=4) :: print_format_1, print_format_2
1565 : CHARACTER(len=default_string_length) :: fname
1566 : INTEGER :: handle, i_E, i_x, i_x_end, i_x_start, &
1567 : i_y, i_y_end, i_y_start, iunit, n_E, &
1568 : n_x, n_y
1569 : REAL(KIND=dp) :: energy
1570 : REAL(KIND=dp), DIMENSION(3) :: coord, idx
1571 :
1572 6 : CALL timeset(routineN, handle)
1573 :
1574 6 : i_x_start = LBOUND(LDOS_2d_bins, 1)
1575 6 : i_x_end = UBOUND(LDOS_2d_bins, 1)
1576 6 : i_y_start = LBOUND(LDOS_2d_bins, 2)
1577 6 : i_y_end = UBOUND(LDOS_2d_bins, 2)
1578 6 : n_E = SIZE(LDOS_2d_bins, 3)
1579 :
1580 6 : n_x = i_x_end - i_x_start + 1
1581 6 : n_y = i_y_end - i_y_start + 1
1582 :
1583 6 : IF (bs_env%para_env%is_source()) THEN
1584 :
1585 15 : DO i_y = i_y_start, i_y_end
1586 63 : DO i_x = i_x_start, i_x_end
1587 :
1588 48 : idx(1) = (REAL(i_x, KIND=dp) - 0.5_dp)/REAL(n_x, KIND=dp)
1589 48 : idx(2) = (REAL(i_y, KIND=dp) - 0.5_dp)/REAL(n_y, KIND=dp)
1590 48 : idx(3) = 0.0_dp
1591 624 : coord(1:3) = MATMUL(bs_env%hmat, idx)
1592 :
1593 48 : CALL get_print_format(coord(1), print_format_1)
1594 48 : CALL get_print_format(coord(2), print_format_2)
1595 :
1596 48 : print_format = "(3A,"//print_format_1//",A,"//print_format_2//",A)"
1597 :
1598 48 : WRITE (fname, print_format) "LDOS_", scf_gw_soc, &
1599 96 : "_at_x_", coord(1)*angstrom, '_A_and_y_', coord(2)*angstrom, '_A'
1600 :
1601 : CALL open_file(TRIM(fname), unit_number=iunit, file_status="REPLACE", &
1602 48 : file_action="WRITE")
1603 :
1604 48 : WRITE (iunit, "(2A)") " Energy E (eV) average LDOS(x,y,E) (1/(eV*Å^2), ", &
1605 96 : "integrated over z, averaged inside bin)"
1606 :
1607 13424 : DO i_E = 1, n_E
1608 13376 : energy = E_min + i_E*bs_env%energy_step_DOS
1609 13376 : WRITE (iunit, "(2F17.3)") energy*evolt, &
1610 : LDOS_2d_bins(i_x, i_y, i_E)* &
1611 26800 : bs_env%unit_ldos_int_z_inv_Ang2_eV
1612 : END DO
1613 :
1614 60 : CALL close_file(iunit)
1615 :
1616 : END DO
1617 : END DO
1618 :
1619 : END IF
1620 :
1621 6 : CALL timestop(handle)
1622 :
1623 6 : END SUBROUTINE print_LDOS_2d_bins
1624 :
1625 : ! **************************************************************************************************
1626 : !> \brief ...
1627 : !> \param coord ...
1628 : !> \param print_format ...
1629 : ! **************************************************************************************************
1630 96 : SUBROUTINE get_print_format(coord, print_format)
1631 : REAL(KIND=dp) :: coord
1632 : CHARACTER(LEN=4) :: print_format
1633 :
1634 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_print_format'
1635 :
1636 : INTEGER :: handle
1637 :
1638 96 : CALL timeset(routineN, handle)
1639 :
1640 96 : IF (coord < -10000/angstrom) THEN
1641 0 : print_format = "F9.2"
1642 96 : ELSE IF (coord < -1000/angstrom) THEN
1643 0 : print_format = "F8.2"
1644 96 : ELSE IF (coord < -100/angstrom) THEN
1645 0 : print_format = "F7.2"
1646 96 : ELSE IF (coord < -10/angstrom) THEN
1647 0 : print_format = "F6.2"
1648 96 : ELSE IF (coord < -1/angstrom) THEN
1649 0 : print_format = "F5.2"
1650 96 : ELSE IF (coord < 10/angstrom) THEN
1651 96 : print_format = "F4.2"
1652 0 : ELSE IF (coord < 100/angstrom) THEN
1653 0 : print_format = "F5.2"
1654 0 : ELSE IF (coord < 1000/angstrom) THEN
1655 0 : print_format = "F6.2"
1656 0 : ELSE IF (coord < 10000/angstrom) THEN
1657 0 : print_format = "F7.2"
1658 : ELSE
1659 0 : print_format = "F8.2"
1660 : END IF
1661 :
1662 96 : CALL timestop(handle)
1663 :
1664 96 : END SUBROUTINE get_print_format
1665 :
1666 : ! **************************************************************************************************
1667 : !> \brief ...
1668 : !> \param bs_env ...
1669 : !> \param qs_env ...
1670 : !> \param ikp ...
1671 : !> \param eigenval_no_SOC ...
1672 : !> \param E_min ...
1673 : !> \param cfm_mos_ikp ...
1674 : !> \param DOS ...
1675 : !> \param PDOS ...
1676 : !> \param band_edges ...
1677 : !> \param eigenval_spinor ...
1678 : !> \param cfm_spinor_wf_ikp ...
1679 : ! **************************************************************************************************
1680 656 : SUBROUTINE SOC_ev(bs_env, qs_env, ikp, eigenval_no_SOC, E_min, cfm_mos_ikp, &
1681 : DOS, PDOS, band_edges, eigenval_spinor, cfm_spinor_wf_ikp)
1682 :
1683 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1684 : TYPE(qs_environment_type), POINTER :: qs_env
1685 : INTEGER :: ikp
1686 : REAL(KIND=dp), DIMENSION(:, :, :) :: eigenval_no_SOC
1687 : REAL(KIND=dp) :: E_min
1688 : TYPE(cp_cfm_type), DIMENSION(2) :: cfm_mos_ikp
1689 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: DOS
1690 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: PDOS
1691 : TYPE(band_edges_type) :: band_edges
1692 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_spinor
1693 : TYPE(cp_cfm_type) :: cfm_spinor_wf_ikp
1694 :
1695 : CHARACTER(LEN=*), PARAMETER :: routineN = 'SOC_ev'
1696 :
1697 : INTEGER :: handle, homo_spinor, n_ao, n_E, nkind
1698 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_spinor_no_SOC
1699 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: proj_mo_on_kind_spinor
1700 : TYPE(cp_cfm_type) :: cfm_eigenvec_ikp_spinor, &
1701 : cfm_ks_ikp_spinor, cfm_mos_ikp_spinor, &
1702 : cfm_SOC_ikp_spinor, cfm_work_ikp_spinor
1703 :
1704 : !TYPE(band_edges_type) :: band_edges_no_SOC
1705 :
1706 656 : CALL timeset(routineN, handle)
1707 :
1708 656 : n_ao = bs_env%n_ao
1709 656 : homo_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin)
1710 656 : CALL get_qs_env(qs_env, nkind=nkind)
1711 :
1712 656 : CALL cp_cfm_create(cfm_ks_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1713 656 : CALL cp_cfm_create(cfm_SOC_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1714 656 : CALL cp_cfm_create(cfm_mos_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1715 656 : CALL cp_cfm_create(cfm_work_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1716 656 : CALL cp_cfm_create(cfm_eigenvec_ikp_spinor, bs_env%cfm_SOC_spinor_ao(1)%matrix_struct)
1717 :
1718 1968 : ALLOCATE (eigenval_spinor_no_SOC(2*n_ao))
1719 2624 : ALLOCATE (proj_mo_on_kind_spinor(2*n_ao, nkind))
1720 : ! PDOS not yet implemented -> projection is just zero -> PDOS is zero
1721 656 : proj_mo_on_kind_spinor(:, :) = 0.0_dp
1722 :
1723 : ! 1. get V^SOC_µν,σσ'(k_i)
1724 676 : SELECT CASE (bs_env%small_cell_full_kp_or_large_cell_Gamma)
1725 : CASE (large_cell_Gamma, large_cell_Gamma_ri_rs, non_periodic_ri_rs)
1726 :
1727 : ! 1. get V^SOC_µν,σσ'(k_i) from V^SOC_µν,σσ'(k=0)
1728 : CALL cfm_ikp_from_cfm_spinor_Gamma(cfm_SOC_ikp_spinor, &
1729 : bs_env%cfm_SOC_spinor_ao(1), &
1730 : bs_env%fm_s_Gamma%matrix_struct, &
1731 20 : ikp, qs_env, bs_env%kpoints_DOS, "ORB")
1732 :
1733 : CASE (small_cell_full_kp)
1734 :
1735 : ! 1. V^SOC_µν,σσ'(k_i) already there
1736 656 : CALL cp_cfm_to_cfm(bs_env%cfm_SOC_spinor_ao(ikp), cfm_SOC_ikp_spinor)
1737 :
1738 : END SELECT
1739 :
1740 : ! 2. V^SOC_nn',σσ'(k_i) = sum_µν C^*_µn,σ(k_i) V^SOC_µν,σσ'(k_i) C_νn'(k_i),
1741 : ! C_µn,σ(k_i): MO coefficiencts from diagonalizing KS-matrix h^KS_nn',σσ'(k_i)
1742 :
1743 : ! 2.1 build matrix C_µn,σ(k_i)
1744 656 : CALL cp_cfm_set_all(cfm_mos_ikp_spinor, z_zero)
1745 656 : CALL add_cfm_submat(cfm_mos_ikp_spinor, cfm_mos_ikp(1), 1, 1)
1746 656 : CALL add_cfm_submat(cfm_mos_ikp_spinor, cfm_mos_ikp(bs_env%n_spin), n_ao + 1, n_ao + 1)
1747 :
1748 : ! 2.2 work_nν,σσ' = sum_µ C^*_µn,σ(k_i) V^SOC_µν,σσ'(k_i)
1749 : CALL parallel_gemm('C', 'N', 2*n_ao, 2*n_ao, 2*n_ao, z_one, &
1750 : cfm_mos_ikp_spinor, cfm_SOC_ikp_spinor, &
1751 656 : z_zero, cfm_work_ikp_spinor)
1752 :
1753 : ! 2.3 V^SOC_nn',σσ'(k_i) = sum_ν work_nν,σσ' C_νn'(k_i)
1754 : CALL parallel_gemm('N', 'N', 2*n_ao, 2*n_ao, 2*n_ao, z_one, &
1755 : cfm_work_ikp_spinor, cfm_mos_ikp_spinor, &
1756 656 : z_zero, cfm_ks_ikp_spinor)
1757 :
1758 : ! 3. remove SOC outside of energy window (otherwise, numerical problems arise
1759 : ! because energetically low semicore states and energetically very high
1760 : ! unbound states couple to the states around the Fermi level)
1761 8800 : eigenval_spinor_no_SOC(1:n_ao) = eigenval_no_SOC(1:n_ao, ikp, 1)
1762 8800 : eigenval_spinor_no_SOC(n_ao + 1:) = eigenval_no_SOC(1:n_ao, ikp, bs_env%n_spin)
1763 656 : IF (bs_env%soc_window_occ > 0.0_dp .OR. bs_env%soc_window_virt > 0.0_dp) THEN
1764 : CALL remove_soc_outside_energy_window_mo(cfm_ks_ikp_spinor, &
1765 : bs_env%soc_window_virt, &
1766 : bs_env%soc_window_smearing, &
1767 : eigenval_spinor_no_SOC, &
1768 256 : bs_env%e_fermi(1))
1769 :
1770 : END IF
1771 :
1772 : ! 4. h^G0W0+SOC_nn',σσ'(k_i) = ε_nσ^G0W0(k_i) δ_nn' δ_σσ' + V^SOC_nn',σσ'(k_i)
1773 656 : CALL cfm_add_on_diag(cfm_ks_ikp_spinor, eigenval_spinor_no_SOC)
1774 :
1775 : ! 5. diagonalize h^G0W0+SOC_nn',σσ'(k_i) to get eigenvalues
1776 656 : CALL cp_cfm_heevd(cfm_ks_ikp_spinor, cfm_eigenvec_ikp_spinor, eigenval_spinor)
1777 :
1778 : ! 6. DOS from spinors, no PDOS
1779 656 : IF (bs_env%do_dos_pdos) THEN
1780 452 : n_E = SIZE(DOS)
1781 : CALL add_to_DOS_PDOS(DOS, PDOS, eigenval_spinor, &
1782 452 : ikp, bs_env, n_E, E_min, proj_mo_on_kind_spinor)
1783 : END IF
1784 :
1785 : ! 7. valence band max. (VBM), conduction band min. (CBM) and direct bandgap (DBG)
1786 656 : band_edges%VBM = MAX(band_edges%VBM, eigenval_spinor(homo_spinor))
1787 656 : band_edges%CBM = MIN(band_edges%CBM, eigenval_spinor(homo_spinor + 1))
1788 : band_edges%DBG = MIN(band_edges%DBG, eigenval_spinor(homo_spinor + 1) &
1789 656 : - eigenval_spinor(homo_spinor))
1790 :
1791 : ! 8. spinor wavefunctions:
1792 : CALL parallel_gemm('N', 'N', 2*n_ao, 2*n_ao, 2*n_ao, z_one, &
1793 : cfm_mos_ikp_spinor, cfm_eigenvec_ikp_spinor, &
1794 656 : z_zero, cfm_spinor_wf_ikp)
1795 :
1796 656 : CALL cp_cfm_release(cfm_ks_ikp_spinor)
1797 656 : CALL cp_cfm_release(cfm_SOC_ikp_spinor)
1798 656 : CALL cp_cfm_release(cfm_work_ikp_spinor)
1799 656 : CALL cp_cfm_release(cfm_eigenvec_ikp_spinor)
1800 656 : CALL cp_cfm_release(cfm_mos_ikp_spinor)
1801 :
1802 656 : CALL timestop(handle)
1803 :
1804 1968 : END SUBROUTINE SOC_ev
1805 :
1806 : ! **************************************************************************************************
1807 : !> \brief ...
1808 : !> \param DOS ...
1809 : !> \param PDOS ...
1810 : !> \param eigenval ...
1811 : !> \param ikp ...
1812 : !> \param bs_env ...
1813 : !> \param n_E ...
1814 : !> \param E_min ...
1815 : !> \param proj_mo_on_kind ...
1816 : ! **************************************************************************************************
1817 920 : SUBROUTINE add_to_DOS_PDOS(DOS, PDOS, eigenval, ikp, bs_env, n_E, E_min, proj_mo_on_kind)
1818 :
1819 : REAL(KIND=dp), DIMENSION(:) :: DOS
1820 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: PDOS
1821 : REAL(KIND=dp), DIMENSION(:) :: eigenval
1822 : INTEGER :: ikp
1823 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1824 : INTEGER :: n_E
1825 : REAL(KIND=dp) :: E_min
1826 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: proj_mo_on_kind
1827 :
1828 : CHARACTER(LEN=*), PARAMETER :: routineN = 'add_to_DOS_PDOS'
1829 :
1830 : INTEGER :: handle, i_E, i_kind, i_mo, n_mo, nkind
1831 : REAL(KIND=dp) :: broadening, energy, energy_step_DOS, wkp
1832 :
1833 920 : CALL timeset(routineN, handle)
1834 :
1835 920 : energy_step_DOS = bs_env%energy_step_DOS
1836 920 : broadening = bs_env%broadening_DOS
1837 :
1838 920 : n_mo = SIZE(eigenval)
1839 920 : nkind = SIZE(proj_mo_on_kind, 2)
1840 :
1841 : ! normalize to closed-shell / open-shell
1842 920 : wkp = bs_env%kpoints_DOS%wkp(ikp)*bs_env%spin_degeneracy
1843 2072664 : DO i_E = 1, n_E
1844 2071744 : energy = E_min + i_E*energy_step_DOS
1845 45624244 : DO i_mo = 1, n_mo
1846 : ! DOS
1847 43551580 : DOS(i_E) = DOS(i_E) + wkp*Gaussian(energy - eigenval(i_mo), broadening)
1848 :
1849 : ! PDOS
1850 132726484 : DO i_kind = 1, nkind
1851 130654740 : IF (proj_mo_on_kind(i_mo, i_kind) > 0.0_dp) THEN
1852 : PDOS(i_E, i_kind) = PDOS(i_E, i_kind) + &
1853 : proj_mo_on_kind(i_mo, i_kind)*wkp* &
1854 27026364 : Gaussian(energy - eigenval(i_mo), broadening)
1855 : END IF
1856 : END DO
1857 : END DO
1858 : END DO
1859 :
1860 920 : CALL timestop(handle)
1861 :
1862 920 : END SUBROUTINE add_to_DOS_PDOS
1863 :
1864 : ! **************************************************************************************************
1865 : !> \brief ...
1866 : !> \param LDOS_2d ...
1867 : !> \param qs_env ...
1868 : !> \param ikp ...
1869 : !> \param bs_env ...
1870 : !> \param cfm_mos_ikp ...
1871 : !> \param eigenval ...
1872 : !> \param band_edges ...
1873 : !> \param do_spinor ...
1874 : !> \param cfm_non_spinor ...
1875 : ! **************************************************************************************************
1876 6 : SUBROUTINE add_to_LDOS_2d(LDOS_2d, qs_env, ikp, bs_env, cfm_mos_ikp, eigenval, &
1877 : band_edges, do_spinor, cfm_non_spinor)
1878 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: LDOS_2d
1879 : TYPE(qs_environment_type), POINTER :: qs_env
1880 : INTEGER :: ikp
1881 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1882 : TYPE(cp_cfm_type) :: cfm_mos_ikp
1883 : REAL(KIND=dp), DIMENSION(:) :: eigenval
1884 : TYPE(band_edges_type) :: band_edges
1885 : LOGICAL, OPTIONAL :: do_spinor
1886 : TYPE(cp_cfm_type), OPTIONAL :: cfm_non_spinor
1887 :
1888 : CHARACTER(LEN=*), PARAMETER :: routineN = 'add_to_LDOS_2d'
1889 :
1890 : INTEGER :: handle, i_E, i_x_end, i_x_start, i_y_end, i_y_start, i_z, i_z_end, i_z_start, &
1891 : j_col, j_mo, n_E, n_mo, n_z, ncol_local, nimages, z_end_global, z_start_global
1892 6 : INTEGER, DIMENSION(:), POINTER :: col_indices
1893 : LOGICAL :: is_any_weight_non_zero, my_do_spinor
1894 : REAL(KIND=dp) :: broadening, E_max, E_min, &
1895 : E_total_window, energy, energy_step, &
1896 : energy_window, spin_degeneracy, weight
1897 : TYPE(cp_cfm_type) :: cfm_weighted_dm_ikp, cfm_work
1898 : TYPE(cp_fm_type) :: fm_non_spinor, fm_weighted_dm_MIC
1899 6 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: weighted_dm_MIC
1900 : TYPE(dft_control_type), POINTER :: dft_control
1901 : TYPE(pw_c1d_gs_type) :: rho_g
1902 : TYPE(pw_env_type), POINTER :: pw_env
1903 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1904 : TYPE(pw_r3d_rs_type) :: LDOS_3d
1905 : TYPE(qs_ks_env_type), POINTER :: ks_env
1906 :
1907 6 : CALL timeset(routineN, handle)
1908 :
1909 6 : my_do_spinor = .FALSE.
1910 6 : IF (PRESENT(do_spinor)) my_do_spinor = do_spinor
1911 :
1912 6 : CALL get_qs_env(qs_env, ks_env=ks_env, pw_env=pw_env, dft_control=dft_control)
1913 :
1914 : ! previously, dft_control%nimages set to # neighbor cells, revert for Γ-only KS matrix
1915 6 : nimages = dft_control%nimages
1916 6 : dft_control%nimages = bs_env%nimages_scf
1917 :
1918 6 : energy_window = bs_env%energy_window_DOS
1919 6 : energy_step = bs_env%energy_step_DOS
1920 6 : broadening = bs_env%broadening_DOS
1921 :
1922 6 : E_min = band_edges%VBM - 0.5_dp*energy_window
1923 6 : E_max = band_edges%CBM + 0.5_dp*energy_window
1924 6 : E_total_window = E_max - E_min
1925 :
1926 6 : n_E = INT(E_total_window/energy_step)
1927 :
1928 6 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1929 :
1930 6 : CALL auxbas_pw_pool%create_pw(LDOS_3d)
1931 6 : CALL auxbas_pw_pool%create_pw(rho_g)
1932 :
1933 6 : i_x_start = LBOUND(LDOS_3d%array, 1)
1934 6 : i_x_end = UBOUND(LDOS_3d%array, 1)
1935 6 : i_y_start = LBOUND(LDOS_3d%array, 2)
1936 6 : i_y_end = UBOUND(LDOS_3d%array, 2)
1937 6 : i_z_start = LBOUND(LDOS_3d%array, 3)
1938 6 : i_z_end = UBOUND(LDOS_3d%array, 3)
1939 :
1940 6 : z_start_global = i_z_start
1941 6 : z_end_global = i_z_end
1942 :
1943 6 : CALL bs_env%para_env%min(z_start_global)
1944 6 : CALL bs_env%para_env%max(z_end_global)
1945 6 : n_z = z_end_global - z_start_global + 1
1946 :
1947 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
1948 0 : CPABORT("Please choose a cell that has 90° angles to the z-direction.")
1949 : END IF
1950 : ! for integration, we need the dz and the conversion from H -> eV and a_Bohr -> Å
1951 6 : bs_env%unit_ldos_int_z_inv_Ang2_eV = bs_env%hmat(3, 3)/REAL(n_z, KIND=dp)/evolt/angstrom**2
1952 :
1953 6 : IF (ikp == 1) THEN
1954 30 : ALLOCATE (LDOS_2d(i_x_start:i_x_end, i_y_start:i_y_end, n_E))
1955 6 : LDOS_2d(:, :, :) = 0.0_dp
1956 : END IF
1957 :
1958 6 : CALL cp_cfm_create(cfm_work, cfm_mos_ikp%matrix_struct)
1959 6 : CALL cp_cfm_create(cfm_weighted_dm_ikp, cfm_mos_ikp%matrix_struct)
1960 6 : CALL cp_fm_create(fm_weighted_dm_MIC, cfm_mos_ikp%matrix_struct)
1961 6 : IF (my_do_spinor) THEN
1962 2 : CALL cp_fm_create(fm_non_spinor, cfm_non_spinor%matrix_struct)
1963 : END IF
1964 :
1965 : CALL cp_cfm_get_info(matrix=cfm_mos_ikp, &
1966 : ncol_global=n_mo, &
1967 : ncol_local=ncol_local, &
1968 6 : col_indices=col_indices)
1969 :
1970 6 : NULLIFY (weighted_dm_MIC)
1971 6 : CALL dbcsr_allocate_matrix_set(weighted_dm_MIC, 1)
1972 6 : ALLOCATE (weighted_dm_MIC(1)%matrix)
1973 : CALL dbcsr_create(weighted_dm_MIC(1)%matrix, template=bs_env%mat_ao_ao%matrix, &
1974 6 : matrix_type=dbcsr_type_symmetric)
1975 :
1976 1678 : DO i_E = 1, n_E
1977 :
1978 1672 : energy = E_min + i_E*energy_step
1979 :
1980 1672 : is_any_weight_non_zero = .FALSE.
1981 :
1982 20950 : DO j_col = 1, ncol_local
1983 :
1984 19278 : j_mo = col_indices(j_col)
1985 :
1986 19278 : IF (my_do_spinor) THEN
1987 : spin_degeneracy = 1.0_dp
1988 : ELSE
1989 10818 : spin_degeneracy = bs_env%spin_degeneracy
1990 : END IF
1991 :
1992 19278 : weight = Gaussian(energy - eigenval(j_mo), broadening)*spin_degeneracy
1993 :
1994 144099 : cfm_work%local_data(:, j_col) = cfm_mos_ikp%local_data(:, j_col)*weight
1995 :
1996 20950 : IF (weight > 1.0E-5_dp) is_any_weight_non_zero = .TRUE.
1997 :
1998 : END DO
1999 :
2000 1672 : CALL bs_env%para_env%sync()
2001 1672 : CALL bs_env%para_env%sum(is_any_weight_non_zero)
2002 1672 : CALL bs_env%para_env%sync()
2003 :
2004 : ! cycle if there are no states at the energy i_E
2005 1678 : IF (is_any_weight_non_zero) THEN
2006 :
2007 : CALL parallel_gemm('N', 'C', n_mo, n_mo, n_mo, z_one, &
2008 24 : cfm_mos_ikp, cfm_work, z_zero, cfm_weighted_dm_ikp)
2009 :
2010 24 : IF (my_do_spinor) THEN
2011 :
2012 : ! contribution from up,up to fm_non_spinor
2013 8 : CALL get_cfm_submat(cfm_non_spinor, cfm_weighted_dm_ikp, 1, 1)
2014 8 : CALL cp_fm_set_all(fm_non_spinor, 0.0_dp)
2015 : CALL MIC_contribution_from_ikp(bs_env, qs_env, fm_non_spinor, &
2016 : cfm_non_spinor, ikp, bs_env%kpoints_DOS, &
2017 8 : "ORB", bs_env%kpoints_DOS%wkp(ikp))
2018 :
2019 : ! add contribution from down,down to fm_non_spinor
2020 8 : CALL get_cfm_submat(cfm_non_spinor, cfm_weighted_dm_ikp, n_mo/2, n_mo/2)
2021 : CALL MIC_contribution_from_ikp(bs_env, qs_env, fm_non_spinor, &
2022 : cfm_non_spinor, ikp, bs_env%kpoints_DOS, &
2023 8 : "ORB", bs_env%kpoints_DOS%wkp(ikp))
2024 : CALL copy_fm_to_dbcsr(fm_non_spinor, weighted_dm_MIC(1)%matrix, &
2025 8 : keep_sparsity=.FALSE.)
2026 : ELSE
2027 16 : CALL cp_fm_set_all(fm_weighted_dm_MIC, 0.0_dp)
2028 : CALL MIC_contribution_from_ikp(bs_env, qs_env, fm_weighted_dm_MIC, &
2029 : cfm_weighted_dm_ikp, ikp, bs_env%kpoints_DOS, &
2030 16 : "ORB", bs_env%kpoints_DOS%wkp(ikp))
2031 : CALL copy_fm_to_dbcsr(fm_weighted_dm_MIC, weighted_dm_MIC(1)%matrix, &
2032 16 : keep_sparsity=.FALSE.)
2033 : END IF
2034 :
2035 338424 : LDOS_3d%array(:, :, :) = 0.0_dp
2036 :
2037 : CALL calculate_rho_elec(matrix_p_kp=weighted_dm_MIC, &
2038 : rho=LDOS_3d, &
2039 : rho_gspace=rho_g, &
2040 24 : ks_env=ks_env)
2041 :
2042 504 : DO i_z = i_z_start, i_z_end
2043 338424 : LDOS_2d(:, :, i_E) = LDOS_2d(:, :, i_E) + LDOS_3d%array(:, :, i_z)
2044 : END DO
2045 :
2046 : END IF
2047 :
2048 : END DO
2049 :
2050 : ! set back nimages
2051 6 : dft_control%nimages = nimages
2052 :
2053 6 : CALL auxbas_pw_pool%give_back_pw(LDOS_3d)
2054 6 : CALL auxbas_pw_pool%give_back_pw(rho_g)
2055 :
2056 6 : CALL cp_cfm_release(cfm_work)
2057 6 : CALL cp_cfm_release(cfm_weighted_dm_ikp)
2058 :
2059 6 : CALL cp_fm_release(fm_weighted_dm_MIC)
2060 :
2061 6 : CALL dbcsr_deallocate_matrix_set(weighted_dm_MIC)
2062 :
2063 6 : IF (my_do_spinor) THEN
2064 2 : CALL cp_fm_release(fm_non_spinor)
2065 : END IF
2066 :
2067 6 : CALL timestop(handle)
2068 :
2069 6 : END SUBROUTINE add_to_LDOS_2d
2070 :
2071 : ! **************************************************************************************************
2072 : !> \brief ...
2073 : !> \param eigenval_spinor ...
2074 : !> \param ikp_for_file ...
2075 : !> \param ikp ...
2076 : !> \param bs_env ...
2077 : !> \param eigenval_spinor_G0W0 ...
2078 : ! **************************************************************************************************
2079 296 : SUBROUTINE write_SOC_eigenvalues(eigenval_spinor, ikp_for_file, ikp, bs_env, eigenval_spinor_G0W0)
2080 :
2081 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_spinor
2082 : INTEGER :: ikp_for_file, ikp
2083 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2084 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: eigenval_spinor_G0W0
2085 :
2086 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_SOC_eigenvalues'
2087 :
2088 : CHARACTER(len=3) :: occ_vir
2089 : CHARACTER(LEN=default_string_length) :: fname
2090 : INTEGER :: handle, i_mo, iunit, n_occ_spinor
2091 :
2092 296 : CALL timeset(routineN, handle)
2093 :
2094 296 : fname = "bandstructure_SCF_and_G0W0_plus_SOC"
2095 :
2096 296 : IF (bs_env%para_env%is_source()) THEN
2097 :
2098 148 : IF (ikp_for_file == 1) THEN
2099 : CALL open_file(TRIM(fname), unit_number=iunit, file_status="REPLACE", &
2100 11 : file_action="WRITE")
2101 : ELSE
2102 : CALL open_file(TRIM(fname), unit_number=iunit, file_status="OLD", &
2103 137 : file_action="WRITE", file_position="APPEND")
2104 : END IF
2105 :
2106 148 : WRITE (iunit, "(A)") " "
2107 148 : WRITE (iunit, "(A10,I7,A25,3F10.4)") "kpoint: ", ikp_for_file, "coordinate: ", &
2108 740 : bs_env%kpoints_DOS%xkp(:, ikp)
2109 148 : WRITE (iunit, "(A)") " "
2110 :
2111 148 : IF (PRESENT(eigenval_spinor_G0W0)) THEN
2112 : ! SCF+SOC and G0W0+SOC eigenvalues
2113 148 : WRITE (iunit, "(A5,A12,2A22)") "n", "k", "ϵ_nk^DFT+SOC (eV)", "ϵ_nk^G0W0+SOC (eV)"
2114 : ELSE
2115 : ! SCF+SOC eigenvalues only
2116 0 : WRITE (iunit, "(A5,A12,A22)") "n", "k", "ϵ_nk^DFT+SOC (eV)"
2117 : END IF
2118 :
2119 148 : n_occ_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin)
2120 :
2121 3932 : DO i_mo = 1, SIZE(eigenval_spinor)
2122 3784 : IF (i_mo <= n_occ_spinor) occ_vir = 'occ'
2123 3784 : IF (i_mo > n_occ_spinor) occ_vir = 'vir'
2124 3932 : IF (PRESENT(eigenval_spinor_G0W0)) THEN
2125 : ! SCF+SOC and G0W0+SOC eigenvalues
2126 3784 : WRITE (iunit, "(I5,3A,I5,4F16.3,2F17.3)") i_mo, ' (', occ_vir, ') ', &
2127 7568 : ikp_for_file, eigenval_spinor(i_mo)*evolt, eigenval_spinor_G0W0(i_mo)*evolt
2128 : ELSE
2129 : ! SCF+SOC eigenvalues only
2130 0 : WRITE (iunit, "(I5,3A,I5,4F16.3,F17.3)") i_mo, ' (', occ_vir, ') ', &
2131 0 : ikp_for_file, eigenval_spinor(i_mo)*evolt
2132 : END IF
2133 : END DO
2134 :
2135 148 : CALL close_file(iunit)
2136 :
2137 : END IF
2138 :
2139 296 : CALL timestop(handle)
2140 :
2141 296 : END SUBROUTINE write_SOC_eigenvalues
2142 :
2143 : ! **************************************************************************************************
2144 : !> \brief ...
2145 : !> \param int_number ...
2146 : !> \return ...
2147 : ! **************************************************************************************************
2148 0 : PURE FUNCTION count_digits(int_number)
2149 :
2150 : INTEGER, INTENT(IN) :: int_number
2151 : INTEGER :: count_digits
2152 :
2153 : INTEGER :: digitCount, tempInt
2154 :
2155 0 : digitCount = 0
2156 :
2157 0 : tempInt = int_number
2158 :
2159 0 : DO WHILE (tempInt /= 0)
2160 0 : tempInt = tempInt/10
2161 0 : digitCount = digitCount + 1
2162 : END DO
2163 :
2164 0 : count_digits = digitCount
2165 :
2166 0 : END FUNCTION count_digits
2167 :
2168 : ! **************************************************************************************************
2169 : !> \brief ...
2170 : !> \param band_edges ...
2171 : !> \param scf_gw_soc ...
2172 : !> \param bs_env ...
2173 : ! **************************************************************************************************
2174 358 : SUBROUTINE write_band_edges(band_edges, scf_gw_soc, bs_env)
2175 :
2176 : TYPE(band_edges_type) :: band_edges
2177 : CHARACTER(LEN=*) :: scf_gw_soc
2178 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2179 :
2180 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_band_edges'
2181 :
2182 : CHARACTER(LEN=17) :: print_format
2183 : INTEGER :: handle, u
2184 :
2185 358 : CALL timeset(routineN, handle)
2186 :
2187 : ! print format
2188 358 : print_format = "(T2,2A,T61,F20.3)"
2189 :
2190 358 : u = bs_env%unit_nr
2191 358 : IF (u > 0) THEN
2192 179 : WRITE (u, '(T2,A)') ''
2193 179 : WRITE (u, print_format) scf_gw_soc, ' valence band maximum (eV):', band_edges%VBM*evolt
2194 179 : WRITE (u, print_format) scf_gw_soc, ' conduction band minimum (eV):', band_edges%CBM*evolt
2195 179 : WRITE (u, print_format) scf_gw_soc, ' indirect band gap (eV):', band_edges%IDBG*evolt
2196 179 : WRITE (u, print_format) scf_gw_soc, ' direct band gap (eV):', band_edges%DBG*evolt
2197 : END IF
2198 :
2199 358 : CALL timestop(handle)
2200 :
2201 358 : END SUBROUTINE write_band_edges
2202 :
2203 : ! **************************************************************************************************
2204 : !> \brief ...
2205 : !> \param DOS ...
2206 : !> \param PDOS ...
2207 : !> \param bs_env ...
2208 : !> \param qs_env ...
2209 : !> \param scf_gw_soc ...
2210 : !> \param E_min ...
2211 : !> \param E_VBM ...
2212 : ! **************************************************************************************************
2213 72 : SUBROUTINE write_dos_pdos(DOS, PDOS, bs_env, qs_env, scf_gw_soc, E_min, E_VBM)
2214 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: DOS
2215 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: PDOS
2216 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2217 : TYPE(qs_environment_type), POINTER :: qs_env
2218 : CHARACTER(LEN=*) :: scf_gw_soc
2219 : REAL(KIND=dp) :: E_min, E_VBM
2220 :
2221 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_dos_pdos'
2222 :
2223 : CHARACTER(LEN=3), DIMENSION(100) :: elements
2224 : CHARACTER(LEN=default_string_length) :: atom_name, fname, output_string
2225 : INTEGER :: handle, i_E, i_kind, iatom, iunit, n_A, &
2226 : n_E, nkind
2227 : REAL(KIND=dp) :: energy
2228 72 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2229 :
2230 72 : CALL timeset(routineN, handle)
2231 :
2232 72 : WRITE (fname, "(3A)") "DOS_PDOS_", scf_gw_soc, ".out"
2233 :
2234 72 : n_E = SIZE(PDOS, 1)
2235 72 : nkind = SIZE(PDOS, 2)
2236 72 : CALL get_qs_env(qs_env, particle_set=particle_set)
2237 :
2238 72 : IF (bs_env%para_env%is_source()) THEN
2239 :
2240 36 : CALL open_file(TRIM(fname), unit_number=iunit, file_status="REPLACE", file_action="WRITE")
2241 :
2242 36 : n_A = 2 + nkind
2243 :
2244 140 : DO iatom = 1, bs_env%n_atom
2245 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
2246 104 : kind_number=i_kind, name=atom_name)
2247 140 : elements(i_kind) = atom_name(1:3)
2248 : END DO
2249 :
2250 36 : WRITE (output_string, "(A,I1,A)") "(", n_A, "A)"
2251 :
2252 36 : WRITE (iunit, TRIM(output_string)) "Energy-E_F (eV) DOS (1/eV) PDOS (1/eV) ", &
2253 72 : " of atom type ", elements(1:nkind)
2254 :
2255 36 : WRITE (output_string, "(A,I1,A)") "(", n_A, "F13.5)"
2256 :
2257 73686 : DO i_E = 1, n_E
2258 : ! energy is relative to valence band maximum => - E_VBM
2259 73650 : energy = E_min + i_E*bs_env%energy_step_DOS - E_VBM
2260 220986 : WRITE (iunit, TRIM(output_string)) energy*evolt, DOS(i_E)/evolt, PDOS(i_E, :)/evolt
2261 : END DO
2262 :
2263 36 : CALL close_file(iunit)
2264 :
2265 : END IF
2266 :
2267 72 : CALL timestop(handle)
2268 :
2269 72 : END SUBROUTINE write_dos_pdos
2270 :
2271 : ! **************************************************************************************************
2272 : !> \brief ...
2273 : !> \param energy ...
2274 : !> \param broadening ...
2275 : !> \return ...
2276 : ! **************************************************************************************************
2277 70597222 : PURE FUNCTION Gaussian(energy, broadening)
2278 :
2279 : REAL(KIND=dp), INTENT(IN) :: energy, broadening
2280 : REAL(KIND=dp) :: Gaussian
2281 :
2282 70597222 : IF (ABS(energy) < 5*broadening) THEN
2283 111792 : Gaussian = 1.0_dp/broadening/SQRT(twopi)*EXP(-0.5_dp*energy**2/broadening**2)
2284 : ELSE
2285 : Gaussian = 0.0_dp
2286 : END IF
2287 :
2288 70597222 : END FUNCTION Gaussian
2289 :
2290 : ! **************************************************************************************************
2291 : !> \brief ...
2292 : !> \param proj_mo_on_kind ...
2293 : !> \param qs_env ...
2294 : !> \param cfm_mos ...
2295 : !> \param cfm_s ...
2296 : ! **************************************************************************************************
2297 466 : SUBROUTINE compute_proj_mo_on_kind(proj_mo_on_kind, qs_env, cfm_mos, cfm_s)
2298 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: proj_mo_on_kind
2299 : TYPE(qs_environment_type), POINTER :: qs_env
2300 : TYPE(cp_cfm_type) :: cfm_mos, cfm_s
2301 :
2302 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_proj_mo_on_kind'
2303 :
2304 : INTEGER :: handle, i_atom, i_global, i_kind, i_row, &
2305 : j_col, n_ao, n_mo, ncol_local, nkind, &
2306 : nrow_local
2307 466 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_bf, kind_of
2308 466 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2309 466 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2310 : TYPE(cp_cfm_type) :: cfm_proj, cfm_s_i_kind, cfm_work
2311 : TYPE(cp_fm_type) :: fm_proj_im, fm_proj_re
2312 :
2313 466 : CALL timeset(routineN, handle)
2314 :
2315 466 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, nkind=nkind)
2316 466 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
2317 :
2318 : CALL cp_cfm_get_info(matrix=cfm_mos, &
2319 : nrow_global=n_mo, &
2320 : nrow_local=nrow_local, &
2321 : ncol_local=ncol_local, &
2322 : row_indices=row_indices, &
2323 466 : col_indices=col_indices)
2324 :
2325 466 : n_ao = qs_env%bs_env%n_ao
2326 :
2327 1398 : ALLOCATE (atom_from_bf(n_ao))
2328 466 : CALL get_atom_index_from_basis_function_index(qs_env, atom_from_bf, n_ao, "ORB")
2329 :
2330 466 : proj_mo_on_kind(:, :) = 0.0_dp
2331 :
2332 466 : CALL cp_cfm_create(cfm_s_i_kind, cfm_s%matrix_struct)
2333 466 : CALL cp_cfm_create(cfm_work, cfm_s%matrix_struct)
2334 466 : CALL cp_cfm_create(cfm_proj, cfm_s%matrix_struct)
2335 466 : CALL cp_fm_create(fm_proj_re, cfm_s%matrix_struct)
2336 466 : CALL cp_fm_create(fm_proj_im, cfm_s%matrix_struct)
2337 :
2338 1286 : DO i_kind = 1, nkind
2339 :
2340 820 : CALL cp_cfm_to_cfm(cfm_s, cfm_s_i_kind)
2341 :
2342 : ! set entries in overlap matrix to zero which do not belong to atoms of i_kind
2343 10676 : DO j_col = 1, ncol_local
2344 74010 : DO i_row = 1, nrow_local
2345 :
2346 63334 : i_global = row_indices(i_row)
2347 :
2348 63334 : IF (i_global <= n_ao) THEN
2349 63334 : i_atom = atom_from_bf(i_global)
2350 0 : ELSE IF (i_global <= 2*n_ao) THEN
2351 0 : i_atom = atom_from_bf(i_global - n_ao)
2352 : ELSE
2353 0 : CPABORT("Wrong indices.")
2354 : END IF
2355 :
2356 73190 : IF (i_kind /= kind_of(i_atom)) THEN
2357 28849 : cfm_s_i_kind%local_data(i_row, j_col) = z_zero
2358 : END IF
2359 :
2360 : END DO
2361 : END DO
2362 :
2363 : CALL parallel_gemm('N', 'N', n_mo, n_mo, n_mo, z_one, &
2364 820 : cfm_s_i_kind, cfm_mos, z_zero, cfm_work)
2365 : CALL parallel_gemm('C', 'N', n_mo, n_mo, n_mo, z_one, &
2366 820 : cfm_mos, cfm_work, z_zero, cfm_proj)
2367 :
2368 820 : CALL cp_cfm_to_fm(cfm_proj, fm_proj_re, fm_proj_im)
2369 :
2370 820 : CALL cp_fm_get_diag(fm_proj_im, proj_mo_on_kind(:, i_kind))
2371 1286 : CALL cp_fm_get_diag(fm_proj_re, proj_mo_on_kind(:, i_kind))
2372 :
2373 : END DO ! i_kind
2374 :
2375 466 : CALL cp_cfm_release(cfm_s_i_kind)
2376 466 : CALL cp_cfm_release(cfm_work)
2377 466 : CALL cp_cfm_release(cfm_proj)
2378 466 : CALL cp_fm_release(fm_proj_re)
2379 466 : CALL cp_fm_release(fm_proj_im)
2380 :
2381 466 : CALL timestop(handle)
2382 :
2383 1864 : END SUBROUTINE compute_proj_mo_on_kind
2384 :
2385 : ! **************************************************************************************************
2386 : !> \brief ...
2387 : !> \param cfm_spinor_ikp ...
2388 : !> \param cfm_spinor_Gamma ...
2389 : !> \param fm_struct_non_spinor ...
2390 : !> \param ikp ...
2391 : !> \param qs_env ...
2392 : !> \param kpoints ...
2393 : !> \param basis_type ...
2394 : ! **************************************************************************************************
2395 120 : SUBROUTINE cfm_ikp_from_cfm_spinor_Gamma(cfm_spinor_ikp, cfm_spinor_Gamma, fm_struct_non_spinor, &
2396 : ikp, qs_env, kpoints, basis_type)
2397 : TYPE(cp_cfm_type) :: cfm_spinor_ikp, cfm_spinor_Gamma
2398 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_non_spinor
2399 : INTEGER :: ikp
2400 : TYPE(qs_environment_type), POINTER :: qs_env
2401 : TYPE(kpoint_type), POINTER :: kpoints
2402 : CHARACTER(LEN=*) :: basis_type
2403 :
2404 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cfm_ikp_from_cfm_spinor_Gamma'
2405 :
2406 : INTEGER :: handle, i_block, i_offset, j_block, &
2407 : j_offset, n_ao
2408 : TYPE(cp_cfm_type) :: cfm_non_spinor_Gamma, cfm_non_spinor_ikp
2409 : TYPE(cp_fm_type) :: fm_non_spinor_Gamma_im, &
2410 : fm_non_spinor_Gamma_re
2411 :
2412 20 : CALL timeset(routineN, handle)
2413 :
2414 20 : CALL cp_cfm_create(cfm_non_spinor_Gamma, fm_struct_non_spinor)
2415 20 : CALL cp_cfm_create(cfm_non_spinor_ikp, fm_struct_non_spinor)
2416 20 : CALL cp_fm_create(fm_non_spinor_Gamma_re, fm_struct_non_spinor)
2417 20 : CALL cp_fm_create(fm_non_spinor_Gamma_im, fm_struct_non_spinor)
2418 :
2419 20 : CALL cp_cfm_get_info(cfm_non_spinor_Gamma, nrow_global=n_ao)
2420 :
2421 20 : CALL cp_cfm_set_all(cfm_spinor_ikp, z_zero)
2422 :
2423 60 : DO i_block = 0, 1
2424 140 : DO j_block = 0, 1
2425 80 : i_offset = i_block*n_ao + 1
2426 80 : j_offset = j_block*n_ao + 1
2427 80 : CALL get_cfm_submat(cfm_non_spinor_Gamma, cfm_spinor_Gamma, i_offset, j_offset)
2428 80 : CALL cp_cfm_to_fm(cfm_non_spinor_Gamma, fm_non_spinor_Gamma_re, fm_non_spinor_Gamma_im)
2429 :
2430 : ! transform real part of Gamma-point matrix to ikp
2431 : CALL cfm_ikp_from_fm_Gamma(cfm_non_spinor_ikp, fm_non_spinor_Gamma_re, &
2432 80 : ikp, qs_env, kpoints, basis_type)
2433 80 : CALL add_cfm_submat(cfm_spinor_ikp, cfm_non_spinor_ikp, i_offset, j_offset)
2434 :
2435 : ! transform imag part of Gamma-point matrix to ikp
2436 : CALL cfm_ikp_from_fm_Gamma(cfm_non_spinor_ikp, fm_non_spinor_Gamma_im, &
2437 80 : ikp, qs_env, kpoints, basis_type)
2438 120 : CALL add_cfm_submat(cfm_spinor_ikp, cfm_non_spinor_ikp, i_offset, j_offset, gaussi)
2439 :
2440 : END DO
2441 : END DO
2442 :
2443 20 : CALL cp_cfm_release(cfm_non_spinor_Gamma)
2444 20 : CALL cp_cfm_release(cfm_non_spinor_ikp)
2445 20 : CALL cp_fm_release(fm_non_spinor_Gamma_re)
2446 20 : CALL cp_fm_release(fm_non_spinor_Gamma_im)
2447 :
2448 20 : CALL timestop(handle)
2449 :
2450 20 : END SUBROUTINE cfm_ikp_from_cfm_spinor_Gamma
2451 :
2452 : ! **************************************************************************************************
2453 : !> \brief ...
2454 : !> \param cfm_ikp ...
2455 : !> \param fm_Gamma ...
2456 : !> \param ikp ...
2457 : !> \param qs_env ...
2458 : !> \param kpoints ...
2459 : !> \param basis_type ...
2460 : ! **************************************************************************************************
2461 7076 : SUBROUTINE cfm_ikp_from_fm_Gamma(cfm_ikp, fm_Gamma, ikp, qs_env, kpoints, basis_type)
2462 : TYPE(cp_cfm_type) :: cfm_ikp
2463 : TYPE(cp_fm_type) :: fm_Gamma
2464 : INTEGER :: ikp
2465 : TYPE(qs_environment_type), POINTER :: qs_env
2466 : TYPE(kpoint_type), POINTER :: kpoints
2467 : CHARACTER(LEN=*) :: basis_type
2468 :
2469 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cfm_ikp_from_fm_Gamma'
2470 :
2471 : INTEGER :: col_global, handle, i_atom, i_atom_old, i_cell, i_mic_cell, i_row, j_atom, &
2472 : j_atom_old, j_cell, j_col, n_bf, ncol_local, nrow_local, num_cells, row_global
2473 7076 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_bf
2474 7076 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2475 7076 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
2476 : LOGICAL :: i_cell_is_the_minimum_image_cell
2477 : REAL(KIND=dp) :: abs_rab_cell_i, abs_rab_cell_j, arg
2478 : REAL(KIND=dp), DIMENSION(3) :: cell_vector, cell_vector_j, rab_cell_i, &
2479 : rab_cell_j
2480 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
2481 : TYPE(cell_type), POINTER :: cell
2482 7076 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2483 :
2484 7076 : CALL timeset(routineN, handle)
2485 :
2486 7076 : IF (.NOT. ASSOCIATED(cfm_ikp%local_data)) THEN
2487 2860 : CALL cp_cfm_create(cfm_ikp, fm_Gamma%matrix_struct)
2488 : END IF
2489 7076 : CALL cp_cfm_set_all(cfm_ikp, z_zero)
2490 :
2491 : CALL cp_fm_get_info(matrix=fm_Gamma, &
2492 : nrow_local=nrow_local, &
2493 : ncol_local=ncol_local, &
2494 : row_indices=row_indices, &
2495 7076 : col_indices=col_indices)
2496 :
2497 : ! get number of basis functions (bf) for different basis sets
2498 7076 : IF (basis_type == "ORB") THEN
2499 4656 : n_bf = qs_env%bs_env%n_ao
2500 2420 : ELSE IF (basis_type == "RI_AUX") THEN
2501 2420 : n_bf = qs_env%bs_env%n_RI
2502 : ELSE
2503 0 : CPABORT("Only ORB and RI_AUX basis implemented.")
2504 : END IF
2505 :
2506 21228 : ALLOCATE (atom_from_bf(n_bf))
2507 7076 : CALL get_atom_index_from_basis_function_index(qs_env, atom_from_bf, n_bf, basis_type)
2508 :
2509 7076 : NULLIFY (cell, particle_set)
2510 7076 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
2511 7076 : CALL get_cell(cell=cell, h=hmat)
2512 :
2513 7076 : index_to_cell => kpoints%index_to_cell
2514 :
2515 7076 : num_cells = SIZE(index_to_cell, 2)
2516 7076 : i_atom_old = 0
2517 7076 : j_atom_old = 0
2518 :
2519 107916 : DO j_col = 1, ncol_local
2520 1398502 : DO i_row = 1, nrow_local
2521 :
2522 1290586 : row_global = row_indices(i_row)
2523 1290586 : col_global = col_indices(j_col)
2524 :
2525 1290586 : i_atom = atom_from_bf(row_global)
2526 1290586 : j_atom = atom_from_bf(col_global)
2527 :
2528 : ! we only need to check for new MIC cell for new i_atom-j_atom pair
2529 1290586 : IF (i_atom /= i_atom_old .OR. j_atom /= j_atom_old) THEN
2530 756984 : DO i_cell = 1, num_cells
2531 :
2532 : ! only check nearest neigbors
2533 1871888 : IF (ANY(ABS(index_to_cell(1:3, i_cell)) > 1)) CYCLE
2534 :
2535 6033280 : cell_vector(1:3) = MATMUL(hmat, REAL(index_to_cell(1:3, i_cell), dp))
2536 :
2537 : rab_cell_i(1:3) = pbc(particle_set(i_atom)%r(1:3), cell) - &
2538 1508320 : (pbc(particle_set(j_atom)%r(1:3), cell) + cell_vector(1:3))
2539 377080 : abs_rab_cell_i = SQRT(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2)
2540 :
2541 : ! minimum image convention
2542 377080 : i_cell_is_the_minimum_image_cell = .TRUE.
2543 3796088 : DO j_cell = 1, num_cells
2544 54704128 : cell_vector_j(1:3) = MATMUL(hmat, REAL(index_to_cell(1:3, j_cell), dp))
2545 : rab_cell_j(1:3) = pbc(particle_set(i_atom)%r(1:3), cell) - &
2546 13676032 : (pbc(particle_set(j_atom)%r(1:3), cell) + cell_vector_j(1:3))
2547 3419008 : abs_rab_cell_j = SQRT(rab_cell_j(1)**2 + rab_cell_j(2)**2 + rab_cell_j(3)**2)
2548 :
2549 3796088 : IF (abs_rab_cell_i > abs_rab_cell_j + 1.0E-6_dp) THEN
2550 676826 : i_cell_is_the_minimum_image_cell = .FALSE.
2551 : END IF
2552 : END DO
2553 :
2554 573312 : IF (i_cell_is_the_minimum_image_cell) THEN
2555 196232 : i_mic_cell = i_cell
2556 : END IF
2557 :
2558 : END DO ! i_cell
2559 : END IF
2560 :
2561 : arg = REAL(index_to_cell(1, i_mic_cell), dp)*kpoints%xkp(1, ikp) + &
2562 : REAL(index_to_cell(2, i_mic_cell), dp)*kpoints%xkp(2, ikp) + &
2563 1290586 : REAL(index_to_cell(3, i_mic_cell), dp)*kpoints%xkp(3, ikp)
2564 :
2565 : cfm_ikp%local_data(i_row, j_col) = COS(twopi*arg)*fm_Gamma%local_data(i_row, j_col)*z_one + &
2566 1290586 : SIN(twopi*arg)*fm_Gamma%local_data(i_row, j_col)*gaussi
2567 :
2568 1290586 : j_atom_old = j_atom
2569 1391426 : i_atom_old = i_atom
2570 :
2571 : END DO ! j_col
2572 : END DO ! i_row
2573 :
2574 7076 : CALL timestop(handle)
2575 :
2576 21228 : END SUBROUTINE cfm_ikp_from_fm_Gamma
2577 :
2578 : ! **************************************************************************************************
2579 : !> \brief ...
2580 : !> \param bs_env ...
2581 : !> \param qs_env ...
2582 : !> \param fm_W_MIC_freq_j ...
2583 : !> \param cfm_W_ikp_freq_j ...
2584 : !> \param ikp ...
2585 : !> \param kpoints ...
2586 : !> \param basis_type ...
2587 : !> \param wkp_ext ...
2588 : ! **************************************************************************************************
2589 2484 : SUBROUTINE MIC_contribution_from_ikp(bs_env, qs_env, fm_W_MIC_freq_j, &
2590 : cfm_W_ikp_freq_j, ikp, kpoints, basis_type, wkp_ext)
2591 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2592 : TYPE(qs_environment_type), POINTER :: qs_env
2593 : TYPE(cp_fm_type) :: fm_W_MIC_freq_j
2594 : TYPE(cp_cfm_type) :: cfm_W_ikp_freq_j
2595 : INTEGER, INTENT(IN) :: ikp
2596 : TYPE(kpoint_type), POINTER :: kpoints
2597 : CHARACTER(LEN=*) :: basis_type
2598 : REAL(KIND=dp), OPTIONAL :: wkp_ext
2599 :
2600 : CHARACTER(LEN=*), PARAMETER :: routineN = 'MIC_contribution_from_ikp'
2601 :
2602 : INTEGER :: handle, i_bf, iatom, iatom_old, irow, &
2603 : j_bf, jatom, jatom_old, jcol, n_bf, &
2604 : ncol_local, nrow_local, num_cells
2605 2484 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_bf_index
2606 2484 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2607 2484 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
2608 : REAL(KIND=dp) :: contribution, weight_im, weight_re, &
2609 : wkp_of_ikp
2610 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
2611 2484 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp
2612 2484 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
2613 : TYPE(cell_type), POINTER :: cell
2614 2484 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2615 :
2616 2484 : CALL timeset(routineN, handle)
2617 :
2618 : ! get number of basis functions (bf) for different basis sets
2619 2484 : IF (basis_type == "ORB") THEN
2620 32 : n_bf = qs_env%bs_env%n_ao
2621 2452 : ELSE IF (basis_type == "RI_AUX") THEN
2622 2452 : n_bf = qs_env%bs_env%n_RI
2623 : ELSE
2624 0 : CPABORT("Only ORB and RI_AUX basis implemented.")
2625 : END IF
2626 :
2627 7452 : ALLOCATE (atom_from_bf_index(n_bf))
2628 2484 : CALL get_atom_index_from_basis_function_index(qs_env, atom_from_bf_index, n_bf, basis_type)
2629 :
2630 2484 : NULLIFY (cell, particle_set)
2631 2484 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
2632 2484 : CALL get_cell(cell=cell, h=hmat)
2633 :
2634 : CALL cp_cfm_get_info(matrix=cfm_W_ikp_freq_j, &
2635 : nrow_local=nrow_local, &
2636 : ncol_local=ncol_local, &
2637 : row_indices=row_indices, &
2638 2484 : col_indices=col_indices)
2639 :
2640 2484 : CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp)
2641 2484 : index_to_cell => kpoints%index_to_cell
2642 2484 : num_cells = SIZE(index_to_cell, 2)
2643 :
2644 2484 : iatom_old = 0
2645 2484 : jatom_old = 0
2646 :
2647 56384 : DO jcol = 1, ncol_local
2648 1066438 : DO irow = 1, nrow_local
2649 :
2650 1010054 : i_bf = row_indices(irow)
2651 1010054 : j_bf = col_indices(jcol)
2652 :
2653 1010054 : iatom = atom_from_bf_index(i_bf)
2654 1010054 : jatom = atom_from_bf_index(j_bf)
2655 :
2656 1010054 : IF (PRESENT(wkp_ext)) THEN
2657 3496 : wkp_of_ikp = wkp_ext
2658 : ELSE
2659 1040270 : SELECT CASE (bs_env%l_RI(i_bf) + bs_env%l_RI(j_bf))
2660 : CASE (0)
2661 : ! both RI functions are s-functions, k-extrapolation for 2D and 3D
2662 33712 : wkp_of_ikp = wkp(ikp)
2663 : CASE (1)
2664 : ! one function is an s-function, the other a p-function, k-extrapolation for 3D
2665 144552 : wkp_of_ikp = bs_env%wkp_s_p(ikp)
2666 : CASE DEFAULT
2667 : ! for any other matrix element of W, there is no need for extrapolation
2668 1006558 : wkp_of_ikp = bs_env%wkp_no_extra(ikp)
2669 : END SELECT
2670 : END IF
2671 :
2672 1010054 : IF (iatom /= iatom_old .OR. jatom /= jatom_old) THEN
2673 :
2674 : CALL compute_weight_re_im(weight_re, weight_im, &
2675 : num_cells, iatom, jatom, xkp(1:3, ikp), wkp_of_ikp, &
2676 104552 : cell, index_to_cell, hmat, particle_set)
2677 :
2678 104552 : iatom_old = iatom
2679 104552 : jatom_old = jatom
2680 :
2681 : END IF
2682 :
2683 : contribution = weight_re*REAL(cfm_W_ikp_freq_j%local_data(irow, jcol)) + &
2684 1010054 : weight_im*AIMAG(cfm_W_ikp_freq_j%local_data(irow, jcol))
2685 :
2686 : fm_W_MIC_freq_j%local_data(irow, jcol) = fm_W_MIC_freq_j%local_data(irow, jcol) &
2687 1063954 : + contribution
2688 :
2689 : END DO
2690 : END DO
2691 :
2692 2484 : CALL timestop(handle)
2693 :
2694 7452 : END SUBROUTINE MIC_contribution_from_ikp
2695 :
2696 : ! **************************************************************************************************
2697 : !> \brief ...
2698 : !> \param xkp ...
2699 : !> \param ikp_start ...
2700 : !> \param ikp_end ...
2701 : !> \param grid ...
2702 : ! **************************************************************************************************
2703 124 : SUBROUTINE compute_xkp(xkp, ikp_start, ikp_end, grid)
2704 :
2705 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
2706 : INTEGER :: ikp_start, ikp_end
2707 : INTEGER, DIMENSION(3) :: grid
2708 :
2709 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_xkp'
2710 :
2711 : INTEGER :: handle, i, ix, iy, iz
2712 :
2713 124 : CALL timeset(routineN, handle)
2714 :
2715 124 : i = ikp_start
2716 276 : DO ix = 1, grid(1)
2717 620 : DO iy = 1, grid(2)
2718 1186 : DO iz = 1, grid(3)
2719 :
2720 690 : IF (i > ikp_end) CYCLE
2721 :
2722 672 : xkp(1, i) = REAL(2*ix - grid(1) - 1, KIND=dp)/(2._dp*REAL(grid(1), KIND=dp))
2723 672 : xkp(2, i) = REAL(2*iy - grid(2) - 1, KIND=dp)/(2._dp*REAL(grid(2), KIND=dp))
2724 672 : xkp(3, i) = REAL(2*iz - grid(3) - 1, KIND=dp)/(2._dp*REAL(grid(3), KIND=dp))
2725 1034 : i = i + 1
2726 :
2727 : END DO
2728 : END DO
2729 : END DO
2730 :
2731 124 : CALL timestop(handle)
2732 :
2733 124 : END SUBROUTINE compute_xkp
2734 :
2735 : ! **************************************************************************************************
2736 : !> \brief ...
2737 : !> \param kpoints ...
2738 : !> \param qs_env ...
2739 : ! **************************************************************************************************
2740 176 : SUBROUTINE kpoint_init_cell_index_simple(kpoints, qs_env)
2741 :
2742 : TYPE(kpoint_type), POINTER :: kpoints
2743 : TYPE(qs_environment_type), POINTER :: qs_env
2744 :
2745 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_init_cell_index_simple'
2746 :
2747 : INTEGER :: handle, nimages
2748 : TYPE(mp_para_env_type), POINTER :: para_env
2749 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2750 88 : POINTER :: sab_orb
2751 :
2752 88 : CALL timeset(routineN, handle)
2753 :
2754 88 : NULLIFY (para_env, sab_orb)
2755 88 : CALL get_qs_env(qs_env=qs_env, para_env=para_env, sab_orb=sab_orb)
2756 88 : CALL kpoint_init_cell_index(kpoints, sab_orb, para_env, nimages)
2757 :
2758 88 : CALL timestop(handle)
2759 :
2760 88 : END SUBROUTINE kpoint_init_cell_index_simple
2761 :
2762 : ! **************************************************************************************************
2763 : !> \brief ...
2764 : !> \param qs_env ...
2765 : !> \param bs_env ...
2766 : ! **************************************************************************************************
2767 22 : SUBROUTINE soc(qs_env, bs_env)
2768 : TYPE(qs_environment_type), POINTER :: qs_env
2769 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2770 :
2771 : CHARACTER(LEN=*), PARAMETER :: routineN = 'soc'
2772 :
2773 : INTEGER :: handle
2774 :
2775 22 : CALL timeset(routineN, handle)
2776 :
2777 : ! V^SOC_µν^(α),R = ħ/2 < ϕ_µ cell O | sum_ℓ ΔV_ℓ^SO(r,r') L^(α) | ϕ_ν cell R>, α = x,y,z
2778 : ! see Hartwigsen, Goedecker, Hutter, Eq.(18), (19) (doi.org/10.1103/PhysRevB.58.3641)
2779 22 : CALL V_SOC_xyz_from_pseudopotential(qs_env, bs_env%mat_V_SOC_xyz)
2780 :
2781 : ! Calculate H^SOC_µν,σσ'(k) = sum_α V^SOC_µν^(α)(k)*Pauli-matrix^(α)_σσ'
2782 : ! see Hartwigsen, Goedecker, Hutter, Eq.(18) (doi.org/10.1103/PhysRevB.58.3641)
2783 28 : SELECT CASE (bs_env%small_cell_full_kp_or_large_cell_Gamma)
2784 : CASE (large_cell_Gamma, large_cell_Gamma_ri_rs, non_periodic_ri_rs)
2785 :
2786 : ! H^SOC_µν,σσ' = sum_α V^SOC_µν^(α)*Pauli-matrix^(α)_σσ'
2787 6 : CALL H_KS_spinor_Gamma(bs_env)
2788 :
2789 : CASE (small_cell_full_kp)
2790 :
2791 : ! V^SOC_µν^(α),R -> V^SOC_µν^(α)(k); then calculate spinor H^SOC_µν,σσ'(k) (see above)
2792 22 : CALL H_KS_spinor_kp(qs_env, bs_env)
2793 :
2794 : END SELECT
2795 :
2796 22 : CALL timestop(handle)
2797 :
2798 22 : END SUBROUTINE soc
2799 :
2800 : ! **************************************************************************************************
2801 : !> \brief ...
2802 : !> \param bs_env ...
2803 : ! **************************************************************************************************
2804 6 : SUBROUTINE H_KS_spinor_Gamma(bs_env)
2805 :
2806 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2807 :
2808 : CHARACTER(LEN=*), PARAMETER :: routineN = 'H_KS_spinor_Gamma'
2809 :
2810 : INTEGER :: handle, nao, s
2811 : TYPE(cp_fm_struct_type), POINTER :: str
2812 :
2813 6 : CALL timeset(routineN, handle)
2814 :
2815 6 : CALL cp_fm_get_info(bs_env%fm_ks_Gamma(1), nrow_global=nao)
2816 :
2817 12 : ALLOCATE (bs_env%cfm_SOC_spinor_ao(1))
2818 6 : CALL create_cfm_double(bs_env%cfm_SOC_spinor_ao(1), fm_orig=bs_env%fm_ks_Gamma(1))
2819 6 : CALL cp_cfm_set_all(bs_env%cfm_SOC_spinor_ao(1), z_zero)
2820 :
2821 6 : str => bs_env%fm_ks_Gamma(1)%matrix_struct
2822 :
2823 6 : s = nao + 1
2824 :
2825 : ! careful: inside add_dbcsr_submat, mat_V_SOC_xyz is multiplied by i because the real matrix
2826 : ! mat_V_SOC_xyz is antisymmetric as V_SOC matrix is purely imaginary and Hermitian
2827 : ! V_x * sigma_x: sigma_x = ((0,1),(1,0))
2828 : ! ud block (1,s): +i*V_x
2829 : CALL add_dbcsr_submat(bs_env%cfm_SOC_spinor_ao(1), bs_env%mat_V_SOC_xyz(1, 1)%matrix, &
2830 6 : str, 1, s, z_one, .FALSE.)
2831 : ! du block (s,1): +i*V_x
2832 : CALL add_dbcsr_submat(bs_env%cfm_SOC_spinor_ao(1), bs_env%mat_V_SOC_xyz(1, 1)%matrix, &
2833 6 : str, s, 1, z_one, .FALSE.)
2834 :
2835 : ! V_y * sigma_y: sigma_y = ((0,-i),(i,0))
2836 : ! ud block (1,s): i*(i*V_y) = -V_y (extra gaussi factor)
2837 : CALL add_dbcsr_submat(bs_env%cfm_SOC_spinor_ao(1), bs_env%mat_V_SOC_xyz(2, 1)%matrix, &
2838 6 : str, 1, s, gaussi, .FALSE.)
2839 : ! du block (s,1): -i*(i*V_y) = +V_y (extra -gaussi factor)
2840 : CALL add_dbcsr_submat(bs_env%cfm_SOC_spinor_ao(1), bs_env%mat_V_SOC_xyz(2, 1)%matrix, &
2841 6 : str, s, 1, -gaussi, .FALSE.)
2842 :
2843 : ! V_z * sigma_z: sigma_z = ((1,0),(0,-1))
2844 : ! uu block (1,1): +i*V_z
2845 : CALL add_dbcsr_submat(bs_env%cfm_SOC_spinor_ao(1), bs_env%mat_V_SOC_xyz(3, 1)%matrix, &
2846 6 : str, 1, 1, z_one, .FALSE.)
2847 : ! dd block (s,s): -i*V_z
2848 : CALL add_dbcsr_submat(bs_env%cfm_SOC_spinor_ao(1), bs_env%mat_V_SOC_xyz(3, 1)%matrix, &
2849 6 : str, s, s, -z_one, .FALSE.)
2850 :
2851 6 : CALL timestop(handle)
2852 :
2853 6 : END SUBROUTINE H_KS_spinor_Gamma
2854 :
2855 : ! **************************************************************************************************
2856 : !> \brief ...
2857 : !> \param qs_env ...
2858 : !> \param bs_env ...
2859 : ! **************************************************************************************************
2860 32 : SUBROUTINE H_KS_spinor_kp(qs_env, bs_env)
2861 : TYPE(qs_environment_type), POINTER :: qs_env
2862 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2863 :
2864 : CHARACTER(LEN=*), PARAMETER :: routineN = 'H_KS_spinor_kp'
2865 :
2866 : INTEGER :: handle, i_dim, ikp, n_spin, &
2867 : nkp_bs_and_DOS, s
2868 16 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_scf
2869 : REAL(KIND=dp), DIMENSION(3) :: xkp
2870 : TYPE(cp_cfm_type) :: cfm_V_SOC_xyz_ikp
2871 : TYPE(cp_fm_struct_type), POINTER :: str
2872 : TYPE(kpoint_type), POINTER :: kpoints_scf
2873 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2874 16 : POINTER :: sab_nl
2875 :
2876 16 : CALL timeset(routineN, handle)
2877 :
2878 16 : nkp_bs_and_DOS = bs_env%nkp_bs_and_DOS
2879 16 : n_spin = bs_env%n_spin
2880 16 : s = bs_env%n_ao + 1
2881 16 : str => bs_env%cfm_ks_kp(1, 1)%matrix_struct
2882 :
2883 16 : CALL cp_cfm_create(cfm_V_SOC_xyz_ikp, bs_env%cfm_work_mo%matrix_struct)
2884 :
2885 16 : CALL alloc_cfm_double_array_1d(bs_env%cfm_SOC_spinor_ao, bs_env%cfm_ks_kp(1, 1), nkp_bs_and_DOS)
2886 :
2887 16 : CALL get_qs_env(qs_env, kpoints=kpoints_scf)
2888 :
2889 16 : NULLIFY (sab_nl)
2890 16 : CALL get_kpoint_info(kpoints_scf, sab_nl=sab_nl, cell_to_index=cell_to_index_scf)
2891 :
2892 64 : DO i_dim = 1, 3
2893 :
2894 1018 : DO ikp = 1, nkp_bs_and_DOS
2895 :
2896 3816 : xkp(1:3) = bs_env%kpoints_DOS%xkp(1:3, ikp)
2897 :
2898 954 : CALL cp_cfm_set_all(cfm_V_SOC_xyz_ikp, z_zero)
2899 :
2900 : CALL rsmat_to_kp(bs_env%mat_V_SOC_xyz, i_dim, xkp, cell_to_index_scf, &
2901 954 : sab_nl, bs_env, cfm_V_SOC_xyz_ikp, imag_rs_mat=.TRUE.)
2902 :
2903 : ! multiply V_SOC with i because bs_env%mat_V_SOC_xyz stores imag. part (real part = 0)
2904 954 : CALL cp_cfm_scale(gaussi, cfm_V_SOC_xyz_ikp)
2905 :
2906 48 : SELECT CASE (i_dim)
2907 : CASE (1)
2908 : ! add V^SOC_x * σ_x for σ_x = ( (0,1) (1,0) )
2909 318 : CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_V_SOC_xyz_ikp, 1, s)
2910 318 : CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_V_SOC_xyz_ikp, s, 1)
2911 : CASE (2)
2912 : ! add V^SOC_y * σ_y for σ_y = ( (0,-i) (i,0) )
2913 318 : CALL cp_cfm_scale(gaussi, cfm_V_SOC_xyz_ikp)
2914 318 : CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_V_SOC_xyz_ikp, 1, s)
2915 318 : CALL cp_cfm_scale(-z_one, cfm_V_SOC_xyz_ikp)
2916 318 : CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_V_SOC_xyz_ikp, s, 1)
2917 : CASE (3)
2918 : ! add V^SOC_z * σ_z for σ_z = ( (1,0) (0,1) )
2919 318 : CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_V_SOC_xyz_ikp, 1, 1)
2920 318 : CALL cp_cfm_scale(-z_one, cfm_V_SOC_xyz_ikp)
2921 1272 : CALL add_cfm_submat(bs_env%cfm_SOC_spinor_ao(ikp), cfm_V_SOC_xyz_ikp, s, s)
2922 : END SELECT
2923 :
2924 : END DO
2925 :
2926 : END DO ! ikp
2927 :
2928 16 : CALL cp_cfm_release(cfm_V_SOC_xyz_ikp)
2929 :
2930 16 : CALL timestop(handle)
2931 :
2932 16 : END SUBROUTINE H_KS_spinor_kp
2933 :
2934 : ! **************************************************************************************************
2935 : !> \brief ...
2936 : !> \param cfm_array ...
2937 : !> \param cfm_template ...
2938 : !> \param n ...
2939 : ! **************************************************************************************************
2940 16 : SUBROUTINE alloc_cfm_double_array_1d(cfm_array, cfm_template, n)
2941 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_array
2942 : TYPE(cp_cfm_type) :: cfm_template
2943 : INTEGER :: n
2944 :
2945 : CHARACTER(LEN=*), PARAMETER :: routineN = 'alloc_cfm_double_array_1d'
2946 :
2947 : INTEGER :: handle, i
2948 :
2949 16 : CALL timeset(routineN, handle)
2950 :
2951 366 : ALLOCATE (cfm_array(n))
2952 334 : DO i = 1, n
2953 318 : CALL create_cfm_double(cfm_array(i), cfm_orig=cfm_template)
2954 334 : CALL cp_cfm_set_all(cfm_array(i), z_zero)
2955 : END DO
2956 :
2957 16 : CALL timestop(handle)
2958 :
2959 16 : END SUBROUTINE alloc_cfm_double_array_1d
2960 :
2961 : ! **************************************************************************************************
2962 : !> \brief ...
2963 : !> \param bs_env ...
2964 : ! **************************************************************************************************
2965 104 : SUBROUTINE get_all_VBM_CBM_bandgaps(bs_env)
2966 :
2967 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2968 :
2969 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_all_VBM_CBM_bandgaps'
2970 :
2971 : INTEGER :: handle
2972 :
2973 104 : CALL timeset(routineN, handle)
2974 :
2975 104 : CALL get_VBM_CBM_bandgaps(bs_env%band_edges_scf, bs_env%eigenval_scf, bs_env)
2976 104 : CALL get_VBM_CBM_bandgaps(bs_env%band_edges_G0W0, bs_env%eigenval_G0W0, bs_env)
2977 104 : CALL get_VBM_CBM_bandgaps(bs_env%band_edges_HF, bs_env%eigenval_HF, bs_env)
2978 :
2979 104 : CALL check_qp_gap_sanity(bs_env)
2980 :
2981 104 : CALL timestop(handle)
2982 :
2983 104 : END SUBROUTINE get_all_VBM_CBM_bandgaps
2984 :
2985 : ! **************************************************************************************************
2986 : !> \brief Warn if the G0W0 fundamental band gap is inverted or implausibly large, i.e. if the
2987 : !> quasiparticle solve has produced a spectrum that cannot be physical.
2988 : !> \param bs_env ...
2989 : ! **************************************************************************************************
2990 104 : SUBROUTINE check_qp_gap_sanity(bs_env)
2991 :
2992 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2993 :
2994 : REAL(KIND=dp) :: gap, gap_scf
2995 :
2996 104 : gap = bs_env%band_edges_G0W0%IDBG
2997 104 : gap_scf = bs_env%band_edges_scf%IDBG
2998 :
2999 : ! requiring a healthy SCF gap keeps the inversion test from firing on a genuine metal
3000 104 : IF (gap < -eps_qp_gap .AND. gap_scf > eps_qp_gap) THEN
3001 : CALL cp_warn(__LOCATION__, &
3002 : "G0W0 band gap is negative ("// &
3003 : TRIM(ADJUSTL(cp_to_string(gap*evolt, '(F12.3)')))//" eV): the quasiparticle "// &
3004 0 : "spectrum is inverted. Check numerical parameters.")
3005 104 : ELSE IF (ABS(gap) > max_qp_gap) THEN
3006 : CALL cp_warn(__LOCATION__, &
3007 : "G0W0 band gap is implausibly large ("// &
3008 : TRIM(ADJUSTL(cp_to_string(gap*evolt, '(F12.3)')))//" eV): the quasiparticle "// &
3009 0 : "solve has likely diverged. Check numerical parameters.")
3010 : END IF
3011 :
3012 104 : END SUBROUTINE check_qp_gap_sanity
3013 :
3014 : ! **************************************************************************************************
3015 : !> \brief ...
3016 : !> \param band_edges ...
3017 : !> \param ev ...
3018 : !> \param bs_env ...
3019 : ! **************************************************************************************************
3020 330 : SUBROUTINE get_VBM_CBM_bandgaps(band_edges, ev, bs_env)
3021 : TYPE(band_edges_type) :: band_edges
3022 : REAL(KIND=dp), DIMENSION(:, :, :) :: ev
3023 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3024 :
3025 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_VBM_CBM_bandgaps'
3026 :
3027 : INTEGER :: handle, homo, homo_1, homo_2, ikp, &
3028 : ispin, lumo, lumo_1, lumo_2, n_mo
3029 : REAL(KIND=dp) :: E_DBG_at_ikp
3030 :
3031 330 : CALL timeset(routineN, handle)
3032 :
3033 330 : n_mo = bs_env%n_ao
3034 :
3035 330 : band_edges%DBG = 1000.0_dp
3036 :
3037 618 : SELECT CASE (bs_env%n_spin)
3038 : CASE (1)
3039 288 : homo = bs_env%n_occ(1)
3040 288 : lumo = homo + 1
3041 7490 : band_edges%VBM = MAXVAL(ev(1:homo, :, 1))
3042 14850 : band_edges%CBM = MINVAL(ev(homo + 1:n_mo, :, 1))
3043 : CASE (2)
3044 42 : homo_1 = bs_env%n_occ(1)
3045 42 : lumo_1 = homo_1 + 1
3046 42 : homo_2 = bs_env%n_occ(2)
3047 42 : lumo_2 = homo_2 + 1
3048 462 : band_edges%VBM = MAX(MAXVAL(ev(1:homo_1, :, 1)), MAXVAL(ev(1:homo_2, :, 2)))
3049 966 : band_edges%CBM = MIN(MINVAL(ev(homo_1 + 1:n_mo, :, 1)), MINVAL(ev(homo_2 + 1:n_mo, :, 2)))
3050 : CASE DEFAULT
3051 330 : CPABORT("Error with number of spins.")
3052 : END SELECT
3053 :
3054 330 : band_edges%IDBG = band_edges%CBM - band_edges%VBM
3055 :
3056 702 : DO ispin = 1, bs_env%n_spin
3057 :
3058 372 : homo = bs_env%n_occ(ispin)
3059 :
3060 2350 : DO ikp = 1, bs_env%nkp_bs_and_DOS
3061 :
3062 21460 : E_DBG_at_ikp = -MAXVAL(ev(1:homo, ikp, ispin)) + MINVAL(ev(homo + 1:n_mo, ikp, ispin))
3063 :
3064 2020 : IF (E_DBG_at_ikp < band_edges%DBG) band_edges%DBG = E_DBG_at_ikp
3065 :
3066 : END DO
3067 :
3068 : END DO
3069 :
3070 330 : CALL timestop(handle)
3071 :
3072 330 : END SUBROUTINE get_VBM_CBM_bandgaps
3073 :
3074 : END MODULE post_scf_bandstructure_utils
|
