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