Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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_blacs_env, ONLY: cp_blacs_env_type
24 : USE cp_control_types, ONLY: dft_control_type
25 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
26 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
27 : USE cp_fm_diag, ONLY: cp_fm_power
28 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
29 : cp_fm_struct_release,&
30 : cp_fm_struct_type
31 : USE cp_fm_types, ONLY: cp_fm_create,&
32 : cp_fm_get_info,&
33 : cp_fm_get_submatrix,&
34 : cp_fm_init_random,&
35 : cp_fm_release,&
36 : cp_fm_type
37 : USE cp_log_handling, ONLY: cp_get_default_logger,&
38 : cp_logger_get_default_io_unit,&
39 : cp_logger_type,&
40 : cp_to_string
41 : USE cp_output_handling, ONLY: cp_p_file,&
42 : cp_print_key_finished_output,&
43 : cp_print_key_should_output,&
44 : cp_print_key_unit_nr
45 : USE input_section_types, ONLY: section_vals_get,&
46 : section_vals_get_subs_vals,&
47 : section_vals_type,&
48 : section_vals_val_get
49 : USE kinds, ONLY: default_string_length,&
50 : dp
51 : USE memory_utilities, ONLY: reallocate
52 : USE message_passing, ONLY: mp_para_env_type
53 : USE orbital_pointers, ONLY: nso,&
54 : nsoset
55 : USE orbital_symbols, ONLY: l_sym,&
56 : sgf_symbol
57 : USE parallel_gemm_api, ONLY: parallel_gemm
58 : USE particle_types, ONLY: particle_type
59 : USE preconditioner_types, ONLY: preconditioner_type
60 : USE pw_env_types, ONLY: pw_env_get,&
61 : pw_env_type
62 : USE pw_pool_types, ONLY: pw_pool_p_type,&
63 : pw_pool_type
64 : USE pw_types, ONLY: pw_c1d_gs_type,&
65 : pw_r3d_rs_type
66 : USE qs_collocate_density, ONLY: calculate_wavefunction
67 : USE qs_environment_types, ONLY: get_qs_env,&
68 : qs_environment_type
69 : USE qs_kind_types, ONLY: get_qs_kind,&
70 : get_qs_kind_set,&
71 : qs_kind_type
72 : USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
73 : USE qs_mo_types, ONLY: get_mo_set,&
74 : mo_set_type
75 : USE qs_ot_eigensolver, ONLY: ot_eigensolver
76 : USE scf_control_types, ONLY: scf_control_type
77 : #include "./base/base_uses.f90"
78 :
79 : IMPLICIT NONE
80 :
81 : PRIVATE
82 :
83 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_pdos'
84 :
85 : ! **************************************************************************************************
86 : ! *** Public subroutines ***
87 :
88 : PUBLIC :: calculate_projected_dos
89 :
90 : TYPE ldos_type
91 : INTEGER :: maxl = -1, nlist = -1
92 : LOGICAL :: separate_components = .FALSE.
93 : INTEGER, DIMENSION(:), POINTER :: list_index => NULL()
94 : REAL(KIND=dp), DIMENSION(:, :), &
95 : POINTER :: pdos_array => NULL()
96 : END TYPE ldos_type
97 :
98 : TYPE r_ldos_type
99 : INTEGER :: nlist = -1, npoints = -1
100 : INTEGER, DIMENSION(:, :), POINTER :: index_grid_local => NULL()
101 : INTEGER, DIMENSION(:), POINTER :: list_index => NULL()
102 : REAL(KIND=dp), DIMENSION(:), POINTER :: x_range => NULL(), y_range => NULL(), z_range => NULL()
103 : REAL(KIND=dp), DIMENSION(:), POINTER :: eval_range => NULL()
104 : REAL(KIND=dp), DIMENSION(:), &
105 : POINTER :: pdos_array => NULL()
106 : END TYPE r_ldos_type
107 :
108 : TYPE ldos_p_type
109 : TYPE(ldos_type), POINTER :: ldos => NULL()
110 : END TYPE ldos_p_type
111 :
112 : TYPE r_ldos_p_type
113 : TYPE(r_ldos_type), POINTER :: ldos => NULL()
114 : END TYPE r_ldos_p_type
115 : CONTAINS
116 :
117 : ! **************************************************************************************************
118 : !> \brief Compute and write projected density of states
119 : !> \param mo_set ...
120 : !> \param atomic_kind_set ...
121 : !> \param qs_kind_set ...
122 : !> \param particle_set ...
123 : !> \param qs_env ...
124 : !> \param dft_section ...
125 : !> \param ispin ...
126 : !> \param xas_mittle ...
127 : !> \param external_matrix_shalf ...
128 : !> \date 26.02.2008
129 : !> \par History:
130 : !> - Added optional external matrix_shalf to avoid recomputing it (A. Bussy, 09.2019)
131 : !> \par Variables
132 : !> -
133 : !> -
134 : !> \author MI
135 : !> \version 1.0
136 : ! **************************************************************************************************
137 60 : SUBROUTINE calculate_projected_dos(mo_set, atomic_kind_set, qs_kind_set, particle_set, qs_env, &
138 : dft_section, ispin, xas_mittle, external_matrix_shalf)
139 :
140 : TYPE(mo_set_type), INTENT(IN) :: mo_set
141 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
142 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
143 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
144 : TYPE(qs_environment_type), POINTER :: qs_env
145 : TYPE(section_vals_type), POINTER :: dft_section
146 : INTEGER, INTENT(IN), OPTIONAL :: ispin
147 : CHARACTER(LEN=default_string_length), INTENT(IN), &
148 : OPTIONAL :: xas_mittle
149 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL, TARGET :: external_matrix_shalf
150 :
151 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_projected_dos'
152 :
153 : CHARACTER(LEN=16) :: fmtstr2
154 : CHARACTER(LEN=27) :: fmtstr1
155 60 : CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:, :, :) :: tmp_str
156 : CHARACTER(LEN=default_string_length) :: kind_name, my_act, my_mittle, my_pos, &
157 : spin(2)
158 : CHARACTER(LEN=default_string_length), &
159 60 : ALLOCATABLE, DIMENSION(:) :: ldos_index, r_ldos_index
160 : INTEGER :: handle, homo, i, iatom, ikind, il, ildos, im, imo, in_x, in_y, in_z, ir, irow, &
161 : iset, isgf, ishell, iso, iterstep, iw, j, jx, jy, jz, k, lcomponent, lshell, maxl, &
162 : maxlgto, my_spin, n_dependent, n_r_ldos, n_rep, nao, natom, ncol_global, nkind, nldos, &
163 : nlumo, nmo, np_tot, npoints, nrow_global, nset, nsgf, nvirt, out_each, output_unit
164 60 : INTEGER, ALLOCATABLE, DIMENSION(:) :: firstrow
165 60 : INTEGER, DIMENSION(:), POINTER :: list, nshell
166 60 : INTEGER, DIMENSION(:, :), POINTER :: bo, l
167 : LOGICAL :: append, calc_matsh, do_ldos, do_r_ldos, &
168 : do_virt, ionode, separate_components, &
169 : should_output
170 60 : LOGICAL, DIMENSION(:, :), POINTER :: read_r
171 : REAL(KIND=dp) :: dh(3, 3), dvol, e_fermi, r(3), r_vec(3), &
172 : ratom(3)
173 60 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, evals_virt, &
174 60 : occupation_numbers
175 60 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
176 60 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pdos_array
177 : TYPE(cell_type), POINTER :: cell
178 : TYPE(cp_blacs_env_type), POINTER :: context
179 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
180 : TYPE(cp_fm_type) :: matrix_shalfc, matrix_work, mo_virt
181 : TYPE(cp_fm_type), POINTER :: matrix_shalf, mo_coeff
182 : TYPE(cp_logger_type), POINTER :: logger
183 60 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: s_matrix
184 : TYPE(dft_control_type), POINTER :: dft_control
185 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
186 60 : TYPE(ldos_p_type), DIMENSION(:), POINTER :: ldos_p
187 : TYPE(mp_para_env_type), POINTER :: para_env
188 : TYPE(pw_c1d_gs_type) :: wf_g
189 : TYPE(pw_env_type), POINTER :: pw_env
190 60 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
191 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
192 : TYPE(pw_r3d_rs_type) :: wf_r
193 60 : TYPE(r_ldos_p_type), DIMENSION(:), POINTER :: r_ldos_p
194 : TYPE(section_vals_type), POINTER :: ldos_section
195 :
196 60 : NULLIFY (logger)
197 120 : logger => cp_get_default_logger()
198 : ionode = logger%para_env%is_source()
199 : should_output = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
200 60 : "PRINT%PDOS"), cp_p_file)
201 60 : output_unit = cp_logger_get_default_io_unit(logger)
202 :
203 60 : spin(1) = "ALPHA"
204 60 : spin(2) = "BETA"
205 60 : IF ((.NOT. should_output)) RETURN
206 :
207 60 : NULLIFY (context, s_matrix, orb_basis_set, para_env, pdos_array)
208 60 : NULLIFY (eigenvalues, fm_struct_tmp, mo_coeff, vecbuffer)
209 60 : NULLIFY (ldos_section, list, cell, pw_env, auxbas_pw_pool, evals_virt)
210 60 : NULLIFY (occupation_numbers, ldos_p, r_ldos_p, dft_control, occupation_numbers)
211 :
212 60 : CALL timeset(routineN, handle)
213 60 : iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
214 :
215 60 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,I10))') &
216 30 : " Calculate PDOS at iteration step ", iterstep
217 : CALL get_qs_env(qs_env=qs_env, &
218 60 : matrix_s=s_matrix)
219 :
220 60 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
221 60 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf, maxlgto=maxlgto)
222 60 : nkind = SIZE(atomic_kind_set)
223 :
224 : CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, homo=homo, nao=nao, nmo=nmo, &
225 60 : mu=e_fermi)
226 : CALL cp_fm_get_info(mo_coeff, &
227 : context=context, para_env=para_env, &
228 : nrow_global=nrow_global, &
229 60 : ncol_global=ncol_global)
230 :
231 60 : CALL section_vals_val_get(dft_section, "PRINT%PDOS%OUT_EACH_MO", i_val=out_each)
232 60 : IF (out_each == -1) out_each = nao + 1
233 60 : CALL section_vals_val_get(dft_section, "PRINT%PDOS%NLUMO", i_val=nlumo)
234 60 : IF (nlumo == -1) nlumo = nao - homo
235 60 : do_virt = (nlumo > (nmo - homo))
236 60 : nvirt = nlumo - (nmo - homo)
237 : ! Generate virtual orbitals
238 60 : IF (do_virt) THEN
239 14 : IF (PRESENT(ispin)) THEN
240 8 : my_spin = ispin
241 : ELSE
242 6 : my_spin = 1
243 : END IF
244 :
245 14 : CALL generate_virtual_mo(qs_env, mo_set, evals_virt, mo_virt, nvirt, ispin=my_spin)
246 : ELSE
247 46 : NULLIFY (evals_virt)
248 46 : nvirt = 0
249 : END IF
250 :
251 60 : calc_matsh = .TRUE.
252 60 : IF (PRESENT(external_matrix_shalf)) calc_matsh = .FALSE.
253 :
254 : ! Create S^1/2 : from sparse to full matrix, if no external available
255 116 : IF (calc_matsh) THEN
256 58 : NULLIFY (matrix_shalf)
257 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
258 58 : nrow_global=nrow_global, ncol_global=nrow_global)
259 58 : ALLOCATE (matrix_shalf)
260 58 : CALL cp_fm_create(matrix_shalf, fm_struct_tmp, name="matrix_shalf")
261 58 : CALL cp_fm_create(matrix_work, fm_struct_tmp, name="matrix_work")
262 58 : CALL cp_fm_struct_release(fm_struct_tmp)
263 58 : CALL copy_dbcsr_to_fm(s_matrix(1)%matrix, matrix_shalf)
264 58 : CALL cp_fm_power(matrix_shalf, matrix_work, 0.5_dp, EPSILON(0.0_dp), n_dependent)
265 58 : CALL cp_fm_release(matrix_work)
266 : ELSE
267 : matrix_shalf => external_matrix_shalf
268 : END IF
269 :
270 : ! Multiply S^(1/2) time the mOS coefficients to get orthonormalized MOS
271 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
272 60 : nrow_global=nrow_global, ncol_global=ncol_global)
273 60 : CALL cp_fm_create(matrix_shalfc, fm_struct_tmp, name="matrix_shalfc")
274 : CALL parallel_gemm("N", "N", nrow_global, ncol_global, nrow_global, &
275 60 : 1.0_dp, matrix_shalf, mo_coeff, 0.0_dp, matrix_shalfc)
276 60 : CALL cp_fm_struct_release(fm_struct_tmp)
277 :
278 60 : IF (do_virt) THEN
279 14 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T14,I10,T27,A))') &
280 7 : " Compute ", nvirt, " additional unoccupied KS orbitals"
281 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
282 14 : nrow_global=nrow_global, ncol_global=nvirt)
283 14 : CALL cp_fm_create(matrix_work, fm_struct_tmp, name="matrix_shalfc")
284 : CALL parallel_gemm("N", "N", nrow_global, nvirt, nrow_global, &
285 14 : 1.0_dp, matrix_shalf, mo_virt, 0.0_dp, matrix_work)
286 14 : CALL cp_fm_struct_release(fm_struct_tmp)
287 : END IF
288 :
289 60 : IF (calc_matsh) THEN
290 58 : CALL cp_fm_release(matrix_shalf)
291 58 : DEALLOCATE (matrix_shalf)
292 : END IF
293 : ! Array to store the PDOS per kind and angular momentum
294 60 : do_ldos = .FALSE.
295 60 : ldos_section => section_vals_get_subs_vals(dft_section, "PRINT%PDOS%LDOS")
296 :
297 60 : CALL section_vals_get(ldos_section, n_repetition=nldos)
298 60 : IF (nldos > 0) THEN
299 8 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,I10))') &
300 4 : " Prepare the list of atoms for LDOS. Number of lists: ", nldos
301 8 : do_ldos = .TRUE.
302 44 : ALLOCATE (ldos_p(nldos))
303 24 : ALLOCATE (ldos_index(nldos))
304 28 : DO ildos = 1, nldos
305 20 : WRITE (ldos_index(ildos), '(I0)') ildos
306 20 : ALLOCATE (ldos_p(ildos)%ldos)
307 : NULLIFY (ldos_p(ildos)%ldos%pdos_array)
308 : NULLIFY (ldos_p(ildos)%ldos%list_index)
309 :
310 20 : CALL section_vals_val_get(ldos_section, "LIST", i_rep_section=ildos, n_rep_val=n_rep)
311 20 : IF (n_rep > 0) THEN
312 20 : ldos_p(ildos)%ldos%nlist = 0
313 40 : DO ir = 1, n_rep
314 20 : NULLIFY (list)
315 : CALL section_vals_val_get(ldos_section, "LIST", i_rep_section=ildos, i_rep_val=ir, &
316 20 : i_vals=list)
317 40 : IF (ASSOCIATED(list)) THEN
318 20 : CALL reallocate(ldos_p(ildos)%ldos%list_index, 1, ldos_p(ildos)%ldos%nlist + SIZE(list))
319 76 : DO i = 1, SIZE(list)
320 76 : ldos_p(ildos)%ldos%list_index(i + ldos_p(ildos)%ldos%nlist) = list(i)
321 : END DO
322 20 : ldos_p(ildos)%ldos%nlist = ldos_p(ildos)%ldos%nlist + SIZE(list)
323 : END IF
324 : END DO
325 : ELSE
326 : ! stop, LDOS without list of atoms is not implemented
327 : END IF
328 :
329 20 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='((T10,A,T18,I6,T25,A,T36,I10,A))') &
330 10 : " List ", ildos, " contains ", ldos_p(ildos)%ldos%nlist, " atoms"
331 : CALL section_vals_val_get(ldos_section, "COMPONENTS", i_rep_section=ildos, &
332 20 : l_val=ldos_p(ildos)%ldos%separate_components)
333 20 : IF (ldos_p(ildos)%ldos%separate_components) THEN
334 16 : ALLOCATE (ldos_p(ildos)%ldos%pdos_array(nsoset(maxlgto), nmo + nvirt))
335 : ELSE
336 64 : ALLOCATE (ldos_p(ildos)%ldos%pdos_array(0:maxlgto, nmo + nvirt))
337 : END IF
338 716 : ldos_p(ildos)%ldos%pdos_array = 0.0_dp
339 48 : ldos_p(ildos)%ldos%maxl = -1
340 :
341 : END DO
342 : END IF
343 :
344 60 : do_r_ldos = .FALSE.
345 60 : ldos_section => section_vals_get_subs_vals(dft_section, "PRINT%PDOS%R_LDOS")
346 60 : CALL section_vals_get(ldos_section, n_repetition=n_r_ldos)
347 60 : IF (n_r_ldos > 0) THEN
348 0 : do_r_ldos = .TRUE.
349 0 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,I10))') &
350 0 : " Prepare the list of points for R_LDOS. Number of lists: ", n_r_ldos
351 0 : ALLOCATE (r_ldos_p(n_r_ldos))
352 0 : ALLOCATE (r_ldos_index(n_r_ldos))
353 : CALL get_qs_env(qs_env=qs_env, &
354 : cell=cell, &
355 : dft_control=dft_control, &
356 0 : pw_env=pw_env)
357 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
358 0 : pw_pools=pw_pools)
359 :
360 0 : CALL auxbas_pw_pool%create_pw(wf_r)
361 0 : CALL auxbas_pw_pool%create_pw(wf_g)
362 0 : ALLOCATE (read_r(4, n_r_ldos))
363 0 : DO ildos = 1, n_r_ldos
364 0 : WRITE (r_ldos_index(ildos), '(I0)') ildos
365 0 : ALLOCATE (r_ldos_p(ildos)%ldos)
366 : NULLIFY (r_ldos_p(ildos)%ldos%pdos_array)
367 : NULLIFY (r_ldos_p(ildos)%ldos%list_index)
368 :
369 0 : CALL section_vals_val_get(ldos_section, "LIST", i_rep_section=ildos, n_rep_val=n_rep)
370 0 : IF (n_rep > 0) THEN
371 0 : r_ldos_p(ildos)%ldos%nlist = 0
372 0 : DO ir = 1, n_rep
373 0 : NULLIFY (list)
374 : CALL section_vals_val_get(ldos_section, "LIST", i_rep_section=ildos, i_rep_val=ir, &
375 0 : i_vals=list)
376 0 : IF (ASSOCIATED(list)) THEN
377 0 : CALL reallocate(r_ldos_p(ildos)%ldos%list_index, 1, r_ldos_p(ildos)%ldos%nlist + SIZE(list))
378 0 : DO i = 1, SIZE(list)
379 0 : r_ldos_p(ildos)%ldos%list_index(i + r_ldos_p(ildos)%ldos%nlist) = list(i)
380 : END DO
381 0 : r_ldos_p(ildos)%ldos%nlist = r_ldos_p(ildos)%ldos%nlist + SIZE(list)
382 : END IF
383 : END DO
384 : ELSE
385 : ! stop, LDOS without list of atoms is not implemented
386 : END IF
387 :
388 0 : ALLOCATE (r_ldos_p(ildos)%ldos%pdos_array(nmo + nvirt))
389 0 : r_ldos_p(ildos)%ldos%pdos_array = 0.0_dp
390 0 : read_r(1:3, ildos) = .FALSE.
391 0 : CALL section_vals_val_get(ldos_section, "XRANGE", i_rep_section=ildos, explicit=read_r(1, ildos))
392 0 : IF (read_r(1, ildos)) THEN
393 : CALL section_vals_val_get(ldos_section, "XRANGE", i_rep_section=ildos, r_vals= &
394 0 : r_ldos_p(ildos)%ldos%x_range)
395 : ELSE
396 0 : ALLOCATE (r_ldos_p(ildos)%ldos%x_range(2))
397 0 : r_ldos_p(ildos)%ldos%x_range = 0.0_dp
398 : END IF
399 0 : CALL section_vals_val_get(ldos_section, "YRANGE", i_rep_section=ildos, explicit=read_r(2, ildos))
400 0 : IF (read_r(2, ildos)) THEN
401 : CALL section_vals_val_get(ldos_section, "YRANGE", i_rep_section=ildos, r_vals= &
402 0 : r_ldos_p(ildos)%ldos%y_range)
403 : ELSE
404 0 : ALLOCATE (r_ldos_p(ildos)%ldos%y_range(2))
405 0 : r_ldos_p(ildos)%ldos%y_range = 0.0_dp
406 : END IF
407 0 : CALL section_vals_val_get(ldos_section, "ZRANGE", i_rep_section=ildos, explicit=read_r(3, ildos))
408 0 : IF (read_r(3, ildos)) THEN
409 : CALL section_vals_val_get(ldos_section, "ZRANGE", i_rep_section=ildos, r_vals= &
410 0 : r_ldos_p(ildos)%ldos%z_range)
411 : ELSE
412 0 : ALLOCATE (r_ldos_p(ildos)%ldos%z_range(2))
413 0 : r_ldos_p(ildos)%ldos%z_range = 0.0_dp
414 : END IF
415 :
416 0 : CALL section_vals_val_get(ldos_section, "ERANGE", i_rep_section=ildos, explicit=read_r(4, ildos))
417 0 : IF (read_r(4, ildos)) THEN
418 : CALL section_vals_val_get(ldos_section, "ERANGE", i_rep_section=ildos, &
419 0 : r_vals=r_ldos_p(ildos)%ldos%eval_range)
420 : ELSE
421 0 : ALLOCATE (r_ldos_p(ildos)%ldos%eval_range(2))
422 0 : r_ldos_p(ildos)%ldos%eval_range(1) = -HUGE(0.0_dp)
423 0 : r_ldos_p(ildos)%ldos%eval_range(2) = +HUGE(0.0_dp)
424 : END IF
425 :
426 0 : bo => wf_r%pw_grid%bounds_local
427 0 : dh = wf_r%pw_grid%dh
428 0 : dvol = wf_r%pw_grid%dvol
429 0 : np_tot = wf_r%pw_grid%npts(1)*wf_r%pw_grid%npts(2)*wf_r%pw_grid%npts(3)
430 0 : ALLOCATE (r_ldos_p(ildos)%ldos%index_grid_local(3, np_tot))
431 :
432 0 : r_ldos_p(ildos)%ldos%npoints = 0
433 0 : DO jz = bo(1, 3), bo(2, 3)
434 0 : DO jy = bo(1, 2), bo(2, 2)
435 0 : DO jx = bo(1, 1), bo(2, 1)
436 : !compute the position of the grid point
437 0 : i = jx - wf_r%pw_grid%bounds(1, 1)
438 0 : j = jy - wf_r%pw_grid%bounds(1, 2)
439 0 : k = jz - wf_r%pw_grid%bounds(1, 3)
440 0 : r(3) = k*dh(3, 3) + j*dh(3, 2) + i*dh(3, 1)
441 0 : r(2) = k*dh(2, 3) + j*dh(2, 2) + i*dh(2, 1)
442 0 : r(1) = k*dh(1, 3) + j*dh(1, 2) + i*dh(1, 1)
443 :
444 0 : DO il = 1, r_ldos_p(ildos)%ldos%nlist
445 0 : iatom = r_ldos_p(ildos)%ldos%list_index(il)
446 0 : ratom = particle_set(iatom)%r
447 0 : r_vec = pbc(ratom, r, cell)
448 0 : IF (cell%orthorhombic) THEN
449 0 : IF (cell%perd(1) == 0) r_vec(1) = MODULO(r_vec(1), cell%hmat(1, 1))
450 0 : IF (cell%perd(2) == 0) r_vec(2) = MODULO(r_vec(2), cell%hmat(2, 2))
451 0 : IF (cell%perd(3) == 0) r_vec(3) = MODULO(r_vec(3), cell%hmat(3, 3))
452 : END IF
453 :
454 0 : in_x = 0
455 0 : in_y = 0
456 0 : in_z = 0
457 0 : IF (r_ldos_p(ildos)%ldos%x_range(1) /= 0.0_dp) THEN
458 0 : IF (r_vec(1) > r_ldos_p(ildos)%ldos%x_range(1) .AND. &
459 : r_vec(1) < r_ldos_p(ildos)%ldos%x_range(2)) THEN
460 0 : in_x = 1
461 : END IF
462 : ELSE
463 : in_x = 1
464 : END IF
465 0 : IF (r_ldos_p(ildos)%ldos%y_range(1) /= 0.0_dp) THEN
466 0 : IF (r_vec(2) > r_ldos_p(ildos)%ldos%y_range(1) .AND. &
467 : r_vec(2) < r_ldos_p(ildos)%ldos%y_range(2)) THEN
468 0 : in_y = 1
469 : END IF
470 : ELSE
471 : in_y = 1
472 : END IF
473 0 : IF (r_ldos_p(ildos)%ldos%z_range(1) /= 0.0_dp) THEN
474 0 : IF (r_vec(3) > r_ldos_p(ildos)%ldos%z_range(1) .AND. &
475 : r_vec(3) < r_ldos_p(ildos)%ldos%z_range(2)) THEN
476 0 : in_z = 1
477 : END IF
478 : ELSE
479 : in_z = 1
480 : END IF
481 0 : IF (in_x*in_y*in_z > 0) THEN
482 0 : r_ldos_p(ildos)%ldos%npoints = r_ldos_p(ildos)%ldos%npoints + 1
483 0 : r_ldos_p(ildos)%ldos%index_grid_local(1, r_ldos_p(ildos)%ldos%npoints) = jx
484 0 : r_ldos_p(ildos)%ldos%index_grid_local(2, r_ldos_p(ildos)%ldos%npoints) = jy
485 0 : r_ldos_p(ildos)%ldos%index_grid_local(3, r_ldos_p(ildos)%ldos%npoints) = jz
486 0 : EXIT
487 : END IF
488 : END DO
489 : END DO
490 : END DO
491 : END DO
492 0 : CALL reallocate(r_ldos_p(ildos)%ldos%index_grid_local, 1, 3, 1, r_ldos_p(ildos)%ldos%npoints)
493 0 : npoints = r_ldos_p(ildos)%ldos%npoints
494 0 : CALL para_env%sum(npoints)
495 0 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='((T10,A,T18,I6,T25,A,T36,I10,A))') &
496 0 : " List ", ildos, " contains ", npoints, " grid points"
497 : END DO
498 : END IF
499 :
500 60 : CALL section_vals_val_get(dft_section, "PRINT%PDOS%COMPONENTS", l_val=separate_components)
501 60 : IF (separate_components) THEN
502 80 : ALLOCATE (pdos_array(nsoset(maxlgto), nkind, nmo + nvirt))
503 : ELSE
504 220 : ALLOCATE (pdos_array(0:maxlgto, nkind, nmo + nvirt))
505 : END IF
506 60 : IF (do_virt) THEN
507 42 : ALLOCATE (eigenvalues(nmo + nvirt))
508 122 : eigenvalues(1:nmo) = mo_set%eigenvalues(1:nmo)
509 346 : eigenvalues(nmo + 1:nmo + nvirt) = evals_virt(1:nvirt)
510 28 : ALLOCATE (occupation_numbers(nmo + nvirt))
511 288 : occupation_numbers(:) = 0.0_dp
512 122 : occupation_numbers(1:nmo) = mo_set%occupation_numbers(1:nmo)
513 : ELSE
514 46 : eigenvalues => mo_set%eigenvalues
515 46 : occupation_numbers => mo_set%occupation_numbers
516 : END IF
517 :
518 8036 : pdos_array = 0.0_dp
519 60 : nao = mo_set%nao
520 180 : ALLOCATE (vecbuffer(1, nao))
521 2432 : vecbuffer = 0.0_dp
522 180 : ALLOCATE (firstrow(natom))
523 224 : firstrow = 0
524 :
525 : !Adjust energy range for r_ldos
526 60 : DO ildos = 1, n_r_ldos
527 0 : IF (eigenvalues(1) > r_ldos_p(ildos)%ldos%eval_range(1)) &
528 0 : r_ldos_p(ildos)%ldos%eval_range(1) = eigenvalues(1)
529 0 : IF (eigenvalues(nmo + nvirt) < r_ldos_p(ildos)%ldos%eval_range(2)) &
530 60 : r_ldos_p(ildos)%ldos%eval_range(2) = eigenvalues(nmo + nvirt)
531 : END DO
532 :
533 60 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T15,A))') &
534 30 : "---- PDOS: start iteration on the KS states --- "
535 :
536 548 : DO imo = 1, nmo + nvirt
537 :
538 488 : IF (output_unit > 0 .AND. MOD(imo, out_each) == 0) WRITE (UNIT=output_unit, FMT='((T20,A,I10))') &
539 0 : " KS state index : ", imo
540 : ! Extract the eigenvector from the distributed full matrix
541 488 : IF (imo > nmo) THEN
542 : CALL cp_fm_get_submatrix(matrix_work, vecbuffer, 1, imo - nmo, &
543 166 : nao, 1, transpose=.TRUE.)
544 : ELSE
545 : CALL cp_fm_get_submatrix(matrix_shalfc, vecbuffer, 1, imo, &
546 322 : nao, 1, transpose=.TRUE.)
547 : END IF
548 :
549 : ! Calculate the pdos for all the kinds
550 488 : irow = 1
551 1700 : DO iatom = 1, natom
552 1212 : firstrow(iatom) = irow
553 1212 : NULLIFY (orb_basis_set)
554 1212 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
555 1212 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
556 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
557 : nset=nset, &
558 : nshell=nshell, &
559 1212 : l=l, maxl=maxl)
560 2912 : IF (separate_components) THEN
561 716 : isgf = 1
562 3256 : DO iset = 1, nset
563 6946 : DO ishell = 1, nshell(iset)
564 3690 : lshell = l(ishell, iset)
565 14544 : DO iso = 1, nso(lshell)
566 8314 : lcomponent = nsoset(lshell - 1) + iso
567 : pdos_array(lcomponent, ikind, imo) = &
568 : pdos_array(lcomponent, ikind, imo) + &
569 8314 : vecbuffer(1, irow)*vecbuffer(1, irow)
570 12004 : irow = irow + 1
571 : END DO ! iso
572 : END DO ! ishell
573 : END DO ! iset
574 : ELSE
575 496 : isgf = 1
576 1392 : DO iset = 1, nset
577 2928 : DO ishell = 1, nshell(iset)
578 1536 : lshell = l(ishell, iset)
579 5536 : DO iso = 1, nso(lshell)
580 : pdos_array(lshell, ikind, imo) = &
581 : pdos_array(lshell, ikind, imo) + &
582 3104 : vecbuffer(1, irow)*vecbuffer(1, irow)
583 4640 : irow = irow + 1
584 : END DO ! iso
585 : END DO ! ishell
586 : END DO ! iset
587 : END IF
588 : END DO ! iatom
589 :
590 : ! Calculate the pdos for all the lists
591 608 : DO ildos = 1, nldos
592 932 : DO il = 1, ldos_p(ildos)%ldos%nlist
593 324 : iatom = ldos_p(ildos)%ldos%list_index(il)
594 :
595 324 : irow = firstrow(iatom)
596 324 : NULLIFY (orb_basis_set)
597 324 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
598 324 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
599 :
600 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
601 : nset=nset, &
602 : nshell=nshell, &
603 324 : l=l, maxl=maxl)
604 324 : ldos_p(ildos)%ldos%maxl = MAX(ldos_p(ildos)%ldos%maxl, maxl)
605 768 : IF (ldos_p(ildos)%ldos%separate_components) THEN
606 108 : isgf = 1
607 324 : DO iset = 1, nset
608 720 : DO ishell = 1, nshell(iset)
609 396 : lshell = l(ishell, iset)
610 1440 : DO iso = 1, nso(lshell)
611 828 : lcomponent = nsoset(lshell - 1) + iso
612 : ldos_p(ildos)%ldos%pdos_array(lcomponent, imo) = &
613 : ldos_p(ildos)%ldos%pdos_array(lcomponent, imo) + &
614 828 : vecbuffer(1, irow)*vecbuffer(1, irow)
615 1224 : irow = irow + 1
616 : END DO ! iso
617 : END DO ! ishell
618 : END DO ! iset
619 : ELSE
620 216 : isgf = 1
621 648 : DO iset = 1, nset
622 1428 : DO ishell = 1, nshell(iset)
623 780 : lshell = l(ishell, iset)
624 2820 : DO iso = 1, nso(lshell)
625 : ldos_p(ildos)%ldos%pdos_array(lshell, imo) = &
626 : ldos_p(ildos)%ldos%pdos_array(lshell, imo) + &
627 1608 : vecbuffer(1, irow)*vecbuffer(1, irow)
628 2388 : irow = irow + 1
629 : END DO ! iso
630 : END DO ! ishell
631 : END DO ! iset
632 : END IF
633 : END DO !il
634 : END DO !ildos
635 :
636 : ! Calculate the DOS projected in a given volume in real space
637 548 : DO ildos = 1, n_r_ldos
638 0 : IF (r_ldos_p(ildos)%ldos%eval_range(1) <= eigenvalues(imo) .AND. &
639 488 : r_ldos_p(ildos)%ldos%eval_range(2) >= eigenvalues(imo)) THEN
640 :
641 0 : IF (imo > nmo) THEN
642 : CALL calculate_wavefunction(mo_virt, imo - nmo, &
643 : wf_r, wf_g, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
644 0 : pw_env)
645 : ELSE
646 : CALL calculate_wavefunction(mo_coeff, imo, &
647 : wf_r, wf_g, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
648 0 : pw_env)
649 : END IF
650 0 : r_ldos_p(ildos)%ldos%pdos_array(imo) = 0.0_dp
651 0 : DO il = 1, r_ldos_p(ildos)%ldos%npoints
652 0 : j = j + 1
653 0 : jx = r_ldos_p(ildos)%ldos%index_grid_local(1, il)
654 0 : jy = r_ldos_p(ildos)%ldos%index_grid_local(2, il)
655 0 : jz = r_ldos_p(ildos)%ldos%index_grid_local(3, il)
656 : r_ldos_p(ildos)%ldos%pdos_array(imo) = r_ldos_p(ildos)%ldos%pdos_array(imo) + &
657 0 : wf_r%array(jx, jy, jz)*wf_r%array(jx, jy, jz)
658 : END DO
659 0 : r_ldos_p(ildos)%ldos%pdos_array(imo) = r_ldos_p(ildos)%ldos%pdos_array(imo)*dvol
660 : END IF
661 : END DO
662 : END DO ! imo
663 :
664 60 : CALL cp_fm_release(matrix_shalfc)
665 60 : DEALLOCATE (vecbuffer)
666 :
667 60 : CALL section_vals_val_get(dft_section, "PRINT%PDOS%APPEND", l_val=append)
668 60 : IF (append .AND. iterstep > 1) THEN
669 6 : my_pos = "APPEND"
670 : ELSE
671 54 : my_pos = "REWIND"
672 : END IF
673 60 : my_act = "WRITE"
674 176 : DO ikind = 1, nkind
675 :
676 116 : NULLIFY (orb_basis_set)
677 116 : CALL get_atomic_kind(atomic_kind_set(ikind), name=kind_name)
678 116 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
679 116 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, maxl=maxl)
680 :
681 : ! basis none has no associated maxl, and no pdos
682 116 : IF (maxl < 0) CYCLE
683 :
684 116 : IF (PRESENT(ispin)) THEN
685 20 : IF (PRESENT(xas_mittle)) THEN
686 20 : my_mittle = TRIM(xas_mittle)//TRIM(spin(ispin))//"_k"//TRIM(ADJUSTL(cp_to_string(ikind)))
687 : ELSE
688 0 : my_mittle = TRIM(spin(ispin))//"_k"//TRIM(ADJUSTL(cp_to_string(ikind)))
689 : END IF
690 20 : my_spin = ispin
691 : ELSE
692 96 : my_mittle = "k"//TRIM(ADJUSTL(cp_to_string(ikind)))
693 96 : my_spin = 1
694 : END IF
695 :
696 : iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%PDOS", &
697 : extension=".pdos", file_position=my_pos, file_action=my_act, &
698 116 : file_form="FORMATTED", middle_name=TRIM(my_mittle))
699 116 : IF (iw > 0) THEN
700 :
701 58 : fmtstr1 = "(I8,2X,2F16.6, (2X,F16.8))"
702 58 : fmtstr2 = "(A42, (10X,A8))"
703 58 : IF (separate_components) THEN
704 16 : WRITE (UNIT=fmtstr1(15:16), FMT="(I2)") nsoset(maxl)
705 16 : WRITE (UNIT=fmtstr2(6:7), FMT="(I2)") nsoset(maxl)
706 : ELSE
707 42 : WRITE (UNIT=fmtstr1(15:16), FMT="(I2)") maxl + 1
708 42 : WRITE (UNIT=fmtstr2(6:7), FMT="(I2)") maxl + 1
709 : END IF
710 :
711 : WRITE (UNIT=iw, FMT="(A,I0,A,F12.6,A)") &
712 58 : "# Projected DOS for atomic kind "//TRIM(kind_name)//" at iteration step i = ", &
713 116 : iterstep, ", E(Fermi) = ", e_fermi, " a.u."
714 58 : IF (separate_components) THEN
715 64 : ALLOCATE (tmp_str(0:0, 0:maxl, -maxl:maxl))
716 496 : tmp_str = ""
717 60 : DO j = 0, maxl
718 184 : DO i = -j, j
719 168 : tmp_str(0, j, i) = sgf_symbol(0, j, i)
720 : END DO
721 : END DO
722 :
723 : WRITE (UNIT=iw, FMT=fmtstr2) &
724 16 : "# MO Eigenvalue [a.u.] Occupation", &
725 200 : ((TRIM(tmp_str(0, il, im)), im=-il, il), il=0, maxl)
726 328 : DO imo = 1, nmo + nvirt
727 312 : WRITE (UNIT=iw, FMT=fmtstr1) imo, eigenvalues(imo), occupation_numbers(imo), &
728 640 : (pdos_array(lshell, ikind, imo), lshell=1, nsoset(maxl))
729 : END DO
730 16 : DEALLOCATE (tmp_str)
731 : ELSE
732 : WRITE (UNIT=iw, FMT=fmtstr2) &
733 42 : "# MO Eigenvalue [a.u.] Occupation", &
734 178 : (TRIM(l_sym(il)), il=0, maxl)
735 210 : DO imo = 1, nmo + nvirt
736 168 : WRITE (UNIT=iw, FMT=fmtstr1) imo, eigenvalues(imo), occupation_numbers(imo), &
737 378 : (pdos_array(lshell, ikind, imo), lshell=0, maxl)
738 : END DO
739 : END IF
740 : END IF
741 : CALL cp_print_key_finished_output(iw, logger, dft_section, &
742 292 : "PRINT%PDOS")
743 :
744 : END DO ! ikind
745 :
746 : ! write the pdos for the lists, each ona different file,
747 : ! the filenames are indexed with the list number
748 80 : DO ildos = 1, nldos
749 : ! basis none has no associated maxl, and no pdos
750 80 : IF (ldos_p(ildos)%ldos%maxl > 0) THEN
751 :
752 20 : IF (PRESENT(ispin)) THEN
753 0 : IF (PRESENT(xas_mittle)) THEN
754 0 : my_mittle = TRIM(xas_mittle)//TRIM(spin(ispin))//"_list"//TRIM(ldos_index(ildos))
755 : ELSE
756 0 : my_mittle = TRIM(spin(ispin))//"_list"//TRIM(ldos_index(ildos))
757 : END IF
758 0 : my_spin = ispin
759 : ELSE
760 20 : my_mittle = "list"//TRIM(ldos_index(ildos))
761 20 : my_spin = 1
762 : END IF
763 :
764 : iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%PDOS", &
765 : extension=".pdos", file_position=my_pos, file_action=my_act, &
766 20 : file_form="FORMATTED", middle_name=TRIM(my_mittle))
767 20 : IF (iw > 0) THEN
768 :
769 10 : fmtstr1 = "(I8,2X,2F16.6, (2X,F16.8))"
770 10 : fmtstr2 = "(A42, (10X,A8))"
771 10 : IF (ldos_p(ildos)%ldos%separate_components) THEN
772 2 : WRITE (UNIT=fmtstr1(15:16), FMT="(I2)") nsoset(ldos_p(ildos)%ldos%maxl)
773 2 : WRITE (UNIT=fmtstr2(6:7), FMT="(I2)") nsoset(ldos_p(ildos)%ldos%maxl)
774 : ELSE
775 8 : WRITE (UNIT=fmtstr1(15:16), FMT="(I2)") ldos_p(ildos)%ldos%maxl + 1
776 8 : WRITE (UNIT=fmtstr2(6:7), FMT="(I2)") ldos_p(ildos)%ldos%maxl + 1
777 : END IF
778 :
779 : WRITE (UNIT=iw, FMT="(A,I0,A,I0,A,I0,A,F12.6,A)") &
780 10 : "# Projected DOS for list ", ildos, " of ", ldos_p(ildos)%ldos%nlist, &
781 10 : " atoms, at iteration step i = ", iterstep, &
782 20 : ", E(Fermi) = ", e_fermi, " a.u."
783 10 : IF (ldos_p(ildos)%ldos%separate_components) THEN
784 8 : ALLOCATE (tmp_str(0:0, 0:ldos_p(ildos)%ldos%maxl, -ldos_p(ildos)%ldos%maxl:ldos_p(ildos)%ldos%maxl))
785 72 : tmp_str = ""
786 8 : DO j = 0, ldos_p(ildos)%ldos%maxl
787 26 : DO i = -j, j
788 24 : tmp_str(0, j, i) = sgf_symbol(0, j, i)
789 : END DO
790 : END DO
791 :
792 : WRITE (UNIT=iw, FMT=fmtstr2) &
793 2 : "# MO Eigenvalue [a.u.] Occupation", &
794 28 : ((TRIM(tmp_str(0, il, im)), im=-il, il), il=0, ldos_p(ildos)%ldos%maxl)
795 20 : DO imo = 1, nmo + nvirt
796 18 : WRITE (UNIT=iw, FMT=fmtstr1) imo, eigenvalues(imo), occupation_numbers(imo), &
797 200 : (ldos_p(ildos)%ldos%pdos_array(lshell, imo), lshell=1, nsoset(ldos_p(ildos)%ldos%maxl))
798 : END DO
799 2 : DEALLOCATE (tmp_str)
800 : ELSE
801 : WRITE (UNIT=iw, FMT=fmtstr2) &
802 8 : "# MO Eigenvalue [a.u.] Occupation", &
803 39 : (TRIM(l_sym(il)), il=0, ldos_p(ildos)%ldos%maxl)
804 50 : DO imo = 1, nmo + nvirt
805 42 : WRITE (UNIT=iw, FMT=fmtstr1) imo, eigenvalues(imo), occupation_numbers(imo), &
806 209 : (ldos_p(ildos)%ldos%pdos_array(lshell, imo), lshell=0, ldos_p(ildos)%ldos%maxl)
807 : END DO
808 : END IF
809 : END IF
810 : CALL cp_print_key_finished_output(iw, logger, dft_section, &
811 20 : "PRINT%PDOS")
812 : END IF ! maxl>0
813 : END DO ! ildos
814 :
815 : ! write the pdos for the lists, each ona different file,
816 : ! the filenames are indexed with the list number
817 60 : DO ildos = 1, n_r_ldos
818 :
819 0 : npoints = r_ldos_p(ildos)%ldos%npoints
820 0 : CALL para_env%sum(npoints)
821 0 : CALL para_env%sum(np_tot)
822 0 : CALL para_env%sum(r_ldos_p(ildos)%ldos%pdos_array)
823 0 : IF (PRESENT(ispin)) THEN
824 0 : IF (PRESENT(xas_mittle)) THEN
825 0 : my_mittle = TRIM(xas_mittle)//TRIM(spin(ispin))//"_r_list"//TRIM(r_ldos_index(ildos))
826 : ELSE
827 0 : my_mittle = TRIM(spin(ispin))//"_r_list"//TRIM(r_ldos_index(ildos))
828 : END IF
829 0 : my_spin = ispin
830 : ELSE
831 0 : my_mittle = "r_list"//TRIM(r_ldos_index(ildos))
832 0 : my_spin = 1
833 : END IF
834 :
835 : iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%PDOS", &
836 : extension=".pdos", file_position=my_pos, file_action=my_act, &
837 0 : file_form="FORMATTED", middle_name=TRIM(my_mittle))
838 0 : IF (iw > 0) THEN
839 0 : fmtstr1 = "(I8,2X,2F16.6, (2X,F16.8))"
840 0 : fmtstr2 = "(A42, (10X,A8))"
841 :
842 : WRITE (UNIT=iw, FMT="(A,I0,A,F12.6,F12.6,A,F12.6,A)") &
843 0 : "# Projected DOS in real space, using ", npoints, &
844 0 : " points of the grid, and eval in the range", r_ldos_p(ildos)%ldos%eval_range(1:2), &
845 0 : " Hartree, E(Fermi) = ", e_fermi, " a.u."
846 : WRITE (UNIT=iw, FMT="(A)") &
847 0 : "# MO Eigenvalue [a.u.] Occupation LDOS"
848 0 : DO imo = 1, nmo + nvirt
849 0 : IF (r_ldos_p(ildos)%ldos%eval_range(1) <= eigenvalues(imo) .AND. &
850 0 : r_ldos_p(ildos)%ldos%eval_range(2) >= eigenvalues(imo)) THEN
851 0 : WRITE (UNIT=iw, FMT="(I8,2X,2F16.6,E20.10,E20.10)") imo, eigenvalues(imo), occupation_numbers(imo), &
852 0 : r_ldos_p(ildos)%ldos%pdos_array(imo), r_ldos_p(ildos)%ldos%pdos_array(imo)*np_tot
853 : END IF
854 : END DO
855 :
856 : END IF
857 : CALL cp_print_key_finished_output(iw, logger, dft_section, &
858 60 : "PRINT%PDOS")
859 : END DO
860 :
861 : ! deallocate local variables
862 60 : DEALLOCATE (pdos_array)
863 60 : DEALLOCATE (firstrow)
864 60 : IF (do_ldos) THEN
865 28 : DO ildos = 1, nldos
866 20 : DEALLOCATE (ldos_p(ildos)%ldos%pdos_array)
867 20 : DEALLOCATE (ldos_p(ildos)%ldos%list_index)
868 28 : DEALLOCATE (ldos_p(ildos)%ldos)
869 : END DO
870 8 : DEALLOCATE (ldos_p)
871 8 : DEALLOCATE (ldos_index)
872 : END IF
873 60 : IF (do_r_ldos) THEN
874 0 : DO ildos = 1, n_r_ldos
875 0 : DEALLOCATE (r_ldos_p(ildos)%ldos%index_grid_local)
876 0 : DEALLOCATE (r_ldos_p(ildos)%ldos%pdos_array)
877 0 : DEALLOCATE (r_ldos_p(ildos)%ldos%list_index)
878 0 : IF (.NOT. read_r(1, ildos)) THEN
879 0 : DEALLOCATE (r_ldos_p(ildos)%ldos%x_range)
880 : END IF
881 0 : IF (.NOT. read_r(2, ildos)) THEN
882 0 : DEALLOCATE (r_ldos_p(ildos)%ldos%y_range)
883 : END IF
884 0 : IF (.NOT. read_r(3, ildos)) THEN
885 0 : DEALLOCATE (r_ldos_p(ildos)%ldos%z_range)
886 : END IF
887 0 : IF (.NOT. read_r(4, ildos)) THEN
888 0 : DEALLOCATE (r_ldos_p(ildos)%ldos%eval_range)
889 : END IF
890 0 : DEALLOCATE (r_ldos_p(ildos)%ldos)
891 : END DO
892 0 : DEALLOCATE (read_r)
893 0 : DEALLOCATE (r_ldos_p)
894 0 : DEALLOCATE (r_ldos_index)
895 0 : CALL auxbas_pw_pool%give_back_pw(wf_r)
896 0 : CALL auxbas_pw_pool%give_back_pw(wf_g)
897 : END IF
898 60 : IF (do_virt) THEN
899 14 : DEALLOCATE (evals_virt)
900 14 : CALL cp_fm_release(mo_virt)
901 14 : CALL cp_fm_release(matrix_work)
902 14 : DEALLOCATE (eigenvalues)
903 14 : DEALLOCATE (occupation_numbers)
904 : END IF
905 :
906 60 : CALL timestop(handle)
907 :
908 600 : END SUBROUTINE calculate_projected_dos
909 :
910 : ! **************************************************************************************************
911 : !> \brief Compute additional virtual states starting from the available MOS
912 : !> \param qs_env ...
913 : !> \param mo_set ...
914 : !> \param evals_virt ...
915 : !> \param mo_virt ...
916 : !> \param nvirt ...
917 : !> \param ispin ...
918 : !> \date 08.03.2008
919 : !> \par History:
920 : !> -
921 : !> \par Variables
922 : !> -
923 : !> -
924 : !> \author MI
925 : !> \version 1.0
926 : ! **************************************************************************************************
927 :
928 14 : SUBROUTINE generate_virtual_mo(qs_env, mo_set, evals_virt, mo_virt, &
929 : nvirt, ispin)
930 :
931 : TYPE(qs_environment_type), POINTER :: qs_env
932 : TYPE(mo_set_type), INTENT(IN) :: mo_set
933 : REAL(KIND=dp), DIMENSION(:), POINTER :: evals_virt
934 : TYPE(cp_fm_type), INTENT(OUT) :: mo_virt
935 : INTEGER, INTENT(IN) :: nvirt, ispin
936 :
937 : INTEGER :: nmo, nrow_global
938 : TYPE(cp_blacs_env_type), POINTER :: context
939 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
940 : TYPE(cp_fm_type), POINTER :: mo_coeff
941 14 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, s_matrix
942 : TYPE(mp_para_env_type), POINTER :: para_env
943 : TYPE(preconditioner_type), POINTER :: local_preconditioner
944 : TYPE(scf_control_type), POINTER :: scf_control
945 :
946 14 : NULLIFY (evals_virt)
947 42 : ALLOCATE (evals_virt(nvirt))
948 :
949 : CALL get_qs_env(qs_env, matrix_ks=ks_matrix, matrix_s=s_matrix, &
950 14 : scf_control=scf_control)
951 14 : CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, nmo=nmo)
952 : CALL cp_fm_get_info(mo_coeff, context=context, para_env=para_env, &
953 14 : nrow_global=nrow_global)
954 :
955 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
956 14 : nrow_global=nrow_global, ncol_global=nvirt)
957 14 : CALL cp_fm_create(mo_virt, fm_struct_tmp, name="virtual")
958 14 : CALL cp_fm_struct_release(fm_struct_tmp)
959 14 : CALL cp_fm_init_random(mo_virt, nvirt)
960 :
961 14 : NULLIFY (local_preconditioner)
962 :
963 : CALL ot_eigensolver(matrix_h=ks_matrix(ispin)%matrix, matrix_s=s_matrix(1)%matrix, &
964 : matrix_c_fm=mo_virt, matrix_orthogonal_space_fm=mo_coeff, &
965 : eps_gradient=scf_control%eps_lumos, &
966 : preconditioner=local_preconditioner, &
967 : iter_max=scf_control%max_iter_lumos, &
968 14 : size_ortho_space=nmo)
969 :
970 : CALL calculate_subspace_eigenvalues(mo_virt, ks_matrix(ispin)%matrix, &
971 14 : evals_virt)
972 :
973 42 : END SUBROUTINE generate_virtual_mo
974 :
975 0 : END MODULE qs_pdos
976 :
|
