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