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