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