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