Octopus
dos.F90
Go to the documentation of this file.
1!! Copyright (C) 2017 N. Tancogne-Dejean
2!! Copyright (C) 2025 N. Tancogne-Dejean
3!!
4!! This program is free software; you can redistribute it and/or modify
5!! it under the terms of the GNU General Public License as published by
6!! the Free Software Foundation; either version 2, or (at your option)
7!! any later version.
8!!
9!! This program is distributed in the hope that it will be useful,
10!! but WITHOUT ANY WARRANTY; without even the implied warranty of
11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12!! GNU General Public License for more details.
13!!
14!! You should have received a copy of the GNU General Public License
15!! along with this program; if not, write to the Free Software
16!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17!! 02110-1301, USA.
18!!
19
20#include "global.h"
21
23module dos_oct_m
24 use accel_oct_m
26 use box_oct_m
27 use comm_oct_m
28 use debug_oct_m
30 use global_oct_m
32 use io_oct_m
34 use ions_oct_m
35 use, intrinsic :: iso_fortran_env
38 use math_oct_m
39 use mesh_oct_m
41 use mpi_oct_m
46 use parser_oct_m
52 use types_oct_m
53 use unit_oct_m
56
57 implicit none
58
59 private
60 public :: &
61 dos_t, &
62 dos_init, &
66
67 type dos_t
68 private
69 real(real64) :: emin
70 real(real64) :: emax
71 integer :: epoints
72 real(real64) :: de
73
74 real(real64) :: gamma
75
76 logical :: computepdos
77
78 integer :: ldos_nenergies = -1
79 real(real64), allocatable :: ldos_energies(:)
80 contains
81 final :: dos_end
82 end type dos_t
83
84contains
85
87 subroutine dos_init(this, namespace, st, kpoints)
88 type(dos_t), intent(out) :: this
89 type(namespace_t), intent(in) :: namespace
90 type(states_elec_t), intent(in) :: st
91 type(kpoints_t), intent(in) :: kpoints
92
93 real(real64) :: evalmin, evalmax, eextend
94 integer :: npath, ie
95 type(block_t) :: blk
96
97 push_sub(dos_init)
98
99 !The range of the dos is only calculated for physical points,
100 !without the one from a k-point path
101 npath = kpoints%nkpt_in_path()
102 if (st%nik > npath) then
103 evalmin = minval(st%eigenval(1:st%nst, 1:(st%nik-npath)))
104 evalmax = maxval(st%eigenval(1:st%nst, 1:(st%nik-npath)))
105 else !In case we only have a path, e.g., a bandstructure calculation
106 evalmin = minval(st%eigenval(1:st%nst, 1:st%nik))
107 evalmax = maxval(st%eigenval(1:st%nst, 1:st%nik))
108 end if
109 ! we extend the energy mesh by this amount
110 eextend = (evalmax - evalmin) / m_four
111
112 !%Variable DOSEnergyMin
113 !%Type float
114 !%Section Output
115 !%Description
116 !% Lower bound for the energy mesh of the DOS.
117 !% The default is the lowest eigenvalue, minus a quarter of the total range of eigenvalues.
118 !% This is ignored for the joint density of states, and the minimal energy is always set to zero.
119 !%End
120 call parse_variable(namespace, 'DOSEnergyMin', evalmin - eextend, this%emin, units_inp%energy)
121
122 !%Variable DOSEnergyMax
123 !%Type float
124 !%Section Output
125 !%Description
126 !% Upper bound for the energy mesh of the DOS.
127 !% The default is the highest eigenvalue, plus a quarter of the total range of eigenvalues.
128 !%End
129 call parse_variable(namespace, 'DOSEnergyMax', evalmax + eextend, this%emax, units_inp%energy)
130
131 !%Variable DOSEnergyPoints
132 !%Type integer
133 !%Default 500
134 !%Section Output
135 !%Description
136 !% Determines how many energy points <tt>Octopus</tt> should use for
137 !% the DOS energy grid.
138 !%End
139 call parse_variable(namespace, 'DOSEnergyPoints', 500, this%epoints)
140
141 !%Variable DOSGamma
142 !%Type float
143 !%Default 0.008 Ha
144 !%Section Output
145 !%Description
146 !% Determines the width of the Lorentzian which is used for the DOS sum.
147 !%End
148 call parse_variable(namespace, 'DOSGamma', 0.008_real64, this%gamma)
149
150 !%Variable DOSComputePDOS
151 !%Type logical
152 !%Default false
153 !%Section Output
154 !%Description
155 !% Determines if projected dos are computed or not.
156 !% At the moment, the PDOS is computed from the bare pseudo-atomic orbitals, directly taken from
157 !% the pseudopotentials. The orbitals are not orthonormalized, in order to preserve their
158 !% atomic orbitals character. As a consequence, the sum of the different PDOS does not integrate
159 !% to the total DOS.
160 !%
161 !% The radii of the orbitals are controled by the threshold defined by <tt>AOThreshold</tt>,
162 !% and the fact that they are normalized or not by <tt>AONormalize</tt>.
163 !%End
164 call parse_variable(namespace, 'DOSComputePDOS', .false., this%computepdos)
166 ! spacing for energy mesh
167 this%de = (this%emax - this%emin) / (this%epoints - 1)
168
169 !%Variable LDOSEnergies
170 !%Type block
171 !%Section Output
172 !%Description
173 !% Specifies the energies at which the LDOS is computed.
174 !%End
175 if (parse_block(global_namespace, 'LDOSEnergies', blk) == 0) then
176 ! There is one high symmetry k-point per line
177 this%ldos_nenergies = parse_block_cols(blk, 0)
178
179 safe_allocate(this%ldos_energies(1:this%ldos_nenergies))
180 do ie = 1, this%ldos_nenergies
181 call parse_block_float(blk, 0, ie-1, this%ldos_energies(ie))
182 end do
183 call parse_block_end(blk)
184 else
185 this%ldos_nenergies = -1
186 end if
187
188 pop_sub(dos_init)
189 end subroutine dos_init
190
192 subroutine dos_end(this)
193 type(dos_t), intent(inout) :: this
194
195 push_sub(dos_end)
196
197 safe_deallocate_a(this%ldos_energies)
198 this%ldos_nenergies = -1
199
200 pop_sub(dos_end)
201 end subroutine
202
203 ! ---------------------------------------------------------
205 subroutine dos_write_dos(this, dir, st, box, ions, mesh, hm, namespace)
206 type(dos_t), intent(in) :: this
207 character(len=*), intent(in) :: dir
208 type(states_elec_t), target, intent(in) :: st
209 class(box_t), intent(in) :: box
210 type(ions_t), target, intent(in) :: ions
211 class(mesh_t), intent(in) :: mesh
212 type(hamiltonian_elec_t), intent(in) :: hm
213 type(namespace_t), intent(in) :: namespace
214
215 integer :: ie, ik, ist, is, ns, maxdos, ib, ind
216 integer, allocatable :: iunit(:)
217 real(real64) :: energy
218 real(real64), allocatable :: tdos(:)
219 real(real64), allocatable :: dos(:,:,:)
220 character(len=64) :: filename,format_str
221 logical :: normalize
222
223 integer :: ii, ll, mm, nn, work, norb, work2
224 integer :: ia, iorb, idim
225 real(real64) :: threshold
226 real(real64), allocatable :: ddot(:,:,:)
227 complex(real64), allocatable :: zdot(:,:,:)
228 real(real64), allocatable :: weight(:,:,:)
229 type(orbitalset_t) :: os
230 type(wfs_elec_t), pointer :: epsib
231
232 push_sub(dos_write_dos)
233
234 ! shortcuts
235 ns = 1
236 if (st%d%nspin == 2) ns = 2
237
238 if (mpi_world%is_root()) then
239 ! space for state-dependent DOS
240 safe_allocate(dos(1:this%epoints, 1:st%nst, 0:ns-1))
241 safe_allocate(iunit(0:ns-1))
242
243 ! compute band/spin-resolved density of states
244 do ist = 1, st%nst
245
246 do is = 0, ns-1
247 if (ns > 1) then
248 write(filename, '(a,i4.4,a,i1.1,a)') 'dos-', ist, '-', is+1,'.dat'
249 else
250 write(filename, '(a,i4.4,a)') 'dos-', ist, '.dat'
251 end if
252 iunit(is) = io_open(trim(dir)//'/'//trim(filename), namespace, action='write')
253 ! write header
254 write(iunit(is), '(3a)') '# energy [', trim(units_abbrev(units_out%energy)), '], band-resolved DOS'
255 end do
256
257 do ie = 1, this%epoints
258 energy = this%emin + (ie - 1) * this%de
259 dos(ie, ist, :) = m_zero
260 ! sum up Lorentzians
261 do ik = 1, st%nik, ns
262 do is = 0, ns-1
263 dos(ie, ist, is) = dos(ie, ist, is) + st%kweights(ik+is) * m_one/m_pi * &
264 this%gamma / ((energy - st%eigenval(ist, ik+is))**2 + this%gamma**2)
265 end do
266 end do
267 do is = 0, ns-1
268 write(message(1), '(2f12.6)') units_from_atomic(units_out%energy, energy), &
269 units_from_atomic(unit_one / units_out%energy, dos(ie, ist, is))
270 call messages_info(1, iunit(is))
271 end do
272 end do
273
274 do is = 0, ns-1
275 call io_close(iunit(is))
276 end do
277 end do
278
279 safe_allocate(tdos(1))
280
281 ! for spin-polarized calculations also output spin-resolved tDOS
282 if (st%d%nspin == spin_polarized) then
283 do is = 0, ns-1
284 write(filename, '(a,i1.1,a)') 'total-dos-', is+1,'.dat'
285 iunit(is) = io_open(trim(dir)//'/'//trim(filename), namespace, action='write')
286 ! write header
287 write(iunit(is), '(3a)') '# energy [', trim(units_abbrev(units_out%energy)), '], total DOS (spin-resolved)'
288
289 do ie = 1, this%epoints
290 energy = this%emin + (ie - 1) * this%de
291 tdos(1) = m_zero
292 do ist = 1, st%nst
293 tdos(1) = tdos(1) + dos(ie, ist, is)
294 end do
295 write(message(1), '(2f12.6)') units_from_atomic(units_out%energy, energy), &
296 units_from_atomic(unit_one / units_out%energy, tdos(1))
297 call messages_info(1, iunit(is))
298 end do
299
300 call io_close(iunit(is))
301 end do
302 end if
303
304
305 iunit(0) = io_open(trim(dir)//'/'//'total-dos.dat', namespace, action='write')
306 write(iunit(0), '(3a)') '# energy [', trim(units_abbrev(units_out%energy)), '], total DOS'
307
308 ! compute total density of states
309 do ie = 1, this%epoints
310 energy = this%emin + (ie - 1) * this%de
311 tdos(1) = m_zero
312 do ist = 1, st%nst
313 do is = 0, ns-1
314 tdos(1) = tdos(1) + dos(ie, ist, is)
315 end do
316 end do
317 write(message(1), '(2f12.6)') units_from_atomic(units_out%energy, energy), &
318 units_from_atomic(unit_one / units_out%energy, tdos(1))
319 call messages_info(1, iunit(0))
320 end do
321
322 call io_close(iunit(0))
323
324 safe_deallocate_a(tdos)
325
326
327 ! write Fermi file
328 iunit(0) = io_open(trim(dir)//'/'//'total-dos-efermi.dat', namespace, action='write')
329 write(message(1), '(3a)') '# Fermi energy [', trim(units_abbrev(units_out%energy)), &
330 '] in a format compatible with total-dos.dat'
331
332 ! this is the maximum that tdos can reach
333 maxdos = st%smear%el_per_state * st%nst
334
335 write(message(2), '(2f12.6)') units_from_atomic(units_out%energy, st%smear%e_fermi), m_zero
336 write(message(3), '(f12.6,i6)') units_from_atomic(units_out%energy, st%smear%e_fermi), maxdos
337
338 call messages_info(3, iunit(0))
339 call io_close(iunit(0))
340
341 end if
342
343 if (this%computepdos) then
344
345 !These variables are defined in basis_set/orbitalbasis.F90
346 call parse_variable(namespace, 'AOThreshold', 0.01_real64, threshold)
347 call parse_variable(namespace, 'AONormalize', .true., normalize)
348
349 do ia = 1, ions%natoms
350 !We first count how many orbital set we have
351 work = orbitalset_utils_count(ions%atom(ia)%species)
352
353 !We loop over the orbital sets of the atom ia
354 do norb = 1, work
355 call orbitalset_init(os)
356 os%spec => ions%atom(ia)%species
357
358 !We count the orbitals
359 work2 = 0
360 do iorb = 1, os%spec%get_niwfs()
361 call os%spec%get_iwf_ilm(iorb, 1, ii, ll, mm)
362 call os%spec%get_iwf_n(iorb, 1, nn)
363 if (ii == norb) then
364 os%ll = ll
365 os%nn = nn
366 os%ii = ii
367 os%radius = atomic_orbital_get_radius(os%spec, mesh, iorb, 1, &
368 option__aotruncation__ao_full, threshold)
369 work2 = work2 + 1
370 end if
371 end do
372 os%norbs = work2
373 os%ndim = 1
374 os%use_submesh = .false.
375 os%allocated_on_mesh = .true.
376 os%spec => ions%atom(ia)%species
377
378 do work = 1, os%norbs
379 ! We obtain the orbital
380 if (states_are_real(st)) then
381 call dget_atomic_orbital(namespace, ions%space, ions%latt, ions%pos(:,ia), &
382 ions%atom(ia)%species, mesh, os%sphere, os%ii, os%ll, os%jj, &
383 os, work, os%radius, os%ndim, use_mesh=.not.os%use_submesh, &
384 normalize = normalize)
385 else
386 call zget_atomic_orbital(namespace, ions%space, ions%latt, ions%pos(:,ia), &
387 ions%atom(ia)%species, mesh, os%sphere, os%ii, os%ll, os%jj, &
388 os, work, os%radius, os%ndim, &
389 use_mesh = .not. hm%phase%is_allocated() .and. .not. os%use_submesh, &
390 normalize = normalize)
391 end if
392 end do
393
394 if (hm%phase%is_allocated()) then
395 ! In case of complex wavefunction, we allocate the array for the phase correction
396 safe_allocate(os%phase(1:os%sphere%np, st%d%kpt%start:st%d%kpt%end))
397 os%phase(:,:) = m_zero
398
399 if (.not. os%use_submesh) then
400 safe_allocate(os%eorb_mesh(1:mesh%np, 1:os%norbs, 1:os%ndim, st%d%kpt%start:st%d%kpt%end))
401 os%eorb_mesh(:,:,:,:) = m_zero
402 else
403 safe_allocate(os%eorb_submesh(1:os%sphere%np, 1:os%ndim, 1:os%norbs, st%d%kpt%start:st%d%kpt%end))
404 os%eorb_submesh(:,:,:,:) = m_zero
405 end if
406
407 if (accel_is_enabled() .and. st%d%dim == 1) then
408 os%ldorbs = max(pad_pow2(os%sphere%np), 1)
409 if(.not. os%use_submesh) os%ldorbs = max(pad_pow2(os%sphere%mesh%np), 1)
410
411 safe_allocate(os%buff_eorb(st%d%kpt%start:st%d%kpt%end))
412 do ik= st%d%kpt%start, st%d%kpt%end
413 call accel_create_buffer(os%buff_eorb(ik), accel_mem_read_only, type_cmplx, os%ldorbs*os%norbs)
414 end do
415 end if
416
417 call orbitalset_update_phase(os, box%dim, st%d%kpt, hm%kpoints, (st%d%ispin==spin_polarized), &
418 vec_pot = hm%hm_base%uniform_vector_potential, &
419 vec_pot_var = hm%hm_base%vector_potential)
420 end if
421
422 if (mpi_world%is_root()) then
423 if (os%nn /= 0) then
424 write(filename,'(a, i3.3, a1, a, i1.1, a1,a)') 'pdos-at', ia, '-', trim(os%spec%get_label()), &
425 os%nn, l_notation(os%ll), '.dat'
426 else
427 write(filename,'(a, i3.3, a1, a, a1,a)') 'pdos-at', ia, '-', trim(os%spec%get_label()), &
428 l_notation(os%ll), '.dat'
429 end if
430
431 iunit(0) = io_open(trim(dir)//'/'//trim(filename), namespace, action='write')
432 ! write header
433 write(iunit(0), '(3a)') '# energy [', trim(units_abbrev(units_out%energy)), &
434 '], projected DOS (total and orbital resolved)'
435 end if
436
437 if (states_are_real(st)) then
438 safe_allocate(ddot(1:st%d%dim, 1:os%norbs, 1:st%block_size))
439 else
440 safe_allocate(zdot(1:st%d%dim, 1:os%norbs, 1:st%block_size))
441 end if
442
443 safe_allocate(weight(1:os%norbs,1:st%nik,1:st%nst))
444 weight(1:os%norbs,1:st%nik,1:st%nst) = m_zero
445
446 do ik = st%d%kpt%start, st%d%kpt%end
447 do ib = st%group%block_start, st%group%block_end
448
449 if (hm%phase%is_allocated()) then
450 safe_allocate(epsib)
451 call st%group%psib(ib, ik)%copy_to(epsib)
452 call hm%phase%apply_to(mesh, mesh%np, .false., epsib, src = st%group%psib(ib, ik))
453 else
454 epsib => st%group%psib(ib, ik)
455 end if
456
457 if (states_are_real(st)) then
458 call dorbitalset_get_coeff_batch(os, st%d%dim, epsib, ddot(:, :, :))
459 do ist = 1, st%group%psib(ib, ik)%nst
460 ind = st%group%psib(ib, ik)%ist(ist)
461 do iorb = 1, os%norbs
462 do idim = 1, st%d%dim
463 weight(iorb, ik, ind) = weight(iorb, ik, ind) + st%kweights(ik) * abs(ddot(idim, iorb, ist))**2
464 end do
465 end do
466 end do
467 else
468 call zorbitalset_get_coeff_batch(os, st%d%dim, epsib, zdot(:, :, :))
469 do ist = 1, st%group%psib(ib, ik)%nst
470 ind = st%group%psib(ib, ik)%ist(ist)
471 do iorb = 1, os%norbs
472 do idim = 1, st%d%dim
473 weight(iorb, ik, ind) = weight(iorb, ik, ind) + st%kweights(ik) * abs(zdot(idim, iorb, ist))**2
474 end do
475 end do
476 end do
477 end if
478
479 if (hm%phase%is_allocated()) then
480 call epsib%end(copy=.false.)
481 safe_deallocate_p(epsib)
482 end if
483
484 end do
485 end do
486
487 if (st%parallel_in_states .or. st%d%kpt%parallel) then
488 call comm_allreduce(st%st_kpt_mpi_grp, weight)
489 end if
490
491 safe_deallocate_a(ddot)
492 safe_deallocate_a(zdot)
493
494 if (mpi_world%is_root()) then
495 write(format_str,'(a,i2,a)') '(', os%norbs+2, 'f12.6)'
496 safe_allocate(tdos(1:os%norbs))
497 do ie = 1, this%epoints
498 energy = this%emin + (ie - 1) * this%de
499 do iorb = 1, os%norbs
500 tdos(iorb) = m_zero
501 do ist = 1, st%nst
502 do ik = 1, st%nik
503 tdos(iorb) = tdos(iorb) + weight(iorb,ik,ist) * m_one/m_pi * &
504 this%gamma / ((energy - st%eigenval(ist, ik))**2 + this%gamma**2)
505 end do
506 end do
507 end do
508
509 write(iunit(0), trim(format_str)) units_from_atomic(units_out%energy, energy), &
510 units_from_atomic(unit_one / units_out%energy, sum(tdos)), &
511 (units_from_atomic(unit_one / units_out%energy, tdos(iorb)), iorb=1,os%norbs)
512 end do
513 safe_deallocate_a(tdos)
514 call io_close(iunit(0))
515 end if
516
517 call orbitalset_end(os)
518 safe_deallocate_a(weight)
519
520 end do
521
522 end do
523
524 end if
525
526 safe_deallocate_a(iunit)
527 safe_deallocate_a(dos)
528
529 pop_sub(dos_write_dos)
530 end subroutine dos_write_dos
531
532 ! ---------------------------------------------------------
534 subroutine dos_write_jdos(this, dir, st, ions, hm, namespace)
535 type(dos_t), intent(in) :: this
536 character(len=*), intent(in) :: dir
537 type(states_elec_t), intent(in) :: st
538 type(ions_t), target, intent(in) :: ions
539 type(hamiltonian_elec_t), intent(in) :: hm
540 type(namespace_t), intent(in) :: namespace
541
542 integer :: ie, ik, val, cond, is, ns
543 integer, allocatable :: iunit(:)
544 real(real64) :: energy
545 real(real64) :: tjdos(1)
546 real(real64), allocatable :: jdos(:,:)
547 character(len=64) :: filename
548
549 push_sub(dos_write_jdos)
550
551 ! shortcuts
552 ns = 1
553 if (st%d%nspin == 2) ns = 2
554
555 if (mpi_world%is_root()) then
556 ! space for state-dependent DOS
557 safe_allocate(jdos(1:this%epoints, 0:ns-1))
558 safe_allocate(iunit(0:ns-1))
559 jdos = m_zero
560
561 ! compute band/spin-resolved density of states
562 do val = 1, st%nst
563 do cond = val, st%nst
564 do ik = 1, st%nik, ns
565 do is = 0, ns-1
566 if(st%occ(val, ik+is) < m_epsilon) cycle
567 if(st%occ(cond, ik+is) > m_epsilon) cycle
568 do ie = 1, this%epoints
569 energy = (ie - 1) * this%de
570 ! sum up Lorentzians
571 jdos(ie, is) = jdos(ie, is) + st%kweights(ik+is) * m_one/m_pi * &
572 this%gamma / ((energy - (st%eigenval(cond, ik+is)-st%eigenval(val, ik+is)))**2 + this%gamma**2)
573 end do
574 end do
575 end do
576 end do
577 end do
578
579 ! for spin-polarized calculations also output spin-resolved tDOS
580 if (st%d%nspin > 1) then
581 do is = 0, ns-1
582 write(filename, '(a,i1.1,a)') 'total-jdos-', is+1,'.dat'
583 iunit(is) = io_open(trim(dir)//'/'//trim(filename), namespace, action='write')
584 ! write header
585 write(iunit(is), '(3a)') '# energy [', trim(units_abbrev(units_out%energy)), '], total JDOS (spin-resolved)'
586
587 do ie = 1, this%epoints
588 energy = (ie - 1) * this%de
589 write(message(1), '(2f12.6)') units_from_atomic(units_out%energy, energy), &
590 units_from_atomic(unit_one / units_out%energy, jdos(ie, is))
591 call messages_info(1, iunit(is))
592 end do
593
594 call io_close(iunit(is))
595 end do
596 end if
597
598
599 iunit(0) = io_open(trim(dir)//'/'//'total-jdos.dat', namespace, action='write')
600 write(iunit(0), '(3a)') '# energy [', trim(units_abbrev(units_out%energy)), '], total JDOS'
601
602 ! compute total joint density of states
603 do ie = 1, this%epoints
604 energy = (ie - 1) * this%de
605 tjdos(1) = m_zero
606 do is = 0, ns-1
607 tjdos(1) = tjdos(1) + jdos(ie, is)
608 end do
609 write(message(1), '(2f12.6)') units_from_atomic(units_out%energy, energy), &
610 units_from_atomic(unit_one / units_out%energy, tjdos(1))
611 call messages_info(1, iunit(0))
612 end do
613
614 call io_close(iunit(0))
615 end if
616
617 safe_deallocate_a(iunit)
618 safe_deallocate_a(jdos)
619
620 pop_sub(dos_write_jdos)
621 end subroutine dos_write_jdos
622
623
624 ! ---------------------------------------------------------
626 subroutine dos_write_ldos(this, dir, st, ions, mesh, how, namespace)
627 type(dos_t), intent(in) :: this
628 character(len=*), intent(in) :: dir
629 type(states_elec_t), intent(in) :: st
630 type(ions_t), target, intent(in) :: ions
631 class(mesh_t), intent(in) :: mesh
632 integer(int64), intent(in) :: how
633 type(namespace_t), intent(in) :: namespace
634
635 integer :: ie, ik, ist, is, ns, ip, ierr
636 character(len=MAX_PATH_LEN) :: fname, name
637 real(real64) :: weight
638 real(real64), allocatable :: ldos(:,:,:), dpsi(:,:), abs_psi2(:)
639 complex(real64), allocatable :: zpsi(:,:)
640 type(unit_t) :: fn_unit
641
642 push_sub(dos_write_ldos)
643
644 if (this%ldos_nenergies < 1) then
645 message(1) = "LDOSEnergies must be defined for Output=ldos"
646 call messages_fatal(1, namespace=namespace)
647 end if
648
649 ! shortcuts
650 ns = 1
651 if (st%d%nspin == 2) ns = 2
652
653 fn_unit = units_out%length**(-ions%space%dim) / units_out%energy
654
655 ! space for state-dependent DOS
656 safe_allocate(ldos(1:mesh%np, 1:this%ldos_nenergies, 1:ns))
657 ldos = m_zero
658
659 safe_allocate(abs_psi2(1:mesh%np))
660 if (states_are_real(st)) then
661 safe_allocate(dpsi(1:mesh%np, 1:st%d%dim))
662 else
663 safe_allocate(zpsi(1:mesh%np, 1:st%d%dim))
664 end if
665
666 ! compute band/spin-resolved density of states
667 do ik = st%d%kpt%start, st%d%kpt%end
668 is = st%d%get_spin_index(ik)
669 do ist = st%st_start, st%st_end
670
671 ! TODO: This will be replaced by a poin-wise batch multiplication
672 if (states_are_real(st)) then
673 call states_elec_get_state(st, mesh, ist, ik, dpsi)
674 do ip = 1, mesh%np
675 abs_psi2(ip) = dpsi(ip, 1)**2
676 end do
677 if (st%d%dim > 1) then
678 abs_psi2(ip) = abs_psi2(ip) + dpsi(ip, 2)**2
679 end if
680 else
681 call states_elec_get_state(st, mesh, ist, ik, zpsi)
682 do ip = 1, mesh%np
683 abs_psi2(ip) = real(conjg(zpsi(ip, 1)) * zpsi(ip, 1), real64)
684 end do
685 if (st%d%dim > 1) then
686 abs_psi2(ip) = abs_psi2(ip) + real(conjg(zpsi(ip, 2)) * zpsi(ip, 2), real64)
687 end if
688 end if
689
690 ! sum up Lorentzians
691 do ie = 1, this%ldos_nenergies
692 weight = st%kweights(ik) * m_one/m_pi * &
693 this%gamma / ((this%ldos_energies(ie) - st%eigenval(ist, ik))**2 + this%gamma**2)
694
695 call lalg_axpy(mesh%np, weight, abs_psi2, ldos(:, ie, is))
696 end do
697 end do
698 end do !ist
699
700 safe_deallocate_a(dpsi)
701 safe_deallocate_a(zpsi)
702 safe_deallocate_a(abs_psi2)
703
704 if (st%parallel_in_states .or. st%d%kpt%parallel) then
705 call comm_allreduce(st%st_kpt_mpi_grp, ldos)
706 end if
707
708 do is = 1, ns
709 do ie = 1, this%ldos_nenergies
710 write(name, '(a,i3.3)') 'ldos_en-', ie
711 fname = get_filename_with_spin(name, st%d%nspin, is)
712
713 call dio_function_output(how, dir, fname, namespace, ions%space, mesh, &
714 ldos(:, ie, is), fn_unit, ierr, pos=ions%pos, atoms=ions%atom, grp = st%dom_st_kpt_mpi_grp)
715 end do
716 end do
717
718 safe_deallocate_a(ldos)
719
720 pop_sub(dos_write_ldos)
721 end subroutine dos_write_ldos
722
723
724end module dos_oct_m
725
726!! Local Variables:
727!! mode: f90
728!! coding: utf-8
729!! End:
constant times a vector plus a vector
Definition: lalg_basic.F90:173
pure logical function, public accel_is_enabled()
Definition: accel.F90:419
integer, parameter, public accel_mem_read_only
Definition: accel.F90:196
subroutine, public zget_atomic_orbital(namespace, space, latt, pos, species, mesh, sm, ii, ll, jj, os, orbind, radius, d_dim, use_mesh, normalize, index_shift)
This routine returns the atomic orbital basis – provided by the pseudopotential structure in geo.
character(len=1), dimension(0:3), parameter, public l_notation
real(real64) function, public atomic_orbital_get_radius(species, mesh, iorb, ispin, truncation, threshold)
subroutine, public dget_atomic_orbital(namespace, space, latt, pos, species, mesh, sm, ii, ll, jj, os, orbind, radius, d_dim, use_mesh, normalize, index_shift)
This routine returns the atomic orbital basis – provided by the pseudopotential structure in geo.
Module that handles computing and output of various density of states.
Definition: dos.F90:118
subroutine dos_end(this)
Finalizer for the dos_t object.
Definition: dos.F90:288
subroutine, public dos_write_dos(this, dir, st, box, ions, mesh, hm, namespace)
Computes and output the DOS and the projected DOS (PDOS)
Definition: dos.F90:301
subroutine, public dos_write_jdos(this, dir, st, ions, hm, namespace)
Computes and output the joint DOS (LDOS)
Definition: dos.F90:630
subroutine, public dos_init(this, namespace, st, kpoints)
Initializes the dot_t object.
Definition: dos.F90:183
subroutine, public dos_write_ldos(this, dir, st, ions, mesh, how, namespace)
Computes and output the local DOS (LDOS)
Definition: dos.F90:722
integer, parameter, public spin_polarized
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_four
Definition: global.F90:195
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:189
real(real64), parameter, public m_epsilon
Definition: global.F90:207
real(real64), parameter, public m_one
Definition: global.F90:192
subroutine, public dio_function_output(how, dir, fname, namespace, space, mesh, ff, unit, ierr, pos, atoms, grp, root)
Top-level IO routine for functions defined on the mesh.
Definition: io.F90:116
subroutine, public io_close(iunit, grp)
Definition: io.F90:466
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:401
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:117
integer pure function, public pad_pow2(size)
create array size, which is padded to powers of 2
Definition: math.F90:892
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:594
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:272
type(namespace_t), public global_namespace
Definition: namespace.F90:134
subroutine, public orbitalset_init(this)
Definition: orbitalset.F90:208
subroutine, public orbitalset_end(this)
Definition: orbitalset.F90:234
subroutine, public orbitalset_update_phase(os, dim, kpt, kpoints, spin_polarized, vec_pot, vec_pot_var, kpt_max)
Build the phase correction to the global phase in case the orbital crosses the border of the simulato...
Definition: orbitalset.F90:284
subroutine, public dorbitalset_get_coeff_batch(os, ndim, psib, dot, reduce)
Definition: orbitalset.F90:588
subroutine, public zorbitalset_get_coeff_batch(os, ndim, psib, dot, reduce)
integer function, public orbitalset_utils_count(species, iselect)
Count the number of orbital sets we have for a given atom.
this module contains the low-level part of the output system
Definition: output_low.F90:117
character(len=max_path_len) function, public get_filename_with_spin(output, nspin, spin_index)
Returns the filame as output, or output-spX is spin polarized.
Definition: output_low.F90:234
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:615
pure logical function, public states_are_real(st)
type(type_t), public type_cmplx
Definition: types.F90:136
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:134
character(len=20) pure function, public units_abbrev(this)
Definition: unit.F90:225
This module defines the unit system, used for input and output.
type(unit_system_t), public units_out
type(unit_system_t), public units_inp
the units systems for reading and writing
type(unit_t), public unit_one
some special units required for particular quantities
class to tell whether a point is inside or outside
Definition: box.F90:143
Describes mesh distribution to nodes.
Definition: mesh.F90:187
The states_elec_t class contains all electronic wave functions.
batches of electronic states
Definition: wfs_elec.F90:141
int true(void)