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