(git:374b731)
Loading...
Searching...
No Matches
post_scf_bandstructure_utils.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-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! **************************************************************************************************
13MODULE post_scf_bandstructure_utils
14 USE atomic_kind_types, ONLY: atomic_kind_type,&
15 get_atomic_kind,&
16 get_atomic_kind_set
17 USE cell_types, ONLY: cell_type,&
18 get_cell,&
19 pbc
20 USE cp_blacs_env, ONLY: cp_blacs_env_type
21 USE cp_cfm_basic_linalg, ONLY: cp_cfm_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
112CONTAINS
113
114! **************************************************************************************************
115!> \brief ...
116!> \param qs_env ...
117!> \param bs_env ...
118!> \param post_scf_bandstructure_section ...
119! **************************************************************************************************
120 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 CALL timeset(routinen, handle)
130
131 ALLOCATE (bs_env)
132
133 CALL print_header(bs_env)
134
135 CALL read_bandstructure_input_parameters(bs_env, post_scf_bandstructure_section)
136
137 CALL get_parameters_from_qs_env(qs_env, bs_env)
138
139 CALL set_heuristic_parameters(bs_env)
140
141 CALL setup_kpoints_dos(qs_env, bs_env, bs_env%kpoints_DOS)
142
143 CALL allocate_and_fill_fm_ks_fm_s(qs_env, bs_env)
144
145 CALL diagonalize_ks_matrix(bs_env)
146
147 CALL check_positive_definite_overlap_mat(bs_env, qs_env)
148
149 IF (bs_env%do_bs) THEN
150 CALL setup_kpoints_bandstructure(bs_env, bs_env%kpoints_bandstructure)
151 CALL setup_primitive_cell_for_bandstructure(qs_env, bs_env)
152 END IF
153
154 CALL timestop(handle)
155
156 END SUBROUTINE create_and_init_bs_env
157
158! **************************************************************************************************
159!> \brief ...
160!> \param bs_env ...
161!> \param bs_sec ...
162! **************************************************************************************************
163 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 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 CALL timeset(routinen, handle)
176
177 NULLIFY (gw_sec)
178 gw_sec => section_vals_get_subs_vals(bs_sec, "GW")
179 CALL section_vals_get(gw_sec, explicit=bs_env%do_gw)
180
181 NULLIFY (soc_sec)
182 soc_sec => section_vals_get_subs_vals(bs_sec, "SOC")
183 CALL section_vals_get(soc_sec, explicit=bs_env%do_soc)
184
185 CALL section_vals_val_get(soc_sec, "ENERGY_WINDOW", r_val=bs_env%energy_window_soc)
186
187 CALL section_vals_val_get(bs_sec, "DOS%KPOINTS", i_vals=bs_env%nkp_grid_DOS_input)
188 CALL section_vals_val_get(bs_sec, "DOS%ENERGY_WINDOW", r_val=bs_env%energy_window_DOS)
189 CALL section_vals_val_get(bs_sec, "DOS%ENERGY_STEP", r_val=bs_env%energy_step_DOS)
190 CALL section_vals_val_get(bs_sec, "DOS%BROADENING", r_val=bs_env%broadening_DOS)
191
192 NULLIFY (ldos_sec)
193 ldos_sec => section_vals_get_subs_vals(bs_sec, "DOS%LDOS")
194 CALL section_vals_get(ldos_sec, explicit=bs_env%do_ldos)
195
196 CALL section_vals_val_get(ldos_sec, "INTEGRATION", i_val=bs_env%int_ldos_xyz)
197 CALL section_vals_val_get(ldos_sec, "BIN_MESH", i_vals=bs_env%bin_mesh)
198
199 NULLIFY (kp_bs_sec)
200 kp_bs_sec => section_vals_get_subs_vals(bs_sec, "BANDSTRUCTURE_PATH")
201 CALL section_vals_get(kp_bs_sec, explicit=bs_env%do_bs)
202 CALL section_vals_val_get(kp_bs_sec, "NPOINTS", i_val=bs_env%input_kp_bs_npoints)
203 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 ALLOCATE (bs_env%xkp_special(3, bs_env%input_kp_bs_n_sp_pts))
207 ALLOCATE (bs_env%kp_special_name(bs_env%input_kp_bs_n_sp_pts))
208 DO ikp = 1, bs_env%input_kp_bs_n_sp_pts
209 CALL section_vals_val_get(kp_bs_sec, "SPECIAL_POINT", i_rep_val=ikp, c_vals=string_ptr)
210 cpassert(SIZE(string_ptr(:), 1) == 4)
211 bs_env%kp_special_name(ikp) = string_ptr(1)
212 DO i = 1, 3
213 CALL read_float_object(string_ptr(i + 1), bs_env%xkp_special(i, ikp), error_msg)
214 IF (len_trim(error_msg) > 0) cpabort(trim(error_msg))
215 END DO
216 END DO
217
218 CALL timestop(handle)
219
220 END SUBROUTINE read_bandstructure_input_parameters
221
222! **************************************************************************************************
223!> \brief ...
224!> \param bs_env ...
225! **************************************************************************************************
226 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 CALL timeset(routinen, handle)
235
236 bs_env%unit_nr = cp_logger_get_default_io_unit()
237
238 u = bs_env%unit_nr
239
240 IF (u > 0) THEN
241 WRITE (u, *) ' '
242 WRITE (u, '(T2,2A)') '-------------------------------------------------', &
243 '------------------------------'
244 WRITE (u, '(T2,2A)') '- ', &
245 ' -'
246 WRITE (u, '(T2,2A)') '- BANDSTRUCTURE CALCULATION', &
247 ' -'
248 WRITE (u, '(T2,2A)') '- ', &
249 ' -'
250 WRITE (u, '(T2,2A)') '--------------------------------------------------', &
251 '-----------------------------'
252 WRITE (u, '(T2,A)') ' '
253 END IF
254
255 CALL timestop(handle)
256
257 END SUBROUTINE print_header
258
259! **************************************************************************************************
260!> \brief ...
261!> \param qs_env ...
262!> \param bs_env ...
263!> \param kpoints ...
264! **************************************************************************************************
265 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 CALL timeset(routinen, handle)
277
278 ! routine adapted from mp2_integrals.F
279 NULLIFY (kpoints)
280 CALL kpoint_create(kpoints)
281
282 kpoints%kp_scheme = "GENERAL"
283
284 periodic(1:3) = bs_env%periodic(1:3)
285
286 DO i_dim = 1, 3
287
288 cpassert(periodic(i_dim) == 0 .OR. periodic(i_dim) == 1)
289
290 IF (bs_env%nkp_grid_DOS_input(i_dim) < 0) THEN
291 IF (periodic(i_dim) == 1) nkp_grid(i_dim) = 2
292 IF (periodic(i_dim) == 0) nkp_grid(i_dim) = 1
293 ELSE
294 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 IF (nkp_grid(1) > 1) THEN
301 nkp = (nkp_grid(1) + 1)/2*nkp_grid(2)*nkp_grid(3)
302 ELSE IF (nkp_grid(2) > 1) THEN
303 nkp = nkp_grid(1)*(nkp_grid(2) + 1)/2*nkp_grid(3)
304 ELSE IF (nkp_grid(3) > 1) THEN
305 nkp = nkp_grid(1)*nkp_grid(2)*(nkp_grid(3) + 1)/2
306 ELSE
307 nkp = 1
308 END IF
309
310 kpoints%nkp_grid(1:3) = nkp_grid(1:3)
311 kpoints%nkp = nkp
312
313 bs_env%nkp_DOS = nkp
314
315 ALLOCATE (kpoints%xkp(3, nkp), kpoints%wkp(nkp))
316 kpoints%wkp(1:nkp) = 1.0_dp/real(nkp, kind=dp)
317
318 CALL compute_xkp(kpoints%xkp, 1, nkp, nkp_grid)
319
320 CALL kpoint_init_cell_index_simple(kpoints, qs_env)
321
322 u = bs_env%unit_nr
323
324 IF (u > 0) THEN
325 WRITE (unit=u, fmt="(T2,1A,T69,3I4)") "K-point mesh for the density of states (DOS)", &
326 nkp_grid(1:3)
327 END IF
328
329 CALL timestop(handle)
330
331 END SUBROUTINE setup_kpoints_dos
332
333! **************************************************************************************************
334!> \brief ...
335!> \param bs_env ...
336!> \param kpoints ...
337! **************************************************************************************************
338 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 CALL timeset(routinen, handle)
349
350 ! routine adapted from mp2_integrals.F
351 NULLIFY (kpoints)
352 CALL kpoint_create(kpoints)
353
354 n_special_kp = bs_env%input_kp_bs_n_sp_pts
355 n_kp_in_line = bs_env%input_kp_bs_npoints
356
357 nkp = n_kp_in_line*(n_special_kp - 1) + 1
358
359 IF (n_special_kp < 1) &
360 cpabort("Please specify special k-points in the Brillouin zone via SPECIAL_POINT.")
361 IF (n_kp_in_line < 1) &
362 cpabort("Please specify the number of k-points between special k-points.")
363
364 ALLOCATE (kpoints%xkp(3, nkp))
365
366 kpoints%nkp = nkp
367 kpoints%xkp(1:3, 1) = bs_env%xkp_special(1:3, 1)
368
369 bs_env%nkp_bs = nkp
370
371 ikk = 1
372
373 DO i_special_kp = 2, n_special_kp
374 DO i_kp_in_line = 1, n_kp_in_line
375 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 bs_env%xkp_special(1:3, i_special_kp - 1))
380 END DO
381 END DO
382
383 u = bs_env%unit_nr
384
385 IF (u > 0) THEN
386 WRITE (unit=u, fmt="(T2,1A,T77,I4)") "Number of special k-points for the bandstructure", &
387 n_special_kp
388 WRITE (unit=u, fmt="(T2,1A,T77,I4)") "Number of k-points for the bandstructure", nkp
389 END IF
390
391 CALL timestop(handle)
392
393 END SUBROUTINE setup_kpoints_bandstructure
394
395! **************************************************************************************************
396!> \brief ...
397!> \param bs_env ...
398! **************************************************************************************************
399 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 CALL timeset(routinen, handle)
408
409 ALLOCATE (bs_env%eigenval_scf_Gamma(bs_env%n_ao, bs_env%n_spin))
410
411 DO ispin = 1, bs_env%n_spin
412
413 ! use work matrices because the matrices are overwritten in cp_fm_geeig_canon
414 CALL cp_fm_to_fm(bs_env%fm_ks_Gamma(ispin), bs_env%fm_work_mo(1))
415 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 bs_env%eps_eigval_mat_s)
425
426 vbm = bs_env%eigenval_scf_Gamma(bs_env%n_occ(ispin), ispin)
427 cbm = bs_env%eigenval_scf_Gamma(bs_env%n_occ(ispin) + 1, ispin)
428
429 bs_env%band_edges_scf_Gamma(ispin)%VBM = vbm
430 bs_env%band_edges_scf_Gamma(ispin)%CBM = cbm
431 bs_env%e_fermi(ispin) = 0.5_dp*(vbm + cbm)
432
433 END DO
434
435 CALL timestop(handle)
436
437 END SUBROUTINE diagonalize_ks_matrix
438
439! **************************************************************************************************
440!> \brief ...
441!> \param bs_env ...
442!> \param qs_env ...
443! **************************************************************************************************
444 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 CALL timeset(routinen, handle)
454
455 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 ikp, qs_env, bs_env%kpoints_DOS, "ORB")
460
461 ! check whether S_µν(k_i) is positive definite
462 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 IF (info .NE. 0) THEN
467 u = bs_env%unit_nr
468
469 IF (u > 0) THEN
470 WRITE (unit=u, fmt="(T2,A)") ""
471 WRITE (unit=u, fmt="(T2,A)") "ERROR: The Cholesky decomposition "// &
472 "of the k-point overlap matrix failed. This is"
473 WRITE (unit=u, fmt="(T2,A)") "because the algorithm is "// &
474 "only correct in the limit of large cells. The cell of "
475 WRITE (unit=u, fmt="(T2,A)") "the calculation is too small. "// &
476 "Use MULTIPLE_UNIT_CELL to create a larger cell "
477 WRITE (unit=u, fmt="(T2,A)") "and to prevent this error."
478 END IF
479
480 CALL bs_env%para_env%sync()
481 cpabort("Please see information on the error above.")
482
483 END IF ! Cholesky decomposition failed
484
485 END DO ! ikp
486
487 CALL cp_cfm_release(cfm_s_ikp)
488
489 CALL timestop(handle)
490
491 END SUBROUTINE check_positive_definite_overlap_mat
492
493! **************************************************************************************************
494!> \brief ...
495!> \param qs_env ...
496!> \param bs_env ...
497! **************************************************************************************************
498 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 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
510 TYPE(mp_para_env_type), POINTER :: para_env
511 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
512 TYPE(scf_control_type), POINTER :: scf_control
513
514 CALL timeset(routinen, handle)
515
516 CALL get_qs_env(qs_env, &
517 dft_control=dft_control, &
518 scf_control=scf_control, &
519 mos=mos)
520
521 bs_env%n_spin = dft_control%nspins
522 IF (bs_env%n_spin == 1) bs_env%spin_degeneracy = 2.0_dp
523 IF (bs_env%n_spin == 2) bs_env%spin_degeneracy = 1.0_dp
524
525 CALL get_mo_set(mo_set=mos(1), nao=n_ao, homo=homo)
526 bs_env%n_ao = n_ao
527 bs_env%n_occ(1:2) = homo
528 bs_env%n_vir(1:2) = n_ao - homo
529
530 IF (bs_env%n_spin == 2) THEN
531 CALL get_mo_set(mo_set=mos(2), homo=homo)
532 bs_env%n_occ(2) = homo
533 bs_env%n_vir(2) = n_ao - homo
534 END IF
535
536 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 CALL get_qs_env(qs_env, para_env=para_env)
540 color_sub = 0
541 ALLOCATE (bs_env%para_env)
542 CALL bs_env%para_env%from_split(para_env, color_sub)
543
544 CALL get_qs_env(qs_env, particle_set=particle_set)
545
546 n_atom = SIZE(particle_set)
547 bs_env%n_atom = n_atom
548
549 CALL get_qs_env(qs_env=qs_env, cell=cell)
550 CALL get_cell(cell=cell, periodic=periodic, h=hmat)
551 bs_env%periodic(1:3) = periodic(1:3)
552 bs_env%hmat(1:3, 1:3) = hmat
553
554 u = bs_env%unit_nr
555
556 IF (u > 0) THEN
557 WRITE (unit=u, fmt="(T2,2A,T73,I8)") "Number of occupied molecular orbitals (MOs) ", &
558 "= Number of occupied bands", homo
559 WRITE (unit=u, fmt="(T2,2A,T73,I8)") "Number of unoccupied (= virtual) MOs ", &
560 "= Number of unoccupied bands", n_ao - homo
561 WRITE (unit=u, fmt="(T2,A,T73,I8)") "Number of Gaussian basis functions for MOs", n_ao
562 END IF
563
564 CALL timestop(handle)
565
566 END SUBROUTINE get_parameters_from_qs_env
567
568! **************************************************************************************************
569!> \brief ...
570!> \param bs_env ...
571! **************************************************************************************************
572 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 CALL timeset(routinen, handle)
580
581 bs_env%n_bins_max_for_printing = 5000
582
583 CALL timestop(handle)
584
585 END SUBROUTINE set_heuristic_parameters
586
587! **************************************************************************************************
588!> \brief ...
589!> \param qs_env ...
590!> \param bs_env ...
591! **************************************************************************************************
592 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
602 TYPE(mp_para_env_type), POINTER :: para_env
603
604 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 matrix_s=matrix_s)
611
612 NULLIFY (fm_struct)
613 CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=bs_env%n_ao, &
614 ncol_global=bs_env%n_ao, para_env=para_env)
615
616 DO i_work = 1, SIZE(bs_env%fm_work_mo)
617 CALL cp_fm_create(bs_env%fm_work_mo(i_work), fm_struct)
618 END DO
619
620 CALL cp_fm_create(bs_env%fm_s_Gamma, fm_struct)
621 CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, bs_env%fm_s_Gamma)
622
623 DO ispin = 1, bs_env%n_spin
624 CALL cp_fm_create(bs_env%fm_ks_Gamma(ispin), fm_struct)
625 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, bs_env%fm_ks_Gamma(ispin))
626 CALL cp_fm_create(bs_env%fm_mo_coeff_Gamma(ispin), fm_struct)
627 END DO
628
629 CALL cp_fm_struct_release(fm_struct)
630
631 NULLIFY (bs_env%mat_ao_ao%matrix)
632 ALLOCATE (bs_env%mat_ao_ao%matrix)
633 CALL dbcsr_create(bs_env%mat_ao_ao%matrix, template=matrix_s(1)%matrix, &
634 matrix_type=dbcsr_type_no_symmetry)
635
636 ALLOCATE (bs_env%eigenval_scf(bs_env%n_ao, bs_env%nkp_DOS, bs_env%n_spin))
637
638 CALL timestop(handle)
639
640 END SUBROUTINE allocate_and_fill_fm_ks_fm_s
641
642! **************************************************************************************************
643!> \brief ...
644!> \param qs_env ...
645!> \param bs_env ...
646! **************************************************************************************************
647 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 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 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
665 TYPE(cell_type), POINTER :: cell
666 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
667
668 CALL timeset(routinen, handle)
669
670 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, particle_set=particle_set, cell=cell)
671 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
672 CALL get_cell(cell=cell, h=hmat)
673
674 ! automatically check structure for smallest possible unit cell
675 n_max_check = 20
676 n_max_x = n_max_check*bs_env%periodic(1) + 1
677 n_max_y = n_max_check*bs_env%periodic(2) + 1
678 n_max_z = n_max_check*bs_env%periodic(3) + 1
679
680 ALLOCATE (valid_multiple_unit_cell(n_max_x, n_max_y, n_max_z))
681 valid_multiple_unit_cell(:, :, :) = .true.
682
683 n_atom = bs_env%n_atom
684
685 DO i_atom = 1, n_atom
686
687 IF (.NOT. modulo(i_atom, bs_env%para_env%num_pe) == bs_env%para_env%mepos) cycle
688
689 ra(1:3) = particle_set(i_atom)%r(1:3)
690
691 DO n_x = 1, n_max_x
692 DO n_y = 1, n_max_y
693 DO n_z = 1, n_max_z
694
695 i_atom_has_image_in_every_subcell = .true.
696
697 DO i_x_cell = 0, n_x - 1
698 DO j_y_cell = 0, n_y - 1
699 DO k_z_cell = 0, n_z - 1
700
701 i_atom_has_image_in_subcell_ijk = .false.
702
703 DO j_atom = 1, n_atom
704
705 IF (kind_of(i_atom) .NE. kind_of(j_atom)) cycle
706
707 IF (i_atom_has_image_in_subcell_ijk) cycle
708
709 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 REAL(k_z_cell, dp)/REAL(n_z, dp)/)
714
715 coord_sub_cell_ijk(1:3) = matmul(hmat, index_sub_cell_ijk)
716
717 ra_ijk(1:3) = ra(1:3) + coord_sub_cell_ijk(1:3)
718
719 rb(1:3) = pbc(particle_set(j_atom)%r(1:3), cell)
720
721 rab(1:3) = rb(1:3) - pbc(ra_ijk(1:3), cell)
722
723 dab = sqrt((rab(1))**2 + (rab(2))**2 + (rab(3))**2)
724
725 IF (dab < 1.0e-5) i_atom_has_image_in_subcell_ijk = .true.
726
727 END DO
728
729 IF (.NOT. i_atom_has_image_in_subcell_ijk) THEN
730 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 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 CALL mpi_and(valid_multiple_unit_cell, bs_env%para_env)
748
749 n_mult_unit_cell_x = 1
750 n_mult_unit_cell_y = 1
751 n_mult_unit_cell_z = 1
752
753 DO n_x = 1, n_max_x
754 DO n_y = 1, n_max_y
755 DO n_z = 1, n_max_z
756 IF (valid_multiple_unit_cell(n_x, n_y, n_z)) THEN
757 n_mult_unit_cell_x = max(n_mult_unit_cell_x, n_x)
758 n_mult_unit_cell_y = max(n_mult_unit_cell_y, n_y)
759 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 bs_env%multiple_unit_cell(1) = n_mult_unit_cell_x
766 bs_env%multiple_unit_cell(2) = n_mult_unit_cell_y
767 bs_env%multiple_unit_cell(3) = n_mult_unit_cell_z
768
769 IF (n_mult_unit_cell_x .NE. 1 .OR. &
770 n_mult_unit_cell_y .NE. 1 .OR. &
771 n_mult_unit_cell_z .NE. 1) THEN
772 bs_env%calculate_bandstructure_of_primitive_cell = .true.
773 ELSE
774 bs_env%calculate_bandstructure_of_primitive_cell = .false.
775 END IF
776
777 n_atom_in_primitive_cell = n_atom/n_mult_unit_cell_x/n_mult_unit_cell_y/n_mult_unit_cell_z
778 bs_env%n_atom_in_primitive_cell = n_atom_in_primitive_cell
779
780 n_primitive_cells = n_atom/n_atom_in_primitive_cell
781 bs_env%n_primitive_cells = n_primitive_cells
782
783 bs_env%hmat_primitive_cell(1, 1:3) = hmat(1, 1:3)/real(n_mult_unit_cell_x)
784 bs_env%hmat_primitive_cell(2, 1:3) = hmat(2, 1:3)/real(n_mult_unit_cell_y)
785 bs_env%hmat_primitive_cell(3, 1:3) = hmat(3, 1:3)/real(n_mult_unit_cell_z)
786
787 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 .AND. n_primitive_cells > 1
791
792 ALLOCATE (bs_env%atoms_i_primitive_cell(n_atom_in_primitive_cell))
793 bs_env%atoms_i_primitive_cell(:) = 0
794
795 ! just a small offset to avoid that atoms are precisely on egdes or faces
796 offset(1:3) = matmul(bs_env%hmat_primitive_cell, (/0.001_dp, 0.001_dp, 0.001_dp/))
797 center_primitive_cell(1:3) = pbc(particle_set(1)%r(1:3), cell) - offset(1:3)
798
799 index_j = 0
800 DO i_atom = 1, n_atom
801
802 rb(1:3) = pbc(particle_set(i_atom)%r(1:3), cell) - center_primitive_cell(1:3)
803
804 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 index_atom_i(3) > -0.5_dp .AND. index_atom_i(3) < 0.5_dp) THEN
809
810 index_j = index_j + 1
811 cpassert(index_j .LE. n_atom_in_primitive_cell)
812 bs_env%atoms_i_primitive_cell(index_j) = i_atom
813
814 END IF
815
816 END DO
817
818 ALLOCATE (bs_env%ref_atom_primitive_cell(n_atom))
819 ALLOCATE (bs_env%cell_of_i_atom(n_atom, 3))
820
821 DO i_atom = 1, n_atom
822
823 ra(1:3) = pbc(particle_set(i_atom)%r(1:3), cell)
824
825 DO j_atom_prim_cell = 1, n_atom_in_primitive_cell
826
827 j_atom = bs_env%atoms_i_primitive_cell(j_atom_prim_cell)
828
829 rb(1:3) = pbc(particle_set(j_atom)%r(1:3), cell)
830
831 DO i_x_cell = -n_mult_unit_cell_x/2, n_mult_unit_cell_x/2
832 DO j_y_cell = -n_mult_unit_cell_y/2, n_mult_unit_cell_y/2
833 DO k_z_cell = -n_mult_unit_cell_z/2, n_mult_unit_cell_z/2
834
835 index_ijk(1:3) = (/real(i_x_cell, dp), real(j_y_cell, dp), real(k_z_cell, dp)/)
836 coord_ijk(1:3) = matmul(bs_env%hmat_primitive_cell, index_ijk)
837
838 ra_ijk(1:3) = ra(1:3) + coord_ijk(1:3)
839
840 rab(1:3) = rb(1:3) - pbc(ra_ijk(1:3), cell)
841
842 dab = sqrt((rab(1))**2 + (rab(2))**2 + (rab(3))**2)
843
844 IF (dab < 1.0e-5) THEN
845 bs_env%ref_atom_primitive_cell(i_atom) = j_atom
846 bs_env%cell_of_i_atom(i_atom, 1) = i_x_cell
847 bs_env%cell_of_i_atom(i_atom, 2) = j_y_cell
848 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 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 'Detected a multiple unit cell (will be used for band structure) ', &
859 bs_env%multiple_unit_cell
860 WRITE (bs_env%unit_nr, '(T2,A,I28)') &
861 'Number of occupied bands in the primitive unit cell', &
862 bs_env%n_occ(1)/bs_env%n_primitive_cells
863 WRITE (bs_env%unit_nr, '(T2,A,I26)') &
864 'Number of unoccupied bands in the primitive unit cell', &
865 bs_env%n_vir(1)/bs_env%n_primitive_cells
866 END IF
867
868 CALL timestop(handle)
869
870 END SUBROUTINE setup_primitive_cell_for_bandstructure
871
872! **************************************************************************************************
873!> \brief ...
874!> \param logical_array_3d ...
875!> \param para_env ...
876! **************************************************************************************************
877 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 CALL timeset(routinen, handle)
887
888 n_1 = SIZE(logical_array_3d, 1)
889 n_2 = SIZE(logical_array_3d, 2)
890 n_3 = SIZE(logical_array_3d, 3)
891
892 ALLOCATE (integer_array_3d(n_1, n_2, n_3))
893 integer_array_3d(:, :, :) = 0
894
895 DO i = 1, n_1
896 DO j = 1, n_2
897 DO k = 1, n_3
898 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 CALL para_env%sync()
904 CALL para_env%sum(integer_array_3d)
905 CALL para_env%sync()
906
907 logical_array_3d(:, :, :) = .false.
908
909 DO i = 1, n_1
910 DO j = 1, n_2
911 DO k = 1, n_3
912 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 CALL timestop(handle)
918
919 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 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 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
945 REAL(kind=dp) :: arg
946
947 CALL timeset(routinen, handle)
948
949 n_ao = bs_env%n_ao
950 n_ao_primitive_cell = n_ao/bs_env%n_primitive_cells
951 n_atom = bs_env%n_atom
952
953 ALLOCATE (h_munu_k(n_ao_primitive_cell, n_ao_primitive_cell))
954 ALLOCATE (s_munu_k(n_ao_primitive_cell, n_ao_primitive_cell))
955 ALLOCATE (mo_coeff_k(n_ao_primitive_cell, n_ao_primitive_cell))
956 ALLOCATE (eigenvalues(n_ao_primitive_cell, bs_env%nkp_bs))
957
958 ALLOCATE (atom_from_bf(n_ao))
959 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 first_bf_from_atom)
962
963 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 first_bf_from_atom)
966
967 IF (bs_env%para_env%is_source()) THEN
968 CALL open_file(filename, unit_number=iunit, file_status="REPLACE", &
969 file_action="WRITE", file_position="APPEND")
970 ELSE
971 iunit = -1
972 END IF
973
974 IF (iunit > 0) THEN
975
976 WRITE (unit=iunit, fmt="(2(A,I0),A)") "# ", &
977 bs_env%input_kp_bs_n_sp_pts, " special points, ", bs_env%nkp_bs, " k-points"
978 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 "# Special point ", ip, bs_env%xkp_special(1:3, ip), &
981 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 col_indices=col_indices)
991
992 DO i_dim = 1, 3
993 min_max_cell(i_dim) = min(maxval(bs_env%cell_of_i_atom(:, i_dim)), &
994 maxval(-bs_env%cell_of_i_atom(:, i_dim)))
995 END DO
996
997 DO ikp = 1, bs_env%nkp_bs
998
999 h_munu_k = z_zero
1000 s_munu_k = z_zero
1001
1002 DO i_row = 1, nrow_local
1003 DO j_col = 1, ncol_local
1004
1005 row_global = row_indices(i_row)
1006 col_global = col_indices(j_col)
1007
1008 i_atom = atom_from_bf(row_global)
1009 j_atom = atom_from_bf(col_global)
1010
1011 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 IF (any(cell_atom_i(1:3) .NE. 0)) cycle
1016
1017 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 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 REAL(cell_atom_j(3), dp)*bs_env%kpoints_bandstructure%xkp(3, ikp))*twopi
1026
1027 ref_atom_j = bs_env%ref_atom_primitive_cell(j_atom)
1028
1029 i = row_global - first_bf_from_atom(i_atom) + first_bf_of_primit_atom(i_atom)
1030 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 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 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 CALL bs_env%para_env%sync()
1043 CALL bs_env%para_env%sum(h_munu_k)
1044 CALL bs_env%para_env%sum(s_munu_k)
1045 CALL bs_env%para_env%sync()
1046
1047 CALL complex_geeig(h_munu_k, s_munu_k, mo_coeff_k, eigenvalues(:, ikp))
1048
1049 IF (iunit > 0) THEN
1050
1051 WRITE (unit=iunit, fmt="(A,I0,T15,A,T24,3(1X,F14.8))") &
1052 "# Point ", ikp, ":", bs_env%kpoints_bandstructure%xkp(1:3, ikp)
1053 WRITE (unit=iunit, fmt="(A)") "# Band Energy [eV]"
1054 DO imo = 1, n_ao_primitive_cell
1055 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 IF (bs_env%para_env%is_source()) CALL close_file(unit_number=iunit)
1063
1064 CALL timestop(handle)
1065
1066 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 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 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 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 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1097
1098 CALL timeset(routinen, handle)
1099
1100 n_ao = bs_env%n_ao
1101 n_ao_primitive_cell = n_ao/bs_env%n_primitive_cells
1102 n_atom = bs_env%n_atom
1103 n_atom_primitive_cell = n_atom/bs_env%n_primitive_cells
1104
1105 ALLOCATE (h_munu_k(2*n_ao_primitive_cell, 2*n_ao_primitive_cell))
1106 ALLOCATE (s_munu_k(2*n_ao_primitive_cell, 2*n_ao_primitive_cell))
1107 ALLOCATE (s_munu_k_single(n_ao_primitive_cell, n_ao_primitive_cell))
1108 ALLOCATE (mo_coeff_k(2*n_ao_primitive_cell, 2*n_ao_primitive_cell))
1109 ALLOCATE (eigenvalues(2*n_ao_primitive_cell, bs_env%nkp_bs))
1110 ALLOCATE (s_dot_mo_coeff_up(n_ao_primitive_cell))
1111 ALLOCATE (s_dot_mo_coeff_down(n_ao_primitive_cell))
1112
1113 ALLOCATE (atom_from_bf(2*n_ao))
1114 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 first_bf_from_atom)
1117 atom_from_bf(n_ao + 1:2*n_ao) = atom_from_bf(1:n_ao) + n_atom
1118 first_bf_from_atom(n_atom + 1:2*n_atom) = first_bf_from_atom(1:n_atom) + n_ao
1119
1120 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 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 + n_ao_primitive_cell
1125
1126 IF (bs_env%para_env%is_source()) THEN
1127 CALL open_file(filename, unit_number=iunit, file_status="REPLACE", &
1128 file_action="WRITE", file_position="APPEND")
1129 ELSE
1130 iunit = -1
1131 END IF
1132
1133 IF (iunit > 0) THEN
1134
1135 WRITE (unit=iunit, fmt="(2(A,I0),A)") "# ", &
1136 bs_env%input_kp_bs_n_sp_pts, " special points, ", bs_env%nkp_bs, " k-points"
1137 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 "# Special point ", ip, bs_env%xkp_special(1:3, ip), &
1140 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 col_indices=col_indices)
1150
1151 DO i_dim = 1, 3
1152 min_max_cell(i_dim) = min(maxval(bs_env%cell_of_i_atom(:, i_dim)), &
1153 maxval(-bs_env%cell_of_i_atom(:, i_dim)))
1154 END DO
1155
1156 DO ikp = 1, bs_env%nkp_bs
1157
1158 h_munu_k = z_zero
1159 s_munu_k = z_zero
1160
1161 DO i_row = 1, nrow_local
1162 DO j_col = 1, ncol_local
1163
1164 row_global = row_indices(i_row)
1165 col_global = col_indices(j_col)
1166
1167 i_atom = atom_from_bf(row_global)
1168 j_atom = atom_from_bf(col_global)
1169
1170 IF (i_atom > n_atom) THEN
1171 i_atom_non_spinor = i_atom - n_atom
1172 ELSE
1173 i_atom_non_spinor = i_atom
1174 END IF
1175
1176 IF (j_atom > n_atom) THEN
1177 j_atom_non_spinor = j_atom - n_atom
1178 ELSE
1179 j_atom_non_spinor = j_atom
1180 END IF
1181
1182 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 IF (any(cell_atom_i(1:3) .NE. 0)) cycle
1187
1188 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 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 REAL(cell_atom_j(3), dp)*bs_env%kpoints_bandstructure%xkp(3, ikp))*twopi
1197
1198 IF (j_atom > n_atom) THEN
1199 ref_atom_j = bs_env%ref_atom_primitive_cell(j_atom_non_spinor) + n_atom
1200 ELSE
1201 ref_atom_j = bs_env%ref_atom_primitive_cell(j_atom)
1202 END IF
1203
1204 i = row_global - first_bf_from_atom(i_atom) + first_bf_of_primit_atom(i_atom)
1205 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 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 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 CALL bs_env%para_env%sync()
1219 CALL bs_env%para_env%sum(h_munu_k)
1220 CALL bs_env%para_env%sum(s_munu_k)
1221 CALL bs_env%para_env%sync()
1222
1223 CALL complex_geeig(h_munu_k, s_munu_k, mo_coeff_k, eigenvalues(:, ikp))
1224
1225 IF (iunit > 0) THEN
1226
1227 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 "# Point ", ikp, ":", bs_env%kpoints_bandstructure%xkp(1:3, ikp)
1231 WRITE (unit=iunit, fmt="(A)") ħ"# Band Energy [eV] <S_z> / (/2) "
1232 DO imo = 1, 2*n_ao_primitive_cell
1233 s_dot_mo_coeff_up(:) = matmul(s_munu_k_single, &
1234 mo_coeff_k(1:n_ao_primitive_cell, imo))
1235 s_dot_mo_coeff_down(:) = matmul(s_munu_k_single, &
1236 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 sum(conjg(mo_coeff_k(n_ao_primitive_cell + 1:, imo))*s_dot_mo_coeff_down)
1239 WRITE (unit=iunit, fmt="(T2,I7,1X,2F14.8)") imo, eigenvalues(imo, ikp)*evolt, &
1240 REAL(s_z, kind=dp)
1241 END DO
1242
1243 END IF
1244
1245 END DO ! ikp
1246
1247 IF (bs_env%para_env%is_source()) CALL close_file(unit_number=iunit)
1248
1249 CALL timestop(handle)
1250
1251 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 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 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: overlap_sqrt_inv, work_1, work_2
1266 INTEGER :: i, n
1267 LOGICAL :: check_size
1268
1269 n = SIZE(matrix, 1)
1270
1271 check_size = SIZE(matrix, 2) == n .AND. SIZE(overlap, 1) == n .AND. &
1272 SIZE(eigenvalues) == n .AND. &
1273 SIZE(eigenvectors, 1) == n .AND. SIZE(eigenvectors, 2) == n
1274 cpassert(check_size)
1275
1276 ALLOCATE (work_1(n, n), work_2(n, n), overlap_sqrt_inv(n, n))
1277
1278 CALL complex_diag(overlap, eigenvectors, eigenvalues)
1279
1280 work_1(:, :) = z_zero
1281 DO i = 1, n
1282 IF (eigenvalues(i) > 1.0e-5_dp) THEN
1283 work_1(i, i) = eigenvalues(i)**(-0.5_dp)
1284 END IF
1285 END DO
1286 work_2(:, :) = matmul(work_1, transpose(conjg(eigenvectors)))
1287 overlap_sqrt_inv(:, :) = matmul(eigenvectors, work_2)
1288
1289 work_1(:, :) = matmul(matrix, overlap_sqrt_inv)
1290 work_2(:, :) = matmul(overlap_sqrt_inv, work_1)
1291
1292 CALL complex_diag(work_2, eigenvectors, eigenvalues)
1293
1294 work_2(:, :) = matmul(overlap_sqrt_inv, eigenvectors)
1295 eigenvectors(:, :) = work_2(:, :)
1296
1297 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 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 CALL timeset(routinen, handle)
1316
1317 first_bf_of_primit_atom(:) = 1
1318
1319 n_atom = bs_env%n_atom
1320
1321 DO i_atom = 2, n_atom
1322 IF (any(bs_env%atoms_i_primitive_cell(:) == i_atom)) THEN
1323 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 + n_bf_of_atom_i
1326 END IF
1327 END DO
1328
1329 CALL timestop(handle)
1330
1331 END SUBROUTINE get_basis_function_index_of_primitive_atoms
1332
1333! **************************************************************************************************
1334!> \brief ...
1335!> \param qs_env ...
1336!> \param bs_env ...
1337! **************************************************************************************************
1338 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dos_g0w0, dos_g0w0_soc, dos_scf, dos_scf_soc, &
1351 eigenval, eigenval_spinor, eigenval_spinor_g0w0, eigenval_spinor_no_soc
1352 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: pdos_g0w0, pdos_g0w0_soc, pdos_scf, &
1353 pdos_scf_soc, proj_mo_on_kind
1354 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ldos_g0w0_2d, ldos_scf_2d, &
1355 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 TYPE(cp_cfm_type), DIMENSION(2) :: cfm_mos_ikp
1363
1364 CALL timeset(routinen, handle)
1365
1366 n_ao = bs_env%n_ao
1367
1368 energy_window_dos = bs_env%energy_window_DOS
1369 energy_step_dos = bs_env%energy_step_DOS
1370 broadening = bs_env%broadening_DOS
1371
1372 ! if we have done GW, we already have the band edges
1373 IF (bs_env%do_gw) THEN
1374 band_edges_scf = bs_env%band_edges_scf
1375 band_edges_scf_guess = band_edges_scf
1376 ELSE
1377
1378 IF (bs_env%n_spin == 1) THEN
1379 homo = bs_env%n_occ(1)
1380 band_edges_scf_guess%VBM = bs_env%eigenval_scf_Gamma(homo, 1)
1381 band_edges_scf_guess%CBM = bs_env%eigenval_scf_Gamma(homo + 1, 1)
1382 ELSE
1383 homo_1 = bs_env%n_occ(1)
1384 homo_2 = bs_env%n_occ(2)
1385 band_edges_scf_guess%VBM = max(bs_env%eigenval_scf_Gamma(homo_1, 1), &
1386 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 bs_env%eigenval_scf_Gamma(homo_2 + 1, 2))
1389 END IF
1390
1391 ! initialization
1392 band_edges_scf%VBM = -1000.0_dp
1393 band_edges_scf%CBM = 1000.0_dp
1394 band_edges_scf%DBG = 1000.0_dp
1395 END IF
1396
1397 e_min = band_edges_scf_guess%VBM - 0.5_dp*energy_window_dos
1398 e_max = band_edges_scf_guess%CBM + 0.5_dp*energy_window_dos
1399
1400 IF (bs_env%do_gw) THEN
1401 band_edges_g0w0 = bs_env%band_edges_G0W0
1402 e_min_g0w0 = band_edges_g0w0%VBM - 0.5_dp*energy_window_dos
1403 e_max_g0w0 = band_edges_g0w0%CBM + 0.5_dp*energy_window_dos
1404 e_min = min(e_min, e_min_g0w0)
1405 e_max = max(e_max, e_max_g0w0)
1406 END IF
1407
1408 e_total_window = e_max - e_min
1409
1410 n_e = int(e_total_window/energy_step_dos)
1411
1412 ALLOCATE (dos_scf(n_e))
1413 dos_scf(:) = 0.0_dp
1414 ALLOCATE (dos_scf_soc(n_e))
1415 dos_scf_soc(:) = 0.0_dp
1416
1417 CALL get_qs_env(qs_env, nkind=nkind)
1418
1419 ALLOCATE (pdos_scf(n_e, nkind))
1420 pdos_scf(:, :) = 0.0_dp
1421 ALLOCATE (pdos_scf_soc(n_e, nkind))
1422 pdos_scf_soc(:, :) = 0.0_dp
1423
1424 ALLOCATE (proj_mo_on_kind(n_ao, nkind))
1425 proj_mo_on_kind(:, :) = 0.0_dp
1426
1427 ALLOCATE (eigenval(n_ao))
1428 ALLOCATE (eigenval_spinor(2*n_ao))
1429 ALLOCATE (eigenval_spinor_no_soc(2*n_ao))
1430 ALLOCATE (eigenval_spinor_g0w0(2*n_ao))
1431
1432 IF (bs_env%do_gw) THEN
1433
1434 ALLOCATE (dos_g0w0(n_e))
1435 dos_g0w0(:) = 0.0_dp
1436 ALLOCATE (dos_g0w0_soc(n_e))
1437 dos_g0w0_soc(:) = 0.0_dp
1438
1439 ALLOCATE (pdos_g0w0(n_e, nkind))
1440 pdos_g0w0(:, :) = 0.0_dp
1441 ALLOCATE (pdos_g0w0_soc(n_e, nkind))
1442 pdos_g0w0_soc(:, :) = 0.0_dp
1443
1444 END IF
1445
1446 CALL cp_cfm_create(cfm_mos_ikp(1), bs_env%fm_ks_Gamma(1)%matrix_struct)
1447 CALL cp_cfm_create(cfm_mos_ikp(2), bs_env%fm_ks_Gamma(1)%matrix_struct)
1448 CALL cp_cfm_create(cfm_work_ikp, bs_env%fm_ks_Gamma(1)%matrix_struct)
1449 CALL cp_cfm_create(cfm_s_ikp_copy, bs_env%fm_ks_Gamma(1)%matrix_struct)
1450
1451 IF (bs_env%do_soc) THEN
1452
1453 CALL cp_cfm_create(cfm_mos_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct)
1454 CALL cp_cfm_create(cfm_work_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct)
1455 CALL cp_cfm_create(cfm_s_ikp_spinor_copy, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct)
1456 CALL cp_cfm_create(cfm_ks_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct)
1457 CALL cp_cfm_create(cfm_soc_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct)
1458 CALL cp_cfm_create(cfm_s_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct)
1459 CALL cp_cfm_create(cfm_spinor_wf_ikp, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct)
1460
1461 homo_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin)
1462
1463 band_edges_scf_soc%VBM = -1000.0_dp
1464 band_edges_scf_soc%CBM = 1000.0_dp
1465 band_edges_scf_soc%DBG = 1000.0_dp
1466
1467 IF (bs_env%do_gw) THEN
1468 band_edges_g0w0_soc%VBM = -1000.0_dp
1469 band_edges_g0w0_soc%CBM = 1000.0_dp
1470 band_edges_g0w0_soc%DBG = 1000.0_dp
1471 END IF
1472
1473 END IF
1474
1475 IF (bs_env%do_ldos) THEN
1476 cpassert(bs_env%int_ldos_xyz == int_ldos_z)
1477 END IF
1478
1479 IF (bs_env%unit_nr > 0) THEN
1480 WRITE (bs_env%unit_nr, '(A)') ''
1481 END IF
1482
1483 DO ikp = 1, bs_env%kpoints_DOS%nkp
1484
1485 bs_env%t1 = m_walltime()
1486
1487 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 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 ikp, qs_env, bs_env%kpoints_DOS, "ORB")
1496 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 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 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 proj_mo_on_kind)
1509 IF (bs_env%do_gw) THEN
1510 CALL add_to_dos_pdos(dos_g0w0, pdos_g0w0, bs_env%eigenval_G0W0(:, ikp, ispin), &
1511 ikp, bs_env, n_e, e_min, proj_mo_on_kind)
1512 END IF
1513
1514 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 eigenval(:), band_edges_scf_guess)
1517
1518 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 bs_env%eigenval_G0W0(:, ikp, 1), band_edges_g0w0)
1521 END IF
1522
1523 END IF
1524
1525 homo = bs_env%n_occ(ispin)
1526
1527 band_edges_scf%VBM = max(band_edges_scf%VBM, eigenval(homo))
1528 band_edges_scf%CBM = min(band_edges_scf%CBM, eigenval(homo + 1))
1529 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 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 cfm_spinor_wf_ikp)
1540
1541 IF (.NOT. bs_env%do_gw) CALL write_soc_eigenvalues(eigenval_spinor, "SCF", ikp, bs_env)
1542
1543 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 eigenval_spinor, band_edges_scf_guess, .true., cfm_work_ikp)
1546 END IF
1547
1548 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 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 eigenval_spinor_g0w0)
1559
1560 END IF ! do_gw
1561
1562 END IF ! do_soc
1563
1564 IF (bs_env%unit_nr > 0) THEN
1565 WRITE (bs_env%unit_nr, '(T2,A,T43,I5,A,I3,A,F7.1,A)') &
1566 'Compute DOS, LDOS for k-point ', ikp, ' /', bs_env%kpoints_DOS%nkp, &
1567 ', Execution time', m_walltime() - bs_env%t1, ' s'
1568 END IF
1569
1570 END DO ! ikp_DOS
1571
1572 band_edges_scf%IDBG = band_edges_scf%CBM - band_edges_scf%VBM
1573 IF (bs_env%do_soc) THEN
1574 band_edges_scf_soc%IDBG = band_edges_scf_soc%CBM - band_edges_scf_soc%VBM
1575 IF (bs_env%do_gw) THEN
1576 band_edges_g0w0_soc%IDBG = band_edges_g0w0_soc%CBM - band_edges_g0w0_soc%VBM
1577 END IF
1578 END IF
1579
1580 CALL write_band_edges(band_edges_scf, "SCF", bs_env)
1581 CALL write_dos_pdos(dos_scf, pdos_scf, bs_env, qs_env, "SCF", e_min, band_edges_scf%VBM)
1582 IF (bs_env%do_ldos) THEN
1583 CALL print_ldos_main(ldos_scf_2d, bs_env, band_edges_scf, "SCF")
1584 END IF
1585
1586 IF (bs_env%do_soc) THEN
1587 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 e_min, band_edges_scf_soc%VBM)
1590 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 "SCF_SOC")
1595 END IF
1596 END IF
1597
1598 IF (bs_env%do_gw) THEN
1599 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 band_edges_g0w0%VBM)
1602 IF (bs_env%do_ldos) THEN
1603 CALL print_ldos_main(ldos_g0w0_2d, bs_env, band_edges_g0w0, "G0W0")
1604 END IF
1605 END IF
1606
1607 IF (bs_env%do_soc .AND. bs_env%do_gw) THEN
1608 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 band_edges_g0w0_soc%VBM)
1611 END IF
1612
1613 CALL cp_cfm_release(cfm_s_ikp)
1614 CALL cp_cfm_release(cfm_ks_ikp)
1615 CALL cp_cfm_release(cfm_mos_ikp(1))
1616 CALL cp_cfm_release(cfm_mos_ikp(2))
1617 CALL cp_cfm_release(cfm_work_ikp)
1618 CALL cp_cfm_release(cfm_s_ikp_copy)
1619
1620 CALL cp_cfm_release(cfm_s_ikp_spinor)
1621 CALL cp_cfm_release(cfm_ks_ikp_spinor)
1622 CALL cp_cfm_release(cfm_soc_ikp_spinor)
1623 CALL cp_cfm_release(cfm_mos_ikp_spinor)
1624 CALL cp_cfm_release(cfm_work_ikp_spinor)
1625 CALL cp_cfm_release(cfm_s_ikp_spinor_copy)
1626 CALL cp_cfm_release(cfm_spinor_wf_ikp)
1627
1628 CALL timestop(handle)
1629
1630 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 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 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 CALL timeset(routinen, handle)
1657
1658 n_e = SIZE(ldos_2d, 3)
1659
1660 energy_window = bs_env%energy_window_DOS
1661 energy_step = bs_env%energy_step_DOS
1662 e_min = band_edges%VBM - 0.5_dp*energy_window
1663
1664 bin_mesh(1:2) = bs_env%bin_mesh(1:2)
1665 do_xy_bins = (bin_mesh(1) > 0 .AND. bin_mesh(2) > 0)
1666
1667 i_x_start = lbound(ldos_2d, 1)
1668 i_x_end = ubound(ldos_2d, 1)
1669 i_y_start = lbound(ldos_2d, 2)
1670 i_y_end = ubound(ldos_2d, 2)
1671
1672 IF (do_xy_bins) THEN
1673 i_x_start_bin = 1
1674 i_x_end_bin = bin_mesh(1)
1675 i_y_start_bin = 1
1676 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 ALLOCATE (ldos_2d_bins(i_x_start_bin:i_x_end_bin, i_y_start_bin:i_y_end_bin, n_e))
1685 ldos_2d_bins(:, :, :) = 0.0_dp
1686
1687 IF (do_xy_bins) THEN
1688
1689 i_x_start_glob = i_x_start
1690 i_x_end_glob = i_x_end
1691 i_y_start_glob = i_y_start
1692 i_y_end_glob = i_y_end
1693
1694 CALL bs_env%para_env%min(i_x_start_glob)
1695 CALL bs_env%para_env%max(i_x_end_glob)
1696 CALL bs_env%para_env%min(i_y_start_glob)
1697 CALL bs_env%para_env%max(i_y_end_glob)
1698
1699 ALLOCATE (n_sum_for_bins(bin_mesh(1), bin_mesh(2)))
1700 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 DO i_x = i_x_start, i_x_end
1704 DO i_y = i_y_start, i_y_end
1705 i_x_bin = bin_mesh(1)*(i_x - i_x_start_glob)/(i_x_end_glob - i_x_start_glob + 1) + 1
1706 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 ldos_2d(i_x, i_y, :)
1709 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 CALL bs_env%para_env%sum(ldos_2d_bins)
1714 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 DO i_x_bin = 1, bin_mesh(1)
1718 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 REAL(n_sum_for_bins(i_x_bin, i_y_bin), kind=dp)
1721 END DO
1722 END DO
1723
1724 ELSE
1725
1726 ldos_2d_bins(:, :, :) = ldos_2d(:, :, :)
1727
1728 END IF
1729
1730 IF (bin_mesh(1)*bin_mesh(2) < bs_env%n_bins_max_for_printing) THEN
1731 CALL print_ldos_2d_bins(ldos_2d_bins, bs_env, e_min, scf_gw_soc)
1732 ELSE
1733 cpwarn("The number of bins for the LDOS is too large. Decrease BIN_MESH.")
1734 END IF
1735
1736 CALL timestop(handle)
1737
1738 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 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 CALL timeset(routinen, handle)
1765
1766 i_x_start = lbound(ldos_2d_bins, 1)
1767 i_x_end = ubound(ldos_2d_bins, 1)
1768 i_y_start = lbound(ldos_2d_bins, 2)
1769 i_y_end = ubound(ldos_2d_bins, 2)
1770 n_e = SIZE(ldos_2d_bins, 3)
1771
1772 n_x = i_x_end - i_x_start + 1
1773 n_y = i_y_end - i_y_start + 1
1774
1775 IF (bs_env%para_env%is_source()) THEN
1776
1777 DO i_x = i_x_start, i_x_end
1778 DO i_y = i_y_start, i_y_end
1779
1780 idx(1) = (real(i_x, kind=dp) - 0.5_dp)/real(n_x, kind=dp)
1781 idx(2) = (real(i_y, kind=dp) - 0.5_dp)/real(n_y, kind=dp)
1782 idx(3) = 0.0_dp
1783 coord(1:3) = matmul(bs_env%hmat, idx)
1784
1785 CALL get_print_format(coord(1), print_format_1)
1786 CALL get_print_format(coord(2), print_format_2)
1787
1788 print_format = "(3A,"//print_format_1//",A,"//print_format_2//",A)"
1789
1790 WRITE (fname, print_format) "LDOS_", scf_gw_soc, &
1791 "_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 file_action="WRITE")
1795
1796 WRITE (iunit, "(2A)") Ã…" Energy E (eV) average LDOS(x,y,E) (1/(eV*^2), ", &
1797 "integrated over z, averaged inside bin)"
1798
1799 DO i_e = 1, n_e
1800 energy = e_min + i_e*bs_env%energy_step_DOS
1801 WRITE (iunit, "(2F17.3)") energy*evolt, &
1802 ldos_2d_bins(i_x, i_y, i_e)* &
1803 bs_env%unit_ldos_int_z_inv_Ang2_eV
1804 END DO
1805
1806 CALL close_file(iunit)
1807
1808 END DO
1809 END DO
1810
1811 END IF
1812
1813 CALL timestop(handle)
1814
1815 END SUBROUTINE print_ldos_2d_bins
1816
1817! **************************************************************************************************
1818!> \brief ...
1819!> \param coord ...
1820!> \param print_format ...
1821! **************************************************************************************************
1822 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 CALL timeset(routinen, handle)
1831
1832 IF (coord < -10000/angstrom) THEN
1833 print_format = "F9.2"
1834 ELSE IF (coord < -1000/angstrom) THEN
1835 print_format = "F8.2"
1836 ELSE IF (coord < -100/angstrom) THEN
1837 print_format = "F7.2"
1838 ELSE IF (coord < -10/angstrom) THEN
1839 print_format = "F6.2"
1840 ELSE IF (coord < -1/angstrom) THEN
1841 print_format = "F5.2"
1842 ELSE IF (coord < 10/angstrom) THEN
1843 print_format = "F4.2"
1844 ELSE IF (coord < 100/angstrom) THEN
1845 print_format = "F5.2"
1846 ELSE IF (coord < 1000/angstrom) THEN
1847 print_format = "F6.2"
1848 ELSE IF (coord < 10000/angstrom) THEN
1849 print_format = "F7.2"
1850 ELSE
1851 print_format = "F8.2"
1852 END IF
1853
1854 CALL timestop(handle)
1855
1856 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 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 CALL timeset(routinen, handle)
1899
1900 n_ao = bs_env%n_ao
1901 homo_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin)
1902 n_e = SIZE(dos)
1903 nkind = SIZE(pdos, 2)
1904
1905 CALL cp_cfm_create(cfm_ks_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct)
1906 CALL cp_cfm_create(cfm_soc_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct)
1907 CALL cp_cfm_create(cfm_mos_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct)
1908 CALL cp_cfm_create(cfm_work_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct)
1909 CALL cp_cfm_create(cfm_eigenvec_ikp_spinor, bs_env%cfm_ks_spinor_ao_Gamma%matrix_struct)
1910
1911 ALLOCATE (eigenval_spinor_no_soc(2*n_ao))
1912 ALLOCATE (proj_mo_on_kind_spinor(2*n_ao, nkind))
1913 ! PDOS not yet implemented -> projection is just zero -> PDOS is zero
1914 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 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 CALL cp_cfm_set_all(cfm_mos_ikp_spinor, z_zero)
1927 CALL add_cfm_submat(cfm_mos_ikp_spinor, cfm_mos_ikp(1), 1, 1)
1928 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 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 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 eigenval_spinor_no_soc(1:n_ao) = eigenval_no_soc(1:n_ao, ikp, 1)
1944 eigenval_spinor_no_soc(n_ao + 1:) = eigenval_no_soc(1:n_ao, ikp, bs_env%n_spin)
1945 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 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 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 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 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 band_edges%VBM = max(band_edges%VBM, eigenval_spinor(homo_spinor))
1965 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 - 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 z_zero, cfm_spinor_wf_ikp)
1973
1974 CALL cp_cfm_release(cfm_ks_ikp_spinor)
1975 CALL cp_cfm_release(cfm_soc_ikp_spinor)
1976 CALL cp_cfm_release(cfm_work_ikp_spinor)
1977 CALL cp_cfm_release(cfm_eigenvec_ikp_spinor)
1978 CALL cp_cfm_release(cfm_mos_ikp_spinor)
1979
1980 CALL timestop(handle)
1981
1982 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 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 CALL timeset(routinen, handle)
2013
2014 energy_step_dos = bs_env%energy_step_DOS
2015 broadening = bs_env%broadening_DOS
2016
2017 n_primitive_cells = 1
2018 IF (bs_env%do_bs) n_primitive_cells = bs_env%n_primitive_cells
2019
2020 n_mo = SIZE(eigenval)
2021 nkind = SIZE(proj_mo_on_kind, 2)
2022
2023 ! normalize to the primitive cell and to closed-shell / open-shell
2024 wkp = bs_env%kpoints_DOS%wkp(ikp)/n_primitive_cells*bs_env%spin_degeneracy
2025 DO i_e = 1, n_e
2026 energy = e_min + i_e*energy_step_dos
2027 DO i_mo = 1, n_mo
2028 ! DOS
2029 dos(i_e) = dos(i_e) + wkp*gaussian(energy - eigenval(i_mo), broadening)
2030
2031 ! PDOS
2032 DO i_kind = 1, nkind
2033 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 gaussian(energy - eigenval(i_mo), broadening)
2037 END IF
2038 END DO
2039 END DO
2040 END DO
2041
2042 CALL timestop(handle)
2043
2044 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 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 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 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 CALL timeset(routinen, handle)
2090
2091 my_do_spinor = .false.
2092 IF (PRESENT(do_spinor)) my_do_spinor = do_spinor
2093
2094 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 nimages = dft_control%nimages
2098 dft_control%nimages = 1
2099
2100 energy_window = bs_env%energy_window_DOS
2101 energy_step = bs_env%energy_step_DOS
2102 broadening = bs_env%broadening_DOS
2103
2104 e_min = band_edges%VBM - 0.5_dp*energy_window
2105 e_max = band_edges%CBM + 0.5_dp*energy_window
2106 e_total_window = e_max - e_min
2107
2108 n_e = int(e_total_window/energy_step)
2109
2110 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
2111
2112 CALL auxbas_pw_pool%create_pw(ldos_3d)
2113 CALL auxbas_pw_pool%create_pw(rho_g)
2114
2115 i_x_start = lbound(ldos_3d%array, 1)
2116 i_x_end = ubound(ldos_3d%array, 1)
2117 i_y_start = lbound(ldos_3d%array, 2)
2118 i_y_end = ubound(ldos_3d%array, 2)
2119 i_z_start = lbound(ldos_3d%array, 3)
2120 i_z_end = ubound(ldos_3d%array, 3)
2121
2122 z_start_global = i_z_start
2123 z_end_global = i_z_end
2124
2125 CALL bs_env%para_env%min(z_start_global)
2126 CALL bs_env%para_env%max(z_end_global)
2127 n_z = z_end_global - z_start_global + 1
2128
2129 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 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 bs_env%unit_ldos_int_z_inv_Ang2_eV = bs_env%hmat(3, 3)/real(n_z, kind=dp)/evolt/angstrom**2
2133
2134 IF (ikp == 1) THEN
2135 ALLOCATE (ldos_2d(i_x_start:i_x_end, i_y_start:i_y_end, n_e))
2136 ldos_2d(:, :, :) = 0.0_dp
2137 END IF
2138
2139 CALL cp_cfm_create(cfm_work, cfm_mos_ikp%matrix_struct)
2140 CALL cp_cfm_create(cfm_weighted_dm_ikp, cfm_mos_ikp%matrix_struct)
2141 CALL cp_fm_create(fm_weighted_dm_mic, cfm_mos_ikp%matrix_struct)
2142 IF (my_do_spinor) THEN
2143 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 col_indices=col_indices)
2150
2151 NULLIFY (weighted_dm_mic)
2152 CALL dbcsr_allocate_matrix_set(weighted_dm_mic, 1)
2153 ALLOCATE (weighted_dm_mic(1)%matrix)
2154 CALL dbcsr_create(weighted_dm_mic(1)%matrix, template=bs_env%mat_ao_ao%matrix, &
2155 matrix_type=dbcsr_type_symmetric)
2156
2157 DO i_e = 1, n_e
2158
2159 energy = e_min + i_e*energy_step
2160
2161 is_any_weight_non_zero = .false.
2162
2163 DO j_col = 1, ncol_local
2164
2165 j_mo = col_indices(j_col)
2166
2167 IF (my_do_spinor) THEN
2168 spin_degeneracy = 1.0_dp
2169 ELSE
2170 spin_degeneracy = bs_env%spin_degeneracy
2171 END IF
2172
2173 weight = gaussian(energy - eigenval(j_mo), broadening)*spin_degeneracy
2174
2175 cfm_work%local_data(:, j_col) = cfm_mos_ikp%local_data(:, j_col)*weight
2176
2177 IF (weight > 1.0e-5_dp) is_any_weight_non_zero = .true.
2178
2179 END DO
2180
2181 CALL bs_env%para_env%sync()
2182 CALL bs_env%para_env%sum(is_any_weight_non_zero)
2183 CALL bs_env%para_env%sync()
2184
2185 ! cycle if there are no states at the energy i_E
2186 IF (is_any_weight_non_zero) THEN
2187
2188 CALL parallel_gemm('N', 'C', n_mo, n_mo, n_mo, z_one, &
2189 cfm_mos_ikp, cfm_work, z_zero, cfm_weighted_dm_ikp)
2190
2191 IF (my_do_spinor) THEN
2192
2193 ! contribution from up,up to fm_non_spinor
2194 CALL get_cfm_submat(cfm_non_spinor, cfm_weighted_dm_ikp, 1, 1)
2195 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 "ORB", bs_env%kpoints_DOS%wkp(ikp))
2199
2200 ! add contribution from down,down to fm_non_spinor
2201 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 "ORB", bs_env%kpoints_DOS%wkp(ikp))
2205 CALL copy_fm_to_dbcsr(fm_non_spinor, weighted_dm_mic(1)%matrix, &
2206 keep_sparsity=.false.)
2207 ELSE
2208 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 "ORB", bs_env%kpoints_DOS%wkp(ikp))
2212 CALL copy_fm_to_dbcsr(fm_weighted_dm_mic, weighted_dm_mic(1)%matrix, &
2213 keep_sparsity=.false.)
2214 END IF
2215
2216 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 ks_env=ks_env)
2222
2223 DO i_z = i_z_start, i_z_end
2224 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 dft_control%nimages = nimages
2233
2234 CALL auxbas_pw_pool%give_back_pw(ldos_3d)
2235 CALL auxbas_pw_pool%give_back_pw(rho_g)
2236
2237 CALL cp_cfm_release(cfm_work)
2238 CALL cp_cfm_release(cfm_weighted_dm_ikp)
2239
2240 CALL cp_fm_release(fm_weighted_dm_mic)
2241
2242 CALL dbcsr_deallocate_matrix_set(weighted_dm_mic)
2243
2244 IF (my_do_spinor) THEN
2245 CALL cp_fm_release(fm_non_spinor)
2246 END IF
2247
2248 CALL timestop(handle)
2249
2250 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 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 CALL timeset(routinen, handle)
2275
2276 CALL get_fname(fname, bs_env, ikp, scf_gw, soc=.true.)
2277
2278 IF (bs_env%para_env%is_source()) THEN
2279
2280 CALL open_file(trim(fname), unit_number=iunit, file_status="REPLACE", file_action="WRITE")
2281
2282 WRITE (iunit, "(A)") " "
2283 WRITE (iunit, "(A10,3F10.4)") "kpoint: ", bs_env%kpoints_DOS%xkp(:, ikp)
2284 WRITE (iunit, "(A)") " "
2285
2286 n_occ_spinor = bs_env%n_occ(1) + bs_env%n_occ(bs_env%n_spin)
2287
2288 DO i_mo = 1, SIZE(eigenval_spinor)
2289 IF (i_mo .LE. n_occ_spinor) occ_vir = 'occ'
2290 IF (i_mo > n_occ_spinor) occ_vir = 'vir'
2291 IF (PRESENT(eigenval_spinor_g0w0)) THEN
2292 ! SCF+SOC and G0W0+SOC eigenvalues
2293 WRITE (iunit, "(I5,3A,4F16.3,2F17.3)") i_mo, ' (', occ_vir, ') ', &
2294 eigenval_spinor(i_mo)*evolt, eigenval_spinor_g0w0(i_mo)*evolt
2295 ELSE
2296 ! SCF+SOC eigenvalues only
2297 WRITE (iunit, "(I5,3A,4F16.3,F17.3)") i_mo, ' (', occ_vir, ') ', &
2298 eigenval_spinor(i_mo)*evolt
2299 END IF
2300 END DO
2301
2302 CALL close_file(iunit)
2303
2304 END IF
2305
2306 CALL timestop(handle)
2307
2308 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 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 CALL timeset(routinen, handle)
2335
2336 my_soc = .false.
2337 IF (PRESENT(soc)) my_soc = soc
2338
2339 n_zeros = count_digits(bs_env%kpoints_DOS%nkp) - count_digits(ikp)
2340
2341 WRITE (digits_ikp, '(I1)') count_digits(ikp)
2342
2343 IF (bs_env%kpoints_DOS%nkp == 1) THEN
2344 ! for molecules
2345 core_name = trim(scf_gw)//"_eigenvalues"
2346 ELSE
2347 ! for periodic systems
2348 core_name = trim(scf_gw)//"_band_struct_kp"
2349 END IF
2350
2351 SELECT CASE (n_zeros)
2352 CASE (0)
2353 WRITE (fname, "(2A,I"//digits_ikp//")") trim(core_name), "_", ikp
2354 CASE (1)
2355 WRITE (fname, "(2A,I"//digits_ikp//")") trim(core_name), "_0", ikp
2356 CASE (2)
2357 WRITE (fname, "(2A,I"//digits_ikp//")") trim(core_name), "_00", ikp
2358 CASE (3)
2359 WRITE (fname, "(2A,I"//digits_ikp//")") trim(core_name), "_000", ikp
2360 CASE DEFAULT
2361 cpabort("Too many zeros.")
2362 END SELECT
2363
2364 ! for molecules, we don't have any k-points; overwrite fname
2365 IF (ikp == 1 .AND. bs_env%kpoints_DOS%nkp == 1) THEN
2366 WRITE (fname, "(A)") trim(core_name)
2367 END IF
2368
2369 IF (my_soc) THEN
2370 WRITE (fname, "(2A)") trim(fname), "_+_SOC"
2371 END IF
2372
2373 IF (bs_env%n_spin == 2 .AND. .NOT. my_soc) THEN
2374 cpassert(PRESENT(ispin))
2375 WRITE (fname, "(2A,I1)") trim(fname), "_spin_", ispin
2376 END IF
2377
2378 CALL timestop(handle)
2379
2380 END SUBROUTINE get_fname
2381
2382! **************************************************************************************************
2383!> \brief ...
2384!> \param int_number ...
2385!> \return ...
2386! **************************************************************************************************
2387 PURE FUNCTION count_digits(int_number)
2388
2389 INTEGER, INTENT(IN) :: int_number
2390 INTEGER :: count_digits
2391
2392 INTEGER :: digitcount, tempint
2393
2394 digitcount = 0
2395
2396 tempint = int_number
2397
2398 DO WHILE (tempint /= 0)
2399 tempint = tempint/10
2400 digitcount = digitcount + 1
2401 END DO
2402
2403 count_digits = digitcount
2404
2405 END FUNCTION count_digits
2406
2407! **************************************************************************************************
2408!> \brief ...
2409!> \param band_edges ...
2410!> \param scf_gw_soc ...
2411!> \param bs_env ...
2412! **************************************************************************************************
2413 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 CALL timeset(routinen, handle)
2425
2426 ! print format
2427 print_format = "(T2,2A,T61,F20.3)"
2428
2429 u = bs_env%unit_nr
2430 IF (u > 0) THEN
2431 WRITE (u, '(T2,A)') ''
2432 WRITE (u, print_format) scf_gw_soc, ' valence band maximum (eV):', band_edges%VBM*evolt
2433 WRITE (u, print_format) scf_gw_soc, ' conduction band minimum (eV):', band_edges%CBM*evolt
2434 WRITE (u, print_format) scf_gw_soc, ' indirect band gap (eV):', band_edges%IDBG*evolt
2435 WRITE (u, print_format) scf_gw_soc, ' direct band gap (eV):', band_edges%DBG*evolt
2436 END IF
2437
2438 CALL timestop(handle)
2439
2440 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 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2468
2469 CALL timeset(routinen, handle)
2470
2471 WRITE (fname, "(3A)") "DOS_PDOS_", scf_gw_soc, ".out"
2472
2473 n_e = SIZE(pdos, 1)
2474 nkind = SIZE(pdos, 2)
2475 CALL get_qs_env(qs_env, particle_set=particle_set)
2476
2477 IF (bs_env%para_env%is_source()) THEN
2478
2479 CALL open_file(trim(fname), unit_number=iunit, file_status="REPLACE", file_action="WRITE")
2480
2481 n_a = 2 + nkind
2482
2483 DO iatom = 1, bs_env%n_atom
2484 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
2485 kind_number=i_kind, name=atom_name)
2486 elements(i_kind) = atom_name
2487 END DO
2488
2489 WRITE (output_string, "(A,I1,A)") "(", n_a, "A)"
2490
2491 WRITE (iunit, trim(output_string)) "Energy-E_F (eV) DOS (1/eV) PDOS (1/eV) ", &
2492 " of atom type ", elements(1:nkind)
2493
2494 WRITE (output_string, "(A,I1,A)") "(", n_a, "F13.5)"
2495
2496 DO i_e = 1, n_e
2497 ! energy is relative to valence band maximum => - E_VBM
2498 energy = e_min + i_e*bs_env%energy_step_DOS - e_vbm
2499 WRITE (iunit, trim(output_string)) energy*evolt, dos(i_e)/evolt, pdos(i_e, :)/evolt
2500 END DO
2501
2502 CALL close_file(iunit)
2503
2504 END IF
2505
2506 CALL timestop(handle)
2507
2508 END SUBROUTINE write_dos_pdos
2509
2510! **************************************************************************************************
2511!> \brief ...
2512!> \param energy ...
2513!> \param broadening ...
2514!> \return ...
2515! **************************************************************************************************
2516 PURE FUNCTION gaussian(energy, broadening)
2517
2518 REAL(kind=dp), INTENT(IN) :: energy, broadening
2519 REAL(kind=dp) :: gaussian
2520
2521 IF (abs(energy) < 5*broadening) THEN
2522 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 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_bf, kind_of
2547 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2548 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 CALL timeset(routinen, handle)
2553
2554 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, nkind=nkind)
2555 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 col_indices=col_indices)
2563
2564 n_ao = qs_env%bs_env%n_ao
2565
2566 ALLOCATE (atom_from_bf(n_ao))
2567 CALL get_atom_index_from_basis_function_index(qs_env, atom_from_bf, n_ao, "ORB")
2568
2569 proj_mo_on_kind(:, :) = 0.0_dp
2570
2571 CALL cp_cfm_create(cfm_s_i_kind, cfm_s%matrix_struct)
2572 CALL cp_cfm_create(cfm_work, cfm_s%matrix_struct)
2573 CALL cp_cfm_create(cfm_proj, cfm_s%matrix_struct)
2574 CALL cp_fm_create(fm_proj_re, cfm_s%matrix_struct)
2575 CALL cp_fm_create(fm_proj_im, cfm_s%matrix_struct)
2576
2577 DO i_kind = 1, nkind
2578
2579 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 DO i_row = 1, nrow_local
2583 DO j_col = 1, ncol_local
2584
2585 i_global = row_indices(i_row)
2586
2587 IF (i_global .LE. n_ao) THEN
2588 i_atom = atom_from_bf(i_global)
2589 ELSE IF (i_global .LE. 2*n_ao) THEN
2590 i_atom = atom_from_bf(i_global - n_ao)
2591 ELSE
2592 cpabort("Wrong indices.")
2593 END IF
2594
2595 IF (i_kind .NE. kind_of(i_atom)) THEN
2596 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 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 cfm_mos, cfm_work, z_zero, cfm_proj)
2606
2607 CALL cp_cfm_to_fm(cfm_proj, fm_proj_re, fm_proj_im)
2608
2609 CALL cp_fm_get_diag(fm_proj_im, proj_mo_on_kind(:, i_kind))
2610 CALL cp_fm_get_diag(fm_proj_re, proj_mo_on_kind(:, i_kind))
2611
2612 END DO ! i_kind
2613
2614 CALL cp_cfm_release(cfm_s_i_kind)
2615 CALL cp_cfm_release(cfm_work)
2616 CALL cp_cfm_release(cfm_proj)
2617 CALL cp_fm_release(fm_proj_re)
2618 CALL cp_fm_release(fm_proj_im)
2619
2620 CALL timestop(handle)
2621
2622 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 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 CALL timeset(routinen, handle)
2652
2653 CALL cp_cfm_create(cfm_non_spinor_gamma, fm_struct_non_spinor)
2654 CALL cp_cfm_create(cfm_non_spinor_ikp, fm_struct_non_spinor)
2655 CALL cp_fm_create(fm_non_spinor_gamma_re, fm_struct_non_spinor)
2656 CALL cp_fm_create(fm_non_spinor_gamma_im, fm_struct_non_spinor)
2657
2658 CALL cp_cfm_get_info(cfm_non_spinor_gamma, nrow_global=n_ao)
2659
2660 CALL cp_cfm_set_all(cfm_spinor_ikp, z_zero)
2661
2662 DO i_block = 0, 1
2663 DO j_block = 0, 1
2664 i_offset = i_block*n_ao + 1
2665 j_offset = j_block*n_ao + 1
2666 CALL get_cfm_submat(cfm_non_spinor_gamma, cfm_spinor_gamma, i_offset, j_offset)
2667 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 ikp, qs_env, kpoints, basis_type)
2672 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 ikp, qs_env, kpoints, basis_type)
2677 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 CALL cp_cfm_release(cfm_non_spinor_gamma)
2683 CALL cp_cfm_release(cfm_non_spinor_ikp)
2684 CALL cp_fm_release(fm_non_spinor_gamma_re)
2685 CALL cp_fm_release(fm_non_spinor_gamma_im)
2686
2687 CALL timestop(handle)
2688
2689 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_bf
2713 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2714 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2722
2723 CALL timeset(routinen, handle)
2724
2725 IF (.NOT. ASSOCIATED(cfm_ikp%local_data)) THEN
2726 CALL cp_cfm_create(cfm_ikp, fm_gamma%matrix_struct)
2727 END IF
2728 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 col_indices=col_indices)
2735
2736 ! get number of basis functions (bf) for different basis sets
2737 IF (basis_type == "ORB") THEN
2738 n_bf = qs_env%bs_env%n_ao
2739 ELSE IF (basis_type == "RI_AUX") THEN
2740 n_bf = qs_env%bs_env%n_RI
2741 ELSE
2742 cpabort("Only ORB and RI_AUX basis implemented.")
2743 END IF
2744
2745 ALLOCATE (atom_from_bf(n_bf))
2746 CALL get_atom_index_from_basis_function_index(qs_env, atom_from_bf, n_bf, basis_type)
2747
2748 NULLIFY (cell, particle_set)
2749 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
2750 CALL get_cell(cell=cell, h=hmat)
2751
2752 index_to_cell => kpoints%index_to_cell
2753
2754 num_cells = SIZE(index_to_cell, 2)
2755 i_atom_old = 0
2756 j_atom_old = 0
2757
2758 DO i_row = 1, nrow_local
2759 DO j_col = 1, ncol_local
2760
2761 row_global = row_indices(i_row)
2762 col_global = col_indices(j_col)
2763
2764 i_atom = atom_from_bf(row_global)
2765 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 IF (i_atom .NE. i_atom_old .OR. j_atom .NE. j_atom_old) THEN
2769 DO i_cell = 1, num_cells
2770
2771 ! only check nearest neigbors
2772 IF (any(abs(index_to_cell(1:3, i_cell)) > 1)) cycle
2773
2774 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 (pbc(particle_set(j_atom)%r(1:3), cell) + cell_vector(1:3))
2778 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 i_cell_is_the_minimum_image_cell = .true.
2782 DO j_cell = 1, num_cells
2783 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 (pbc(particle_set(j_atom)%r(1:3), cell) + cell_vector_j(1:3))
2786 abs_rab_cell_j = sqrt(rab_cell_j(1)**2 + rab_cell_j(2)**2 + rab_cell_j(3)**2)
2787
2788 IF (abs_rab_cell_i > abs_rab_cell_j + 1.0e-6_dp) THEN
2789 i_cell_is_the_minimum_image_cell = .false.
2790 END IF
2791 END DO
2792
2793 IF (i_cell_is_the_minimum_image_cell) THEN
2794 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 REAL(index_to_cell(3, i_mic_cell), dp)*kpoints%xkp(3, ikp)
2803
2804 i = i_row
2805 j = j_col
2806
2807 cfm_ikp%local_data(i, j) = cos(twopi*arg)*fm_gamma%local_data(i, j)*z_one + &
2808 sin(twopi*arg)*fm_gamma%local_data(i, j)*gaussi
2809
2810 j_atom_old = j_atom
2811 i_atom_old = i_atom
2812
2813 END DO ! j_col
2814 END DO ! i_row
2815
2816 CALL timestop(handle)
2817
2818 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_bf_index
2848 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2849 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 REAL(kind=dp), DIMENSION(:), POINTER :: wkp
2854 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
2855 TYPE(cell_type), POINTER :: cell
2856 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2857
2858 CALL timeset(routinen, handle)
2859
2860 ! get number of basis functions (bf) for different basis sets
2861 IF (basis_type == "ORB") THEN
2862 n_bf = qs_env%bs_env%n_ao
2863 ELSE IF (basis_type == "RI_AUX") THEN
2864 n_bf = qs_env%bs_env%n_RI
2865 ELSE
2866 cpabort("Only ORB and RI_AUX basis implemented.")
2867 END IF
2868
2869 ALLOCATE (atom_from_bf_index(n_bf))
2870 CALL get_atom_index_from_basis_function_index(qs_env, atom_from_bf_index, n_bf, basis_type)
2871
2872 NULLIFY (cell, particle_set)
2873 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
2874 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 col_indices=col_indices)
2881
2882 CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp)
2883 index_to_cell => kpoints%index_to_cell
2884 num_cells = SIZE(index_to_cell, 2)
2885
2886 iatom_old = 0
2887 jatom_old = 0
2888
2889 DO irow = 1, nrow_local
2890 DO jcol = 1, ncol_local
2891
2892 i_bf = row_indices(irow)
2893 j_bf = col_indices(jcol)
2894
2895 iatom = atom_from_bf_index(i_bf)
2896 jatom = atom_from_bf_index(j_bf)
2897
2898 IF (PRESENT(wkp_ext)) THEN
2899 wkp_of_ikp = wkp_ext
2900 ELSE
2901 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 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 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 wkp_of_ikp = bs_env%wkp_no_extra(ikp)
2911 END SELECT
2912 END IF
2913
2914 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 cell, index_to_cell, hmat, particle_set)
2919
2920 iatom_old = iatom
2921 jatom_old = jatom
2922
2923 END IF
2924
2925 contribution = weight_re*real(cfm_w_ikp_freq_j%local_data(irow, jcol)) + &
2926 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 + contribution
2930
2931 END DO
2932 END DO
2933
2934 CALL timestop(handle)
2935
2936 END SUBROUTINE mic_contribution_from_ikp
2937
2938END MODULE post_scf_bandstructure_utils
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:117
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition grid_common.h:153
cell_types::pbc
Definition cell_types.F:92
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:55
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:89
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:195
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_cholesky_decompose
subroutine, public cp_cfm_cholesky_decompose(matrix, n, info_out)
Used to replace a symmetric positive definite matrix M with its Cholesky decomposition U: M = U^T * U...
Definition cp_cfm_basic_linalg.F:926
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_geeig
subroutine, public cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition cp_cfm_diag.F:145
cp_cfm_diag::cp_cfm_heevd
subroutine, public cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
Perform a diagonalisation of a complex matrix.
Definition cp_cfm_diag.F:52
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:121
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:159
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:607
cp_cfm_types::cp_cfm_set_all
subroutine, public cp_cfm_set_all(matrix, alpha, beta)
Set all elements of the full matrix to alpha. Besides, set all diagonal matrix elements to beta (if g...
Definition cp_cfm_types.F:179
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition cp_cfm_types.F:765
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:208
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:118
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:119
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::cp_fm_geeig_canon
subroutine, public cp_fm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
General Eigenvalue Problem AX = BXE Use canonical diagonalization : U*s**(-1/2)
Definition cp_fm_diag.F:1316
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
Definition cp_fm_types.F:570
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1016
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:515
cp_parser_methods
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition cp_parser_methods.F:20
cp_parser_methods::read_float_object
elemental subroutine, public read_float_object(string, object, error_message)
Returns a floating point number read from a string including fraction like z1/z2.
Definition cp_parser_methods.F:1258
gw_utils
Definition gw_utils.F:13
gw_utils::compute_xkp
subroutine, public compute_xkp(xkp, ikp_start, ikp_end, grid)
...
Definition gw_utils.F:519
gw_utils::kpoint_init_cell_index_simple
subroutine, public kpoint_init_cell_index_simple(kpoints, qs_env)
...
Definition gw_utils.F:489
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::int_ldos_z
integer, parameter, public int_ldos_z
Definition input_constants.F:1063
input_constants::gaussian
integer, parameter, public gaussian
Definition input_constants.F:153
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:724
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:697
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1040
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::max_line_length
integer, parameter, public max_line_length
Definition kinds.F:59
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::kpoint_create
subroutine, public kpoint_create(kpoint)
Create a kpoint environment.
Definition kpoint_types.F:188
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:333
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:123
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:143
mathconstants::gaussi
complex(kind=dp), parameter, public gaussi
Definition mathconstants.F:142
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:143
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::complex_diag
subroutine, public complex_diag(matrix, eigenvectors, eigenvalues)
Diagonalizes a local complex Hermitian matrix using LAPACK. Based on cp_cfm_heevd.
Definition mathlib.F:1746
mathlib::inv_3x3
pure real(kind=dp) function, dimension(3, 3), public inv_3x3(a)
Returns the inverse of the 3 x 3 matrix a.
Definition mathlib.F:516
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
post_scf_bandstructure_types
Definition post_scf_bandstructure_types.F:13
post_scf_bandstructure_utils
Definition post_scf_bandstructure_utils.F:13
post_scf_bandstructure_utils::dos_pdos_ldos
subroutine, public dos_pdos_ldos(qs_env, bs_env)
...
Definition post_scf_bandstructure_utils.F:1339
post_scf_bandstructure_utils::bandstructure_primitive_cell_spinor
subroutine, public bandstructure_primitive_cell_spinor(qs_env, bs_env, eigenvalues, filename, cfm_h_gamma_spinor)
...
Definition post_scf_bandstructure_utils.F:1078
post_scf_bandstructure_utils::cfm_ikp_from_fm_gamma
subroutine, public cfm_ikp_from_fm_gamma(cfm_ikp, fm_gamma, ikp, qs_env, kpoints, basis_type)
...
Definition post_scf_bandstructure_utils.F:2701
post_scf_bandstructure_utils::mic_contribution_from_ikp
subroutine, public mic_contribution_from_ikp(bs_env, qs_env, fm_w_mic_freq_j, cfm_w_ikp_freq_j, ikp, kpoints, basis_type, wkp_ext)
...
Definition post_scf_bandstructure_utils.F:2833
post_scf_bandstructure_utils::get_fname
subroutine, public get_fname(fname, bs_env, ikp, scf_gw, ispin, soc)
...
Definition post_scf_bandstructure_utils.F:2320
post_scf_bandstructure_utils::bandstructure_primitive_cell
subroutine, public bandstructure_primitive_cell(qs_env, bs_env, eigenvalues, filename, fm_h_gamma)
...
Definition post_scf_bandstructure_utils.F:930
post_scf_bandstructure_utils::create_and_init_bs_env
subroutine, public create_and_init_bs_env(qs_env, bs_env, post_scf_bandstructure_section)
...
Definition post_scf_bandstructure_utils.F:121
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition qs_collocate_density.F:38
qs_collocate_density::calculate_rho_elec
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
Definition qs_collocate_density.F:1710
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:502
qs_ks_types
Definition qs_ks_types.F:15
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:397
rpa_gw_im_time_util
Utility routines for GW with imaginary time.
Definition rpa_gw_im_time_util.F:13
rpa_gw_im_time_util::compute_weight_re_im
subroutine, public compute_weight_re_im(weight_re, weight_im, num_cells, iatom, jatom, xkp, wkp_w, cell, index_to_cell, hmat, particle_set)
...
Definition rpa_gw_im_time_util.F:802
rpa_gw_im_time_util::get_atom_index_from_basis_function_index
subroutine, public get_atom_index_from_basis_function_index(qs_env, atom_from_basis_index, basis_size, basis_type, first_bf_from_atom)
...
Definition rpa_gw_im_time_util.F:747
scf_control_types
parameters that control an scf iteration
Definition scf_control_types.F:17
soc_pseudopotential_methods
Definition soc_pseudopotential_methods.F:8
soc_pseudopotential_methods::remove_soc_outside_energy_window_mo
subroutine, public remove_soc_outside_energy_window_mo(cfm_ks_spinor, e_win, eigenval, e_homo, e_lumo)
...
Definition soc_pseudopotential_methods.F:338
soc_pseudopotential_utils
Definition soc_pseudopotential_utils.F:8
soc_pseudopotential_utils::add_cfm_submat
subroutine, public add_cfm_submat(cfm_mat_target, cfm_mat_source, nstart_row, nstart_col, factor)
...
Definition soc_pseudopotential_utils.F:206
soc_pseudopotential_utils::get_cfm_submat
subroutine, public get_cfm_submat(cfm_mat_target, cfm_mat_source, nstart_row, nstart_col)
...
Definition soc_pseudopotential_utils.F:270
soc_pseudopotential_utils::cfm_add_on_diag
subroutine, public cfm_add_on_diag(cfm, alpha)
...
Definition soc_pseudopotential_utils.F:121
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:68
cp_control_types::dft_control_type
Definition cp_control_types.F:559
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:138
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
post_scf_bandstructure_types::band_edges_type
Definition post_scf_bandstructure_types.F:45
post_scf_bandstructure_types::post_scf_bandstructure_type
Definition post_scf_bandstructure_types.F:52
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:135
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
scf_control_types::scf_control_type
Definition scf_control_types.F:122