Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation and writing of projected density of states
10 : !> The DOS is computed per angular momentum and per kind
11 : !> \par History
12 : !> -
13 : !> \author Marcella (29.02.2008,MK)
14 : ! **************************************************************************************************
15 : MODULE qs_pdos
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind,&
18 : get_atomic_kind_set
19 : USE basis_set_types, ONLY: get_gto_basis_set,&
20 : gto_basis_set_type
21 : USE cell_types, ONLY: cell_type,&
22 : pbc
23 : USE cp_array_utils, ONLY: cp_1d_r_p_type
24 : USE cp_blacs_env, ONLY: cp_blacs_env_type
25 : USE cp_cfm_types, ONLY: cp_cfm_create,&
26 : cp_cfm_get_submatrix,&
27 : cp_cfm_release,&
28 : cp_cfm_type
29 : USE cp_control_types, ONLY: dft_control_type
30 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
31 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
32 : USE cp_fm_diag, ONLY: cp_fm_power
33 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
34 : cp_fm_struct_release,&
35 : cp_fm_struct_type
36 : USE cp_fm_types, ONLY: cp_fm_create,&
37 : cp_fm_get_info,&
38 : cp_fm_get_submatrix,&
39 : cp_fm_release,&
40 : cp_fm_type
41 : USE cp_log_handling, ONLY: cp_get_default_logger,&
42 : cp_logger_get_default_io_unit,&
43 : cp_logger_type,&
44 : cp_to_string
45 : USE cp_output_handling, ONLY: cp_p_file,&
46 : cp_print_key_finished_output,&
47 : cp_print_key_should_output,&
48 : cp_print_key_unit_nr
49 : USE input_section_types, ONLY: section_vals_get,&
50 : section_vals_get_subs_vals,&
51 : section_vals_type,&
52 : section_vals_val_get
53 : USE kinds, ONLY: default_string_length,&
54 : dp
55 : USE kpoint_methods, ONLY: lowdin_kp_mo_coeff
56 : USE kpoint_types, ONLY: kpoint_env_type,&
57 : kpoint_type
58 : USE memory_utilities, ONLY: reallocate
59 : USE message_passing, ONLY: mp_para_env_type
60 : USE orbital_pointers, ONLY: nso,&
61 : nsoset
62 : USE orbital_symbols, ONLY: l_sym,&
63 : sgf_symbol
64 : USE parallel_gemm_api, ONLY: parallel_gemm
65 : USE particle_types, ONLY: particle_type
66 : USE pw_env_types, ONLY: pw_env_get,&
67 : pw_env_type
68 : USE pw_pool_types, ONLY: pw_pool_p_type,&
69 : pw_pool_type
70 : USE pw_types, ONLY: pw_c1d_gs_type,&
71 : pw_r3d_rs_type
72 : USE qs_collocate_density, ONLY: calculate_wavefunction
73 : USE qs_dos_utils, ONLY: &
74 : add_broadened_value, broadening_cutoff, broadening_function, dos_density_scale, &
75 : dos_energy_label, dos_energy_scale, dos_energy_unit_ev, dos_energy_zero_absolute, &
76 : dos_energy_zero_auto, dos_energy_zero_hoco, dos_energy_zero_label, &
77 : dos_resolve_energy_zero, write_broadening_info
78 : USE qs_environment_types, ONLY: get_qs_env,&
79 : qs_environment_type
80 : USE qs_kind_types, ONLY: get_qs_kind,&
81 : get_qs_kind_set,&
82 : qs_kind_type
83 : USE qs_mo_types, ONLY: get_mo_set,&
84 : mo_set_type
85 : USE qs_scf_diagonalization, ONLY: diag_kp_smat
86 : USE qs_scf_types, ONLY: qs_scf_env_type
87 : #include "./base/base_uses.f90"
88 :
89 : IMPLICIT NONE
90 :
91 : PRIVATE
92 :
93 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_pdos'
94 :
95 : ! **************************************************************************************************
96 : ! *** Public subroutines ***
97 :
98 : PUBLIC :: calculate_projected_dos, calculate_projected_dos_kp
99 :
100 : TYPE ldos_type
101 : INTEGER :: maxl = -1, nlist = -1
102 : LOGICAL :: separate_components = .FALSE.
103 : INTEGER, DIMENSION(:), POINTER :: list_index => NULL()
104 : REAL(KIND=dp), DIMENSION(:, :), &
105 : POINTER :: pdos_array => NULL()
106 : END TYPE ldos_type
107 :
108 : TYPE r_ldos_type
109 : INTEGER :: nlist = -1, npoints = -1
110 : INTEGER, DIMENSION(:, :), POINTER :: index_grid_local => NULL()
111 : INTEGER, DIMENSION(:), POINTER :: list_index => NULL()
112 : REAL(KIND=dp), DIMENSION(:), POINTER :: x_range => NULL(), y_range => NULL(), z_range => NULL()
113 : REAL(KIND=dp), DIMENSION(:), POINTER :: eval_range => NULL()
114 : REAL(KIND=dp), DIMENSION(:), &
115 : POINTER :: pdos_array => NULL()
116 : END TYPE r_ldos_type
117 :
118 : TYPE ldos_p_type
119 : TYPE(ldos_type), POINTER :: ldos => NULL()
120 : END TYPE ldos_p_type
121 :
122 : TYPE r_ldos_p_type
123 : TYPE(r_ldos_type), POINTER :: ldos => NULL()
124 : END TYPE r_ldos_p_type
125 : CONTAINS
126 :
127 : ! **************************************************************************************************
128 : !> \brief Compute and write projected density of states
129 : !> \param mo_set ...
130 : !> \param atomic_kind_set ...
131 : !> \param qs_kind_set ...
132 : !> \param particle_set ...
133 : !> \param qs_env ...
134 : !> \param dft_section ...
135 : !> \param ispin ...
136 : !> \param xas_mittle ...
137 : !> \param external_matrix_shalf ...
138 : !> \param unoccupied_orbs ...
139 : !> \param unoccupied_evals ...
140 : !> \param pdos_print_key ...
141 : !> \param write_pdos ...
142 : !> \param write_pdos_raw ...
143 : !> \date 26.02.2008
144 : !> \par History:
145 : !> - Added optional external matrix_shalf to avoid recomputing it (A. Bussy, 09.2019)
146 : !> \par Variables
147 : !> -
148 : !> -
149 : !> \author MI
150 : !> \version 1.0
151 : ! **************************************************************************************************
152 34 : SUBROUTINE calculate_projected_dos(mo_set, atomic_kind_set, qs_kind_set, particle_set, qs_env, &
153 : dft_section, ispin, xas_mittle, external_matrix_shalf, &
154 : unoccupied_orbs, unoccupied_evals, pdos_print_key, write_pdos, write_pdos_raw)
155 :
156 : TYPE(mo_set_type), INTENT(IN) :: mo_set
157 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
158 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
159 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
160 : TYPE(qs_environment_type), POINTER :: qs_env
161 : TYPE(section_vals_type), POINTER :: dft_section
162 : INTEGER, INTENT(IN), OPTIONAL :: ispin
163 : CHARACTER(LEN=default_string_length), INTENT(IN), &
164 : OPTIONAL :: xas_mittle
165 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL, TARGET :: external_matrix_shalf, unoccupied_orbs
166 : TYPE(cp_1d_r_p_type), INTENT(IN), OPTIONAL, TARGET :: unoccupied_evals
167 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: pdos_print_key
168 : LOGICAL, INTENT(IN), OPTIONAL :: write_pdos, write_pdos_raw
169 :
170 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_projected_dos'
171 :
172 : CHARACTER(LEN=16) :: energy_label, fmtstr2
173 : CHARACTER(LEN=27) :: fmtstr1
174 : CHARACTER(LEN=32) :: zero_label
175 34 : CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:, :, :) :: tmp_str
176 : CHARACTER(LEN=default_string_length) :: kind_name, my_act, my_mittle, my_pos, &
177 : my_print_key, spin(2)
178 : CHARACTER(LEN=default_string_length), &
179 34 : ALLOCATABLE, DIMENSION(:) :: ldos_index, r_ldos_index
180 : INTEGER :: broaden_type, energy_unit, energy_zero, handle, homo, i, iatom, ikind, il, ildos, &
181 : im, imo, in_x, in_y, in_z, ir, irow, iset, isgf, ishell, iso, iterstep, iw, j, jx, jy, &
182 : jz, k, lcomponent, lshell, maxl, maxlgto, my_spin, n_dependent, n_r_ldos, n_rep, nao, &
183 : natom, ncol_global, ndigits, nkind, nldos, nmo, np_tot, npoints, nrow_global, nset, nsgf, &
184 : nvirt, out_each, output_unit, resolved_energy_zero
185 34 : INTEGER, ALLOCATABLE, DIMENSION(:) :: firstrow
186 34 : INTEGER, DIMENSION(:), POINTER :: list, nshell
187 34 : INTEGER, DIMENSION(:, :), POINTER :: bo, l
188 : LOGICAL :: append, calc_matsh, do_broaden, do_ldos, do_r_ldos, do_virt, &
189 : fractional_occupation, ionode, separate_components, should_output, write_raw, &
190 : write_standard
191 34 : LOGICAL, DIMENSION(:, :), POINTER :: read_r
192 : REAL(KIND=dp) :: broaden_width, de, dh(3, 3), dvol, e_fermi, energy_factor, energy_ref, &
193 : ev_factor, hoco, r(3), r_vec(3), ratom(3), voigt_mixing
194 34 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, evals_virt, &
195 34 : occupation_numbers
196 34 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
197 34 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pdos_array
198 : TYPE(cell_type), POINTER :: cell
199 : TYPE(cp_blacs_env_type), POINTER :: context
200 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
201 : TYPE(cp_fm_type) :: matrix_shalfc, matrix_work
202 : TYPE(cp_fm_type), POINTER :: matrix_shalf, mo_coeff, mo_virt
203 : TYPE(cp_logger_type), POINTER :: logger
204 34 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: s_matrix
205 : TYPE(dft_control_type), POINTER :: dft_control
206 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
207 34 : TYPE(ldos_p_type), DIMENSION(:), POINTER :: ldos_p
208 : TYPE(mp_para_env_type), POINTER :: para_env
209 : TYPE(pw_c1d_gs_type) :: wf_g
210 : TYPE(pw_env_type), POINTER :: pw_env
211 34 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
212 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
213 : TYPE(pw_r3d_rs_type) :: wf_r
214 34 : TYPE(r_ldos_p_type), DIMENSION(:), POINTER :: r_ldos_p
215 : TYPE(section_vals_type), POINTER :: ldos_section
216 :
217 34 : NULLIFY (logger)
218 68 : logger => cp_get_default_logger()
219 : ionode = logger%para_env%is_source()
220 34 : my_print_key = "PRINT%PDOS"
221 34 : IF (PRESENT(pdos_print_key)) my_print_key = TRIM(pdos_print_key)
222 34 : write_standard = .TRUE.
223 34 : IF (PRESENT(write_pdos)) write_standard = write_pdos
224 34 : write_raw = .TRUE.
225 34 : IF (PRESENT(write_pdos_raw)) write_raw = write_pdos_raw
226 34 : write_raw = write_raw .AND. write_standard
227 : should_output = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
228 34 : TRIM(my_print_key)), cp_p_file)
229 34 : output_unit = cp_logger_get_default_io_unit(logger)
230 :
231 34 : spin(1) = "ALPHA"
232 34 : spin(2) = "BETA"
233 34 : IF ((.NOT. should_output)) RETURN
234 :
235 34 : NULLIFY (context, s_matrix, orb_basis_set, para_env, pdos_array)
236 34 : NULLIFY (eigenvalues, fm_struct_tmp, mo_coeff, vecbuffer, mo_virt)
237 34 : NULLIFY (ldos_section, list, cell, pw_env, auxbas_pw_pool, evals_virt)
238 34 : NULLIFY (occupation_numbers, ldos_p, r_ldos_p, dft_control, occupation_numbers)
239 :
240 34 : CALL timeset(routineN, handle)
241 34 : iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
242 :
243 34 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,I10))') &
244 17 : " Calculate PDOS at iteration step ", iterstep
245 : CALL get_qs_env(qs_env=qs_env, &
246 : matrix_s=s_matrix, &
247 34 : dft_control=dft_control)
248 :
249 34 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
250 34 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf, maxlgto=maxlgto)
251 34 : nkind = SIZE(atomic_kind_set)
252 :
253 : CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, homo=homo, nao=nao, nmo=nmo, &
254 34 : mu=e_fermi)
255 : CALL cp_fm_get_info(mo_coeff, &
256 : context=context, para_env=para_env, &
257 : nrow_global=nrow_global, &
258 34 : ncol_global=ncol_global)
259 :
260 34 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%OUT_EACH_MO", i_val=out_each)
261 34 : IF (out_each == -1) out_each = nao + 1
262 34 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%DELTA_E", r_val=de)
263 34 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%BROADEN_TYPE", i_val=broaden_type)
264 34 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%BROADEN_WIDTH", r_val=broaden_width)
265 34 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%VOIGT_MIXING", r_val=voigt_mixing)
266 34 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%NDIGITS", i_val=ndigits)
267 34 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%ENERGY_UNIT", i_val=energy_unit)
268 34 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%ENERGY_ZERO", i_val=energy_zero)
269 34 : ndigits = MIN(MAX(ndigits, 1), 10)
270 34 : do_broaden = (broaden_width > 0.0_dp)
271 34 : IF (do_broaden) de = MAX(de, 0.00001_dp)
272 34 : nvirt = 0
273 34 : NULLIFY (evals_virt)
274 34 : IF (PRESENT(unoccupied_orbs) .AND. PRESENT(unoccupied_evals)) THEN
275 2 : IF (ASSOCIATED(unoccupied_evals%array)) THEN
276 2 : nvirt = SIZE(unoccupied_evals%array)
277 2 : IF (nvirt > 0) THEN
278 2 : mo_virt => unoccupied_orbs
279 2 : evals_virt => unoccupied_evals%array
280 : END IF
281 : END IF
282 : END IF
283 34 : do_virt = (nvirt > 0)
284 :
285 34 : calc_matsh = .TRUE.
286 34 : IF (PRESENT(external_matrix_shalf)) calc_matsh = .FALSE.
287 :
288 : ! Create S^1/2 : from sparse to full matrix, if no external available
289 64 : IF (calc_matsh) THEN
290 32 : NULLIFY (matrix_shalf)
291 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
292 32 : nrow_global=nrow_global, ncol_global=nrow_global)
293 32 : ALLOCATE (matrix_shalf)
294 32 : CALL cp_fm_create(matrix_shalf, fm_struct_tmp, name="matrix_shalf")
295 32 : CALL cp_fm_create(matrix_work, fm_struct_tmp, name="matrix_work")
296 32 : CALL cp_fm_struct_release(fm_struct_tmp)
297 32 : CALL copy_dbcsr_to_fm(s_matrix(1)%matrix, matrix_shalf)
298 32 : CALL cp_fm_power(matrix_shalf, matrix_work, 0.5_dp, EPSILON(0.0_dp), n_dependent)
299 32 : CALL cp_fm_release(matrix_work)
300 : ELSE
301 : matrix_shalf => external_matrix_shalf
302 : END IF
303 :
304 : ! Multiply S^(1/2) time the mOS coefficients to get orthonormalized MOS
305 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
306 34 : nrow_global=nrow_global, ncol_global=ncol_global)
307 34 : CALL cp_fm_create(matrix_shalfc, fm_struct_tmp, name="matrix_shalfc")
308 : CALL parallel_gemm("N", "N", nrow_global, ncol_global, nrow_global, &
309 34 : 1.0_dp, matrix_shalf, mo_coeff, 0.0_dp, matrix_shalfc)
310 34 : CALL cp_fm_struct_release(fm_struct_tmp)
311 :
312 34 : IF (do_virt) THEN
313 2 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T14,I10,T27,A))') &
314 1 : " Use ", nvirt, " additional unoccupied KS orbitals"
315 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
316 2 : nrow_global=nrow_global, ncol_global=nvirt)
317 2 : CALL cp_fm_create(matrix_work, fm_struct_tmp, name="matrix_shalfc")
318 : CALL parallel_gemm("N", "N", nrow_global, nvirt, nrow_global, &
319 2 : 1.0_dp, matrix_shalf, mo_virt, 0.0_dp, matrix_work)
320 2 : CALL cp_fm_struct_release(fm_struct_tmp)
321 : END IF
322 :
323 34 : IF (calc_matsh) THEN
324 32 : CALL cp_fm_release(matrix_shalf)
325 32 : DEALLOCATE (matrix_shalf)
326 : END IF
327 : ! Array to store the PDOS per kind and angular momentum
328 34 : do_ldos = .FALSE.
329 34 : ldos_section => section_vals_get_subs_vals(dft_section, TRIM(my_print_key)//"%LDOS")
330 :
331 34 : CALL section_vals_get(ldos_section, n_repetition=nldos)
332 34 : IF (nldos > 0) THEN
333 8 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,I10))') &
334 4 : " Prepare the list of atoms for LDOS. Number of lists: ", nldos
335 8 : do_ldos = .TRUE.
336 44 : ALLOCATE (ldos_p(nldos))
337 24 : ALLOCATE (ldos_index(nldos))
338 28 : DO ildos = 1, nldos
339 20 : WRITE (ldos_index(ildos), '(I0)') ildos
340 20 : ALLOCATE (ldos_p(ildos)%ldos)
341 20 : NULLIFY (ldos_p(ildos)%ldos%pdos_array)
342 20 : NULLIFY (ldos_p(ildos)%ldos%list_index)
343 :
344 20 : CALL section_vals_val_get(ldos_section, "LIST", i_rep_section=ildos, n_rep_val=n_rep)
345 20 : IF (n_rep > 0) THEN
346 20 : ldos_p(ildos)%ldos%nlist = 0
347 40 : DO ir = 1, n_rep
348 20 : NULLIFY (list)
349 : CALL section_vals_val_get(ldos_section, "LIST", i_rep_section=ildos, i_rep_val=ir, &
350 20 : i_vals=list)
351 40 : IF (ASSOCIATED(list)) THEN
352 20 : CALL reallocate(ldos_p(ildos)%ldos%list_index, 1, ldos_p(ildos)%ldos%nlist + SIZE(list))
353 76 : DO i = 1, SIZE(list)
354 76 : ldos_p(ildos)%ldos%list_index(i + ldos_p(ildos)%ldos%nlist) = list(i)
355 : END DO
356 20 : ldos_p(ildos)%ldos%nlist = ldos_p(ildos)%ldos%nlist + SIZE(list)
357 : END IF
358 : END DO
359 : ELSE
360 : ! stop, LDOS without list of atoms is not implemented
361 : END IF
362 :
363 20 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='((T10,A,T18,I6,T25,A,T36,I10,A))') &
364 10 : " List ", ildos, " contains ", ldos_p(ildos)%ldos%nlist, " atoms"
365 : CALL section_vals_val_get(ldos_section, "COMPONENTS", i_rep_section=ildos, &
366 20 : l_val=ldos_p(ildos)%ldos%separate_components)
367 20 : IF (ldos_p(ildos)%ldos%separate_components) THEN
368 16 : ALLOCATE (ldos_p(ildos)%ldos%pdos_array(nsoset(maxlgto), nmo + nvirt))
369 : ELSE
370 64 : ALLOCATE (ldos_p(ildos)%ldos%pdos_array(0:maxlgto, nmo + nvirt))
371 : END IF
372 716 : ldos_p(ildos)%ldos%pdos_array = 0.0_dp
373 48 : ldos_p(ildos)%ldos%maxl = -1
374 :
375 : END DO
376 : END IF
377 :
378 34 : do_r_ldos = .FALSE.
379 34 : ldos_section => section_vals_get_subs_vals(dft_section, TRIM(my_print_key)//"%R_LDOS")
380 34 : CALL section_vals_get(ldos_section, n_repetition=n_r_ldos)
381 34 : IF (n_r_ldos > 0) THEN
382 0 : do_r_ldos = .TRUE.
383 0 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,I10))') &
384 0 : " Prepare the list of points for R_LDOS. Number of lists: ", n_r_ldos
385 0 : ALLOCATE (r_ldos_p(n_r_ldos))
386 0 : ALLOCATE (r_ldos_index(n_r_ldos))
387 : CALL get_qs_env(qs_env=qs_env, &
388 : cell=cell, &
389 : dft_control=dft_control, &
390 0 : pw_env=pw_env)
391 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
392 0 : pw_pools=pw_pools)
393 :
394 0 : CALL auxbas_pw_pool%create_pw(wf_r)
395 0 : CALL auxbas_pw_pool%create_pw(wf_g)
396 0 : ALLOCATE (read_r(4, n_r_ldos))
397 0 : DO ildos = 1, n_r_ldos
398 0 : WRITE (r_ldos_index(ildos), '(I0)') ildos
399 0 : ALLOCATE (r_ldos_p(ildos)%ldos)
400 0 : NULLIFY (r_ldos_p(ildos)%ldos%pdos_array)
401 0 : NULLIFY (r_ldos_p(ildos)%ldos%list_index)
402 :
403 0 : CALL section_vals_val_get(ldos_section, "LIST", i_rep_section=ildos, n_rep_val=n_rep)
404 0 : IF (n_rep > 0) THEN
405 0 : r_ldos_p(ildos)%ldos%nlist = 0
406 0 : DO ir = 1, n_rep
407 0 : NULLIFY (list)
408 : CALL section_vals_val_get(ldos_section, "LIST", i_rep_section=ildos, i_rep_val=ir, &
409 0 : i_vals=list)
410 0 : IF (ASSOCIATED(list)) THEN
411 0 : CALL reallocate(r_ldos_p(ildos)%ldos%list_index, 1, r_ldos_p(ildos)%ldos%nlist + SIZE(list))
412 0 : DO i = 1, SIZE(list)
413 0 : r_ldos_p(ildos)%ldos%list_index(i + r_ldos_p(ildos)%ldos%nlist) = list(i)
414 : END DO
415 0 : r_ldos_p(ildos)%ldos%nlist = r_ldos_p(ildos)%ldos%nlist + SIZE(list)
416 : END IF
417 : END DO
418 : ELSE
419 : ! stop, LDOS without list of atoms is not implemented
420 : END IF
421 :
422 0 : ALLOCATE (r_ldos_p(ildos)%ldos%pdos_array(nmo + nvirt))
423 0 : r_ldos_p(ildos)%ldos%pdos_array = 0.0_dp
424 0 : read_r(1:3, ildos) = .FALSE.
425 0 : CALL section_vals_val_get(ldos_section, "XRANGE", i_rep_section=ildos, explicit=read_r(1, ildos))
426 0 : IF (read_r(1, ildos)) THEN
427 : CALL section_vals_val_get(ldos_section, "XRANGE", i_rep_section=ildos, r_vals= &
428 0 : r_ldos_p(ildos)%ldos%x_range)
429 : ELSE
430 0 : ALLOCATE (r_ldos_p(ildos)%ldos%x_range(2))
431 0 : r_ldos_p(ildos)%ldos%x_range = 0.0_dp
432 : END IF
433 0 : CALL section_vals_val_get(ldos_section, "YRANGE", i_rep_section=ildos, explicit=read_r(2, ildos))
434 0 : IF (read_r(2, ildos)) THEN
435 : CALL section_vals_val_get(ldos_section, "YRANGE", i_rep_section=ildos, r_vals= &
436 0 : r_ldos_p(ildos)%ldos%y_range)
437 : ELSE
438 0 : ALLOCATE (r_ldos_p(ildos)%ldos%y_range(2))
439 0 : r_ldos_p(ildos)%ldos%y_range = 0.0_dp
440 : END IF
441 0 : CALL section_vals_val_get(ldos_section, "ZRANGE", i_rep_section=ildos, explicit=read_r(3, ildos))
442 0 : IF (read_r(3, ildos)) THEN
443 : CALL section_vals_val_get(ldos_section, "ZRANGE", i_rep_section=ildos, r_vals= &
444 0 : r_ldos_p(ildos)%ldos%z_range)
445 : ELSE
446 0 : ALLOCATE (r_ldos_p(ildos)%ldos%z_range(2))
447 0 : r_ldos_p(ildos)%ldos%z_range = 0.0_dp
448 : END IF
449 :
450 0 : CALL section_vals_val_get(ldos_section, "ERANGE", i_rep_section=ildos, explicit=read_r(4, ildos))
451 0 : IF (read_r(4, ildos)) THEN
452 : CALL section_vals_val_get(ldos_section, "ERANGE", i_rep_section=ildos, &
453 0 : r_vals=r_ldos_p(ildos)%ldos%eval_range)
454 : ELSE
455 0 : ALLOCATE (r_ldos_p(ildos)%ldos%eval_range(2))
456 0 : r_ldos_p(ildos)%ldos%eval_range(1) = -HUGE(0.0_dp)
457 0 : r_ldos_p(ildos)%ldos%eval_range(2) = +HUGE(0.0_dp)
458 : END IF
459 :
460 0 : bo => wf_r%pw_grid%bounds_local
461 0 : dh = wf_r%pw_grid%dh
462 0 : dvol = wf_r%pw_grid%dvol
463 0 : np_tot = wf_r%pw_grid%npts(1)*wf_r%pw_grid%npts(2)*wf_r%pw_grid%npts(3)
464 0 : ALLOCATE (r_ldos_p(ildos)%ldos%index_grid_local(3, np_tot))
465 :
466 0 : r_ldos_p(ildos)%ldos%npoints = 0
467 0 : DO jz = bo(1, 3), bo(2, 3)
468 0 : DO jy = bo(1, 2), bo(2, 2)
469 0 : DO jx = bo(1, 1), bo(2, 1)
470 : !compute the position of the grid point
471 0 : i = jx - wf_r%pw_grid%bounds(1, 1)
472 0 : j = jy - wf_r%pw_grid%bounds(1, 2)
473 0 : k = jz - wf_r%pw_grid%bounds(1, 3)
474 0 : r(3) = k*dh(3, 3) + j*dh(3, 2) + i*dh(3, 1)
475 0 : r(2) = k*dh(2, 3) + j*dh(2, 2) + i*dh(2, 1)
476 0 : r(1) = k*dh(1, 3) + j*dh(1, 2) + i*dh(1, 1)
477 :
478 0 : DO il = 1, r_ldos_p(ildos)%ldos%nlist
479 0 : iatom = r_ldos_p(ildos)%ldos%list_index(il)
480 0 : ratom = particle_set(iatom)%r
481 0 : r_vec = pbc(ratom, r, cell)
482 0 : IF (cell%orthorhombic) THEN
483 0 : IF (cell%perd(1) == 0) r_vec(1) = MODULO(r_vec(1), cell%hmat(1, 1))
484 0 : IF (cell%perd(2) == 0) r_vec(2) = MODULO(r_vec(2), cell%hmat(2, 2))
485 0 : IF (cell%perd(3) == 0) r_vec(3) = MODULO(r_vec(3), cell%hmat(3, 3))
486 : END IF
487 :
488 0 : in_x = 0
489 0 : in_y = 0
490 0 : in_z = 0
491 0 : IF (r_ldos_p(ildos)%ldos%x_range(1) /= 0.0_dp) THEN
492 0 : IF (r_vec(1) > r_ldos_p(ildos)%ldos%x_range(1) .AND. &
493 : r_vec(1) < r_ldos_p(ildos)%ldos%x_range(2)) THEN
494 0 : in_x = 1
495 : END IF
496 : ELSE
497 : in_x = 1
498 : END IF
499 0 : IF (r_ldos_p(ildos)%ldos%y_range(1) /= 0.0_dp) THEN
500 0 : IF (r_vec(2) > r_ldos_p(ildos)%ldos%y_range(1) .AND. &
501 : r_vec(2) < r_ldos_p(ildos)%ldos%y_range(2)) THEN
502 0 : in_y = 1
503 : END IF
504 : ELSE
505 : in_y = 1
506 : END IF
507 0 : IF (r_ldos_p(ildos)%ldos%z_range(1) /= 0.0_dp) THEN
508 0 : IF (r_vec(3) > r_ldos_p(ildos)%ldos%z_range(1) .AND. &
509 : r_vec(3) < r_ldos_p(ildos)%ldos%z_range(2)) THEN
510 0 : in_z = 1
511 : END IF
512 : ELSE
513 : in_z = 1
514 : END IF
515 0 : IF (in_x*in_y*in_z > 0) THEN
516 0 : r_ldos_p(ildos)%ldos%npoints = r_ldos_p(ildos)%ldos%npoints + 1
517 0 : r_ldos_p(ildos)%ldos%index_grid_local(1, r_ldos_p(ildos)%ldos%npoints) = jx
518 0 : r_ldos_p(ildos)%ldos%index_grid_local(2, r_ldos_p(ildos)%ldos%npoints) = jy
519 0 : r_ldos_p(ildos)%ldos%index_grid_local(3, r_ldos_p(ildos)%ldos%npoints) = jz
520 0 : EXIT
521 : END IF
522 : END DO
523 : END DO
524 : END DO
525 : END DO
526 0 : CALL reallocate(r_ldos_p(ildos)%ldos%index_grid_local, 1, 3, 1, r_ldos_p(ildos)%ldos%npoints)
527 0 : npoints = r_ldos_p(ildos)%ldos%npoints
528 0 : CALL para_env%sum(npoints)
529 0 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='((T10,A,T18,I6,T25,A,T36,I10,A))') &
530 0 : " List ", ildos, " contains ", npoints, " grid points"
531 : END DO
532 : END IF
533 :
534 34 : IF (TRIM(my_print_key) == "PRINT%DOS") THEN
535 24 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%PDOS%COMPONENTS", l_val=separate_components)
536 : ELSE
537 10 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%COMPONENTS", l_val=separate_components)
538 : END IF
539 34 : IF (separate_components) THEN
540 90 : ALLOCATE (pdos_array(nsoset(maxlgto), nkind, nmo + nvirt))
541 : ELSE
542 80 : ALLOCATE (pdos_array(0:maxlgto, nkind, nmo + nvirt))
543 : END IF
544 34 : IF (do_virt) THEN
545 6 : ALLOCATE (eigenvalues(nmo + nvirt))
546 10 : eigenvalues(1:nmo) = mo_set%eigenvalues(1:nmo)
547 22 : eigenvalues(nmo + 1:nmo + nvirt) = evals_virt(1:nvirt)
548 6 : ALLOCATE (occupation_numbers(nmo + nvirt))
549 20 : occupation_numbers(:) = 0.0_dp
550 10 : occupation_numbers(1:nmo) = mo_set%occupation_numbers(1:nmo)
551 : ELSE
552 32 : eigenvalues => mo_set%eigenvalues
553 32 : occupation_numbers => mo_set%occupation_numbers
554 : END IF
555 :
556 34 : hoco = -HUGE(0.0_dp)
557 34 : fractional_occupation = .FALSE.
558 318 : DO imo = 1, nmo + nvirt
559 284 : IF (occupation_numbers(imo) > 1.0e-10_dp) hoco = MAX(hoco, eigenvalues(imo))
560 284 : IF (ABS(occupation_numbers(imo) - REAL(NINT(occupation_numbers(imo)), KIND=dp)) > &
561 38 : 1.0e-8_dp) fractional_occupation = .TRUE.
562 : END DO
563 34 : IF (hoco < -0.5_dp*HUGE(0.0_dp)) hoco = e_fermi
564 34 : resolved_energy_zero = dos_resolve_energy_zero(energy_zero, dft_control%smear, fractional_occupation)
565 0 : SELECT CASE (resolved_energy_zero)
566 : CASE (dos_energy_zero_absolute)
567 0 : energy_ref = 0.0_dp
568 : CASE (dos_energy_zero_hoco)
569 32 : energy_ref = hoco
570 : CASE DEFAULT
571 34 : energy_ref = e_fermi
572 : END SELECT
573 34 : energy_factor = dos_energy_scale(energy_unit)
574 34 : energy_label = dos_energy_label(energy_unit)
575 34 : zero_label = dos_energy_zero_label(resolved_energy_zero)
576 34 : IF (energy_zero == dos_energy_zero_auto) zero_label = "AUTO -> "//TRIM(zero_label)
577 34 : ev_factor = dos_energy_scale(dos_energy_unit_ev)
578 :
579 4190 : pdos_array = 0.0_dp
580 34 : nao = mo_set%nao
581 102 : ALLOCATE (vecbuffer(1, nao))
582 1582 : vecbuffer = 0.0_dp
583 102 : ALLOCATE (firstrow(natom))
584 34 : firstrow = 0
585 :
586 : !Adjust energy range for r_ldos
587 34 : DO ildos = 1, n_r_ldos
588 0 : IF (eigenvalues(1) > r_ldos_p(ildos)%ldos%eval_range(1)) &
589 0 : r_ldos_p(ildos)%ldos%eval_range(1) = eigenvalues(1)
590 0 : IF (eigenvalues(nmo + nvirt) < r_ldos_p(ildos)%ldos%eval_range(2)) &
591 34 : r_ldos_p(ildos)%ldos%eval_range(2) = eigenvalues(nmo + nvirt)
592 : END DO
593 :
594 34 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T15,A))') &
595 17 : "---- PDOS: start iteration on the KS states --- "
596 :
597 318 : DO imo = 1, nmo + nvirt
598 :
599 284 : IF (output_unit > 0 .AND. MOD(imo, out_each) == 0) WRITE (UNIT=output_unit, FMT='((T20,A,I10))') &
600 0 : " KS state index : ", imo
601 : ! Extract the eigenvector from the distributed full matrix
602 284 : IF (imo > nmo) THEN
603 : CALL cp_fm_get_submatrix(matrix_work, vecbuffer, 1, imo - nmo, &
604 10 : nao, 1, transpose=.TRUE.)
605 : ELSE
606 : CALL cp_fm_get_submatrix(matrix_shalfc, vecbuffer, 1, imo, &
607 274 : nao, 1, transpose=.TRUE.)
608 : END IF
609 :
610 : ! Calculate the pdos for all the kinds
611 284 : irow = 1
612 2464 : DO iatom = 1, natom
613 2180 : firstrow(iatom) = irow
614 2180 : NULLIFY (orb_basis_set)
615 2180 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
616 2180 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
617 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
618 : nset=nset, &
619 : nshell=nshell, &
620 2180 : l=l, maxl=maxl)
621 4644 : IF (separate_components) THEN
622 1988 : isgf = 1
623 4912 : DO iset = 1, nset
624 8394 : DO ishell = 1, nshell(iset)
625 3482 : lshell = l(ishell, iset)
626 12144 : DO iso = 1, nso(lshell)
627 5738 : lcomponent = nsoset(lshell - 1) + iso
628 : pdos_array(lcomponent, ikind, imo) = &
629 : pdos_array(lcomponent, ikind, imo) + &
630 5738 : vecbuffer(1, irow)*vecbuffer(1, irow)
631 9220 : irow = irow + 1
632 : END DO ! iso
633 : END DO ! ishell
634 : END DO ! iset
635 : ELSE
636 192 : isgf = 1
637 576 : DO iset = 1, nset
638 1200 : DO ishell = 1, nshell(iset)
639 624 : lshell = l(ishell, iset)
640 2240 : DO iso = 1, nso(lshell)
641 : pdos_array(lshell, ikind, imo) = &
642 : pdos_array(lshell, ikind, imo) + &
643 1232 : vecbuffer(1, irow)*vecbuffer(1, irow)
644 1856 : irow = irow + 1
645 : END DO ! iso
646 : END DO ! ishell
647 : END DO ! iset
648 : END IF
649 : END DO ! iatom
650 :
651 : ! Calculate the pdos for all the lists
652 404 : DO ildos = 1, nldos
653 728 : DO il = 1, ldos_p(ildos)%ldos%nlist
654 324 : iatom = ldos_p(ildos)%ldos%list_index(il)
655 :
656 324 : irow = firstrow(iatom)
657 324 : NULLIFY (orb_basis_set)
658 324 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
659 324 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
660 :
661 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
662 : nset=nset, &
663 : nshell=nshell, &
664 324 : l=l, maxl=maxl)
665 324 : ldos_p(ildos)%ldos%maxl = MAX(ldos_p(ildos)%ldos%maxl, maxl)
666 768 : IF (ldos_p(ildos)%ldos%separate_components) THEN
667 108 : isgf = 1
668 324 : DO iset = 1, nset
669 720 : DO ishell = 1, nshell(iset)
670 396 : lshell = l(ishell, iset)
671 1440 : DO iso = 1, nso(lshell)
672 828 : lcomponent = nsoset(lshell - 1) + iso
673 : ldos_p(ildos)%ldos%pdos_array(lcomponent, imo) = &
674 : ldos_p(ildos)%ldos%pdos_array(lcomponent, imo) + &
675 828 : vecbuffer(1, irow)*vecbuffer(1, irow)
676 1224 : irow = irow + 1
677 : END DO ! iso
678 : END DO ! ishell
679 : END DO ! iset
680 : ELSE
681 216 : isgf = 1
682 648 : DO iset = 1, nset
683 1428 : DO ishell = 1, nshell(iset)
684 780 : lshell = l(ishell, iset)
685 2820 : DO iso = 1, nso(lshell)
686 : ldos_p(ildos)%ldos%pdos_array(lshell, imo) = &
687 : ldos_p(ildos)%ldos%pdos_array(lshell, imo) + &
688 1608 : vecbuffer(1, irow)*vecbuffer(1, irow)
689 2388 : irow = irow + 1
690 : END DO ! iso
691 : END DO ! ishell
692 : END DO ! iset
693 : END IF
694 : END DO !il
695 : END DO !ildos
696 :
697 : ! Calculate the DOS projected in a given volume in real space
698 318 : DO ildos = 1, n_r_ldos
699 0 : IF (r_ldos_p(ildos)%ldos%eval_range(1) <= eigenvalues(imo) .AND. &
700 284 : r_ldos_p(ildos)%ldos%eval_range(2) >= eigenvalues(imo)) THEN
701 :
702 0 : IF (imo > nmo) THEN
703 : CALL calculate_wavefunction(mo_virt, imo - nmo, &
704 : wf_r, wf_g, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
705 0 : pw_env)
706 : ELSE
707 : CALL calculate_wavefunction(mo_coeff, imo, &
708 : wf_r, wf_g, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
709 0 : pw_env)
710 : END IF
711 0 : r_ldos_p(ildos)%ldos%pdos_array(imo) = 0.0_dp
712 0 : DO il = 1, r_ldos_p(ildos)%ldos%npoints
713 0 : j = j + 1
714 0 : jx = r_ldos_p(ildos)%ldos%index_grid_local(1, il)
715 0 : jy = r_ldos_p(ildos)%ldos%index_grid_local(2, il)
716 0 : jz = r_ldos_p(ildos)%ldos%index_grid_local(3, il)
717 : r_ldos_p(ildos)%ldos%pdos_array(imo) = r_ldos_p(ildos)%ldos%pdos_array(imo) + &
718 0 : wf_r%array(jx, jy, jz)*wf_r%array(jx, jy, jz)
719 : END DO
720 0 : r_ldos_p(ildos)%ldos%pdos_array(imo) = r_ldos_p(ildos)%ldos%pdos_array(imo)*dvol
721 : END IF
722 : END DO
723 : END DO ! imo
724 :
725 34 : CALL cp_fm_release(matrix_shalfc)
726 34 : DEALLOCATE (vecbuffer)
727 :
728 34 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%APPEND", l_val=append)
729 34 : IF (append .AND. iterstep > 1) THEN
730 6 : my_pos = "APPEND"
731 : ELSE
732 28 : my_pos = "REWIND"
733 : END IF
734 34 : my_act = "WRITE"
735 34 : IF (write_standard .OR. write_raw) THEN
736 100 : DO ikind = 1, nkind
737 :
738 66 : NULLIFY (orb_basis_set)
739 66 : CALL get_atomic_kind(atomic_kind_set(ikind), name=kind_name)
740 66 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
741 66 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, maxl=maxl)
742 :
743 : ! basis none has no associated maxl, and no pdos
744 66 : IF (maxl < 0) CYCLE
745 :
746 66 : IF (PRESENT(ispin)) THEN
747 20 : IF (PRESENT(xas_mittle)) THEN
748 20 : my_mittle = TRIM(xas_mittle)//TRIM(spin(ispin))//"_k"//TRIM(ADJUSTL(cp_to_string(ikind)))
749 : ELSE
750 0 : my_mittle = TRIM(spin(ispin))//"_k"//TRIM(ADJUSTL(cp_to_string(ikind)))
751 : END IF
752 20 : my_spin = ispin
753 : ELSE
754 46 : my_mittle = "k"//TRIM(ADJUSTL(cp_to_string(ikind)))
755 46 : my_spin = 1
756 : END IF
757 :
758 66 : IF (write_standard .AND. do_broaden) THEN
759 66 : IF (separate_components) THEN
760 : CALL write_broadened_pdos(logger, dft_section, TRIM(my_mittle), my_pos, my_act, iterstep, &
761 : e_fermi, hoco, energy_ref, TRIM(zero_label), &
762 : "Projected DOS for atomic kind "//TRIM(kind_name), &
763 : maxl, .TRUE., pdos_array(1:nsoset(maxl), ikind, :), &
764 : eigenvalues, nmo + nvirt, de, broaden_type, broaden_width, &
765 34 : voigt_mixing, ndigits, pdos_print_key=TRIM(my_print_key))
766 : ELSE
767 : CALL write_broadened_pdos(logger, dft_section, TRIM(my_mittle), my_pos, my_act, iterstep, &
768 : e_fermi, hoco, energy_ref, TRIM(zero_label), &
769 : "Projected DOS for atomic kind "//TRIM(kind_name), &
770 : maxl, .FALSE., pdos_array(0:maxl, ikind, :), &
771 : eigenvalues, nmo + nvirt, de, broaden_type, broaden_width, &
772 32 : voigt_mixing, ndigits, pdos_print_key=TRIM(my_print_key))
773 : END IF
774 : END IF
775 :
776 166 : IF (write_raw) THEN
777 : iw = cp_print_key_unit_nr(logger, dft_section, TRIM(my_print_key), &
778 : extension=".pdos_raw", file_position=my_pos, file_action=my_act, &
779 20 : file_form="FORMATTED", middle_name=TRIM(my_mittle))
780 20 : IF (iw > 0) THEN
781 :
782 10 : fmtstr1 = "(I8,2X,2F16.6, (2X,F16.8))"
783 10 : fmtstr2 = "(A42, (10X,A8))"
784 10 : IF (separate_components) THEN
785 10 : WRITE (UNIT=fmtstr1(15:16), FMT="(I2)") nsoset(maxl) + 1
786 10 : WRITE (UNIT=fmtstr2(6:7), FMT="(I2)") nsoset(maxl) + 1
787 : ELSE
788 0 : WRITE (UNIT=fmtstr1(15:16), FMT="(I2)") maxl + 2
789 0 : WRITE (UNIT=fmtstr2(6:7), FMT="(I2)") maxl + 2
790 : END IF
791 :
792 : WRITE (UNIT=iw, FMT="(A,I0)") &
793 10 : "# Projected DOS for atomic kind "//TRIM(kind_name)//" at iteration step i = ", iterstep
794 : WRITE (UNIT=iw, FMT="(A,F12.6,A,F12.6,A)") &
795 10 : "# E(Fermi) = ", e_fermi, " a.u. = ", e_fermi*ev_factor, " eV"
796 : WRITE (UNIT=iw, FMT="(A,F12.6,A,F12.6,A)") &
797 10 : "# E(HOCO) = ", hoco, " a.u. = ", hoco*ev_factor, " eV"
798 10 : WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label)
799 10 : IF (separate_components) THEN
800 40 : ALLOCATE (tmp_str(0:0, 0:maxl, -maxl:maxl))
801 340 : tmp_str = ""
802 39 : DO j = 0, maxl
803 124 : DO i = -j, j
804 114 : tmp_str(0, j, i) = sgf_symbol(0, j, i)
805 : END DO
806 : END DO
807 :
808 : WRITE (UNIT=iw, FMT=fmtstr2) &
809 10 : "# MO "//TRIM(energy_label)//" Occupation", "Total", &
810 134 : ((TRIM(tmp_str(0, il, im)), im=-il, il), il=0, maxl)
811 120 : DO imo = 1, nmo + nvirt
812 110 : WRITE (UNIT=iw, FMT=fmtstr1) imo, (eigenvalues(imo) - energy_ref)*energy_factor, &
813 1005 : occupation_numbers(imo), SUM(pdos_array(1:nsoset(maxl), ikind, imo)), &
814 1125 : (pdos_array(lshell, ikind, imo), lshell=1, nsoset(maxl))
815 : END DO
816 10 : DEALLOCATE (tmp_str)
817 : ELSE
818 : WRITE (UNIT=iw, FMT=fmtstr2) &
819 0 : "# MO "//TRIM(energy_label)//" Occupation", "Total", &
820 0 : (TRIM(l_sym(il)), il=0, maxl)
821 0 : DO imo = 1, nmo + nvirt
822 0 : WRITE (UNIT=iw, FMT=fmtstr1) imo, (eigenvalues(imo) - energy_ref)*energy_factor, &
823 0 : occupation_numbers(imo), SUM(pdos_array(0:maxl, ikind, imo)), &
824 0 : (pdos_array(lshell, ikind, imo), lshell=0, maxl)
825 : END DO
826 : END IF
827 : END IF
828 : CALL cp_print_key_finished_output(iw, logger, dft_section, &
829 20 : TRIM(my_print_key))
830 : END IF
831 :
832 : END DO ! ikind
833 : END IF
834 :
835 : ! write the pdos for the lists, each ona different file,
836 : ! the filenames are indexed with the list number
837 54 : DO ildos = 1, nldos
838 : ! basis none has no associated maxl, and no pdos
839 54 : IF (ldos_p(ildos)%ldos%maxl > 0) THEN
840 :
841 20 : IF (PRESENT(ispin)) THEN
842 0 : IF (PRESENT(xas_mittle)) THEN
843 0 : my_mittle = TRIM(xas_mittle)//TRIM(spin(ispin))//"_list"//TRIM(ldos_index(ildos))
844 : ELSE
845 0 : my_mittle = TRIM(spin(ispin))//"_list"//TRIM(ldos_index(ildos))
846 : END IF
847 0 : my_spin = ispin
848 : ELSE
849 20 : my_mittle = "list"//TRIM(ldos_index(ildos))
850 20 : my_spin = 1
851 : END IF
852 :
853 20 : IF (do_broaden) THEN
854 20 : IF (ldos_p(ildos)%ldos%separate_components) THEN
855 : CALL write_broadened_pdos(logger, dft_section, TRIM(my_mittle), my_pos, my_act, iterstep, &
856 : e_fermi, hoco, energy_ref, TRIM(zero_label), &
857 : "Projected DOS for atom list "//TRIM(ldos_index(ildos)), &
858 : ldos_p(ildos)%ldos%maxl, .TRUE., &
859 : ldos_p(ildos)%ldos%pdos_array(1:nsoset(ldos_p(ildos)%ldos%maxl), :), &
860 : eigenvalues, nmo + nvirt, de, broaden_type, broaden_width, &
861 4 : voigt_mixing, ndigits, pdos_print_key=TRIM(my_print_key))
862 : ELSE
863 : CALL write_broadened_pdos(logger, dft_section, TRIM(my_mittle), my_pos, my_act, iterstep, &
864 : e_fermi, hoco, energy_ref, TRIM(zero_label), &
865 : "Projected DOS for atom list "//TRIM(ldos_index(ildos)), &
866 : ldos_p(ildos)%ldos%maxl, .FALSE., &
867 : ldos_p(ildos)%ldos%pdos_array(0:ldos_p(ildos)%ldos%maxl, :), &
868 : eigenvalues, nmo + nvirt, de, broaden_type, broaden_width, &
869 16 : voigt_mixing, ndigits, pdos_print_key=TRIM(my_print_key))
870 : END IF
871 : END IF
872 :
873 : iw = cp_print_key_unit_nr(logger, dft_section, TRIM(my_print_key), &
874 : extension=".pdos_raw", file_position=my_pos, file_action=my_act, &
875 20 : file_form="FORMATTED", middle_name=TRIM(my_mittle))
876 20 : IF (iw > 0) THEN
877 :
878 10 : fmtstr1 = "(I8,2X,2F16.6, (2X,F16.8))"
879 10 : fmtstr2 = "(A42, (10X,A8))"
880 10 : IF (ldos_p(ildos)%ldos%separate_components) THEN
881 2 : WRITE (UNIT=fmtstr1(15:16), FMT="(I2)") nsoset(ldos_p(ildos)%ldos%maxl) + 1
882 2 : WRITE (UNIT=fmtstr2(6:7), FMT="(I2)") nsoset(ldos_p(ildos)%ldos%maxl) + 1
883 : ELSE
884 8 : WRITE (UNIT=fmtstr1(15:16), FMT="(I2)") ldos_p(ildos)%ldos%maxl + 2
885 8 : WRITE (UNIT=fmtstr2(6:7), FMT="(I2)") ldos_p(ildos)%ldos%maxl + 2
886 : END IF
887 :
888 : WRITE (UNIT=iw, FMT="(A,I0,A,I0,A,I0)") &
889 10 : "# Projected DOS for list ", ildos, " of ", ldos_p(ildos)%ldos%nlist, &
890 20 : " atoms, at iteration step i = ", iterstep
891 : WRITE (UNIT=iw, FMT="(A,F12.6,A,F12.6,A)") &
892 10 : "# E(Fermi) = ", e_fermi, " a.u. = ", e_fermi*ev_factor, " eV"
893 : WRITE (UNIT=iw, FMT="(A,F12.6,A,F12.6,A)") &
894 10 : "# E(HOCO) = ", hoco, " a.u. = ", hoco*ev_factor, " eV"
895 10 : WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label)
896 10 : IF (ldos_p(ildos)%ldos%separate_components) THEN
897 8 : ALLOCATE (tmp_str(0:0, 0:ldos_p(ildos)%ldos%maxl, -ldos_p(ildos)%ldos%maxl:ldos_p(ildos)%ldos%maxl))
898 72 : tmp_str = ""
899 8 : DO j = 0, ldos_p(ildos)%ldos%maxl
900 26 : DO i = -j, j
901 24 : tmp_str(0, j, i) = sgf_symbol(0, j, i)
902 : END DO
903 : END DO
904 :
905 : WRITE (UNIT=iw, FMT=fmtstr2) &
906 2 : "# MO "//TRIM(energy_label)//" Occupation", "Total", &
907 28 : ((TRIM(tmp_str(0, il, im)), im=-il, il), il=0, ldos_p(ildos)%ldos%maxl)
908 20 : DO imo = 1, nmo + nvirt
909 18 : WRITE (UNIT=iw, FMT=fmtstr1) imo, (eigenvalues(imo) - energy_ref)*energy_factor, &
910 18 : occupation_numbers(imo), &
911 180 : SUM(ldos_p(ildos)%ldos%pdos_array(1:nsoset(ldos_p(ildos)%ldos%maxl), imo)), &
912 180 : (ldos_p(ildos)%ldos%pdos_array(lshell, imo), &
913 218 : lshell=1, nsoset(ldos_p(ildos)%ldos%maxl))
914 : END DO
915 2 : DEALLOCATE (tmp_str)
916 : ELSE
917 : WRITE (UNIT=iw, FMT=fmtstr2) &
918 8 : "# MO "//TRIM(energy_label)//" Occupation", "Total", &
919 39 : (TRIM(l_sym(il)), il=0, ldos_p(ildos)%ldos%maxl)
920 50 : DO imo = 1, nmo + nvirt
921 42 : WRITE (UNIT=iw, FMT=fmtstr1) imo, (eigenvalues(imo) - energy_ref)*energy_factor, &
922 42 : occupation_numbers(imo), &
923 159 : SUM(ldos_p(ildos)%ldos%pdos_array(0:ldos_p(ildos)%ldos%maxl, imo)), &
924 209 : (ldos_p(ildos)%ldos%pdos_array(lshell, imo), lshell=0, ldos_p(ildos)%ldos%maxl)
925 : END DO
926 : END IF
927 : END IF
928 : CALL cp_print_key_finished_output(iw, logger, dft_section, &
929 20 : TRIM(my_print_key))
930 : END IF ! maxl>0
931 : END DO ! ildos
932 :
933 : ! write the pdos for the lists, each ona different file,
934 : ! the filenames are indexed with the list number
935 34 : DO ildos = 1, n_r_ldos
936 :
937 0 : npoints = r_ldos_p(ildos)%ldos%npoints
938 0 : CALL para_env%sum(npoints)
939 0 : CALL para_env%sum(np_tot)
940 0 : CALL para_env%sum(r_ldos_p(ildos)%ldos%pdos_array)
941 0 : IF (PRESENT(ispin)) THEN
942 0 : IF (PRESENT(xas_mittle)) THEN
943 0 : my_mittle = TRIM(xas_mittle)//TRIM(spin(ispin))//"_r_list"//TRIM(r_ldos_index(ildos))
944 : ELSE
945 0 : my_mittle = TRIM(spin(ispin))//"_r_list"//TRIM(r_ldos_index(ildos))
946 : END IF
947 0 : my_spin = ispin
948 : ELSE
949 0 : my_mittle = "r_list"//TRIM(r_ldos_index(ildos))
950 0 : my_spin = 1
951 : END IF
952 :
953 : iw = cp_print_key_unit_nr(logger, dft_section, TRIM(my_print_key), &
954 : extension=".pdos_raw", file_position=my_pos, file_action=my_act, &
955 0 : file_form="FORMATTED", middle_name=TRIM(my_mittle))
956 0 : IF (iw > 0) THEN
957 0 : fmtstr1 = "(I8,2X,2F16.6, (2X,F16.8))"
958 0 : fmtstr2 = "(A42, (10X,A8))"
959 :
960 : WRITE (UNIT=iw, FMT="(A,I0,A,F12.6,F12.6,A)") &
961 0 : "# Projected DOS in real space, using ", npoints, &
962 0 : " points of the grid, and eval in the range", r_ldos_p(ildos)%ldos%eval_range(1:2), " Hartree"
963 : WRITE (UNIT=iw, FMT="(A,F12.6,A,F12.6,A)") &
964 0 : "# E(Fermi) = ", e_fermi, " a.u. = ", e_fermi*ev_factor, " eV"
965 : WRITE (UNIT=iw, FMT="(A,F12.6,A,F12.6,A)") &
966 0 : "# E(HOCO) = ", hoco, " a.u. = ", hoco*ev_factor, " eV"
967 0 : WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label)
968 : WRITE (UNIT=iw, FMT="(A,A,A)") &
969 0 : "# MO ", TRIM(energy_label), " Occupation LDOS"
970 0 : DO imo = 1, nmo + nvirt
971 0 : IF (r_ldos_p(ildos)%ldos%eval_range(1) <= eigenvalues(imo) .AND. &
972 0 : r_ldos_p(ildos)%ldos%eval_range(2) >= eigenvalues(imo)) THEN
973 0 : WRITE (UNIT=iw, FMT="(I8,2X,2F16.6,E20.10,E20.10)") imo, &
974 0 : (eigenvalues(imo) - energy_ref)*energy_factor, occupation_numbers(imo), &
975 0 : r_ldos_p(ildos)%ldos%pdos_array(imo), r_ldos_p(ildos)%ldos%pdos_array(imo)*np_tot
976 : END IF
977 : END DO
978 :
979 : END IF
980 : CALL cp_print_key_finished_output(iw, logger, dft_section, &
981 34 : TRIM(my_print_key))
982 : END DO
983 :
984 : ! deallocate local variables
985 34 : DEALLOCATE (pdos_array)
986 34 : DEALLOCATE (firstrow)
987 34 : IF (do_ldos) THEN
988 28 : DO ildos = 1, nldos
989 20 : DEALLOCATE (ldos_p(ildos)%ldos%pdos_array)
990 20 : DEALLOCATE (ldos_p(ildos)%ldos%list_index)
991 28 : DEALLOCATE (ldos_p(ildos)%ldos)
992 : END DO
993 8 : DEALLOCATE (ldos_p)
994 8 : DEALLOCATE (ldos_index)
995 : END IF
996 34 : IF (do_r_ldos) THEN
997 0 : DO ildos = 1, n_r_ldos
998 0 : DEALLOCATE (r_ldos_p(ildos)%ldos%index_grid_local)
999 0 : DEALLOCATE (r_ldos_p(ildos)%ldos%pdos_array)
1000 0 : DEALLOCATE (r_ldos_p(ildos)%ldos%list_index)
1001 0 : IF (.NOT. read_r(1, ildos)) THEN
1002 0 : DEALLOCATE (r_ldos_p(ildos)%ldos%x_range)
1003 : END IF
1004 0 : IF (.NOT. read_r(2, ildos)) THEN
1005 0 : DEALLOCATE (r_ldos_p(ildos)%ldos%y_range)
1006 : END IF
1007 0 : IF (.NOT. read_r(3, ildos)) THEN
1008 0 : DEALLOCATE (r_ldos_p(ildos)%ldos%z_range)
1009 : END IF
1010 0 : IF (.NOT. read_r(4, ildos)) THEN
1011 0 : DEALLOCATE (r_ldos_p(ildos)%ldos%eval_range)
1012 : END IF
1013 0 : DEALLOCATE (r_ldos_p(ildos)%ldos)
1014 : END DO
1015 0 : DEALLOCATE (read_r)
1016 0 : DEALLOCATE (r_ldos_p)
1017 0 : DEALLOCATE (r_ldos_index)
1018 0 : CALL auxbas_pw_pool%give_back_pw(wf_r)
1019 0 : CALL auxbas_pw_pool%give_back_pw(wf_g)
1020 : END IF
1021 34 : IF (do_virt) THEN
1022 2 : CALL cp_fm_release(matrix_work)
1023 2 : DEALLOCATE (eigenvalues)
1024 2 : DEALLOCATE (occupation_numbers)
1025 : END IF
1026 :
1027 34 : CALL timestop(handle)
1028 :
1029 238 : END SUBROUTINE calculate_projected_dos
1030 :
1031 : ! **************************************************************************************************
1032 : !> \brief Compute and write broadened projected density of states for k-point calculations.
1033 : !> \param qs_env ...
1034 : !> \param dft_section ...
1035 : !> \param pdos_print_key ...
1036 : !> \param write_pdos ...
1037 : !> \param write_pdos_raw ...
1038 : ! **************************************************************************************************
1039 16 : SUBROUTINE calculate_projected_dos_kp(qs_env, dft_section, pdos_print_key, write_pdos, write_pdos_raw)
1040 :
1041 : TYPE(qs_environment_type), POINTER :: qs_env
1042 : TYPE(section_vals_type), POINTER :: dft_section
1043 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: pdos_print_key
1044 : LOGICAL, INTENT(IN), OPTIONAL :: write_pdos, write_pdos_raw
1045 :
1046 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_projected_dos_kp'
1047 :
1048 : CHARACTER(LEN=32) :: zero_label
1049 : CHARACTER(LEN=default_string_length) :: kind_name, my_act, my_mittle, my_pos, &
1050 : my_print_key, spin(2)
1051 16 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: zvecbuffer
1052 : INTEGER :: broaden_type, energy_zero, fractional_occupation_int, handle, icomp, ik, ikind, &
1053 : imo, ispin, iterstep, maxl, maxlgto, n_r_ldos, nao, ncomp, ndigits, nhist, nkind, nldos, &
1054 : nmo_kp, nspins, output_unit, resolved_energy_zero
1055 16 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ao_comp, ao_kind, ao_l, kind_maxl
1056 : LOGICAL :: append, fractional_occupation, &
1057 : separate_components, should_output, &
1058 : write_raw, write_standard
1059 : REAL(KIND=dp) :: broaden_cutoff, broaden_width, de, e1, &
1060 : e2, e_fermi(2), emax, emin, &
1061 : energy_ref(2), hoco(2), voigt_mixing, &
1062 : wkp
1063 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ao_weight
1064 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: proj_weight, vecbuffer
1065 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: pdos_curve
1066 16 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
1067 16 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1068 : TYPE(cp_cfm_type) :: cshalfc
1069 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
1070 : TYPE(cp_fm_type) :: shalfc
1071 : TYPE(cp_logger_type), POINTER :: logger
1072 16 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp
1073 : TYPE(dft_control_type), POINTER :: dft_control
1074 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1075 : TYPE(kpoint_env_type), POINTER :: kp
1076 : TYPE(kpoint_type), POINTER :: kpoints
1077 : TYPE(mo_set_type), POINTER :: mo_set
1078 : TYPE(mp_para_env_type), POINTER :: para_env
1079 16 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1080 16 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1081 : TYPE(qs_scf_env_type), POINTER :: scf_env
1082 : TYPE(section_vals_type), POINTER :: ldos_section
1083 :
1084 16 : NULLIFY (logger, kpoints, dft_control, para_env, atomic_kind_set, qs_kind_set, particle_set)
1085 16 : NULLIFY (matrix_s_kp, scf_env)
1086 16 : NULLIFY (kp, mo_set, eigenvalues, fm_struct_tmp, orb_basis_set, ldos_section)
1087 32 : logger => cp_get_default_logger()
1088 16 : my_print_key = "PRINT%PDOS"
1089 16 : IF (PRESENT(pdos_print_key)) my_print_key = TRIM(pdos_print_key)
1090 16 : write_standard = .TRUE.
1091 16 : IF (PRESENT(write_pdos)) write_standard = write_pdos
1092 16 : write_raw = .FALSE.
1093 16 : IF (PRESENT(write_pdos_raw)) write_raw = write_pdos_raw
1094 16 : write_raw = write_raw .AND. write_standard
1095 : should_output = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
1096 16 : TRIM(my_print_key)), cp_p_file)
1097 16 : output_unit = cp_logger_get_default_io_unit(logger)
1098 16 : IF ((.NOT. should_output)) RETURN
1099 :
1100 16 : CALL timeset(routineN, handle)
1101 16 : iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
1102 :
1103 16 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,I10))') &
1104 8 : " Calculate k-point PDOS at iteration step ", iterstep
1105 :
1106 : CALL get_qs_env(qs_env=qs_env, &
1107 : kpoints=kpoints, &
1108 : dft_control=dft_control, &
1109 : matrix_s_kp=matrix_s_kp, &
1110 : scf_env=scf_env, &
1111 : atomic_kind_set=atomic_kind_set, &
1112 : qs_kind_set=qs_kind_set, &
1113 16 : particle_set=particle_set)
1114 16 : para_env => kpoints%para_env_inter_kp
1115 16 : nspins = dft_control%nspins
1116 16 : nkind = SIZE(atomic_kind_set)
1117 16 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto)
1118 16 : IF (.NOT. ASSOCIATED(kpoints%kp_env)) THEN
1119 0 : CPWARN("No local k points available for k-point PDOS")
1120 0 : CALL timestop(handle)
1121 0 : RETURN
1122 : END IF
1123 16 : IF (SIZE(kpoints%kp_env) == 0) THEN
1124 0 : CPWARN("No local k points available for k-point PDOS")
1125 0 : CALL timestop(handle)
1126 0 : RETURN
1127 : END IF
1128 :
1129 16 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%DELTA_E", r_val=de)
1130 16 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%APPEND", l_val=append)
1131 16 : IF (TRIM(my_print_key) == "PRINT%DOS") THEN
1132 16 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%PDOS%COMPONENTS", l_val=separate_components)
1133 : ELSE
1134 0 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%COMPONENTS", l_val=separate_components)
1135 : END IF
1136 16 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%NDIGITS", i_val=ndigits)
1137 16 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%BROADEN_TYPE", i_val=broaden_type)
1138 16 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%BROADEN_WIDTH", r_val=broaden_width)
1139 16 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%VOIGT_MIXING", r_val=voigt_mixing)
1140 16 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%ENERGY_ZERO", i_val=energy_zero)
1141 16 : ndigits = MIN(MAX(ndigits, 1), 10)
1142 16 : de = MAX(de, 0.00001_dp)
1143 :
1144 16 : ldos_section => section_vals_get_subs_vals(dft_section, TRIM(my_print_key)//"%LDOS")
1145 16 : CALL section_vals_get(ldos_section, n_repetition=nldos)
1146 16 : ldos_section => section_vals_get_subs_vals(dft_section, TRIM(my_print_key)//"%R_LDOS")
1147 16 : CALL section_vals_get(ldos_section, n_repetition=n_r_ldos)
1148 16 : IF (nldos > 0 .OR. n_r_ldos > 0) THEN
1149 0 : CPWARN("LDOS/R_LDOS are not implemented for k-point PDOS and will be ignored")
1150 : END IF
1151 16 : IF (write_raw) THEN
1152 0 : CPWARN("Raw k-point PDOS output is not implemented and will be ignored")
1153 : END IF
1154 16 : IF (broaden_width <= 0.0_dp) THEN
1155 0 : CPWARN("Broadened k-point PDOS output requires a finite BROADEN_WIDTH")
1156 0 : CALL timestop(handle)
1157 0 : RETURN
1158 : END IF
1159 :
1160 16 : IF (separate_components) THEN
1161 2 : ncomp = nsoset(maxlgto)
1162 : ELSE
1163 14 : ncomp = maxlgto + 1
1164 : END IF
1165 48 : ALLOCATE (kind_maxl(nkind))
1166 32 : kind_maxl = -1
1167 32 : DO ikind = 1, nkind
1168 16 : NULLIFY (orb_basis_set)
1169 16 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1170 16 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, maxl=maxl)
1171 32 : kind_maxl(ikind) = maxl
1172 : END DO
1173 :
1174 16 : emin = HUGE(0.0_dp)
1175 16 : emax = -HUGE(0.0_dp)
1176 16 : e_fermi(:) = 0.0_dp
1177 48 : hoco(:) = -HUGE(0.0_dp)
1178 16 : fractional_occupation = .FALSE.
1179 16 : IF (kpoints%nkp /= 0) THEN
1180 116 : DO ik = 1, SIZE(kpoints%kp_env)
1181 100 : kp => kpoints%kp_env(ik)%kpoint_env
1182 216 : DO ispin = 1, nspins
1183 100 : mo_set => kp%mos(1, ispin)
1184 100 : CALL get_mo_set(mo_set=mo_set, nmo=nmo_kp, mu=e_fermi(ispin))
1185 100 : eigenvalues => mo_set%eigenvalues
1186 100 : occupation_numbers => mo_set%occupation_numbers
1187 1906 : DO imo = 1, nmo_kp
1188 1806 : IF (occupation_numbers(imo) > 1.0e-10_dp) hoco(ispin) = MAX(hoco(ispin), eigenvalues(imo))
1189 1806 : IF (ABS(occupation_numbers(imo) - REAL(NINT(occupation_numbers(imo)), KIND=dp)) > &
1190 282 : 1.0e-8_dp) fractional_occupation = .TRUE.
1191 : END DO
1192 1906 : e1 = MINVAL(eigenvalues(1:nmo_kp))
1193 1906 : e2 = MAXVAL(eigenvalues(1:nmo_kp))
1194 100 : emin = MIN(emin, e1)
1195 300 : emax = MAX(emax, e2)
1196 : END DO
1197 : END DO
1198 : END IF
1199 16 : CALL para_env%min(emin)
1200 16 : CALL para_env%max(emax)
1201 16 : CALL para_env%max(e_fermi)
1202 16 : CALL para_env%max(hoco)
1203 16 : fractional_occupation_int = MERGE(1, 0, fractional_occupation)
1204 16 : CALL para_env%max(fractional_occupation_int)
1205 16 : fractional_occupation = (fractional_occupation_int /= 0)
1206 32 : DO ispin = 1, nspins
1207 32 : IF (hoco(ispin) < -0.5_dp*HUGE(0.0_dp)) hoco(ispin) = e_fermi(ispin)
1208 : END DO
1209 16 : resolved_energy_zero = dos_resolve_energy_zero(energy_zero, dft_control%smear, fractional_occupation)
1210 0 : SELECT CASE (resolved_energy_zero)
1211 : CASE (dos_energy_zero_absolute)
1212 0 : energy_ref(:) = 0.0_dp
1213 : CASE (dos_energy_zero_hoco)
1214 2 : energy_ref(:) = hoco(:)
1215 : CASE DEFAULT
1216 16 : energy_ref(:) = e_fermi(:)
1217 : END SELECT
1218 16 : zero_label = dos_energy_zero_label(resolved_energy_zero)
1219 16 : IF (energy_zero == dos_energy_zero_auto) zero_label = "AUTO -> "//TRIM(zero_label)
1220 16 : broaden_cutoff = broadening_cutoff(broaden_type, broaden_width)
1221 16 : emin = emin - broaden_cutoff
1222 16 : emax = emax + broaden_cutoff
1223 16 : nhist = NINT((emax - emin)/de) + 1
1224 96 : ALLOCATE (pdos_curve(nhist, ncomp, nkind, nspins))
1225 16 : pdos_curve = 0.0_dp
1226 :
1227 : ! Ensure that S(k)^1/2 is available for the Lowdin projection.
1228 : ! This is normally only constructed for Lowdin population/DFT+U paths.
1229 16 : CALL diag_kp_smat(matrix_s_kp, kpoints, scf_env%scf_work1)
1230 :
1231 : ! Use the first local k point to construct the AO -> kind/l/component map.
1232 16 : kp => kpoints%kp_env(1)%kpoint_env
1233 16 : mo_set => kp%mos(1, 1)
1234 16 : CALL get_mo_set(mo_set=mo_set, nao=nao)
1235 16 : CALL build_pdos_ao_map(qs_kind_set, particle_set, nao, ao_kind, ao_l, ao_comp)
1236 96 : ALLOCATE (ao_weight(nao), proj_weight(ncomp, nkind))
1237 64 : ALLOCATE (vecbuffer(1, nao), zvecbuffer(1, nao))
1238 :
1239 16 : IF (kpoints%nkp /= 0) THEN
1240 116 : DO ik = 1, SIZE(kpoints%kp_env)
1241 100 : kp => kpoints%kp_env(ik)%kpoint_env
1242 100 : wkp = kp%wkp
1243 216 : DO ispin = 1, nspins
1244 100 : mo_set => kp%mos(1, ispin)
1245 100 : CALL get_mo_set(mo_set=mo_set, nao=nao, nmo=nmo_kp)
1246 100 : eigenvalues => mo_set%eigenvalues
1247 100 : CALL cp_fm_get_info(mo_set%mo_coeff, matrix_struct=fm_struct_tmp)
1248 100 : IF (kpoints%use_real_wfn) THEN
1249 0 : CALL cp_fm_create(shalfc, fm_struct_tmp, name="shalfc")
1250 0 : CALL lowdin_kp_mo_coeff(kp, ispin, kpoints%use_real_wfn, shalfc=shalfc)
1251 : ELSE
1252 100 : CALL cp_cfm_create(cshalfc, fm_struct_tmp, name="cshalfc")
1253 100 : CALL lowdin_kp_mo_coeff(kp, ispin, kpoints%use_real_wfn, cshalfc=cshalfc)
1254 : END IF
1255 :
1256 1906 : DO imo = 1, nmo_kp
1257 1806 : IF (kpoints%use_real_wfn) THEN
1258 0 : CALL cp_fm_get_submatrix(shalfc, vecbuffer, 1, imo, nao, 1, transpose=.TRUE.)
1259 0 : ao_weight(:) = vecbuffer(1, 1:nao)**2
1260 : ELSE
1261 1806 : CALL cp_cfm_get_submatrix(cshalfc, zvecbuffer, 1, imo, nao, 1, transpose=.TRUE.)
1262 34902 : ao_weight(:) = REAL(CONJG(zvecbuffer(1, 1:nao))*zvecbuffer(1, 1:nao), KIND=dp)
1263 : END IF
1264 1806 : proj_weight = 0.0_dp
1265 : CALL accumulate_pdos_weights(ao_weight, ao_kind, ao_l, ao_comp, &
1266 1806 : separate_components, proj_weight)
1267 3712 : DO ikind = 1, nkind
1268 1806 : IF (kind_maxl(ikind) < 0) CYCLE
1269 3612 : IF (separate_components) THEN
1270 210 : DO icomp = 1, nsoset(kind_maxl(ikind))
1271 : CALL add_broadened_value(pdos_curve(:, icomp, ikind, ispin), &
1272 : emin, de, eigenvalues(imo), &
1273 : wkp*proj_weight(icomp, ikind), &
1274 210 : broaden_type, broaden_width, voigt_mixing)
1275 : END DO
1276 : ELSE
1277 7056 : DO icomp = 1, kind_maxl(ikind) + 1
1278 : CALL add_broadened_value(pdos_curve(:, icomp, ikind, ispin), &
1279 : emin, de, eigenvalues(imo), &
1280 : wkp*proj_weight(icomp, ikind), &
1281 7056 : broaden_type, broaden_width, voigt_mixing)
1282 : END DO
1283 : END IF
1284 : END DO
1285 : END DO
1286 :
1287 300 : IF (kpoints%use_real_wfn) THEN
1288 0 : CALL cp_fm_release(shalfc)
1289 : ELSE
1290 100 : CALL cp_cfm_release(cshalfc)
1291 : END IF
1292 : END DO
1293 : END DO
1294 : END IF
1295 16 : CALL para_env%sum(pdos_curve)
1296 :
1297 16 : IF (append .AND. iterstep > 1) THEN
1298 0 : my_pos = "APPEND"
1299 : ELSE
1300 16 : my_pos = "REWIND"
1301 : END IF
1302 16 : my_act = "WRITE"
1303 16 : spin(1) = "ALPHA"
1304 16 : spin(2) = "BETA"
1305 16 : IF (write_standard) THEN
1306 32 : DO ikind = 1, nkind
1307 16 : IF (kind_maxl(ikind) < 0) CYCLE
1308 16 : CALL get_atomic_kind(atomic_kind_set(ikind), name=kind_name)
1309 48 : DO ispin = 1, nspins
1310 16 : IF (nspins == 2) THEN
1311 0 : my_mittle = TRIM(spin(ispin))//"_k"//TRIM(ADJUSTL(cp_to_string(ikind)))
1312 : ELSE
1313 16 : my_mittle = "k"//TRIM(ADJUSTL(cp_to_string(ikind)))
1314 : END IF
1315 : CALL write_broadened_pdos_curve(logger, dft_section, TRIM(my_mittle), my_pos, my_act, &
1316 : iterstep, e_fermi(ispin), hoco(ispin), energy_ref(ispin), &
1317 : TRIM(zero_label), &
1318 : "K-point projected DOS for atomic kind "//TRIM(kind_name), &
1319 : kind_maxl(ikind), separate_components, &
1320 : pdos_curve(:, :, ikind, ispin), emin, de, &
1321 32 : broaden_type, broaden_width, voigt_mixing, ndigits, pdos_print_key=TRIM(my_print_key))
1322 : END DO
1323 : END DO
1324 : END IF
1325 :
1326 0 : DEALLOCATE (ao_comp, ao_kind, ao_l, ao_weight, kind_maxl, pdos_curve, proj_weight, &
1327 16 : vecbuffer, zvecbuffer)
1328 :
1329 16 : CALL timestop(handle)
1330 :
1331 112 : END SUBROUTINE calculate_projected_dos_kp
1332 :
1333 : ! **************************************************************************************************
1334 : !> \brief Build AO mapping arrays for PDOS accumulation.
1335 : !> \param qs_kind_set ...
1336 : !> \param particle_set ...
1337 : !> \param nao ...
1338 : !> \param ao_kind ...
1339 : !> \param ao_l ...
1340 : !> \param ao_comp ...
1341 : ! **************************************************************************************************
1342 16 : SUBROUTINE build_pdos_ao_map(qs_kind_set, particle_set, nao, ao_kind, ao_l, ao_comp)
1343 :
1344 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1345 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1346 : INTEGER, INTENT(IN) :: nao
1347 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: ao_kind, ao_l, ao_comp
1348 :
1349 : INTEGER :: iatom, ikind, irow, iset, ishell, iso, &
1350 : lshell, maxl, nset
1351 16 : INTEGER, DIMENSION(:), POINTER :: nshell
1352 16 : INTEGER, DIMENSION(:, :), POINTER :: l
1353 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1354 :
1355 80 : ALLOCATE (ao_kind(nao), ao_l(nao), ao_comp(nao))
1356 16 : irow = 0
1357 60 : DO iatom = 1, SIZE(particle_set)
1358 44 : NULLIFY (orb_basis_set)
1359 44 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1360 44 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1361 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1362 : nset=nset, &
1363 : nshell=nshell, &
1364 44 : l=l, maxl=maxl)
1365 204 : DO iset = 1, nset
1366 260 : DO ishell = 1, nshell(iset)
1367 116 : lshell = l(ishell, iset)
1368 532 : DO iso = 1, nso(lshell)
1369 316 : irow = irow + 1
1370 316 : CPASSERT(irow <= nao)
1371 316 : ao_kind(irow) = ikind
1372 316 : ao_l(irow) = lshell
1373 432 : ao_comp(irow) = nsoset(lshell - 1) + iso
1374 : END DO
1375 : END DO
1376 : END DO
1377 : END DO
1378 16 : CPASSERT(irow == nao)
1379 :
1380 16 : END SUBROUTINE build_pdos_ao_map
1381 :
1382 : ! **************************************************************************************************
1383 : !> \brief Accumulate AO weights into kind/l or kind/component projected weights.
1384 : !> \param ao_weight ...
1385 : !> \param ao_kind ...
1386 : !> \param ao_l ...
1387 : !> \param ao_comp ...
1388 : !> \param separate_components ...
1389 : !> \param proj_weight ...
1390 : ! **************************************************************************************************
1391 1806 : SUBROUTINE accumulate_pdos_weights(ao_weight, ao_kind, ao_l, ao_comp, &
1392 1806 : separate_components, proj_weight)
1393 :
1394 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: ao_weight
1395 : INTEGER, DIMENSION(:), INTENT(IN) :: ao_kind, ao_l, ao_comp
1396 : LOGICAL, INTENT(IN) :: separate_components
1397 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: proj_weight
1398 :
1399 : INTEGER :: iao, icomp, ikind
1400 :
1401 34902 : DO iao = 1, SIZE(ao_weight)
1402 33096 : ikind = ao_kind(iao)
1403 33096 : IF (separate_components) THEN
1404 1344 : icomp = ao_comp(iao)
1405 : ELSE
1406 31752 : icomp = ao_l(iao) + 1
1407 : END IF
1408 34902 : proj_weight(icomp, ikind) = proj_weight(icomp, ikind) + ao_weight(iao)
1409 : END DO
1410 :
1411 1806 : END SUBROUTINE accumulate_pdos_weights
1412 :
1413 : ! **************************************************************************************************
1414 : !> \brief Write a broadened k-point PDOS curve.
1415 : !> \param logger ...
1416 : !> \param dft_section ...
1417 : !> \param middle_name ...
1418 : !> \param file_position ...
1419 : !> \param file_action ...
1420 : !> \param iterstep ...
1421 : !> \param e_fermi ...
1422 : !> \param hoco ...
1423 : !> \param energy_ref ...
1424 : !> \param zero_label ...
1425 : !> \param title ...
1426 : !> \param maxl ...
1427 : !> \param separate_components ...
1428 : !> \param pdos_curve ...
1429 : !> \param emin ...
1430 : !> \param de ...
1431 : !> \param broaden_type ...
1432 : !> \param broaden_width ...
1433 : !> \param voigt_mixing ...
1434 : !> \param ndigits ...
1435 : !> \param pdos_print_key ...
1436 : ! **************************************************************************************************
1437 102 : SUBROUTINE write_broadened_pdos_curve(logger, dft_section, middle_name, file_position, &
1438 : file_action, iterstep, e_fermi, hoco, energy_ref, zero_label, &
1439 : title, maxl, &
1440 204 : separate_components, pdos_curve, emin, de, &
1441 : broaden_type, broaden_width, voigt_mixing, ndigits, pdos_print_key)
1442 :
1443 : TYPE(cp_logger_type), POINTER :: logger
1444 : TYPE(section_vals_type), POINTER :: dft_section
1445 : CHARACTER(LEN=*), INTENT(IN) :: middle_name, file_position, file_action
1446 : INTEGER, INTENT(IN) :: iterstep
1447 : REAL(KIND=dp), INTENT(IN) :: e_fermi, hoco, energy_ref
1448 : CHARACTER(LEN=*), INTENT(IN) :: zero_label, title
1449 : INTEGER, INTENT(IN) :: maxl
1450 : LOGICAL, INTENT(IN) :: separate_components
1451 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: pdos_curve
1452 : REAL(KIND=dp), INTENT(IN) :: emin, de
1453 : INTEGER, INTENT(IN) :: broaden_type
1454 : REAL(KIND=dp), INTENT(IN) :: broaden_width, voigt_mixing
1455 : INTEGER, INTENT(IN) :: ndigits
1456 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: pdos_print_key
1457 :
1458 : CHARACTER(LEN=16) :: energy_label
1459 : CHARACTER(LEN=20) :: fmtstr_data
1460 102 : CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:, :, :) :: tmp_str
1461 : CHARACTER(LEN=default_string_length) :: my_print_key
1462 : INTEGER :: energy_unit, i, icomp, il, im, iw, &
1463 : ncomponents, nhist
1464 : REAL(KIND=dp) :: density_factor, energy_factor, &
1465 : ev_factor, eval
1466 :
1467 102 : my_print_key = "PRINT%PDOS"
1468 102 : IF (PRESENT(pdos_print_key)) my_print_key = TRIM(pdos_print_key)
1469 :
1470 102 : nhist = SIZE(pdos_curve, 1)
1471 102 : IF (separate_components) THEN
1472 40 : ncomponents = nsoset(maxl)
1473 : ELSE
1474 62 : ncomponents = maxl + 1
1475 : END IF
1476 102 : CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%ENERGY_UNIT", i_val=energy_unit)
1477 102 : energy_factor = dos_energy_scale(energy_unit)
1478 102 : density_factor = dos_density_scale(energy_unit)
1479 102 : energy_label = dos_energy_label(energy_unit)
1480 102 : ev_factor = dos_energy_scale(dos_energy_unit_ev)
1481 : iw = cp_print_key_unit_nr(logger, dft_section, TRIM(my_print_key), &
1482 : extension=".pdos", file_position=file_position, file_action=file_action, &
1483 102 : file_form="FORMATTED", middle_name=TRIM(middle_name))
1484 102 : IF (iw > 0) THEN
1485 51 : WRITE (UNIT=iw, FMT="(A,I0)") "# "//TRIM(title)//" at iteration step i = ", iterstep
1486 : WRITE (UNIT=iw, FMT="(A,F12.6,A,F12.6,A)") &
1487 51 : "# E(Fermi) = ", e_fermi, " a.u. = ", e_fermi*ev_factor, " eV"
1488 : WRITE (UNIT=iw, FMT="(A,F12.6,A,F12.6,A)") &
1489 51 : "# E(HOCO) = ", hoco, " a.u. = ", hoco*ev_factor, " eV"
1490 51 : WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label)
1491 51 : CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
1492 51 : WRITE (UNIT=iw, FMT="(A)", ADVANCE="NO") "# "//TRIM(energy_label)
1493 51 : WRITE (UNIT=iw, FMT="(2X,A)", ADVANCE="NO") "total"
1494 51 : IF (separate_components) THEN
1495 80 : ALLOCATE (tmp_str(0:0, 0:maxl, -maxl:maxl))
1496 588 : tmp_str = ""
1497 73 : DO il = 0, maxl
1498 220 : DO im = -il, il
1499 147 : tmp_str(0, il, im) = sgf_symbol(0, il, im)
1500 200 : WRITE (UNIT=iw, FMT="(2X,A)", ADVANCE="NO") TRIM(tmp_str(0, il, im))
1501 : END DO
1502 : END DO
1503 20 : DEALLOCATE (tmp_str)
1504 : ELSE
1505 111 : DO il = 0, maxl
1506 111 : WRITE (UNIT=iw, FMT="(2X,A)", ADVANCE="NO") TRIM(l_sym(il))
1507 : END DO
1508 : END IF
1509 51 : WRITE (UNIT=iw, FMT="()")
1510 51 : WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(2X,F20.", ndigits, ")"
1511 105398 : DO i = 1, nhist
1512 105347 : eval = (emin + (i - 1)*de - energy_ref)*energy_factor
1513 105347 : WRITE (UNIT=iw, FMT="(F15.8)", ADVANCE="NO") eval
1514 510849 : WRITE (UNIT=iw, FMT=fmtstr_data, ADVANCE="NO") SUM(pdos_curve(i, 1:ncomponents))*density_factor
1515 510849 : DO icomp = 1, ncomponents
1516 510849 : WRITE (UNIT=iw, FMT=fmtstr_data, ADVANCE="NO") pdos_curve(i, icomp)*density_factor
1517 : END DO
1518 105398 : WRITE (UNIT=iw, FMT="()")
1519 : END DO
1520 : END IF
1521 102 : CALL cp_print_key_finished_output(iw, logger, dft_section, TRIM(my_print_key))
1522 :
1523 102 : END SUBROUTINE write_broadened_pdos_curve
1524 :
1525 : ! **************************************************************************************************
1526 : !> \brief Write a broadened PDOS curve for a projected weight matrix.
1527 : !> \param logger ...
1528 : !> \param dft_section ...
1529 : !> \param middle_name ...
1530 : !> \param file_position ...
1531 : !> \param file_action ...
1532 : !> \param iterstep ...
1533 : !> \param e_fermi ...
1534 : !> \param hoco ...
1535 : !> \param energy_ref ...
1536 : !> \param zero_label ...
1537 : !> \param title ...
1538 : !> \param maxl ...
1539 : !> \param separate_components ...
1540 : !> \param weights ...
1541 : !> \param eigenvalues ...
1542 : !> \param nstates ...
1543 : !> \param de ...
1544 : !> \param broaden_type ...
1545 : !> \param broaden_width ...
1546 : !> \param voigt_mixing ...
1547 : !> \param ndigits ...
1548 : !> \param pdos_print_key ...
1549 : ! **************************************************************************************************
1550 86 : SUBROUTINE write_broadened_pdos(logger, dft_section, middle_name, file_position, file_action, &
1551 86 : iterstep, e_fermi, hoco, energy_ref, zero_label, title, maxl, separate_components, weights, &
1552 86 : eigenvalues, nstates, de, broaden_type, broaden_width, &
1553 : voigt_mixing, ndigits, pdos_print_key)
1554 :
1555 : TYPE(cp_logger_type), POINTER :: logger
1556 : TYPE(section_vals_type), POINTER :: dft_section
1557 : CHARACTER(LEN=*), INTENT(IN) :: middle_name, file_position, file_action
1558 : INTEGER, INTENT(IN) :: iterstep
1559 : REAL(KIND=dp), INTENT(IN) :: e_fermi, hoco, energy_ref
1560 : CHARACTER(LEN=*), INTENT(IN) :: zero_label, title
1561 : INTEGER, INTENT(IN) :: maxl
1562 : LOGICAL, INTENT(IN) :: separate_components
1563 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: weights
1564 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: eigenvalues
1565 : INTEGER, INTENT(IN) :: nstates
1566 : REAL(KIND=dp), INTENT(IN) :: de
1567 : INTEGER, INTENT(IN) :: broaden_type
1568 : REAL(KIND=dp), INTENT(IN) :: broaden_width, voigt_mixing
1569 : INTEGER, INTENT(IN) :: ndigits
1570 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: pdos_print_key
1571 :
1572 : INTEGER :: i, icomp, imo, ncomponents, nhist
1573 : REAL(KIND=dp) :: cutoff, emax, emin, eval, line_shape
1574 86 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pdos_curve
1575 :
1576 86 : IF (broaden_width <= 0.0_dp) RETURN
1577 :
1578 86 : ncomponents = SIZE(weights, 1)
1579 86 : cutoff = broadening_cutoff(broaden_type, broaden_width)
1580 718 : emin = MINVAL(eigenvalues(1:nstates)) - cutoff
1581 718 : emax = MAXVAL(eigenvalues(1:nstates)) + cutoff
1582 86 : nhist = NINT((emax - emin)/de) + 1
1583 344 : ALLOCATE (pdos_curve(nhist, ncomponents))
1584 86 : pdos_curve = 0.0_dp
1585 :
1586 718 : DO imo = 1, nstates
1587 30640 : DO i = MAX(1, FLOOR((eigenvalues(imo) - cutoff - emin)/de) + 1), &
1588 718 : MIN(nhist, CEILING((eigenvalues(imo) + cutoff - emin)/de) + 1)
1589 30008 : eval = emin + (i - 1)*de
1590 : line_shape = broadening_function(eval - eigenvalues(imo), broaden_type, broaden_width, &
1591 30008 : voigt_mixing)
1592 161716 : DO icomp = 1, ncomponents
1593 161084 : pdos_curve(i, icomp) = pdos_curve(i, icomp) + weights(icomp, imo)*line_shape
1594 : END DO
1595 : END DO
1596 : END DO
1597 :
1598 : CALL write_broadened_pdos_curve(logger, dft_section, middle_name, file_position, file_action, &
1599 : iterstep, e_fermi, hoco, energy_ref, zero_label, title, maxl, separate_components, &
1600 : pdos_curve, emin, de, broaden_type, broaden_width, &
1601 86 : voigt_mixing, ndigits, pdos_print_key=pdos_print_key)
1602 :
1603 86 : DEALLOCATE (pdos_curve)
1604 :
1605 : END SUBROUTINE write_broadened_pdos
1606 :
1607 0 : END MODULE qs_pdos
|
