34 use,
intrinsic :: iso_fortran_env
76 integer,
intent(in) :: nst
77 type(states_elec_t),
intent(in) :: st
78 class(space_t),
intent(in) :: space
79 type(kpoints_t),
intent(in) :: kpoints
80 real(real64),
optional,
intent(in) :: error(:,:)
81 integer,
optional,
intent(in) :: st_start
82 logical,
optional,
intent(in) :: compact
83 integer,
optional,
intent(in) :: iunit
84 type(namespace_t),
optional,
intent(in) :: namespace
86 integer :: ik, ikk, ist, ns, is, idir, st_start_, iflat, iqn, homo_index, not_printed
87 logical :: print_eigenval
88 real(real64) :: kpoint(space%dim), max_error
89 character(len=120) :: tmp_str(max(space%dim, 3)), cspin
91 real(real64),
allocatable :: flat_eigenval(:)
92 integer,
allocatable :: flat_indices(:, :)
93 integer,
parameter :: print_range = 8
98 if (
present(st_start)) st_start_ = st_start
110 if (st%d%nspin == 2) ns = 2
115 if (st%d%ispin ==
spinors)
then
116 write(
message(1),
'(a4,1x,a5,1x,a12,1x,a12,2x,a4,4x,a4,4x,a4)') &
117 '#st',
' Spin',
' Eigenvalue',
'Occupation ',
'<Sx>',
'<Sy>',
'<Sz>'
119 write(
message(1),
'(a4,1x,a5,1x,a12,4x,a12)') &
120 '#st',
' Spin',
' Eigenvalue',
'Occupation'
122 if (
present(error))
then
127 do ik = 1, st%nik, ns
128 if (space%is_periodic())
then
129 ikk = st%d%get_kpoint_index(ik)
130 kpoint(1:space%dim) = kpoints%get_point(ikk, absolute_coordinates = .false.)
131 write(
message(1),
'(a,i8,a)')
'#k =', ikk,
', k = ('
132 do idir = 1, space%dim
133 write(tmp_str(1),
'(f10.6)') kpoint(idir)
141 do ist = st_start_, nst
143 if (is == 0) cspin =
'up'
144 if (is == 1) cspin =
'dn'
147 write(tmp_str(1),
'(i4,3x,a2)') ist, trim(cspin)
148 if (st%d%ispin ==
spinors)
then
149 write(tmp_str(2),
'(1x,f12.6,5x,f5.2,3x,3f8.4)') &
151 if (
present(error))
write(tmp_str(3),
'(a3,es8.1,a1)')
' (', error(ist, ik),
')'
153 write(tmp_str(2),
'(1x,f12.6,3x,f12.6)') &
155 if (
present(error))
write(tmp_str(3),
'(a7,es8.1,a1)')
' (', error(ist, ik+is),
')'
157 if (
present(error))
then
158 message(1) = trim(tmp_str(1))//trim(tmp_str(2))//trim(tmp_str(3))
160 message(1) = trim(tmp_str(1))//trim(tmp_str(2))
171 safe_allocate(flat_eigenval(1:st%nik*nst))
172 safe_allocate(flat_indices(1:2, 1:st%nik*nst))
178 flat_eigenval(iflat) = st%eigenval(ist, iqn)
179 flat_indices(1:2, iflat) = (/iqn, ist/)
185 call sort(flat_eigenval, flat_indices(:, :))
187 homo_index = st%nik*nst
188 do iflat = 1, st%nik*nst
189 iqn = flat_indices(1, iflat)
190 ist = flat_indices(2, iflat)
191 if (abs(st%occ(ist, iqn)) < 0.1_real64)
then
192 homo_index = iflat - 1
197 tmp_str(1) =
'# State'
199 if (space%is_periodic()) tmp_str(1) = trim(tmp_str(1))//
' KPoint'
201 if (st%d%ispin ==
spin_polarized) tmp_str(1) = trim(tmp_str(1))//
' Spin'
205 tmp_str(1) = trim(tmp_str(1))//
' Occupation'
207 if (st%d%ispin ==
spinors)
then
208 tmp_str(1) = trim(tmp_str(1))//
' <Sx> <Sy> <Sz>'
211 if (
present(error)) tmp_str(1) = trim(tmp_str(1))//
' Error'
217 max_error = 0.0_real64
218 do iflat = 1, st%nik*nst
219 iqn = flat_indices(1, iflat)
220 ist = flat_indices(2, iflat)
221 ik = st%d%get_kpoint_index(iqn)
222 is = st%d%get_spin_index(iqn)
224 print_eigenval = iflat <= print_range
225 print_eigenval = print_eigenval .or. st%nik*nst - iflat < print_range
226 print_eigenval = print_eigenval .or. abs(iflat - homo_index) <= print_range
228 if (print_eigenval)
then
230 if (not_printed > 0)
then
236 if (
present(error))
then
246 max_error = 0.0_real64
250 write(tmp_str(1),
'(i7)') ist
252 if (space%is_periodic())
then
253 write(tmp_str(1),
'(2a,i7)') trim(tmp_str(1)),
' ', ik
257 if (is == 1) cspin =
' up'
258 if (is == 2) cspin =
' dn'
259 write(tmp_str(1),
'(2a,a5)') trim(tmp_str(1)),
' ', cspin
269 write(tmp_str(1),
'(2a,f11.6)') trim(tmp_str(1)),
' ', st%occ(ist, iqn)
271 if (st%d%ispin ==
spinors)
then
272 write(tmp_str(1),
'(2a,3f9.4)') trim(tmp_str(1)),
' ', st%spin(1:3, ist, iqn)
275 if (
present(error))
then
276 write(tmp_str(1),
'(2a,es8.1,a)') trim(tmp_str(1)),
' (', error(ist, iqn),
')'
284 not_printed = not_printed + 1
286 if (
present(error))
then
287 max_error = max(max_error, error(ist, iqn))
296 safe_deallocate_a(flat_indices)
297 safe_deallocate_a(flat_eigenval)
309 integer,
parameter :: ndiv = 70, height = 10
310 integer :: histogram(1:ndiv), iev, ien, iline, ife
311 real(real64) :: dhistogram(1:ndiv)
312 character(len=ndiv) :: line
313 real(real64) :: emin, emax, de, maxhist
317 emin = flat_eigenval(1)
318 emax = flat_eigenval(st%nik*nst)
319 de = (emax - emin)/(ndiv -
m_one)
326 ife = nint((st%smear%e_fermi - emin)/de) + 1
330 do iev = 1, st%nik*nst
331 ien = nint((flat_eigenval(iev) - emin)/de) + 1
334 dhistogram(ien) = dhistogram(ien) + st%kweights(flat_indices(1, iev))*kpoints%full%npoints
338 maxhist = maxval(dhistogram)
340 histogram(ien) = nint(dhistogram(ien)*height/maxhist)
349 do iline = height, 1, -1
351 if (histogram(ien) >= iline)
then
375 integer,
intent(in) :: iunit
377 class(
space_t),
intent(in) :: space
381 real(real64) :: homo, lumo, egdir, egindir
382 integer :: homok, lumok, egdirk
393 if (ceiling(st%qtot/st%smear%el_per_state) == st%nst)
then
407 if (abs(st%kweights(ik)) <
m_epsilon) cycle
409 if (st%occ(ist, ik) >
m_epsilon .and. st%eigenval(ist, ik) > homo)
then
410 homo = st%eigenval(ist, ik)
414 if (st%occ(ist+1, ik) <=
m_epsilon .and. st%eigenval(ist+1, ik) < lumo)
then
415 lumo = st%eigenval(ist+1, ik)
420 .and. (st%eigenval(ist+1, ik) - st%eigenval(ist, ik)) < egdir)
then
421 egdir = (st%eigenval(ist+1, ik) - st%eigenval(ist, ik))
429 if (lumok == -1 .or. egdir <=
m_epsilon)
then
430 write(
message(1),
'(a)')
'The system seems to have no gap.'
433 write(
message(1),
'(a,i5,a,f7.4,a)')
'Direct gap at ik=', egdirk,
' of ', &
435 write(
message(2),
'(a,i5,a,i5,a,f7.4,a)')
'Indirect gap between ik=', homok,
' and ik=', lumok, &
451 character(len=*),
intent(in) :: dir
453 class(
space_t),
intent(in) :: space
454 class(
mesh_t),
intent(in) :: mesh
458 integer :: ncols, icoord, ist, ik, tpa_initialst, tpa_initialk
461 real(real64),
allocatable :: ff(:)
462 complex(real64),
allocatable :: cff(:)
463 real(real64),
allocatable :: osc(:)
464 real(real64) :: transition_energy, osc_strength, dsf
466 logical :: use_qvector
467 real(real64),
allocatable :: qvector(:), psi_initial(:, :), psi_ist(:, :)
471 use_qvector = .false.
485 if (tpa_initialst == -1)
then
488 call messages_write(
'No orbital with half-occupancy found. TPA output is not written.')
512 if (
parse_block(namespace,
'MomentumTransfer', blk) == 0)
then
517 if (ncols /= space%dim)
then
520 call messages_write(
'Inconsistent size of momentum-transfer vector. It will not be used in the TPA calculation.')
527 safe_allocate(qvector(1:space%dim))
529 do icoord = 1, space%dim
540 safe_allocate(ff(1:mesh%np))
541 if (use_qvector)
then
542 safe_allocate(cff(1:mesh%np))
544 safe_allocate(osc(1:space%dim))
550 iunit =
io_open(trim(dir)//
'/'//trim(
'tpa_xas'), namespace, action=
'write')
553 if (use_qvector)
then
554 write (
message(1),
'(a1,a30,3(es14.5,1x),a1)')
'#',
' momentum-transfer vector : (', &
556 select case (space%dim)
558 write(
message(2),
'(a1,4(a15,1x))')
'#',
'E' ,
'<x>',
'<f>',
'S(q,omega)'
560 write(
message(2),
'(a1,5(a15,1x))')
'#',
'E' ,
'<x>',
'<y>',
'<f>',
'S(q,omega)'
562 write(
message(2),
'(a1,6(a15,1x))')
'#',
'E' ,
'<x>',
'<y>',
'<z>',
'<f>',
'S(q,omega)'
566 select case (space%dim)
568 write(
message(1),
'(a1,3(a15,1x))')
'#',
'E' ,
'<x>',
'<f>'
570 write(
message(1),
'(a1,4(a15,1x))')
'#',
'E' ,
'<x>',
'<y>',
'<f>'
572 write(
message(1),
'(a1,5(a15,1x))')
'#',
'E' ,
'<x>',
'<y>',
'<z>',
'<f>'
579 safe_allocate(psi_initial(1:mesh%np, 1:st%d%dim))
580 safe_allocate(psi_ist(1:mesh%np, 1:st%d%dim))
590 if (abs(st%occ(ist,tpa_initialk)) <
m_min_occ)
then
593 transition_energy = st%eigenval(ist, tpa_initialk) - st%eigenval(tpa_initialst, tpa_initialk)
596 do icoord = 1, space%dim
598 ff(1:mesh%np) = psi_initial(1:mesh%np, 1)*mesh%x(1:mesh%np, icoord)* &
599 psi_ist(1:mesh%np, 1)
601 osc_strength = osc_strength +
m_two/real(space%dim, real64) *transition_energy*abs(osc(icoord))**2
606 if (use_qvector)
then
608 cff(1:mesh%np) = psi_initial(1:mesh%np, 1)*psi_ist(1:mesh%np, 1)
610 do icoord = 1, space%dim
611 cff(1:mesh%np) = cff(1:mesh%np)*
exp(
m_zi*mesh%x(1:mesh%np, icoord)*qvector(icoord))
620 if (use_qvector)
then
640 safe_deallocate_a(psi_initial)
641 safe_deallocate_a(psi_ist)
643 safe_deallocate_a(ff)
644 if (use_qvector)
then
645 safe_deallocate_a(cff)
647 safe_deallocate_a(osc)
648 if (use_qvector)
then
649 safe_deallocate_a(qvector)
661 phase, vec_pot, vec_pot_var)
662 character(len=*),
intent(in) :: dir
664 integer,
intent(in) :: nst
666 type(
ions_t),
target,
intent(in) :: ions
667 class(
mesh_t),
intent(in) :: mesh
669 type(
phase_t),
intent(in) :: phase
670 real(real64),
optional,
allocatable,
intent(in) :: vec_pot(:)
671 real(real64),
optional,
allocatable,
intent(in) :: vec_pot_var(:, :)
674 integer :: idir, ist, ik, ns, is,npath
675 integer,
allocatable :: iunit(:)
676 real(real64) :: red_kpoint(mesh%box%dim)
677 character(len=80) :: filename
679 logical :: projection
680 integer :: ii, ll, mm, nn, work, norb, work2
681 integer :: ia, iorb, idim, maxnorb
682 real(real64),
allocatable :: dpsi(:,:), ddot(:,:)
683 complex(real64),
allocatable :: zpsi(:,:), zdot(:,:)
684 real(real64),
allocatable :: weight(:,:,:,:,:)
692 if (st%d%nspin == 2) ns = 2
702 call parse_variable(namespace,
'BandStructureComputeProjections', .false., projection)
706 safe_allocate(iunit(0:ns-1))
710 write(filename,
'(a,i1.1,a)')
'bandstructure-sp', is+1
712 write(filename,
'(a)')
'bandstructure'
714 iunit(is) =
io_open(trim(dir)//
'/'//trim(filename), namespace, action=
'write')
717 write(iunit(is),
'(a)',advance=
'no')
'# coord. '
718 do idir = 1, mesh%box%dim
719 write(iunit(is),
'(3a)',advance=
'no')
'k',
index2axis(idir),
' '
721 if (.not. projection)
then
722 write(iunit(is),
'(a,i6,3a)')
'(red. coord.), bands:', nst,
' [', trim(
units_abbrev(
units_out%energy)),
']'
724 write(iunit(is),
'(a,i6,3a)',advance=
'no')
'(red. coord.), bands:', nst,
' [', trim(
units_abbrev(
units_out%energy)),
'] '
725 do ia = 1, ions%natoms
729 write(iunit(is),
'(a, i3.3,a,i1.1,a)',advance=
'no')
'w(at=',ia,
',os=',norb,
') '
732 write(iunit(is),
'(a)')
''
737 npath = kpoints%nkpt_in_path()*ns
744 safe_allocate(dpsi(1:mesh%np, 1:st%d%dim))
746 safe_allocate(zpsi(1:mesh%np, 1:st%d%dim))
750 do ia = 1, ions%natoms
754 safe_allocate(weight(1:st%nik,1:st%nst, 1:(2*
max_l+1), 1:maxnorb, 1:ions%natoms))
755 weight(1:st%nik,1:st%nst, 1:(2*
max_l+1), 1:maxnorb, 1:ions%natoms) =
m_zero
757 do ia = 1, ions%natoms
768 do iorb = 1, ions%atom(ia)%species%get_niwfs()
769 call ions%atom(ia)%species%get_iwf_ilm(iorb, 1, ii, ll, mm)
770 call ions%atom(ia)%species%get_iwf_n(iorb, 1, nn)
781 os%use_submesh = .false.
782 os%spec => ions%atom(ia)%species
784 do iorb = 1, os%norbs
788 ions%atom(ia)%species, mesh, os%sphere, os%ii, os%ll, os%jj, &
789 os, iorb, os%radius, os%ndim, use_mesh=.not.os%use_submesh, &
793 ions%atom(ia)%species, mesh, os%sphere, os%ii, os%ll, os%jj, &
794 os, iorb, os%radius, os%ndim, &
795 use_mesh = .not. phase%is_allocated() .and. .not. os%use_submesh, &
800 if (phase%is_allocated())
then
802 safe_allocate(os%phase(1:os%sphere%np, st%d%kpt%start:st%d%kpt%end))
804 if (.not. os%use_submesh)
then
805 safe_allocate(os%eorb_mesh(1:mesh%np, 1:os%norbs, 1:os%ndim, st%d%kpt%start:st%d%kpt%end))
806 os%eorb_mesh(:,:,:,:) =
m_zero
808 safe_allocate(os%eorb_submesh(1:os%sphere%np, 1:os%ndim, 1:os%norbs, st%d%kpt%start:st%d%kpt%end))
809 os%eorb_submesh(:,:,:,:) =
m_zero
813 os%ldorbs = max(
pad_pow2(os%sphere%np), 1)
814 if(.not. os%use_submesh) os%ldorbs = max(
pad_pow2(os%sphere%mesh%np), 1)
816 safe_allocate(os%buff_eorb(st%d%kpt%start:st%d%kpt%end))
817 do ik= st%d%kpt%start, st%d%kpt%end
823 vec_pot, vec_pot_var)
827 safe_allocate(ddot(1:st%d%dim,1:os%norbs))
829 safe_allocate(zdot(1:st%d%dim,1:os%norbs))
832 do ist = st%st_start, st%st_end
833 do ik = st%d%kpt%start, st%d%kpt%end
834 if (ik < st%nik - npath + 1) cycle
838 do iorb = 1, os%norbs
839 do idim = 1, st%d%dim
840 weight(ik,ist,iorb,norb,ia) = weight(ik,ist,iorb,norb,ia) + abs(ddot(idim,iorb))**2
845 if (phase%is_allocated())
then
847 call phase%apply_to_single(zpsi, mesh%np, st%d%dim, ik, .false.)
850 do iorb = 1, os%norbs
851 do idim = 1, st%d%dim
852 weight(ik,ist,iorb,norb,ia) = weight(ik,ist,iorb,norb,ia) + abs(zdot(idim,iorb))**2
859 safe_deallocate_a(ddot)
860 safe_deallocate_a(zdot)
865 if (st%parallel_in_states .or. st%d%kpt%parallel)
then
870 safe_deallocate_a(dpsi)
871 safe_deallocate_a(zpsi)
877 do ik = st%nik-npath+1, st%nik, ns
879 red_kpoint(1:mesh%box%dim) = kpoints%get_point(st%d%get_kpoint_index(ik + is), absolute_coordinates=.false.)
880 write(iunit(is),
'(1x)',advance=
'no')
881 if (st%nik > npath)
then
883 st%d%get_kpoint_index(ik + is)-st%d%get_kpoint_index(st%nik -npath))
886 st%d%get_kpoint_index(ik + is))
888 do idir = 1, mesh%box%dim
889 write(iunit(is),
'(f14.8)',advance=
'no') red_kpoint(idir)
895 do ia = 1, ions%natoms
901 write(iunit(is),
'(es15.8)',advance=
'no') weight(ik+is,ist,iorb,norb,ia)
910 write(iunit(is),
'(a)')
920 safe_deallocate_a(weight)
923 safe_deallocate_a(iunit)
This is the common interface to a sorting routine. It performs the shell algorithm,...
double exp(double __x) __attribute__((__nothrow__
pure logical function, public accel_is_enabled()
integer, parameter, public accel_mem_read_only
integer, parameter, public max_l
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.
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.
integer, parameter, public unpolarized
Parameters...
integer, parameter, public spinors
integer, parameter, public spin_polarized
real(real64), parameter, public m_two
real(real64), parameter, public m_huge
real(real64), parameter, public m_zero
complex(real64), parameter, public m_zi
real(real64), parameter, public m_epsilon
real(real64), parameter, public m_half
real(real64), parameter, public m_one
real(real64), parameter, public m_min_occ
Minimal occupation that is considered to be non-zero.
This module implements the underlying real-space grid.
subroutine, public io_close(iunit, grp)
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
real(real64) pure function, public kpoints_get_path_coord(this, ind)
This module is intended to contain "only mathematical" functions and procedures.
integer pure function, public pad_pow2(size)
create array size, which is padded to powers of 2
This module defines various routines, operating on mesh functions.
complex(real64) function, public zmf_integrate(mesh, ff, mask, reduce)
Integrate a function on the mesh.
real(real64) function, public dmf_integrate(mesh, ff, mask, reduce)
Integrate a function on the mesh.
This module defines the meshes, which are used in Octopus.
subroutine, public messages_warning(no_lines, all_nodes, namespace)
subroutine, public messages_new_line()
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
logical function mpi_grp_is_root(grp)
Is the current MPI process of grpcomm, root.
type(mpi_grp_t), public mpi_world
subroutine, public orbitalset_init(this)
subroutine, public orbitalset_end(this)
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...
subroutine, public dorbitalset_get_coefficients(os, ndim, psi, ik, has_phase, dot, reduce)
subroutine, public zorbitalset_get_coefficients(os, ndim, psi, ik, has_phase, dot, reduce)
integer function, public orbitalset_utils_count(species, iselect)
Count the number of orbital sets we have for a given atom.
integer function, public parse_block(namespace, name, blk, check_varinfo_)
subroutine, public smear_write_info(this, namespace, iunit)
integer, parameter, public smear_semiconductor
This module is intended to contain "only mathematical" functions and procedures.
pure logical function, public states_are_real(st)
This module handles spin dimensions of the states and the k-point distribution.
This module defines routines to write information about states.
subroutine, public states_elec_write_tpa(dir, namespace, space, mesh, st)
calculate the transition potential and write to a file
subroutine, public states_elec_write_eigenvalues(nst, st, space, kpoints, error, st_start, compact, iunit, namespace)
write the eigenvalues for some states to a file.
subroutine, public states_elec_write_gaps(iunit, st, space)
calculate gaps and write to a file.
subroutine, public states_elec_write_bandstructure(dir, namespace, nst, st, ions, mesh, kpoints, phase, vec_pot, vec_pot_var)
calculate and write the bandstructure
type(type_t), public type_cmplx
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
character(len=20) pure function, public units_abbrev(this)
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
This module is intended to contain simple general-purpose utility functions and procedures.
character pure function, public index2axis(idir)
Describes mesh distribution to nodes.
A container for the phase.
The states_elec_t class contains all electronic wave functions.