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