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