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