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