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